QuantumCore

Lesson 29 · 100 min · advanced

Quantum Phase Estimation

Loading progress…

Given a unitary UU and an eigenvector ∣u⟩\lvert u\rangle with U∣u⟩=e2πiφ∣u⟩U\lvert u\rangle=e^{2\pi i\varphi}\lvert u\rangle, quantum phase estimation (QPE) writes a tt-bit approximation of φ∈[0,1)\varphi\in[0,1) into a counting register. This lesson derives the state 2−t/2∑ke2πiφk∣k⟩∣u⟩2^{-t/2}\sum_{k}e^{2\pi i\varphi k}\lvert k\rangle\lvert u\rangle produced by phase kickback from the controlled powers U2jU^{2^j}, shows that QFT†\mathrm{QFT}^\dagger returns ∣m⟩\lvert m\rangle exactly when φ=m/2t\varphi=m/2^t, proves the bound P(b)≥4/π2P(b)\ge4/\pi^2 for the nearest tt-bit estimate in the general case, proves that t=n+⌈log⁡2(2+1/(2ε))⌉t=n+\lceil\log_2(2+1/(2\varepsilon))\rceil counting qubits give nn correct bits with probability at least 1−ε1-\varepsilon, and explains how a superposition of eigenvectors makes QPE sample the spectrum, which is the engine of the HHL algorithm.

Prerequisites: The Quantum Fourier Transform

Objectives and key results

  • State the phase-estimation problem U∣u⟩=e2πiφ∣u⟩U\lvert u\rangle=e^{2\pi i\varphi}\lvert u\rangle and derive the counting state 2−t/2∑k=02t−1e2πiφk∣k⟩∣u⟩2^{-t/2}\sum_{k=0}^{2^t-1}e^{2\pi i\varphi k}\lvert k\rangle\lvert u\rangle by phase kickback, with the weight-2j2^j counting qubit controlling U2jU^{2^j}.

  • Prove that QFT†\mathrm{QFT}^\dagger maps the counting state to ∣m⟩\lvert m\rangle exactly when φ=m/2t\varphi=m/2^t.

  • For a non-dyadic phase, derive the outcome distribution as a geometric sum and prove P(b)≥4/π2P(b)\ge4/\pi^2 for b=round(2tφ)b=\mathrm{round}(2^t\varphi).

  • Prove the precision/success trade-off t=n+⌈log⁡2(2+12ε)⌉t=n+\lceil\log_2(2+\tfrac1{2\varepsilon})\rceil and explain how QPE on a superposition of eigenvectors samples eigenphases.

definition · Phase-estimation problem
U∣u⟩=e2πiφ∣u⟩,φ∈[0,1): output a t-bit estimate φ~=b/2tU\lvert u\rangle=e^{2\pi i\varphi}\lvert u\rangle,\quad\varphi\in[0,1):\ \text{output a }t\text{-bit estimate }\tilde\varphi=b/2^t
identity · Counting state after the controlled powers
12t/2⨂j=t−10(∣0⟩+e2πi2jφ∣1⟩)⊗∣u⟩=12t/2∑k=02t−1e2πiφk∣k⟩∣u⟩\frac1{2^{t/2}}\bigotimes_{j=t-1}^{0}\big(\lvert0\rangle+e^{2\pi i2^j\varphi}\lvert1\rangle\big)\otimes\lvert u\rangle=\frac1{2^{t/2}}\sum_{k=0}^{2^t-1}e^{2\pi i\varphi k}\lvert k\rangle\lvert u\rangle
theorem · Nearest-estimate bound
P(b)=∣sin⁡(π2tδ)2tsin⁡(πδ)∣2≥4π2,b=round(2tφ), δ=φ−b2tP(b)=\left\lvert\frac{\sin(\pi2^t\delta)}{2^t\sin(\pi\delta)}\right\rvert^2\ge\frac4{\pi^2},\qquad b=\mathrm{round}(2^t\varphi),\ \delta=\varphi-\frac b{2^t}
theorem · Precision versus success probability
t=n+⌈log⁡2(2+12ε)⌉ ⇒ P(∣m−b∣≤2t−n−1)≥1−εt=n+\Big\lceil\log_2\Big(2+\frac1{2\varepsilon}\Big)\Big\rceil\ \Rightarrow\ P\big(\lvert m-b\rvert\le2^{t-n}-1\big)\ge1-\varepsilon

The phase-estimation problem#

Definition

Phase estimation

Let UU be a unitary on mm qubits and ∣u⟩\lvert u\rangle an eigenvector of UU. Because UU is unitary, ∥U∣u⟩∥=∥∣u⟩∥\lVert U\lvert u\rangle\rVert=\lVert\lvert u\rangle\rVert, so its eigenvalue has modulus 11 and can be written e2πiφe^{2\pi i\varphi} with a unique φ∈[0,1)\varphi\in[0,1): U∣u⟩=e2πiφ∣u⟩.U\lvert u\rangle=e^{2\pi i\varphi}\lvert u\rangle. Phase estimation is the task: given the ability to apply controlled-U2jU^{2^j} for j=0,…,t−1j=0,\dots,t-1 and a copy of ∣u⟩\lvert u\rangle, output a tt-bit number b∈{0,…,2t−1}b\in\{0,\dots,2^t-1\} such that b/2tb/2^t approximates φ\varphi. The register holding bb is the counting register; it has tt qubits, qubit 00 is the most significant bit, so k=∑j=0t−1kt−1−j 2jk=\sum_{j=0}^{t-1}k_{t-1-j}\,2^j when qubit qq holds bit kqk_q. The counting qubit qq therefore carries the weight 2t−1−q2^{t-1-q}, and the last counting qubit q=t−1q=t-1 has weight 202^0.

Phase kickback and the counting state#

Lemma

Phase kickback

If U∣u⟩=e2πiφ∣u⟩U\lvert u\rangle=e^{2\pi i\varphi}\lvert u\rangle then for every integer p≥0p\ge0 the controlled operation C-Up=∣0⟩⟨0∣⊗I+∣1⟩⟨1∣⊗UpC\text{-}U^{p}=\lvert0\rangle\langle0\rvert\otimes I+\lvert1\rangle\langle1\rvert\otimes U^{p} satisfies C-Up 12(∣0⟩+∣1⟩)∣u⟩=12(∣0⟩+e2πipφ∣1⟩)∣u⟩.C\text{-}U^{p}\,\tfrac1{\sqrt2}\big(\lvert0\rangle+\lvert1\rangle\big)\lvert u\rangle=\tfrac1{\sqrt2}\big(\lvert0\rangle+e^{2\pi ip\varphi}\lvert1\rangle\big)\lvert u\rangle.

Proof (3 steps)
  1. Step 1:
    Up∣u⟩=Up−1e2πiφ∣u⟩=⋯=e2πipφ∣u⟩U^{p}\lvert u\rangle=U^{p-1}e^{2\pi i\varphi}\lvert u\rangle=\cdots=e^{2\pi ip\varphi}\lvert u\rangle
    Apply the eigenvalue equation pp times; scalars commute with UU. For p=0p=0 both sides are ∣u⟩\lvert u\rangle.
  2. Step 2:
    C-Up∣0⟩∣u⟩=∣0⟩∣u⟩,C-Up∣1⟩∣u⟩=∣1⟩Up∣u⟩=e2πipφ∣1⟩∣u⟩C\text{-}U^{p}\lvert0\rangle\lvert u\rangle=\lvert0\rangle\lvert u\rangle,\qquad C\text{-}U^{p}\lvert1\rangle\lvert u\rangle=\lvert1\rangle U^{p}\lvert u\rangle=e^{2\pi ip\varphi}\lvert1\rangle\lvert u\rangle
    Definition of the controlled operation on the two control basis states, then the previous step.
  3. Step 3:
    C-Up12(∣0⟩+∣1⟩)∣u⟩=12(∣0⟩+e2πipφ∣1⟩)∣u⟩C\text{-}U^{p}\tfrac1{\sqrt2}\big(\lvert0\rangle+\lvert1\rangle\big)\lvert u\rangle=\tfrac1{\sqrt2}\big(\lvert0\rangle+e^{2\pi ip\varphi}\lvert1\rangle\big)\lvert u\rangle
    Linearity. The scalar e2πipφe^{2\pi ip\varphi} belongs to the whole tensor product and may be attached to the ∣1⟩\lvert1\rangle of the control: the phase is "kicked back" onto the control while the target stays ∣u⟩\lvert u\rangle.

∎

Derivation

The state after H⊗tH^{\otimes t} and the controlled powers

  1. Step 1:
    ∣0⟩⊗t∣u⟩ →H⊗t 12t/2(∣0⟩+∣1⟩)⊗t∣u⟩\lvert0\rangle^{\otimes t}\lvert u\rangle\ \xrightarrow{H^{\otimes t}}\ \frac1{2^{t/2}}\big(\lvert0\rangle+\lvert1\rangle\big)^{\otimes t}\lvert u\rangle
    H∣0⟩=12(∣0⟩+∣1⟩)H\lvert0\rangle=\tfrac1{\sqrt2}(\lvert0\rangle+\lvert1\rangle) on each of the tt counting qubits.
  2. Step 2:
    qubit q controls U2t−1−q:q=0↦U2t−1, …, q=t−1↦U1\text{qubit }q\text{ controls }U^{2^{t-1-q}}:\quad q=0\mapsto U^{2^{t-1}},\ \dots,\ q=t-1\mapsto U^{1}
    The counting qubit of weight 2j2^j controls U2jU^{2^j}. In big-endian order (qubit 00 = most significant) qubit qq has weight 2t−1−q2^{t-1-q}, so the last counting qubit controls U1U^1 and the first controls U2t−1U^{2^{t-1}}.
  3. Step 3:
    12t/2(∣0⟩+e2πi2t−1φ∣1⟩)⊗(∣0⟩+e2πi2t−2φ∣1⟩)⊗⋯⊗(∣0⟩+e2πi20φ∣1⟩)⊗∣u⟩\frac1{2^{t/2}}\big(\lvert0\rangle+e^{2\pi i2^{t-1}\varphi}\lvert1\rangle\big)\otimes\big(\lvert0\rangle+e^{2\pi i2^{t-2}\varphi}\lvert1\rangle\big)\otimes\cdots\otimes\big(\lvert0\rangle+e^{2\pi i2^{0}\varphi}\lvert1\rangle\big)\otimes\lvert u\rangle
    Apply the kickback lemma with p=2t−1−qp=2^{t-1-q} to each counting qubit. The controlled operations act on different controls and leave ∣u⟩\lvert u\rangle unchanged, so they can be applied one after another and the state stays a product.
  4. Step 4:
    =12t/2⨂j=t−10(∣0⟩+e2πi2jφ∣1⟩)⊗∣u⟩=12t/2⨂j=t−10(∑kj′=01e2πiφ 2jkj′∣kj′⟩)⊗∣u⟩=\frac1{2^{t/2}}\bigotimes_{j=t-1}^{0}\big(\lvert0\rangle+e^{2\pi i2^{j}\varphi}\lvert1\rangle\big)\otimes\lvert u\rangle=\frac1{2^{t/2}}\bigotimes_{j=t-1}^{0}\Big(\sum_{k_j'=0}^{1}e^{2\pi i\varphi\,2^{j}k_j'}\lvert k_j'\rangle\Big)\otimes\lvert u\rangle
    Write each factor as a sum over its bit kj′∈{0,1}k_j'\in\{0,1\}, where kj′k_j' is the bit of weight 2j2^j (stored on qubit t−1−jt-1-j).
  5. Step 5:
    =12t/2∑kt−1′=01⋯∑k0′=01e2πiφ∑j2jkj′∣kt−1′⋯k0′⟩∣u⟩=\frac1{2^{t/2}}\sum_{k_{t-1}'=0}^{1}\cdots\sum_{k_0'=0}^{1}e^{2\pi i\varphi\sum_{j}2^{j}k_j'}\lvert k_{t-1}'\cdots k_0'\rangle\lvert u\rangle
    Expand the tensor product of sums by multilinearity; the product of the exponentials is the exponential of the sum of exponents.
  6. Step 6:
    =12t/2∑k=02t−1e2πiφk∣k⟩∣u⟩=\frac1{2^{t/2}}\sum_{k=0}^{2^t-1}e^{2\pi i\varphi k}\lvert k\rangle\lvert u\rangle
    k=∑j2jkj′k=\sum_j2^jk_j' runs over all integers 0,…,2t−10,\dots,2^t-1 exactly once as the bits range over {0,1}t\{0,1\}^t, and ∣kt−1′⋯k0′⟩=∣k⟩\lvert k_{t-1}'\cdots k_0'\rangle=\lvert k\rangle in big-endian order.

Result

∣Ψφ⟩:=12t/2∑k=02t−1e2πiφk∣k⟩\boxed{\lvert\Psi_\varphi\rangle:=\frac1{2^{t/2}}\sum_{k=0}^{2^t-1}e^{2\pi i\varphi k}\lvert k\rangle}

Exact readout for dyadic phases#

Theorem

Dyadic phases are read exactly

If φ=m/2t\varphi=m/2^t for an integer 0≤m<2t0\le m<2^t, then ∣Ψφ⟩=QFT2t∣m⟩\lvert\Psi_\varphi\rangle=\mathrm{QFT}_{2^t}\lvert m\rangle and therefore QFT2t†∣Ψφ⟩=∣m⟩\mathrm{QFT}_{2^t}^\dagger\lvert\Psi_\varphi\rangle=\lvert m\rangle: measuring the counting register returns mm with probability 11.

Proof (4 steps)
  1. Step 1:
    QFTN∣m⟩=1N∑k=0N−1e2πimk/N∣k⟩,N=2t\mathrm{QFT}_{N}\lvert m\rangle=\frac1{\sqrt N}\sum_{k=0}^{N-1}e^{2\pi imk/N}\lvert k\rangle,\qquad N=2^t
    Definition of the QFT from the previous lesson.
  2. Step 2:
    ∣Ψm/N⟩=12t/2∑k=0N−1e2πi(m/N)k∣k⟩=1N∑k=0N−1e2πimk/N∣k⟩\lvert\Psi_{m/N}\rangle=\frac1{2^{t/2}}\sum_{k=0}^{N-1}e^{2\pi i(m/N)k}\lvert k\rangle=\frac1{\sqrt N}\sum_{k=0}^{N-1}e^{2\pi imk/N}\lvert k\rangle
    Substitute φ=m/N\varphi=m/N into the counting state; 2t/2=N2^{t/2}=\sqrt N.
  3. Step 3:
    ∣Ψm/N⟩=QFTN∣m⟩\lvert\Psi_{m/N}\rangle=\mathrm{QFT}_N\lvert m\rangle
    The two expansions agree coefficient by coefficient.
  4. Step 4:
    QFTN†∣Ψm/N⟩=QFTN†QFTN∣m⟩=∣m⟩\mathrm{QFT}_N^\dagger\lvert\Psi_{m/N}\rangle=\mathrm{QFT}_N^\dagger\mathrm{QFT}_N\lvert m\rangle=\lvert m\rangle
    QFTN\mathrm{QFT}_N is unitary, QFTN†QFTN=I\mathrm{QFT}_N^\dagger\mathrm{QFT}_N=I. The Born rule then gives P(m)=∣⟨m∣m⟩∣2=1P(m)=\lvert\langle m\vert m\rangle\rvert^2=1.

∎

Derivation

Worked example: φ=1/8\varphi=1/8, t=3t=3 (U=TU=T)

  1. Step 1:
    φ=18=123 ⇒ m=1=0012\varphi=\tfrac18=\tfrac{1}{2^3}\ \Rightarrow\ m=1=001_2
    The phase is dyadic with three bits.
  2. Step 2:
    ∣Ψ1/8⟩=18(∣0⟩+e2πi⋅48∣1⟩)(∣0⟩+e2πi⋅28∣1⟩)(∣0⟩+e2πi⋅18∣1⟩)=18(∣0⟩−∣1⟩)(∣0⟩+i∣1⟩)(∣0⟩+eiπ/4∣1⟩)\lvert\Psi_{1/8}\rangle=\tfrac1{\sqrt8}\big(\lvert0\rangle+e^{2\pi i\cdot\frac48}\lvert1\rangle\big)\big(\lvert0\rangle+e^{2\pi i\cdot\frac28}\lvert1\rangle\big)\big(\lvert0\rangle+e^{2\pi i\cdot\frac18}\lvert1\rangle\big)=\tfrac1{\sqrt8}\big(\lvert0\rangle-\lvert1\rangle\big)\big(\lvert0\rangle+i\lvert1\rangle\big)\big(\lvert0\rangle+e^{i\pi/4}\lvert1\rangle\big)
    Qubit 00 controls U4=T4=ZU^4=T^4=Z, qubit 11 controls U2=T2=SU^2=T^2=S, qubit 22 controls U=TU=T; eiπ=−1e^{i\pi}=-1, eiπ/2=ie^{i\pi/2}=i.
  3. Step 3:
    ∣Ψ1/8⟩=18∑k=07e2πik/8∣k⟩=QFT8∣001⟩\lvert\Psi_{1/8}\rangle=\tfrac1{\sqrt8}\sum_{k=0}^{7}e^{2\pi ik/8}\lvert k\rangle=\mathrm{QFT}_8\lvert001\rangle
    The counting-state formula with φ=1/8\varphi=1/8, which is column 11 of the QFT matrix.
  4. Step 4:
    QFT8†∣Ψ1/8⟩=∣001⟩,P(001)=1\mathrm{QFT}_8^\dagger\lvert\Psi_{1/8}\rangle=\lvert001\rangle,\qquad P(001)=1
    The dyadic theorem. The inverse QFT on three qubits is SWAP(0,2)\mathrm{SWAP}(0,2), H(2)H(2), CP(2,1,−π2)\mathrm{CP}(2,1,-\tfrac\pi2), H(1)H(1), CP(2,0,−π4)\mathrm{CP}(2,0,-\tfrac\pi4), CP(1,0,−π2)\mathrm{CP}(1,0,-\tfrac\pi2), H(0)H(0).

Quick check

With t=3t=3 counting qubits and U=P(2π⋅58)U=P(2\pi\cdot\tfrac58), which integer mm (in decimal) does QPE return with probability 11?

tolerance ±1e-9

Non-dyadic phases: the 4/π² bound#

Derivation

Outcome amplitude as a geometric sum

  1. Step 1:
    ⟨ℓ∣QFTN†=(QFTN∣ℓ⟩)†=1N∑k=0N−1e−2πiℓk/N⟨k∣\langle\ell\rvert\mathrm{QFT}_N^\dagger=\Big(\mathrm{QFT}_N\lvert\ell\rangle\Big)^\dagger=\frac1{\sqrt N}\sum_{k=0}^{N-1}e^{-2\pi i\ell k/N}\langle k\rvert
    Take the adjoint of QFTN∣ℓ⟩\mathrm{QFT}_N\lvert\ell\rangle; complex conjugation flips the exponent.
  2. Step 2:
    αℓ:=⟨ℓ∣QFTN†∣Ψφ⟩=1N∑k=0N−1e−2πiℓk/Ne2πiφk=1N∑k=0N−1e2πi(φ−ℓ/N)k\alpha_\ell:=\langle\ell\rvert\mathrm{QFT}_N^\dagger\lvert\Psi_\varphi\rangle=\frac1N\sum_{k=0}^{N-1}e^{-2\pi i\ell k/N}e^{2\pi i\varphi k}=\frac1N\sum_{k=0}^{N-1}e^{2\pi i(\varphi-\ell/N)k}
    Insert ∣Ψφ⟩=1N∑ke2πiφk∣k⟩\lvert\Psi_\varphi\rangle=\tfrac1{\sqrt N}\sum_ke^{2\pi i\varphi k}\lvert k\rangle and use ⟨k∣k′⟩=δkk′\langle k\vert k'\rangle=\delta_{kk'}.
  3. Step 3:
    z:=e2πiδℓ, δℓ:=φ−ℓN∉Z ⇒ αℓ=1N∑k=0N−1zk=1N⋅1−zN1−zz:=e^{2\pi i\delta_\ell},\ \delta_\ell:=\varphi-\frac\ell N\notin\mathbb Z\ \Rightarrow\ \alpha_\ell=\frac1N\sum_{k=0}^{N-1}z^k=\frac1N\cdot\frac{1-z^N}{1-z}
    Geometric series with ratio z≠1z\ne1 (the case z=1z=1 is the dyadic case already treated).
  4. Step 4:
    ∣1−eiθ∣=∣eiθ/2∣ ∣e−iθ/2−eiθ/2∣=2∣sin⁡(θ/2)∣\lvert1-e^{i\theta}\rvert=\lvert e^{i\theta/2}\rvert\,\lvert e^{-i\theta/2}-e^{i\theta/2}\rvert=2\lvert\sin(\theta/2)\rvert
    Factor out eiθ/2e^{i\theta/2}, which has modulus 11, and use eix−e−ix=2isin⁡xe^{i x}-e^{-ix}=2i\sin x.
  5. Step 5:
    P(ℓ)=∣αℓ∣2=1N2⋅4sin⁡2(πNδℓ)4sin⁡2(πδℓ)=∣sin⁡(πNδℓ)Nsin⁡(πδℓ)∣2P(\ell)=\lvert\alpha_\ell\rvert^2=\frac{1}{N^2}\cdot\frac{4\sin^2(\pi N\delta_\ell)}{4\sin^2(\pi\delta_\ell)}=\left\lvert\frac{\sin(\pi N\delta_\ell)}{N\sin(\pi\delta_\ell)}\right\rvert^2
    Apply the previous identity with θ=2πNδℓ\theta=2\pi N\delta_\ell in the numerator and θ=2πδℓ\theta=2\pi\delta_\ell in the denominator.

Result

P(ℓ)=∣sin⁡(π2tδℓ)2tsin⁡(πδℓ)∣2,δℓ=φ−ℓ2t\boxed{P(\ell)=\left\lvert\frac{\sin(\pi2^t\delta_\ell)}{2^t\sin(\pi\delta_\ell)}\right\rvert^2,\qquad\delta_\ell=\varphi-\frac{\ell}{2^t}}
Theorem

The nearest estimate is found with probability at least 4/π²

Let b=round(2tφ)b=\mathrm{round}(2^t\varphi) (taken modulo 2t2^t) and δ=φ−b/2t\delta=\varphi-b/2^t, so ∣δ∣≤2−(t+1)\lvert\delta\rvert\le2^{-(t+1)}. Then P(b)=∣sin⁡(π2tδ)2tsin⁡(πδ)∣2≥4π2≈0.405285P(b)=\left\lvert\frac{\sin(\pi2^t\delta)}{2^t\sin(\pi\delta)}\right\rvert^2\ge\frac4{\pi^2}\approx0.405285.

Proof (6 steps)
  1. Step 1:
    ∣2tφ−b∣≤12 ⇒ ∣δ∣≤12t+1 ⇒ ∣π2tδ∣≤π2\lvert2^t\varphi-b\rvert\le\tfrac12\ \Rightarrow\ \lvert\delta\rvert\le\frac1{2^{t+1}}\ \Rightarrow\ \lvert\pi2^t\delta\rvert\le\frac\pi2
    Rounding to the nearest integer moves by at most 12\tfrac12; divide by 2t2^t, then multiply by π2t\pi2^t. If δ=0\delta=0 the phase is dyadic and P(b)=1P(b)=1, so assume δ≠0\delta\ne0.
  2. Step 2:
    ∣sin⁡x∣≥2π∣x∣for ∣x∣≤π2\lvert\sin x\rvert\ge\frac2\pi\lvert x\rvert\quad\text{for }\lvert x\rvert\le\frac\pi2
    Jordan's inequality: sin⁡\sin is concave on [0,π2][0,\tfrac\pi2], so it lies above the chord from (0,0)(0,0) to (π2,1)(\tfrac\pi2,1), which is 2πx\tfrac2\pi x; ∣sin⁡∣\lvert\sin\rvert is even.
  3. Step 3:
    ∣sin⁡(π2tδ)∣≥2π⋅π2t∣δ∣=2t+1∣δ∣\lvert\sin(\pi2^t\delta)\rvert\ge\frac2\pi\cdot\pi2^t\lvert\delta\rvert=2^{t+1}\lvert\delta\rvert
    Apply Jordan's inequality with x=π2tδx=\pi2^t\delta, admissible by the first step.
  4. Step 4:
    ∣sin⁡(πδ)∣≤π∣δ∣\lvert\sin(\pi\delta)\rvert\le\pi\lvert\delta\rvert
    ∣sin⁡x∣≤∣x∣\lvert\sin x\rvert\le\lvert x\rvert for all real xx (the derivative of x−sin⁡xx-\sin x is 1−cos⁡x≥01-\cos x\ge0).
  5. Step 5:
    P(b)=∣sin⁡(π2tδ)∣2t∣sin⁡(πδ)∣≥2t+1∣δ∣2t⋅π∣δ∣=2π\sqrt{P(b)}=\frac{\lvert\sin(\pi2^t\delta)\rvert}{2^t\lvert\sin(\pi\delta)\rvert}\ge\frac{2^{t+1}\lvert\delta\rvert}{2^t\cdot\pi\lvert\delta\rvert}=\frac2\pi
    Lower-bound the numerator and upper-bound the denominator; ∣δ∣≠0\lvert\delta\rvert\ne0 cancels.
  6. Step 6:
    P(b)≥4π2≈0.405285P(b)\ge\frac4{\pi^2}\approx0.405285
    Square both sides. The bound is approached as t→∞t\to\infty with ∣δ∣=2−(t+1)\lvert\delta\rvert=2^{-(t+1)}, where the phase sits exactly halfway between two grid points.

∎

Derivation

Numerical example: φ=0.3\varphi=0.3, t=3t=3

  1. Step 1:
    2tφ=8⋅0.3=2.4 ⇒ b=round(2.4)=2=01022^t\varphi=8\cdot0.3=2.4\ \Rightarrow\ b=\mathrm{round}(2.4)=2=010_2
    The nearest three-bit grid point is 2/8=0.252/8=0.25.
  2. Step 2:
    δ=0.3−28=0.3−0.25=0.05,π2tδ=0.4π,πδ=0.05π\delta=0.3-\tfrac28=0.3-0.25=0.05,\qquad\pi2^t\delta=0.4\pi,\qquad\pi\delta=0.05\pi
    Definition of δ\delta.
  3. Step 3:
    P(010)=∣sin⁡(0.4π)8sin⁡(0.05π)∣2=(0.9510578⋅0.156434)2=(0.9510571.251476)2≈0.577521P(010)=\left\lvert\frac{\sin(0.4\pi)}{8\sin(0.05\pi)}\right\rvert^2=\left(\frac{0.951057}{8\cdot0.156434}\right)^2=\left(\frac{0.951057}{1.251476}\right)^2\approx0.577521
    Evaluate the geometric-sum formula; sin⁡(0.4π)=sin⁡72∘\sin(0.4\pi)=\sin72^\circ, sin⁡(0.05π)=sin⁡9∘\sin(0.05\pi)=\sin9^\circ.
  4. Step 4:
    0.577521≥4π2≈0.4052850.577521\ge\frac4{\pi^2}\approx0.405285
    Consistent with the theorem. The next most likely outcome is 011011 (38=0.375\tfrac38=0.375, distance 0.0750.075) with probability ≈0.259336\approx0.259336, so the two neighbours of φ\varphi together carry ≈0.836857\approx0.836857.

Quick check

For φ=0.3\varphi=0.3 and t=3t=3, what is the probability that QPE returns the nearest estimate b=010b=010? Give four decimals.

tolerance ±0.0001

Precision versus success probability#

Theorem

Counting qubits for n bits with failure probability ε

Let N=2tN=2^t, let bb be the integer in {0,…,N−1}\{0,\dots,N-1\} with b/N≤φ<(b+1)/Nb/N\le\varphi<(b+1)/N (the best estimate from below), and let e≥2e\ge2 be an integer. Then the probability that the measured outcome mm satisfies ∣m−b∣>e\lvert m-b\rvert>e (distance taken modulo NN) is at most 12(e−1)\frac1{2(e-1)}. Consequently, with t=n+⌈log⁡2(2+12ε)⌉t=n+\big\lceil\log_2\big(2+\frac1{2\varepsilon}\big)\big\rceil and e=2t−n−1e=2^{t-n}-1, QPE returns φ\varphi to nn bits of accuracy, ∣m/2t−φ∣<2−n\lvert m/2^t-\varphi\rvert<2^{-n}, with probability at least 1−ε1-\varepsilon.

Proof (12 steps)
  1. Step 1:
    δ:=φ−bN∈[0,1N),α(b+l) mod N=1N∑k=0N−1e2πi(δ−l/N)k=1N⋅1−e2πi(Nδ−l)1−e2πi(δ−l/N)\delta:=\varphi-\frac bN\in\Big[0,\frac1N\Big),\qquad\alpha_{(b+l)\bmod N}=\frac1N\sum_{k=0}^{N-1}e^{2\pi i(\delta-l/N)k}=\frac1N\cdot\frac{1-e^{2\pi i(N\delta-l)}}{1-e^{2\pi i(\delta-l/N)}}
    The geometric-sum amplitude from the previous section with ℓ=b+l\ell=b+l, −N/2<l≤N/2-N/2<l\le N/2, so δℓ=δ−l/N\delta_\ell=\delta-l/N.
  2. Step 2:
    ∣αb+l∣≤1N⋅2∣1−e2πi(δ−l/N)∣=2N⋅2∣sin⁡(π(δ−l/N))∣\lvert\alpha_{b+l}\rvert\le\frac1N\cdot\frac{2}{\lvert1-e^{2\pi i(\delta-l/N)}\rvert}=\frac{2}{N\cdot2\lvert\sin(\pi(\delta-l/N))\rvert}
    The numerator has modulus at most ∣1∣+∣eiθ∣=2\lvert1\rvert+\lvert e^{i\theta}\rvert=2; the denominator is 2∣sin⁡(θ/2)∣2\lvert\sin(\theta/2)\rvert by the identity ∣1−eiθ∣=2∣sin⁡(θ/2)∣\lvert1-e^{i\theta}\rvert=2\lvert\sin(\theta/2)\rvert.
  3. Step 3:
    −π≤2π(δ−lN)≤π ⇒ ∣sin⁡(π(δ−lN))∣≥2π⋅π∣δ−lN∣=2∣δ−lN∣-\pi\le2\pi\Big(\delta-\frac lN\Big)\le\pi\ \Rightarrow\ \Big\lvert\sin\Big(\pi\Big(\delta-\frac lN\Big)\Big)\Big\rvert\ge\frac2\pi\cdot\pi\Big\lvert\delta-\frac lN\Big\rvert=2\Big\lvert\delta-\frac lN\Big\rvert
    Since 0≤Nδ<10\le N\delta<1 and −N/2<l≤N/2-N/2<l\le N/2, ∣π(δ−l/N)∣≤π/2\lvert\pi(\delta-l/N)\rvert\le\pi/2; apply Jordan's inequality ∣sin⁡x∣≥2π∣x∣\lvert\sin x\rvert\ge\tfrac2\pi\lvert x\rvert.
  4. Step 4:
    ∣αb+l∣≤12N∣δ−l/N∣=12∣Nδ−l∣ ⇒ P(b+l)≤14(l−Nδ)2\lvert\alpha_{b+l}\rvert\le\frac{1}{2N\lvert\delta-l/N\rvert}=\frac{1}{2\lvert N\delta-l\rvert}\ \Rightarrow\ P(b+l)\le\frac1{4(l-N\delta)^2}
    Combine the two bounds and multiply numerator and denominator by NN.
  5. Step 5:
    P(∣m−b∣>e)≤∑l=−N/2+1−(e+1)14(l−Nδ)2+∑l=e+1N/214(l−Nδ)2P(\lvert m-b\rvert>e)\le\sum_{l=-N/2+1}^{-(e+1)}\frac1{4(l-N\delta)^2}+\sum_{l=e+1}^{N/2}\frac1{4(l-N\delta)^2}
    The failure event is the union of the outcomes b+lb+l with l<−el<-e or l>el>e; sum their probabilities.
  6. Step 6:
    ≤14[∑l=−N/2+1−(e+1)1l2+∑l=e+1N/21(l−1)2]\le\frac14\Big[\sum_{l=-N/2+1}^{-(e+1)}\frac1{l^2}+\sum_{l=e+1}^{N/2}\frac1{(l-1)^2}\Big]
    0≤Nδ<10\le N\delta<1. For l<0l<0: ∣l−Nδ∣=∣l∣+Nδ≥∣l∣\lvert l-N\delta\rvert=\lvert l\rvert+N\delta\ge\lvert l\rvert. For l>0l>0: l−Nδ>l−1≥e≥1>0l-N\delta>l-1\ge e\ge1>0.
  7. Step 7:
    ≤14[∑l=e+1∞1l2+∑l=e∞1l2]≤12∑l=e∞1l2\le\frac14\Big[\sum_{l=e+1}^{\infty}\frac1{l^2}+\sum_{l=e}^{\infty}\frac1{l^2}\Big]\le\frac12\sum_{l=e}^{\infty}\frac1{l^2}
    Substitute l→−ll\to-l in the first sum and l−1→ll-1\to l in the second, extend both to infinity (all terms are positive), and bound the first sum by the second, which has one more positive term.
  8. Step 8:
    ∑l=e∞1l2≤1e2+∫e∞dxx2≤∫e−1∞dxx2=1e−1(e≥2)\sum_{l=e}^{\infty}\frac1{l^2}\le\frac1{e^2}+\int_{e}^{\infty}\frac{dx}{x^2}\le\int_{e-1}^{\infty}\frac{dx}{x^2}=\frac1{e-1}\quad(e\ge2)
    1/l2≤∫l−1lx−2 dx1/l^2\le\int_{l-1}^{l}x^{-2}\,dx because x−2x^{-2} is decreasing, so the whole sum is bounded by ∫e−1∞x−2dx=1e−1\int_{e-1}^\infty x^{-2}dx=\frac1{e-1}.
  9. Step 9:
    P(∣m−b∣>e)≤12(e−1)P(\lvert m-b\rvert>e)\le\frac1{2(e-1)}
    Combine the last two steps. This is Nielsen and Chuang equation (5.34).
  10. Step 10:
    e=2t−n−1:∣m−b∣≤e ⇒ ∣m2t−φ∣≤∣m−b∣2t+∣b2t−φ∣<2t−n−12t+12t=2−ne=2^{t-n}-1:\quad\lvert m-b\rvert\le e\ \Rightarrow\ \Big\lvert\frac m{2^t}-\varphi\Big\rvert\le\frac{\lvert m-b\rvert}{2^t}+\Big\lvert\frac b{2^t}-\varphi\Big\rvert<\frac{2^{t-n}-1}{2^t}+\frac1{2^t}=2^{-n}
    Triangle inequality with 0≤φ−b/2t<2−t0\le\varphi-b/2^t<2^{-t}: success within distance ee means nn bits of accuracy.
  11. Step 11:
    12(e−1)=12(2t−n−2)≤ε  ⟺  2t−n≥2+12ε  ⟺  t−n≥log⁡2(2+12ε)\frac1{2(e-1)}=\frac1{2(2^{t-n}-2)}\le\varepsilon\iff2^{t-n}\ge2+\frac1{2\varepsilon}\iff t-n\ge\log_2\Big(2+\frac1{2\varepsilon}\Big)
    Require the failure bound to be at most ε\varepsilon and solve for t−nt-n; 2t−n−2>02^{t-n}-2>0 holds since 2t−n≥2+12ε>22^{t-n}\ge2+\tfrac1{2\varepsilon}>2.
  12. Step 12:
    t=n+⌈log⁡2(2+12ε)⌉ ⇒ P(∣m2t−φ∣<2−n)≥1−εt=n+\Big\lceil\log_2\Big(2+\frac1{2\varepsilon}\Big)\Big\rceil\ \Rightarrow\ P\Big(\Big\lvert\frac m{2^t}-\varphi\Big\rvert<2^{-n}\Big)\ge1-\varepsilon
    The smallest integer tt satisfying the previous condition. For example n=3n=3, ε=0.1\varepsilon=0.1: log⁡2(2+5)=log⁡27≈2.807\log_2(2+5)=\log_2 7\approx2.807, so t=3+3=6t=3+3=6 counting qubits.

∎

Quick check

How many counting qubits tt does the formula t=n+⌈log⁡2(2+12ε)⌉t=n+\lceil\log_2(2+\tfrac1{2\varepsilon})\rceil prescribe for n=4n=4 bits with failure probability ε=0.05\varepsilon=0.05?

tolerance ±1e-9

Superpositions of eigenvectors and the road to HHL#

Derivation

QPE on ∣ψ⟩=∑jcj∣uj⟩\lvert\psi\rangle=\sum_j c_j\lvert u_j\rangle

  1. Step 1:
    ∣ψ⟩=∑jcj∣uj⟩,U∣uj⟩=e2πiφj∣uj⟩,∑j∣cj∣2=1\lvert\psi\rangle=\sum_jc_j\lvert u_j\rangle,\qquad U\lvert u_j\rangle=e^{2\pi i\varphi_j}\lvert u_j\rangle,\qquad\sum_j\lvert c_j\rvert^2=1
    A unitary is normal, so by the spectral theorem its eigenvectors ∣uj⟩\lvert u_j\rangle form an orthonormal basis and every input expands in it.
  2. Step 2:
    QPE(∣0⟩⊗t∣uj⟩)=(∑ℓαℓ(j)∣ℓ⟩)∣uj⟩\mathrm{QPE}\big(\lvert0\rangle^{\otimes t}\lvert u_j\rangle\big)=\Big(\sum_{\ell}\alpha_\ell^{(j)}\lvert\ell\rangle\Big)\lvert u_j\rangle
    QPE on an eigenvector, with the amplitudes αℓ(j)\alpha_\ell^{(j)} of the geometric-sum formula for φ=φj\varphi=\varphi_j.
  3. Step 3:
    QPE(∣0⟩⊗t∣ψ⟩)=∑jcj(∑ℓαℓ(j)∣ℓ⟩)∣uj⟩\mathrm{QPE}\big(\lvert0\rangle^{\otimes t}\lvert\psi\rangle\big)=\sum_jc_j\Big(\sum_\ell\alpha_\ell^{(j)}\lvert\ell\rangle\Big)\lvert u_j\rangle
    Every gate of QPE is linear, so the circuit acts term by term.
  4. Step 4:
    P(ℓ)=∑j∣cj∣2 ∣αℓ(j)∣2P(\ell)=\sum_j\lvert c_j\rvert^2\,\lvert\alpha_\ell^{(j)}\rvert^2
    The vectors ∣uj⟩\lvert u_j\rangle are orthonormal, so ∥∑jcjαℓ(j)∣uj⟩∥2=∑j∣cj∣2∣αℓ(j)∣2\big\lVert\sum_jc_j\alpha_\ell^{(j)}\lvert u_j\rangle\big\rVert^2=\sum_j\lvert c_j\rvert^2\lvert\alpha_\ell^{(j)}\rvert^2: there is no interference between different eigenvectors.
  5. Step 5:
    φj=mj2t ⇒ QPE(∣0⟩⊗t∣ψ⟩)=∑jcj∣mj⟩∣uj⟩,P(mj)=∣cj∣2\varphi_j=\frac{m_j}{2^t}\ \Rightarrow\ \mathrm{QPE}\big(\lvert0\rangle^{\otimes t}\lvert\psi\rangle\big)=\sum_jc_j\lvert m_j\rangle\lvert u_j\rangle,\qquad P(m_j)=\lvert c_j\rvert^2
    With dyadic eigenphases αℓ(j)=δℓmj\alpha_\ell^{(j)}=\delta_{\ell m_j}. Measuring samples the eigenphase φj\varphi_j with probability ∣cj∣2\lvert c_j\rvert^2 and collapses the target onto ∣uj⟩\lvert u_j\rangle.
TypeScriptrunnable
import { runStatevector } from '@/lib/quantum-simulator';

// Phase estimation of U = P(2 pi phi) = diag(1, e^{2 pi i phi}) with eigenvector |1>.
// Counting qubits 0 .. t-1 (qubit 0 = most significant bit), target qubit t.
const phi = 0.3;
const t = 3;
const ops = [];
ops.push({ kind: 'gate', gate: 'X', targets: [t] });
for (let q = 0; q < t; q++) ops.push({ kind: 'gate', gate: 'H', targets: [q] });
// Counting qubit q has weight 2^(t-1-q) and controls U^(2^(t-1-q)) = P(2 pi phi 2^(t-1-q)).
for (let q = 0; q < t; q++) {
  ops.push({ kind: 'gate', gate: 'CP', targets: [q, t], params: [2 * Math.PI * phi * 2 ** (t - 1 - q)] });
}
// Inverse QFT on qubits 0 .. t-1: the QFT circuit reversed, with every angle negated.
for (let q = 0; q < Math.floor(t / 2); q++) ops.push({ kind: 'gate', gate: 'SWAP', targets: [q, t - 1 - q] });
for (let q = t - 1; q >= 0; q--) {
  for (let m = t - q; m >= 2; m--) {
    ops.push({ kind: 'gate', gate: 'CP', targets: [q + m - 1, q], params: [(-2 * Math.PI) / 2 ** m] });
  }
  ops.push({ kind: 'gate', gate: 'H', targets: [q] });
}
const amps = runStatevector({ version: 1, numQubits: t + 1, numClbits: 0, ops }).amplitudes();
// Marginal distribution of the counting register (sum over the target qubit, the least significant bit).
const p = new Array(2 ** t).fill(0);
amps.forEach((a, i) => { p[i >> 1] += a.re * a.re + a.im * a.im; });
p.forEach((v, k) => console.log('P(' + k.toString(2).padStart(t, '0') + ') = ' + v.toFixed(6)));
console.log('round(2^t phi) = ' + Math.round(2 ** t * phi));

Expected output

P(000) = 0.021593
P(001) = 0.051768
P(010) = 0.577521
P(011) = 0.259336
P(100) = 0.040907
P(101) = 0.019440
P(110) = 0.014487
P(111) = 0.014948
round(2^t phi) = 2

References

  1. M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, 10th anniversary ed., Cambridge University Press (2010), §5.2 Phase estimation.doi:10.1017/CBO9780511976667
  2. R. Cleve, A. Ekert, C. Macchiavello and M. Mosca, "Quantum algorithms revisited", Proceedings of the Royal Society A 454, 339–354 (1998).doi:10.1098/rspa.1998.0164
  3. A. Yu. Kitaev, "Quantum measurements and the Abelian stabilizer problem" (1995); arXiv:quant-ph/9511026.arxiv.org

Hands-on labs

  • Three-bit phase estimation of a phase gate

    1 step · auto-graded

    Complete a three-qubit phase-estimation circuit around the controlled powers C-UC\text{-}U, C-U2C\text{-}U^2, C-U4C\text{-}U^4 of a black-box phase gate U=P(2πφ)=diag(1,e2πiφ)U=P(2\pi\varphi)=\mathrm{diag}(1,e^{2\pi i\varphi}), whose eigenvector is ∣1⟩\lvert1\rangle with eigenvalue e2πiφe^{2\pi i\varphi}. For U=TU=T the phase is φ=1/8\varphi=1/8. The grader substitutes φ=1/8\varphi=1/8, 3/83/8, 5/85/8 and 1/21/2 and requires the counting register to read ∣m⟩\lvert m\rangle with φ=m/8\varphi=m/8 with certainty.

    Open lab