Quantum Phase Estimation: The Missing Link Between the QFT and Shor's Algorithm

Our QFT tutorial ended on a deliberately unsatisfying note: the Quantum Fourier Transform is fast, exact, and useless on its own. You can't load its input efficiently and you can't read its output — measure a transformed basis state and you get a uniformly random bitstring.

This post is the sequel that resolves that. Quantum phase estimation (QPE) is the machine that turns the inverse QFT into a readout primitive. It's the reason the QFT matters, and it is the algorithm sitting underneath Shor's — Shor's is QPE applied to modular exponentiation, plus number theory.

Read the QFT post first if you haven't; this one assumes the construction and won't rebuild it.

All code below was executed against Qiskit 2.5.0 and Qiskit Aer 0.17.2 on Python 3.11. The outputs shown are real.

The problem QPE solves

Given a unitary U and an eigenstate |ψ⟩ satisfying

U|ψ⟩ = e^{2πiθ} |ψ⟩

estimate θ ∈ [0, 1).

That sounds like a linear-algebra exercise, not an application. It is in fact one of the most consequential problems in the field, because enormous numbers of questions are secretly eigenvalue questions:

• What is a molecule's ground-state energy? Build U = e^{-iHt} for the molecular Hamiltonian H. Its eigenphases are the energy levels. Estimating θ estimates the energy.
• What is the order of a mod N? Build U|y⟩ = |ay mod N⟩. Its eigenphases are multiples of 1/r where r is the order. Estimating θ gives you r — and factoring falls out.
• How many solutions does an oracle mark? Amplitude estimation is QPE on the Grover operator, whose eigenphase encodes the count. That's the quantum-counting subroutine we flagged as a missing piece in the Grover tutorial.

One primitive, three flagship algorithms. That's why "estimate an eigenphase" is worth a dedicated circuit.

Note the precondition, which people gloss over: you need to have the eigenstate |ψ⟩, or at least a state with decent overlap with it. Preparing good eigenstates is its own hard problem, and in chemistry it's often the binding constraint.

Phase kickback: the key insight

Here's the thing that makes QPE work, and it's worth slowing down for.

A controlled-U applied to a control qubit in state |1⟩ and a target in the eigenstate |ψ⟩ gives:

CU |1⟩|ψ⟩ = |1⟩ ⊗ U|ψ⟩ = |1⟩ ⊗ e^{2πiθ}|ψ⟩ = e^{2πiθ} |1⟩|ψ⟩

The eigenstate came out unchanged. The phase e^{2πiθ} is a scalar, so it's now attached to the whole state — and since it only appeared for the |1⟩ branch of the control, it's really a relative phase on the control qubit. Put the control in superposition with a Hadamard and you get:

CU (|0⟩+|1⟩)/√2 ⊗ |ψ⟩ = (|0⟩ + e^{2πiθ}|1⟩)/√2 ⊗ |ψ⟩

The phase kicked back from the target onto the control. The target register is untouched and unentangled — it's a catalyst. All of the information now lives in the control's phase.

Do this with n control qubits, and give control qubit k a controlled-U^(2^k) instead, so it picks up phase e^{2πi·2^k·θ}. The counting register ends up in:

(1/√N) Σ_{x=0}^{N-1} e^{2πi·x·θ} |x⟩

Stare at that and compare it to the QFT formula from the previous post: QFT|j⟩ = (1/√N) Σ_x e^{2πi·jx/N}|x⟩. They are the same state, with j/N = θ. The counting register is holding the QFT of the number j = Nθ.

Why the inverse QFT

So the controlled-U powers have written θ into the counting register's phases — and, as the QFT post hammered, phases are invisible to measurement. Measure now and you get uniform noise.

But we know exactly what transform produced this state. So we undo it. Apply QFT⁻¹ and the register collapses to |j⟩ where j = Nθ, a plain binary integer you can read off in one shot.

That's the payoff the QFT post set up. The QFT direction moves position information into phases; the inverse direction pulls phase information back out into a readable number. QPE is the canonical example of the pattern we described there — the QFT is never the whole algorithm, it's the readout stage of a larger circuit.

The full recipe:

1. n counting qubits + enough qubits to hold |ψ⟩.
2. Hadamard the counting register.
3. For k = 0…n−1: apply controlled-U^(2^k) from counting qubit k onto the eigenstate register.
4. Apply the inverse QFT to the counting register.
5. Measure the counting register. The result, read as a binary fraction, is θ.

Building it in Qiskit

Simplest possible U with a known answer: the phase gate P(λ), which maps |1⟩ → e^{iλ}|1⟩. So |1⟩ is an eigenstate, and setting λ = 2πθ makes the eigenphase exactly θ. We'll reuse the qft function from the previous post and just call .inverse() on it.

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


def qft_rotations(circuit, n):
    if n == 0:
        return circuit
    n -= 1
    circuit.h(n)
    for qubit in range(n):
        circuit.cp(np.pi / 2 ** (n - qubit), qubit, n)
    qft_rotations(circuit, n)
    return circuit


def qft(n):
    qc = QuantumCircuit(n, name="QFT")
    qft_rotations(qc, n)
    for q in range(n // 2):          # the swaps matter here — see below
        qc.swap(q, n - q - 1)
    return qc


def qpe_phase_gate(n_count, theta):
    """QPE for U = P(2*pi*theta), whose eigenstate is |1>."""
    qc = QuantumCircuit(n_count + 1, n_count)
    qc.h(range(n_count))             # counting register in superposition
    qc.x(n_count)                    # prepare the eigenstate |1>

    for k in range(n_count):         # controlled-U^(2^k) = CP(2*pi*theta*2^k)
        qc.cp(2 * np.pi * theta * 2 ** k, k, n_count)

    qc.compose(qft(n_count).inverse(), qubits=range(n_count), inplace=True)
    qc.measure(range(n_count), range(n_count))
    return qc


sim = AerSimulator()
qc = qpe_phase_gate(4, 0.625)        # theta = 0.625 = 0.1010 in binary
counts = sim.run(transpile(qc, sim), shots=4096).result().get_counts()

for bits, c in sorted(counts.items(), key=lambda kv: -kv[1]):
    est = int(bits, 2) / 2 ** 4
    print(f"{bits} -> theta = {int(bits, 2)}/16 = {est:.4f}   counts = {c}")

Real output:

1010 -> theta = 10/16 = 0.6250   counts = 4096

4096 out of 4096 shots. θ = 0.625 is exactly representable in 4 bits (0.1010₂ = 1/2 + 1/8), so QPE returns it with certainty. Three counting qubits and θ = 0.125 behaves the same way:

[('001', 4096)]

001 = 1, and 1/8 = 0.125. Exact again.

The endianness trap

int(bits, 2) / 2**n worked above, and it is worth understanding why rather than copying it. Qiskit prints bitstrings with classical bit 0 on the right, and our counting qubit 0 — the one that got the smallest rotation, U^1 — is the least significant bit of j. The QFT's final swap network is what lines those two conventions up.

Drop the swaps from qft() — a tempting "optimization", since they look cosmetic — and rerun the exact θ = 0.625 case. Real output over 2048 shots:

{'1101': 1007, '0101': 470, '0001': 251, '0100': 80, '1100': 77,
 '0010': 66, '1110': 64, '0110': 16, '1010': 9, '1001': 8}

The correct answer 1010 now gets 9 shots out of 2048. The distribution smeared across the register with its largest peak at 1101, a confidently wrong 0.8125. Nothing errors. Nothing warns. You just publish the wrong number, and with a θ you didn't already know you would have no way to tell. This is the same silent-failure class as the endianness bug in the Grover oracle, and it's why the QFT post insisted on keeping the swaps in a reusable block.

Always sanity-check a QPE implementation against a θ you already know before pointing it at one you don't.

When θ doesn't fit in n bits

Real eigenphases are not tidy binary fractions. Try θ = 0.3, which has an infinite binary expansion, with 4 counting qubits:

qc = qpe_phase_gate(4, 0.3)
counts = sim.run(transpile(qc, sim), shots=4096).result().get_counts()
for bits, c in sorted(counts.items(), key=lambda kv: -kv[1])[:6]:
    print(f"{bits}  theta_est = {int(bits, 2) / 2**4:.4f}   p = {c / 4096:.4f}")

Real output:

0101  theta_est = 0.3125   p = 0.8823
0100  theta_est = 0.2500   p = 0.0542
0110  theta_est = 0.3750   p = 0.0220
0011  theta_est = 0.1875   p = 0.0093
0111  theta_est = 0.4375   p = 0.0063
0010  theta_est = 0.1250   p = 0.0049

The distribution is peaked, not delta-shaped. 88% of shots land on 5/16 = 0.3125, the closest 4-bit value to 0.3. The rest spread out, decaying fast as you move away from the peak. This is the generic behaviour: you get the nearest n-bit value with high probability, not certainty. The standard guarantee is that the best n-bit estimate appears with probability at least 4/π² ≈ 40.5%; here we're comfortably above that.

Two things improve together when you add counting qubits — finer resolution and a tighter peak. Same θ = 0.3 with 8 counting qubits:

01001101  theta_est = 0.30078   p = 0.8718
01001100  theta_est = 0.29688   p = 0.0596
01001110  theta_est = 0.30469   p = 0.0276
01001111  theta_est = 0.30859   p = 0.0076
01001011  theta_est = 0.29297   p = 0.0073

Resolution went from ±0.031 to ±0.002, and the peak stayed at 87%. In general, to get n bits of precision with success probability 1 − ε you use n + O(log(1/ε)) counting qubits — a handful of extra ancillas buys you a lot of confidence.

The cost is the part people underestimate. Counting qubit k needs U applied 2^k times. Total controlled-U applications across the register: 2ⁿ − 1. Every bit of precision doubles the circuit depth. That's fine when U^(2^k) has an efficient closed form — for a phase gate it's just one gate with a bigger angle, which is why the demo above is cheap. For modular exponentiation in Shor's it's a genuine arithmetic circuit repeated exponentially many times, and the depth is brutal.

Where QPE actually gets used

Shor's algorithm. Take U|y⟩ = |ay mod N⟩. Its eigenphases are s/r for integer s, where r is the multiplicative order of a mod N. Run QPE, measure s/r as a binary fraction, recover r by continued fractions, and factor N with a GCD. That's the whole quantum part — everything else in Shor's 1994 paper is number theory. Order finding is QPE; QPE is controlled powers plus an inverse QFT. The chain from the QFT post now closes. (The algorithm was recently machine-verified in Lean, if you want the formal treatment.)

Quantum chemistry. Set U = e^{-iHt}. QPE returns the ground-state energy to whatever precision you're willing to pay depth for — and crucially, it returns it with a provable error bound rather than a variational estimate. This is the long-term goal of quantum simulation, and it's the application most likely to matter commercially before cryptanalysis does.

Contrast with VQE. VQE attacks the same chemistry problem from the opposite direction. Instead of one deep coherent circuit, it uses a shallow variational circuit plus a classical optimizer, as introduced in Peruzzo et al. 2014. The trade is explicit:

VQE runs today and gives you an answer you can't fully trust. QPE gives you an answer you can trust and can't run today. Our VQE walkthrough in PennyLane covers the near-term side; PennyLane is generally the better fit for variational work while Qiskit suits structured circuits like this one.

The honest caveat

QPE is not a NISQ algorithm and no amount of clever transpilation will make it one.

The circuit must stay coherent through 2ⁿ − 1 controlled-U applications, each of which decomposes into many native two-qubit gates. With two-qubit error around 0.5% on current hardware, a few hundred entangling gates leaves you with noise. Decoherence sets a hard wall long before you reach chemically useful precision. Published hardware demonstrations of QPE are toy instances with 2–4 counting qubits — real, but not computations anyone needed done.

QPE becomes practical when circuits run on logical qubits rather than physical ones. Google's below-threshold result is the first real evidence that path is open, and our fault tolerance explainer covers what still has to happen. Until then: develop QPE on simulators, where you can verify it exactly, and reach for VQE or error mitigation if you need results from today's hardware.

That's not a reason to skip learning it. QPE is the template for what a fault-tolerant algorithm looks like — deep, exact, and structured — and it's the piece that makes the QFT worth building.

Next steps

• Compare the SDKs — build the same QPE circuit in Cirq or Braket
• Courses — structured paths through the algorithm canon
• Glossary — quantum gate, ancilla qubit, and the rest of the vocabulary
• Extend the code above to a 2-qubit U with two eigenstates, feed in a superposition of both, and watch the counting register entangle with the eigenstate register — that's the mechanism behind order finding.