Week 7 | reading the eigenvalue of a unitary to t bits with t control qubits and an inverse QFT | reading: Watrous Ch. 7 (phase estimation); Nielsen–Chuang §5.2
Lecture 13 showed that the QFT turns a pattern of phases \(\omega^{xy}\) into a single basis state. This lecture produces such a pattern from an unknown unitary and reads off its eigenvalue. The resulting algorithm, phase estimation, is the engine inside Shor's algorithm (Lecture 15), quantum chemistry's energy estimation, and many other algorithms. By the end you should be able to state the problem, build the circuit, explain why it is exact when the phase has \(t\) bits, compute its output distribution when it does not, choose the number of control qubits for a target precision and confidence, and say what happens when the input is not an eigenvector.
Every eigenvalue of a unitary has modulus 1 (Lecture 2), so it is a point on the unit circle and can be written \(e^{2\pi i\theta}\) for a unique \(\theta\in[0,1)\): the fraction of a full turn.
Given a unitary \(U\) on \(m\) qubits, which we can apply in controlled form, and a state \(|\psi\rangle\) with
$$U|\psi\rangle=e^{2\pi i\theta}|\psi\rangle,\qquad 0\le\theta\lt1,$$estimate the eigenphase \(\theta\). We call an estimate \(\tilde\theta\) accurate to \(n\) bits if \(|\tilde\theta-\theta|\lt2^{-n}\), with the distance measured around the circle (so 0.98 and 0.01 are 0.03 apart).
A unitary has the eigenvalue \(-i\). What is its eigenphase \(\theta\in[0,1)\)? And how far apart, around the circle, are the estimates 0.95 and 0.02?
Two features make this a quantum problem. First, \(U\) is an operation, not a matrix written down: it may act on hundreds of qubits, so we cannot diagonalise it classically. Second, a global phase is invisible (Lecture 3), so \(U|\psi\rangle\) and \(|\psi\rangle\) give identical measurement statistics. The phase \(e^{2\pi i\theta}\) only becomes visible as a relative phase, which is what control qubits provide.
We hold the eigenvector \(|\psi\rangle\) and apply \(U\) once, with no control qubit. Measured in any basis we like, the state \(U|\psi\rangle\) gives statistics that are .
Lecture 6's phase kickback (Theorem 6.5) says that a controlled-\(U\) acting on \(|1\rangle\otimes|\psi\rangle\) multiplies the state by the eigenvalue, and that this phase lands on the control. If \(U|\psi\rangle=e^{2\pi i\theta}|\psi\rangle\), then \(U^k|\psi\rangle=e^{2\pi ik\theta}|\psi\rangle\) as well, so a controlled power of \(U\) kicks back a multiple of the phase.
Let \(U|\psi\rangle=e^{2\pi i\theta}|\psi\rangle\) and let \(C\text{-}U^{k}\) be \(U^{k}\) controlled by one qubit. Then
$$C\text{-}U^{k}\Big(\tfrac{1}{\sqrt2}\big(|0\rangle+|1\rangle\big)\otimes|\psi\rangle\Big)=\tfrac{1}{\sqrt2}\big(|0\rangle+e^{2\pi i\,k\theta}|1\rangle\big)\otimes|\psi\rangle .$$The target register is unchanged and stays in a product state with the control.
By linearity, \(C\text{-}U^k\) leaves \(|0\rangle|\psi\rangle\) alone and maps \(|1\rangle|\psi\rangle\) to \(|1\rangle U^k|\psi\rangle=e^{2\pi ik\theta}|1\rangle|\psi\rangle\). Adding the two halves and factoring out \(|\psi\rangle\) gives the result.
The phase \(k\theta\) only matters mod 1, and that is where the information is. With \(k=2^j\), multiplying \(\theta=0.b_1b_2b_3\ldots\) (in binary) by \(2^j\) shifts the binary point \(j\) places: \(2^j\theta\equiv0.b_{j+1}b_{j+2}\ldots\pmod1\). So the control of \(C\text{-}U^{2^j}\) carries the digits of \(\theta\) from position \(j+1\) on. Using \(t\) control qubits with \(k=1,2,4,\dots,2^{t-1}\) records \(\theta\) at \(t\) different "zoom levels", just as the turntable in the cartoon is checked after 1, 2, 4 and 8 time steps.
The eigenphase is \(\theta=0.1101\) in binary. The control qubit of \(C\text{-}U^{4}\) ends in \(\tfrac1{\sqrt2}\big(|0\rangle+e^{2\pi i\varphi}|1\rangle\big)\). What is \(\varphi\in[0,1)\)?
Use \(t\) control qubits \(0,\dots,t-1\), each starting in \(|0\rangle\), and a target register holding \(|\psi\rangle\). Write \(T=2^t\).
After step 2, the control register is in the state
$$|\phi\rangle=\frac{1}{\sqrt T}\sum_{x=0}^{T-1}e^{2\pi i\,\theta x}\,|x\rangle ,$$unentangled from the target. If \(\theta=y_0/T\) for an integer \(y_0\) (that is, \(\theta\) has an exact \(t\)-bit binary expansion), then \(|\phi\rangle=\mathrm{QFT}_T|y_0\rangle\), and the algorithm outputs \(y_0\), hence \(\tilde\theta=\theta\), with probability 1.
By Proposition 14.1 control qubit \(j\) ends in \(\tfrac{1}{\sqrt2}\big(|0\rangle+e^{2\pi i\theta 2^j}|1\rangle\big)\). Multiply out the \(t\) factors: the term with control bits \(x_{t-1}\cdots x_0\) picks up the phase \(\prod_je^{2\pi i\theta x_j2^j}=e^{2\pi i\theta x}\), where \(x=\sum_jx_j2^j\). That gives \(|\phi\rangle\). If \(\theta=y_0/T\), then \(e^{2\pi i\theta x}=e^{2\pi i\,y_0x/T}=\omega^{y_0x}\) with \(\omega=e^{2\pi i/T}\), so \(|\phi\rangle=\tfrac{1}{\sqrt T}\sum_x\omega^{y_0x}|x\rangle\), which is \(\mathrm{QFT}_T|y_0\rangle\) by Definition 13.1. Applying \(\mathrm{QFT}_T^{-1}\) gives \(|y_0\rangle\).
Phase estimation with \(t=4\) control qubits is run on an eigenvector with \(\theta=\tfrac{11}{16}\). Which outcome \(y\) does it return, and with what probability?
Suppose \(U|\psi\rangle=e^{2\pi i\cdot5/8}|\psi\rangle\). What does each control qubit hold after the controlled powers, and what does the measurement return?
\(\theta=5/8=0.101\) in binary. Control 0 gets the phase \(e^{2\pi i\cdot5/8}=e^{5\pi i/4}\); control 1 gets \(e^{2\pi i\cdot10/8}=e^{2\pi i/4}=i\) (one whole turn discarded); control 2 gets \(e^{2\pi i\cdot20/8}=e^{\pi i}=-1\). Control 2 is therefore \(|-\rangle\), which records the last bit of \(\theta\), and control 1 is \(|{+i}\rangle\). Since \(\theta=5/8\) is a 3-bit fraction, Theorem 14.2 says the inverse QFT returns \(y=5=101_2\) with certainty, and \(\tilde\theta=5/8\) exactly.
With \(t=1\) the circuit is \(H\), a controlled-\(U\), \(H\): the inverse QFT on one qubit is just \(H\). For \(\theta\in\{0,\tfrac12\}\) it decides which one with certainty, which is Deutsch's algorithm of Lecture 9 in disguise.
Most phases are not multiples of \(1/T\). Then \(|\phi\rangle\) is not a Fourier basis state, and the output is random, but concentrated near \(T\theta\).
For every \(\theta\in[0,1)\), phase estimation with \(t\) control qubits returns \(y\) with probability
$$P(y)=\frac{1}{T^2}\Big|\sum_{x=0}^{T-1}e^{2\pi i\,x(\theta-y/T)}\Big|^2=\frac{\sin^2(\pi T\delta)}{T^2\sin^2(\pi\delta)},\qquad\delta=\theta-\frac yT ,$$with the value 1 when \(\delta\) is an integer.
By Proposition 13.1 the amplitude of \(|y\rangle\) in \(\mathrm{QFT}_T^{-1}|\phi\rangle\) is \(\tfrac{1}{\sqrt T}\sum_x\omega^{-xy}\langle x|\phi\rangle=\tfrac1T\sum_xe^{2\pi ix(\theta-y/T)}\). This is a geometric series with ratio \(q=e^{2\pi i\delta}\). If \(q\neq1\), its sum is \(\frac{q^T-1}{q-1}\), and \(|e^{i\alpha}-1|=2|\sin(\alpha/2)|\) gives \(\big|\frac{q^T-1}{q-1}\big|=\frac{|\sin(\pi T\delta)|}{|\sin(\pi\delta)|}\).
Take \(t=1\) (so \(T=2\)) and \(\theta=0.45\). Use Theorem 14.3 to find the probability of the output \(y=1\). (Three decimals is fine.)
Find the most likely outputs for \(\theta=0.3\) with \(t=3\) and with \(t=6\).
With \(t=3\), \(T\theta=2.4\), between the grid points 2 and 3, closer to 2. Theorem 14.3 gives
$$P(2)=0.578,\qquad P(3)=0.259,\qquad P(1)=0.052,\qquad P(4)=0.041,$$and the remaining four outcomes share about 0.07. The best 3-bit estimate, \(y=2\) (\(\tilde\theta=0.25\)), comes out 58% of the time, and one of the two neighbours of \(2.4\) comes out 84% of the time. With \(t=6\), \(T\theta=19.2\), and \(P(19)=0.875\), \(P(20)=0.055\), \(P(18)=0.024\): the estimate \(19/64\approx0.297\) is much more precise and also more likely, because \(19.2\) sits closer to a grid point than \(2.4\) did.
How bad can it get? The worst case for the nearest grid point is when \(T\theta\) sits exactly half-way between two of them.
Let \(y^*\) be the integer nearest to \(T\theta\) (mod \(T\)). Then phase estimation returns \(y^*\) with probability at least \(4/\pi^2\approx0.405\), for every \(\theta\) and every \(t\). It returns one of the two grid points adjacent to \(T\theta\) with probability at least \(8/\pi^2\approx0.81\).
For \(y^*\), \(|\delta|=|\theta-y^*/T|\le\tfrac1{2T}\), so \(\pi T|\delta|\le\tfrac\pi2\). On \([0,\tfrac\pi2]\) the sine lies above its chord: \(\sin u\ge\tfrac2\pi u\). So \(\sin^2(\pi T\delta)\ge4T^2\delta^2\). Also \(\sin^2(\pi\delta)\le\pi^2\delta^2\). Theorem 14.3 gives
$$P(y^*)\ge\frac{4T^2\delta^2}{T^2\pi^2\delta^2}=\frac4{\pi^2}.$$For the two neighbours \(y_-\le T\theta\le y_+\), with \(T\theta-y_-=\Delta\in[0,1]\): here \(|T\delta|\) is \(\Delta\) for one and \(1-\Delta\) for the other, so \(\sin^2(\pi T\delta)=\sin^2(\pi\Delta)\) for both, and \(\sin^2(\pi\delta)\le\pi^2\delta^2\) gives \(P(y_-)+P(y_+)\ge\tfrac{\sin^2(\pi\Delta)}{\pi^2}\big(\tfrac1{\Delta^2}+\tfrac1{(1-\Delta)^2}\big)\), whose minimum over \(\Delta\) is at \(\Delta=\tfrac12\) and equals \(\tfrac{8}{\pi^2}\).
Tick every statement that holds for every \(\theta\) and every \(t\).
Theorem 14.4 gives the nearest \(t\)-bit estimate only with probability around 40%. In practice we want \(n\) correct bits with high confidence. Adding a few control qubits beyond \(n\) does both: the extra digits are thrown away, and the probability that the estimate is off by more than \(2^{-n}\) becomes small.
Let \(n\ge1\) and \(0\lt\varepsilon\lt1\). Phase estimation with
$$t=n+\Big\lceil\log_2\Big(2+\frac{1}{2\varepsilon}\Big)\Big\rceil$$control qubits returns an estimate with \(|\tilde\theta-\theta|\lt2^{-n}\) with probability at least \(1-\varepsilon\).
Let \(b\) be the integer with \(b\le T\theta\lt b+1\) and write \(T\theta=b+\Delta\), \(0\le\Delta\lt1\). An outcome \(y=b+k\) (with \(k\) taken in \((-T/2,T/2]\)) has \(T\delta=\Delta-k\). Using \(\sin^2(\pi T\delta)\le1\) and \(|\sin(\pi\delta)|\ge2|\delta|\) for \(|\delta|\le\tfrac12\), Theorem 14.3 gives
$$P(b+k)\le\frac{1}{4(\Delta-k)^2}.$$Fix an integer \(e\ge2\). For \(k\ge e+1\), \(|\Delta-k|\ge k-1\ge e\); for \(k\le-e\), \(|\Delta-k|\ge|k|\ge e\). So
$$P(|y-b|\gt e)\le2\sum_{j\ge e}\frac{1}{4j^2}\le\frac12\cdot\frac{1}{e-1},$$using \(\sum_{j\ge e}j^{-2}\le\int_{e-1}^\infty u^{-2}\,du\). If \(|y-b|\le e\) then \(|y/T-\theta|\le(|y-b|+\Delta)/T\lt(e+1)/T\). Choose \(e=2^{t-n}-1\), so that \((e+1)/T=2^{-n}\). The failure probability is at most \(\frac{1}{2(2^{t-n}-2)}\), which is at most \(\varepsilon\) exactly when \(2^{t-n}\ge2+\frac1{2\varepsilon}\). (This is the argument of Nielsen–Chuang §5.2.1.)
How many control qubits does Theorem 14.5 prescribe for a 6-bit estimate with failure probability at most \(\varepsilon=0.02\)?
How many control qubits give a 10-bit estimate with probability at least 0.99, and how many applications of \(U\) does that cost if each \(U^{2^j}\) is done as \(2^j\) copies of \(U\)?
\(2+\frac{1}{2\cdot0.01}=52\) and \(\lceil\log_252\rceil=6\), so \(t=10+6=16\). The controlled powers use \(1+2+4+\dots+2^{15}=2^{16}-1=65{,}535\) controlled applications of \(U\). The bound is conservative: for \(n=4\), \(\varepsilon=0.05\), it prescribes \(t=8\), and the exact worst case over all \(\theta\) (computed in verify_L14.py) is a success probability of 0.9875, not just 0.95.
There is a second way to buy confidence: run phase estimation with fewer control qubits several times and take the median of the estimates. Each run is correct to \(n\) bits with probability above \(\tfrac12\) (say 0.8), so by a Chernoff bound the median of \(O(\log(1/\varepsilon))\) runs is correct with probability \(1-\varepsilon\). Fewer qubits per run, more runs: which is better depends on the hardware.
Preparing an eigenvector is often as hard as the problem itself. Fortunately phase estimation is linear, so it handles superpositions of eigenvectors gracefully.
Let \(U|\psi_k\rangle=e^{2\pi i\theta_k}|\psi_k\rangle\) for an orthonormal family \(|\psi_k\rangle\), and run phase estimation on \(|\psi\rangle=\sum_kc_k|\psi_k\rangle\). Then the output is distributed as follows: choose \(k\) with probability \(|c_k|^2\), then output \(y\) with the distribution of Theorem 14.3 for \(\theta_k\). If the output is exact, the target collapses to an eigenvector with that phase.
By linearity and Theorem 14.2, before measurement the state is \(\sum_kc_k\,\big(\mathrm{QFT}^{-1}|\phi_k\rangle\big)\otimes|\psi_k\rangle\), where \(|\phi_k\rangle\) is the control state for \(\theta_k\). The amplitude of \(|y\rangle|\psi_k\rangle\) is \(c_k\,a_k(y)\), with \(|a_k(y)|^2\) the probability of Theorem 14.3. Because the \(|\psi_k\rangle\) are orthonormal, \(P(y)=\sum_k|c_k|^2|a_k(y)|^2\), which is the two-stage description. The target registers for different \(k\) never interfere, which is why the mixture is classical.
\(U\) has eigenvectors \(|\psi_1\rangle,|\psi_2\rangle\) with eigenphases \(\tfrac14\) and \(\tfrac34\). Phase estimation with \(t=2\) is run on \(0.6\,|\psi_1\rangle+0.8i\,|\psi_2\rangle\). What is the output distribution?
For example, \(U=P(2\pi\theta)=\mathrm{diag}(1,e^{2\pi i\theta})\) has eigenvectors \(|0\rangle\) (phase 0) and \(|1\rangle\) (phase \(\theta\)). Run on \(|+\rangle\), phase estimation returns 0 half the time and an estimate of \(\theta\) half the time; the second code box below shows it. The crucial use is next lecture.
Let \(N\) be the RSA modulus to factor and \(a\) a number coprime to \(N\). The unitary \(U_a|y\rangle=|ay\bmod N\rangle\) permutes the residues, and its eigenvalues are \(e^{2\pi is/r}\), where \(r\) is the order of \(a\) (the least \(r\gt0\) with \(a^r\equiv1\bmod N\)) and \(s=0,\dots,r-1\). We cannot prepare any one eigenvector without knowing \(r\), but the easy state \(|1\rangle\) is an equal superposition of them. By Proposition 14.6, phase estimation on \(|1\rangle\) returns an approximation of \(s/r\) for a uniformly random \(s\), and Lecture 15 recovers \(r\) from it with continued fractions. Knowing \(r\) factors \(N\). Everything that makes Shor's algorithm work is already in this lecture except two things: the efficient powers \(U_a^{2^j}\) (Section 7) and the number theory.
The QFT costs about \(t^2/2\) gates, which is cheap. The expensive part is step 2. Implemented naively, \(U^{2^j}\) is \(2^j\) copies of \(U\), and the circuit uses \(2^t-1\) controlled applications of \(U\) in total:
| control qubits \(t\) | 3 | 8 | 16 | 32 |
|---|---|---|---|---|
| applications of \(U\), naively | 7 | 255 | 65,535 | about \(4.3\times10^9\) |
Since the precision is \(2^{-t}\), naive phase estimation costs time proportional to \(1/\text{precision}\): no better than measuring the phase by brute repetition. Phase estimation is efficient only when \(U^{2^j}\) can be implemented with cost polynomial in \(j\). Two cases matter.
Order finding with \(N=55\) and \(a=2\). The controlled power \(U_2^{2^3}\) is the single multiplication \(U_c\). What is \(c\in\{0,\dots,54\}\)?
For a generic \(U\), such as the time evolution of a molecule, no shortcut exists, and phase estimation pays the full \(2^t\). This is why phase estimation for quantum chemistry is so expensive: every extra bit of the energy doubles the number of applications of \(U\).
import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit.library import QFTGate
from qiskit.quantum_info import Statevector
def qpe(theta, t):
"""phase estimation of U = P(2 pi theta) on its eigenvector |1>; controls are qubits 0..t-1"""
qc = QuantumCircuit(t + 1)
qc.x(t) # target = |1>, eigenvalue e^{2 pi i theta}
qc.h(range(t))
for j in range(t):
qc.cp(2 * np.pi * theta * 2 ** j, j, t) # controlled-U^(2^j) is one controlled phase here
qc.append(QFTGate(t).inverse(), range(t))
return qc
def closed_form(theta, t):
T = 2 ** t; x = np.arange(T)
return np.array([abs(np.exp(2j * np.pi * x * (theta - y / T)).sum() / T) ** 2 for y in range(T)])
theta = 0.3
for t in range(3, 9):
p = Statevector.from_instruction(qpe(theta, t)).probabilities(qargs=list(range(t)))
y = int(np.argmax(p))
print(f"t = {t}: best y = {y:3d}, estimate {y / 2**t:.5f}, P(best) = {p[y]:.3f}, "
f"matches closed form: {np.allclose(p, closed_form(theta, t))}")
t = 3: best y = 2, estimate 0.25000, P(best) = 0.578, matches closed form: True t = 4: best y = 5, estimate 0.31250, P(best) = 0.876, matches closed form: True t = 5: best y = 10, estimate 0.31250, P(best) = 0.573, matches closed form: True t = 6: best y = 19, estimate 0.29688, P(best) = 0.875, matches closed form: True t = 7: best y = 38, estimate 0.29688, P(best) = 0.573, matches closed form: True t = 8: best y = 77, estimate 0.30078, P(best) = 0.875, matches closed form: True
The circuit agrees with Theorem 14.3 for every \(t\). The probability of the best estimate alternates between 0.58 and 0.88 because \(2^t\cdot0.3\) alternates between fractional parts 0.4 (as in \(2.4\), \(9.6\), \(38.4\)) and 0.8 (as in \(4.8\), \(19.2\), \(76.8\)), that is, between 0.4 and 0.2 away from the nearest grid point. QFTGate(t).inverse() is the inverse QFT with the same qubit ordering as Lecture 13.
import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit.library import QFTGate
from qiskit.quantum_info import Statevector
t, theta = 4, 0.3
qc = QuantumCircuit(t + 1)
qc.h(t) # target = |+> = (|0> + |1>)/sqrt2: NOT an eigenvector
qc.h(range(t))
for j in range(t):
qc.cp(2 * np.pi * theta * 2 ** j, j, t)
qc.append(QFTGate(t).inverse(), range(t))
p = Statevector.from_instruction(qc).probabilities(qargs=list(range(t)))
print({y: round(float(p[y]), 3) for y in range(2 ** t) if p[y] > 0.02})
{0: 0.501, 4: 0.028, 5: 0.438}
Proposition 14.6 in action: \(y=0\) (the eigenphase 0 of \(|0\rangle\)) with probability one half, plus a little from the other branch, and otherwise the distribution for \(\theta=0.3\) at \(t=4\), with \(y=5\) (\(\tilde\theta=0.3125\)) most likely.
The midterm (next meeting, 75 minutes) covers Lectures 1–14. A formula sheet with the gate matrices, the Hadamard-transform formula, Grover's angle formulas and the QFT definition is provided. You should be able to do the following without notes.
| lectures | you can |
|---|---|
| 1–3 | compute measurement probabilities in the Z, X and Y bases; tell global from relative phase; place a state on the Bloch sphere; give a 95% interval for a probability from counts |
| 4 | multiply gate matrices; describe gates as rotations; choose the basis change for an X or Y measurement |
| 5–6 | compute tensor products in Qiskit's ordering; test entanglement with the determinant; carry out a partial measurement; build the Bell states |
| 7–8 | build classical functions reversibly with Toffolis; trace a teleportation branch; explain no-cloning; compute a CHSH correlator |
| 9–10 | use \(H^{\otimes n}|x\rangle=2^{-n/2}\sum_y(-1)^{x\cdot y}|y\rangle\); run Deutsch–Jozsa, Bernstein–Vazirani and one step of Simon; solve a small system over \(\mathbb F_2\) |
| 11–12 | compute Grover's success probability \(\sin^2((2k+1)\theta)\) and the iteration count; explain why Grover halves symmetric key length |
| 13–14 | apply the QFT definition and product formula; compute a phase-estimation output distribution; choose \(t\) for a given precision and confidence |
| all | read a short Qiskit snippet and predict its counts, including the bit order |
Next time: the midterm, then order finding and Shor's algorithm, and the launch of Lab 4. HW7 (Lectures 13–14) is due at the start of Week 8.
all numbers and figures in these notes are computed in python (verify_L14.py) and follow the page theme