QuantumCore

Lesson 28 · 90 min · advanced

The Quantum Fourier Transform

Loading progress…

The quantum Fourier transform QFTN∣j⟩=1N∑k=0N−1e2πijk/N∣k⟩\mathrm{QFT}_N\lvert j\rangle=\tfrac1{\sqrt N}\sum_{k=0}^{N-1}e^{2\pi ijk/N}\lvert k\rangle on N=2nN=2^n amplitudes is the discrete Fourier transform written as a unitary operator. This lesson proves its unitarity from the geometric series, derives the product representation QFT∣j⟩=2−n/2⨂l=1n(∣0⟩+e2πi 0.jn−l+1⋯jn∣1⟩)\mathrm{QFT}\lvert j\rangle=2^{-n/2}\bigotimes_{l=1}^{n}\big(\lvert0\rangle+e^{2\pi i\,0.j_{n-l+1}\cdots j_n}\lvert1\rangle\big) by factoring the exponential, turns that formula into a circuit of n(n+1)2\tfrac{n(n+1)}2 Hadamard and controlled-phase gates plus ⌊n/2⌋\lfloor n/2\rfloor SWAPs, bounds the error of the approximate QFT, inverts the circuit, works the case n=2n=2 in full, and explains why the output amplitudes cannot be read out directly.

Prerequisites: The Deutsch–Jozsa Algorithm, Rotation Gates & ZYZ Euler Decomposition

Objectives and key results

  • Define QFTN\mathrm{QFT}_N on N=2nN=2^n basis states and prove QFT†QFT=I\mathrm{QFT}^\dagger\mathrm{QFT}=I by summing the geometric series ∑k=0N−1e2πi(j−j′)k/N=Nδjj′\sum_{k=0}^{N-1}e^{2\pi i(j-j')k/N}=N\delta_{jj'}.

  • Derive the product representation of QFT∣j⟩\mathrm{QFT}\lvert j\rangle with binary fractions 0.jl⋯jn=∑m=lnjm2l−m−10.j_l\cdots j_n=\sum_{m=l}^{n}j_m2^{l-m-1}, explaining why integer parts of the phase drop out.

  • Build the QFT circuit from HH, CP(2π/2m)\mathrm{CP}(2\pi/2^m) and SWAP, count its n(n+1)2+⌊n/2⌋=O(n2)\tfrac{n(n+1)}2+\lfloor n/2\rfloor=O(n^2) gates against the O(n2n)O(n2^n) operations of the classical FFT, and bound the error of the approximate QFT.

  • Construct the inverse QFT by reversing the gate order and negating the angles, with a proof of (U1⋯Uk)†=Uk†⋯U1†(U_1\cdots U_k)^\dagger=U_k^\dagger\cdots U_1^\dagger.

  • Compute the 4×44\times4 QFT matrix and QFT∣01⟩\mathrm{QFT}\lvert01\rangle both from the definition and from the product formula, and explain why Born-rule sampling limits what the QFT can deliver.

definition · Quantum Fourier transform
QFTN∣j⟩=1N∑k=0N−1e2πijk/N∣k⟩,N=2n\mathrm{QFT}_N\lvert j\rangle=\frac1{\sqrt N}\sum_{k=0}^{N-1}e^{2\pi ijk/N}\lvert k\rangle,\qquad N=2^n
identity · Orthogonality of the characters
∑k=0N−1e2πi(j−j′)k/N=N δjj′(0≤j,j′<N)\sum_{k=0}^{N-1}e^{2\pi i(j-j')k/N}=N\,\delta_{jj'}\qquad(0\le j,j'<N)
theorem · Product representation
QFT∣j1⋯jn⟩=12n/2(∣0⟩+e2πi 0.jn∣1⟩)⊗(∣0⟩+e2πi 0.jn−1jn∣1⟩)⊗⋯⊗(∣0⟩+e2πi 0.j1⋯jn∣1⟩)\mathrm{QFT}\lvert j_1\cdots j_n\rangle=\frac1{2^{n/2}}\big(\lvert0\rangle+e^{2\pi i\,0.j_n}\lvert1\rangle\big)\otimes\big(\lvert0\rangle+e^{2\pi i\,0.j_{n-1}j_n}\lvert1\rangle\big)\otimes\cdots\otimes\big(\lvert0\rangle+e^{2\pi i\,0.j_1\cdots j_n}\lvert1\rangle\big)
algorithm · QFT circuit
for q=0,…,n−1: Hq, CPq+m−1→q(2π/2m) (m=2,…,n−q);then SWAPq,n−1−q (q<⌊n/2⌋)\text{for }q=0,\dots,n-1:\ H_q,\ \mathrm{CP}_{q+m-1\to q}\big(2\pi/2^m\big)\ (m=2,\dots,n-q);\quad\text{then }\mathrm{SWAP}_{q,n-1-q}\ (q<\lfloor n/2\rfloor)
formula · Gate count
n⏟H+n(n−1)2⏟CP+⌊n/2⌋⏟SWAP=n(n+1)2+⌊n/2⌋=O(n2)\underbrace{n}_{H}+\underbrace{\tfrac{n(n-1)}2}_{\mathrm{CP}}+\underbrace{\lfloor n/2\rfloor}_{\mathrm{SWAP}}=\tfrac{n(n+1)}2+\lfloor n/2\rfloor=O(n^2)

Definition and unitarity#

Definition

Quantum Fourier transform

Let N=2nN=2^n. The quantum Fourier transform on nn qubits is the linear operator defined on the computational basis by QFTN∣j⟩=1N∑k=0N−1e2πijk/N∣k⟩,j=0,1,…,N−1,\mathrm{QFT}_N\lvert j\rangle=\frac1{\sqrt N}\sum_{k=0}^{N-1}e^{2\pi ijk/N}\lvert k\rangle,\qquad j=0,1,\dots,N-1, and extended by linearity: QFTN∑jxj∣j⟩=∑kyk∣k⟩\mathrm{QFT}_N\sum_j x_j\lvert j\rangle=\sum_k y_k\lvert k\rangle with yk=1N∑je2πijk/Nxjy_k=\tfrac1{\sqrt N}\sum_j e^{2\pi ijk/N}x_j. Its matrix in the computational basis is Fkj=⟨k∣QFTN∣j⟩=ωjk/NF_{kj}=\langle k\rvert\mathrm{QFT}_N\lvert j\rangle=\omega^{jk}/\sqrt N with ω=e2πi/N\omega=e^{2\pi i/N}, a primitive NN-th root of unity. The matrix is symmetric, Fkj=FjkF_{kj}=F_{jk}. Up to the sign convention of the exponent and the normalisation 1/N1/\sqrt N, the vector map x↦yx\mapsto y is the classical discrete Fourier transform.

Lemma

Orthogonality of the characters

For integers 0≤j,j′≤N−10\le j,j'\le N-1: ∑k=0N−1e2πi(j−j′)k/N=N δjj′\sum_{k=0}^{N-1}e^{2\pi i(j-j')k/N}=N\,\delta_{jj'}.

Proof (6 steps)
  1. Step 1:
    z:=e2πi(j−j′)/N,∑k=0N−1e2πi(j−j′)k/N=∑k=0N−1zk=:Sz:=e^{2\pi i(j-j')/N},\qquad\sum_{k=0}^{N-1}e^{2\pi i(j-j')k/N}=\sum_{k=0}^{N-1}z^k=:S
    Each term is the kk-th power of the single complex number zz, because eak=(ea)ke^{a k}=(e^{a})^k.
  2. Step 2:
    j=j′ ⇒ z=1 ⇒ S=∑k=0N−11=Nj=j'\ \Rightarrow\ z=1\ \Rightarrow\ S=\sum_{k=0}^{N-1}1=N
    If j=j′j=j' every term equals 11 and there are NN terms.
  3. Step 3:
    j≠j′ ⇒ 0<∣j−j′∣<N ⇒ z≠1j\neq j'\ \Rightarrow\ 0<\lvert j-j'\rvert<N\ \Rightarrow\ z\neq1
    z=1z=1 would require (j−j′)/N∈Z(j-j')/N\in\mathbb Z, which is impossible for 0<∣j−j′∣≤N−10<\lvert j-j'\rvert\le N-1.
  4. Step 4:
    S−zS=∑k=0N−1zk−∑k=1Nzk=1−zNS-zS=\sum_{k=0}^{N-1}z^k-\sum_{k=1}^{N}z^k=1-z^N
    Multiplying by zz shifts every exponent by one; all terms z1,…,zN−1z^1,\dots,z^{N-1} cancel and only z0z^0 and −zN-z^N remain.
  5. Step 5:
    S=1−zN1−zS=\frac{1-z^N}{1-z}
    Divide by 1−z1-z, which is non-zero because z≠1z\neq1.
  6. Step 6:
    zN=e2πi(j−j′)=1 ⇒ S=0z^N=e^{2\pi i(j-j')}=1\ \Rightarrow\ S=0
    j−j′j-j' is an integer and e2πim=1e^{2\pi i m}=1 for every integer mm. Together with the case j=j′j=j' this is S=Nδjj′S=N\delta_{jj'}.

∎

Derivation

QFTN†QFTN=I\mathrm{QFT}_N^\dagger\mathrm{QFT}_N=I

  1. Step 1:
    (F†F)j′j=∑k=0N−1Fkj′‾ Fkj\big(F^\dagger F\big)_{j'j}=\sum_{k=0}^{N-1}\overline{F_{kj'}}\,F_{kj}
    Matrix product with (F†)j′k=Fkj′‾(F^\dagger)_{j'k}=\overline{F_{kj'}}.
  2. Step 2:
    =∑k=0N−1e−2πij′k/NN⋅e2πijk/NN=1N∑k=0N−1e2πi(j−j′)k/N=\sum_{k=0}^{N-1}\frac{e^{-2\pi ij'k/N}}{\sqrt N}\cdot\frac{e^{2\pi ijk/N}}{\sqrt N}=\frac1N\sum_{k=0}^{N-1}e^{2\pi i(j-j')k/N}
    Insert Fkj=e2πijk/N/NF_{kj}=e^{2\pi ijk/N}/\sqrt N; the complex conjugate flips the sign of the exponent; combine the exponentials.
  3. Step 3:
    =1N⋅Nδjj′=δjj′=\frac1N\cdot N\delta_{jj'}=\delta_{jj'}
    Orthogonality of the characters (the lemma above).
  4. Step 4:
    F†F=I ⇒ FF†=IF^\dagger F=I\ \Rightarrow\ FF^\dagger=I
    For square matrices a left inverse is also a right inverse, so FF is unitary and QFTN\mathrm{QFT}_N preserves norms and inner products.

Result

QFTN∈U(N),QFTN−1=QFTN†: ∣k⟩↦1N∑je−2πijk/N∣j⟩\boxed{\mathrm{QFT}_N\in U(N),\qquad\mathrm{QFT}_N^{-1}=\mathrm{QFT}_N^\dagger:\ \lvert k\rangle\mapsto\frac1{\sqrt N}\sum_{j}e^{-2\pi ijk/N}\lvert j\rangle}

Quick check

For n=3n=3, what is the probability of measuring ∣000⟩\lvert000\rangle after applying QFT8\mathrm{QFT}_8 to the basis state ∣101⟩=∣5⟩\lvert101\rangle=\lvert5\rangle?

tolerance ±0.000001

The product representation#

Definition

Binary expansion and binary fraction

Write j∈{0,…,2n−1}j\in\{0,\dots,2^n-1\} in binary as j=j1j2⋯jnj=j_1j_2\cdots j_n, i.e. j=∑m=1njm2n−mj=\sum_{m=1}^{n}j_m2^{n-m} with j1j_1 the most significant bit. In this platform qubit 00 holds j1j_1 and qubit n−1n-1 holds jnj_n. The binary fraction is 0.jljl+1⋯jn=∑m=lnjm2l−m−1=jl2+jl+14+⋯+jn2n−l+10.j_lj_{l+1}\cdots j_n=\sum_{m=l}^{n}j_m2^{l-m-1}=\tfrac{j_l}2+\tfrac{j_{l+1}}4+\cdots+\tfrac{j_n}{2^{n-l+1}}. For example 0.1=120.1=\tfrac12, 0.01=140.01=\tfrac14, 0.101=580.101=\tfrac58.

Theorem

Product representation of the QFT

For j=j1⋯jnj=j_1\cdots j_n and N=2nN=2^n, QFT∣j⟩=12n/2⨂l=1n(∣0⟩+e2πi 0.jn−l+1⋯jn∣1⟩)=(∣0⟩+e2πi 0.jn∣1⟩)⊗(∣0⟩+e2πi 0.jn−1jn∣1⟩)⊗⋯⊗(∣0⟩+e2πi 0.j1⋯jn∣1⟩)2n/2,\mathrm{QFT}\lvert j\rangle=\frac1{2^{n/2}}\bigotimes_{l=1}^{n}\big(\lvert0\rangle+e^{2\pi i\,0.j_{n-l+1}\cdots j_n}\lvert1\rangle\big)=\frac{\big(\lvert0\rangle+e^{2\pi i\,0.j_n}\lvert1\rangle\big)\otimes\big(\lvert0\rangle+e^{2\pi i\,0.j_{n-1}j_n}\lvert1\rangle\big)\otimes\cdots\otimes\big(\lvert0\rangle+e^{2\pi i\,0.j_1\cdots j_n}\lvert1\rangle\big)}{2^{n/2}}, where the ll-th factor acts on qubit l−1l-1. In particular QFT∣j⟩\mathrm{QFT}\lvert j\rangle is a product state for every basis input.

Proof (7 steps)
  1. Step 1:
    k=∑l=1nkl2n−l ⇒ kN=k2n=∑l=1nkl2−lk=\sum_{l=1}^{n}k_l2^{n-l}\ \Rightarrow\ \frac kN=\frac{k}{2^n}=\sum_{l=1}^{n}k_l2^{-l}
    Binary expansion of the summation index k=k1⋯knk=k_1\cdots k_n, divided by 2n2^n.
  2. Step 2:
    QFT∣j⟩=12n/2∑k1=01⋯∑kn=01e2πij∑lkl2−l∣k1⋯kn⟩\mathrm{QFT}\lvert j\rangle=\frac1{2^{n/2}}\sum_{k_1=0}^{1}\cdots\sum_{k_n=0}^{1}e^{2\pi ij\sum_{l}k_l2^{-l}}\lvert k_1\cdots k_n\rangle
    The sum over k=0,…,2n−1k=0,\dots,2^n-1 is the same as the sum over all nn bits k1,…,kn∈{0,1}k_1,\dots,k_n\in\{0,1\}.
  3. Step 3:
    =12n/2∑k1=01⋯∑kn=01 ⨂l=1ne2πijkl2−l∣kl⟩=\frac1{2^{n/2}}\sum_{k_1=0}^{1}\cdots\sum_{k_n=0}^{1}\ \bigotimes_{l=1}^{n}e^{2\pi ijk_l2^{-l}}\lvert k_l\rangle
    The exponential of a sum is the product ∏le2πijkl2−l\prod_l e^{2\pi ijk_l2^{-l}}, and ∣k1⋯kn⟩=∣k1⟩⊗⋯⊗∣kn⟩\lvert k_1\cdots k_n\rangle=\lvert k_1\rangle\otimes\cdots\otimes\lvert k_n\rangle; each scalar factor is attached to the tensor factor it depends on (multilinearity of ⊗\otimes).
  4. Step 4:
    =12n/2⨂l=1n(∑kl=01e2πijkl2−l∣kl⟩)=12n/2⨂l=1n(∣0⟩+e2πij2−l∣1⟩)=\frac1{2^{n/2}}\bigotimes_{l=1}^{n}\Big(\sum_{k_l=0}^{1}e^{2\pi ijk_l2^{-l}}\lvert k_l\rangle\Big)=\frac1{2^{n/2}}\bigotimes_{l=1}^{n}\big(\lvert0\rangle+e^{2\pi ij2^{-l}}\lvert1\rangle\big)
    The ll-th factor depends only on klk_l, so the nn nested sums distribute over the tensor product (expand the right-hand side to check). Then evaluate the two terms kl=0k_l=0 (phase 11) and kl=1k_l=1.
  5. Step 5:
    j2−l=∑m=1njm2n−m−l=∑m=1n−ljm2n−l−m⏟∈Z≥0+∑m=n−l+1njm2n−l−m⏟= 0.jn−l+1⋯jnj2^{-l}=\sum_{m=1}^{n}j_m2^{n-m-l}=\underbrace{\sum_{m=1}^{n-l}j_m2^{n-l-m}}_{\in\mathbb Z_{\ge0}}+\underbrace{\sum_{m=n-l+1}^{n}j_m2^{n-l-m}}_{=\,0.j_{n-l+1}\cdots j_n}
    Split the binary expansion of jj at the position of the binary point after dividing by 2l2^l: terms with m≤n−lm\le n-l have a non-negative exponent of 22 (an integer), terms with m≥n−l+1m\ge n-l+1 have exponents −1,−2,…,−l-1,-2,\dots,-l, which is the binary fraction 0.jn−l+1⋯jn0.j_{n-l+1}\cdots j_n by the definition with first index n−l+1n-l+1.
  6. Step 6:
    e2πij2−l=e2πi (integer) e2πi 0.jn−l+1⋯jn=e2πi 0.jn−l+1⋯jne^{2\pi ij2^{-l}}=e^{2\pi i\,(\text{integer})}\,e^{2\pi i\,0.j_{n-l+1}\cdots j_n}=e^{2\pi i\,0.j_{n-l+1}\cdots j_n}
    The integer part drops out because e2πim=1e^{2\pi im}=1 for every integer mm: only the fractional part of j/2lj/2^l matters.
  7. Step 7:
    QFT∣j⟩=12n/2⨂l=1n(∣0⟩+e2πi 0.jn−l+1⋯jn∣1⟩)\mathrm{QFT}\lvert j\rangle=\frac1{2^{n/2}}\bigotimes_{l=1}^{n}\big(\lvert0\rangle+e^{2\pi i\,0.j_{n-l+1}\cdots j_n}\lvert1\rangle\big)
    Substitute into the product. For l=1l=1 (qubit 00) the phase is 0.jn0.j_n; for l=nl=n (qubit n−1n-1) it is 0.j1⋯jn0.j_1\cdots j_n.

∎

Quick check

For n=3n=3 and j=5=1012j=5=101_2, the factor on qubit 22 is 12(∣0⟩+e2πi 0.j1j2j3∣1⟩)\tfrac1{\sqrt2}\big(\lvert0\rangle+e^{2\pi i\,0.j_1j_2j_3}\lvert1\rangle\big). What is 0.j1j2j30.j_1j_2j_3 as a decimal?

tolerance ±1e-9

The circuit and its cost#

Definition

Controlled phase gate

CP(φ)=diag(1,1,1,eiφ)\mathrm{CP}(\varphi)=\mathrm{diag}(1,1,1,e^{i\varphi}) multiplies ∣11⟩\lvert11\rangle by eiφe^{i\varphi} and leaves the other basis states unchanged. It is symmetric in its two qubits, so "control" and "target" are interchangeable. With control cc and target qq, it multiplies the ∣1⟩q\lvert1\rangle_q component by eiφjce^{i\varphi j_c} when the control is in the basis state ∣jc⟩\lvert j_c\rangle. We write Rm:=CP(2π/2m)R_m:=\mathrm{CP}(2\pi/2^m), so R2=CP(π/2)R_2=\mathrm{CP}(\pi/2) and R3=CP(π/4)R_3=\mathrm{CP}(\pi/4).

Derivation

Qubit 00 after HH and the controlled phases R2,…,RnR_2,\dots,R_n

  1. Step 1:
    H∣j1⟩=12(∣0⟩+(−1)j1∣1⟩)=12(∣0⟩+e2πi 0.j1∣1⟩)H\lvert j_1\rangle=\tfrac1{\sqrt2}\big(\lvert0\rangle+(-1)^{j_1}\lvert1\rangle\big)=\tfrac1{\sqrt2}\big(\lvert0\rangle+e^{2\pi i\,0.j_1}\lvert1\rangle\big)
    H∣0⟩=12(∣0⟩+∣1⟩)H\lvert0\rangle=\tfrac1{\sqrt2}(\lvert0\rangle+\lvert1\rangle), H∣1⟩=12(∣0⟩−∣1⟩)H\lvert1\rangle=\tfrac1{\sqrt2}(\lvert0\rangle-\lvert1\rangle), and e2πi 0.j1=e2πij1/2=eiπj1=(−1)j1e^{2\pi i\,0.j_1}=e^{2\pi ij_1/2}=e^{i\pi j_1}=(-1)^{j_1}.
  2. Step 2:
    R2 (control 1) ⁣: 12(∣0⟩+e2πi 0.j1e2πij2/4∣1⟩)=12(∣0⟩+e2πi 0.j1j2∣1⟩)R_2\ (\text{control }1)\!:\ \tfrac1{\sqrt2}\big(\lvert0\rangle+e^{2\pi i\,0.j_1}e^{2\pi ij_2/4}\lvert1\rangle\big)=\tfrac1{\sqrt2}\big(\lvert0\rangle+e^{2\pi i\,0.j_1j_2}\lvert1\rangle\big)
    Qubit 11 is in the basis state ∣j2⟩\lvert j_2\rangle, so CP(2π/4)\mathrm{CP}(2\pi/4) multiplies the ∣1⟩\lvert1\rangle component of qubit 00 by e2πij2/4e^{2\pi ij_2/4}; and j12+j24=0.j1j2\tfrac{j_1}2+\tfrac{j_2}4=0.j_1j_2.
  3. Step 3:
    Rm (control m−1) ⁣: e2πi 0.j1⋯jm−1⋅e2πijm/2m=e2πi 0.j1⋯jmR_m\ (\text{control }m-1)\!:\ e^{2\pi i\,0.j_1\cdots j_{m-1}}\cdot e^{2\pi ij_m/2^m}=e^{2\pi i\,0.j_1\cdots j_m}
    By induction on mm: the control qubit m−1m-1 holds jmj_m, and adding jm/2mj_m/2^m appends the bit jmj_m at position mm of the binary fraction.
  4. Step 4:
    Rn⋯R2H∣j1⟩⊗∣j2⋯jn⟩=12(∣0⟩+e2πi 0.j1j2⋯jn∣1⟩)⊗∣j2⋯jn⟩R_n\cdots R_2H\lvert j_1\rangle\otimes\lvert j_2\cdots j_n\rangle=\tfrac1{\sqrt2}\big(\lvert0\rangle+e^{2\pi i\,0.j_1j_2\cdots j_n}\lvert1\rangle\big)\otimes\lvert j_2\cdots j_n\rangle
    After m=nm=n qubit 00 holds exactly the last factor of the product formula; qubits 1,…,n−11,\dots,n-1 are untouched because a controlled phase on a basis-state control returns the control unchanged.
  5. Step 5:
    qubit q ⁣: H, R2,…,Rn−q ⇒ 12(∣0⟩+e2πi 0.jq+1⋯jn∣1⟩)\text{qubit }q\!:\ H,\ R_2,\dots,R_{n-q}\ \Rightarrow\ \tfrac1{\sqrt2}\big(\lvert0\rangle+e^{2\pi i\,0.j_{q+1}\cdots j_n}\lvert1\rangle\big)
    The same argument on qubit qq (holding jq+1j_{q+1}), with RmR_m controlled by qubit q+m−1q+m-1 (holding jq+mj_{q+m}), m=2,…,n−qm=2,\dots,n-q. Qubits q′>qq'>q are still basis states when qubit qq is processed, which is why the qubits are treated in the order 0,1,…,n−10,1,\dots,n-1.
  6. Step 6:
    12n/2(∣0⟩+e2πi 0.j1⋯jn∣1⟩)⊗(∣0⟩+e2πi 0.j2⋯jn∣1⟩)⊗⋯⊗(∣0⟩+e2πi 0.jn∣1⟩)\frac1{2^{n/2}}\big(\lvert0\rangle+e^{2\pi i\,0.j_1\cdots j_n}\lvert1\rangle\big)\otimes\big(\lvert0\rangle+e^{2\pi i\,0.j_2\cdots j_n}\lvert1\rangle\big)\otimes\cdots\otimes\big(\lvert0\rangle+e^{2\pi i\,0.j_n}\lvert1\rangle\big)
    State after all qubits: these are the factors of the product formula in reverse order.
  7. Step 7:
    ∏q=0⌊n/2⌋−1SWAPq, n−1−q ⇒ QFT∣j⟩\prod_{q=0}^{\lfloor n/2\rfloor-1}\mathrm{SWAP}_{q,\,n-1-q}\ \Rightarrow\ \mathrm{QFT}\lvert j\rangle
    ⌊n/2⌋\lfloor n/2\rfloor SWAPs reverse the qubit order, giving exactly the product representation. Since the circuit and QFT\mathrm{QFT} agree on every basis state and both are linear, they are equal as operators.

Circuit · 3 qubits

Quantum circuit on 3 qubitsX on q0; X on q2; H on q0; P(π/2) on q0 controlled by q1; P(π/4) on q0 controlled by q2; H on q1; P(π/2) on q1 controlled by q2; H on q2; SWAP on q0, q2q0q1q2XXHP(π/2)P(π/4)HP(π/2)H

The three-qubit QFT on ∣101⟩=∣5⟩\lvert101\rangle=\lvert5\rangle: X(0),X(2)X(0),X(2) prepare the input, then H(0); CP(1,0,π/2); CP(2,0,π/4); H(1); CP(2,1,π/2); H(2); SWAP(0,2)H(0);\ \mathrm{CP}(1,0,\pi/2);\ \mathrm{CP}(2,0,\pi/4);\ H(1);\ \mathrm{CP}(2,1,\pi/2);\ H(2);\ \mathrm{SWAP}(0,2). Every amplitude has modulus 1/8≈0.35361/\sqrt8\approx0.3536 and amplitude kk has phase 2π⋅5k/82\pi\cdot5k/8, e.g. ⟨1∣QFT∣5⟩=e5πi/4/8≈−0.25−0.25i\langle1\vert\mathrm{QFT}\vert5\rangle=e^{5\pi i/4}/\sqrt8\approx-0.25-0.25i.

∣ψ⟩=122∣000⟩+(−14−14 i)∣001⟩+122 i∣010⟩+(14−14 i)∣011⟩−122∣100⟩+(14+14 i)∣101⟩−122 i∣110⟩+(−14+14 i)∣111⟩\lvert\psi\rangle = \tfrac{1}{2\sqrt{2}}\lvert 000\rangle + \left(-\tfrac{1}{4} - \tfrac{1}{4}\,i\right)\lvert 001\rangle + \tfrac{1}{2\sqrt{2}}\,i\lvert 010\rangle + \left(\tfrac{1}{4} - \tfrac{1}{4}\,i\right)\lvert 011\rangle - \tfrac{1}{2\sqrt{2}}\lvert 100\rangle + \left(\tfrac{1}{4} + \tfrac{1}{4}\,i\right)\lvert 101\rangle - \tfrac{1}{2\sqrt{2}}\,i\lvert 110\rangle + \left(-\tfrac{1}{4} + \tfrac{1}{4}\,i\right)\lvert 111\rangle
Derivation

Gate count: O(n2)O(n^2) against O(n2n)O(n2^n) for the FFT

  1. Step 1:
    #H=n,#CP=∑q=0n−1(n−q−1)=∑r=0n−1r=n(n−1)2\#H=n,\qquad\#\mathrm{CP}=\sum_{q=0}^{n-1}(n-q-1)=\sum_{r=0}^{n-1}r=\frac{n(n-1)}2
    Qubit qq receives one HH and the rotations R2,…,Rn−qR_2,\dots,R_{n-q}, i.e. n−q−1n-q-1 controlled phases; substitute r=n−q−1r=n-q-1.
  2. Step 2:
    #H+#CP=n+n(n−1)2=n(n+1)2,#SWAP=⌊n2⌋\#H+\#\mathrm{CP}=n+\frac{n(n-1)}2=\frac{n(n+1)}2,\qquad\#\mathrm{SWAP}=\Big\lfloor\frac n2\Big\rfloor
    Add the two counts. For n=3n=3: 66 gates plus one SWAP, the seven gates of lab-qft3.
  3. Step 3:
    n(n+1)2+⌊n2⌋≤n2(n≥1) ⇒ O(n2)=O(log⁡2N)\frac{n(n+1)}2+\Big\lfloor\frac n2\Big\rfloor\le n^2\quad(n\ge1)\ \Rightarrow\ O(n^2)=O(\log^2N)
    n2+n2+n2=n22+n≤n2\tfrac{n^2+n}2+\tfrac n2=\tfrac{n^2}2+n\le n^2 for n≥2n\ge2, and the count is 11 for n=1n=1. Each SWAP is three CNOTs, which changes only the constant.
  4. Step 4:
    TFFT(N)=2TFFT(N/2)+O(N) ⇒ TFFT(N)=O(Nlog⁡N)=O(n2n)T_{\mathrm{FFT}}(N)=2T_{\mathrm{FFT}}(N/2)+O(N)\ \Rightarrow\ T_{\mathrm{FFT}}(N)=O(N\log N)=O(n2^n)
    Cooley–Tukey splits the sum into even and odd indices, yk=Ek+ωkOky_k=E_k+\omega^kO_k, two half-size transforms plus O(N)O(N) combinations; unrolling log⁡2N\log_2N levels of O(N)O(N) work each gives O(Nlog⁡N)O(N\log N).

Result

QFT: O(n2) gates acting on 2n amplitudesvs.FFT: O(n2n) arithmetic operations\boxed{\text{QFT: }O(n^2)\text{ gates acting on }2^n\text{ amplitudes}\qquad\text{vs.}\qquad\text{FFT: }O(n2^n)\text{ arithmetic operations}}
Proposition

Approximate QFT

Let QFT~\widetilde{\mathrm{QFT}} be the QFT circuit with every rotation RmR_m, m>m0m>m_0, removed. Then ∥QFT~−QFT∥≤2πn/2m0\lVert\widetilde{\mathrm{QFT}}-\mathrm{QFT}\rVert\le2\pi n/2^{m_0} in operator norm. With m0=⌈log⁡2(n/ε)⌉m_0=\lceil\log_2(n/\varepsilon)\rceil the error is at most 2πε2\pi\varepsilon, and the circuit uses O(nlog⁡(n/ε))O(n\log(n/\varepsilon)) gates.

Proof (4 steps)
  1. Step 1:
    ∥U1⋯Uk−V1⋯Vk∥≤∑i=1k∥Ui−Vi∥(Ui,Vi unitary)\lVert U_1\cdots U_k-V_1\cdots V_k\rVert\le\sum_{i=1}^{k}\lVert U_i-V_i\rVert\qquad(U_i,V_i\text{ unitary})
    Telescoping: U1⋯Uk−V1⋯Vk=∑iU1⋯Ui−1(Ui−Vi)Vi+1⋯VkU_1\cdots U_k-V_1\cdots V_k=\sum_i U_1\cdots U_{i-1}(U_i-V_i)V_{i+1}\cdots V_k; the triangle inequality and ∥AB∥≤∥A∥∥B∥\lVert AB\rVert\le\lVert A\rVert\lVert B\rVert with ∥U∥=1\lVert U\rVert=1 for unitaries bound each term by ∥Ui−Vi∥\lVert U_i-V_i\rVert.
  2. Step 2:
    ∥Rm−I∥=∥diag(0,0,0,eiθ−1)∥=∣eiθ−1∣=2∣sin⁡(θ/2)∣≤θ=2π2m\lVert R_m-I\rVert=\lVert\mathrm{diag}(0,0,0,e^{i\theta}-1)\rVert=\lvert e^{i\theta}-1\rvert=2\lvert\sin(\theta/2)\rvert\le\theta=\frac{2\pi}{2^m}
    Removing RmR_m replaces it by II. ∣eiθ−1∣2=(cos⁡θ−1)2+sin⁡2θ=2−2cos⁡θ=4sin⁡2(θ/2)\lvert e^{i\theta}-1\rvert^2=(\cos\theta-1)^2+\sin^2\theta=2-2\cos\theta=4\sin^2(\theta/2) and ∣sin⁡x∣≤∣x∣\lvert\sin x\rvert\le\lvert x\rvert.
  3. Step 3:
    ∑m=m0+1∞2π2m=2π2m0\sum_{m=m_0+1}^{\infty}\frac{2\pi}{2^m}=\frac{2\pi}{2^{m_0}}
    Each qubit loses at most the rotations RmR_m with m=m0+1,…,nm=m_0+1,\dots,n; the geometric series bounds their total error.
  4. Step 4:
    ∥QFT~−QFT∥≤n⋅2π2m0≤2πεfor 2m0≥n/ε\lVert\widetilde{\mathrm{QFT}}-\mathrm{QFT}\rVert\le n\cdot\frac{2\pi}{2^{m_0}}\le2\pi\varepsilon\quad\text{for }2^{m_0}\ge n/\varepsilon
    Sum over the nn qubits and choose m0=⌈log⁡2(n/ε)⌉m_0=\lceil\log_2(n/\varepsilon)\rceil. Each qubit keeps at most m0m_0 gates, so the gate count is O(nm0)=O(nlog⁡(n/ε))O(nm_0)=O(n\log(n/\varepsilon)).

∎

Quick check

How many gates (counting HH, CP\mathrm{CP} and SWAP) does the exact QFT circuit use on n=5n=5 qubits?

tolerance ±1e-9

The inverse QFT#

Lemma

Adjoint of a product

For matrices U1,…,UkU_1,\dots,U_k of compatible size, (U1U2⋯Uk)†=Uk†⋯U2†U1†(U_1U_2\cdots U_k)^\dagger=U_k^\dagger\cdots U_2^\dagger U_1^\dagger. If every UiU_i is unitary, the inverse of the circuit U1⋯UkU_1\cdots U_k is the circuit Uk†⋯U1†U_k^\dagger\cdots U_1^\dagger.

Proof (3 steps)
  1. Step 1:
    ((AB)†)ij=(AB)ji‾=∑lAjl‾ Bli‾=∑l(B†)il(A†)lj=(B†A†)ij\big((AB)^\dagger\big)_{ij}=\overline{(AB)_{ji}}=\sum_l\overline{A_{jl}}\,\overline{B_{li}}=\sum_l(B^\dagger)_{il}(A^\dagger)_{lj}=(B^\dagger A^\dagger)_{ij}
    Definition of the adjoint, the conjugate of a sum of products is the sum of the conjugated products, and reordering the scalar factors.
  2. Step 2:
    (U1⋯Uk)†=((U1⋯Uk−1)Uk)†=Uk†(U1⋯Uk−1)†=Uk†Uk−1†⋯U1†(U_1\cdots U_k)^\dagger=\big((U_1\cdots U_{k-1})U_k\big)^\dagger=U_k^\dagger(U_1\cdots U_{k-1})^\dagger=U_k^\dagger U_{k-1}^\dagger\cdots U_1^\dagger
    Induction on kk: apply the case k=2k=2 with A=U1⋯Uk−1A=U_1\cdots U_{k-1}, B=UkB=U_k, then the induction hypothesis.
  3. Step 3:
    (U1⋯Uk)(Uk†⋯U1†)=U1⋯Uk−1(UkUk†)Uk−1†⋯U1†=⋯=I(U_1\cdots U_k)(U_k^\dagger\cdots U_1^\dagger)=U_1\cdots U_{k-1}(U_kU_k^\dagger)U_{k-1}^\dagger\cdots U_1^\dagger=\cdots=I
    UiUi†=IU_iU_i^\dagger=I removes the innermost pair repeatedly, so the reversed adjoint circuit is the inverse.

∎

Derivation

Gate-by-gate adjoints: reverse the order and negate the angles

  1. Step 1:
    H†=H,SWAP†=SWAPH^\dagger=H,\qquad\mathrm{SWAP}^\dagger=\mathrm{SWAP}
    Both matrices are real and symmetric.
  2. Step 2:
    CP(φ)†=diag(1,1,1,eiφ)†=diag(1,1,1,e−iφ)=CP(−φ)\mathrm{CP}(\varphi)^\dagger=\mathrm{diag}(1,1,1,e^{i\varphi})^\dagger=\mathrm{diag}(1,1,1,e^{-i\varphi})=\mathrm{CP}(-\varphi)
    The adjoint of a diagonal matrix conjugates the diagonal.
  3. Step 3:
    QFT8†=H0 CP1,0(−π2) CP2,0(−π4) H1 CP2,1(−π2) H2 SWAP0,2\mathrm{QFT}_8^\dagger=H_0\,\mathrm{CP}_{1,0}(-\tfrac\pi2)\,\mathrm{CP}_{2,0}(-\tfrac\pi4)\,H_1\,\mathrm{CP}_{2,1}(-\tfrac\pi2)\,H_2\,\mathrm{SWAP}_{0,2}
    The circuit SWAP(0,2)→H(0)→⋯→H(2)\mathrm{SWAP}(0,2)\to H(0)\to\cdots\to H(2) read as the operator product SWAP02H2 CP2,1(π2)H1 CP2,0(π4) CP1,0(π2)H0\mathrm{SWAP}_{02}H_2\,\mathrm{CP}_{2,1}(\tfrac\pi2)H_1\,\mathrm{CP}_{2,0}(\tfrac\pi4)\,\mathrm{CP}_{1,0}(\tfrac\pi2)H_0; by the lemma its adjoint is the reversed product of the adjoints, written here as an operator product (first gate applied = rightmost factor SWAP0,2\mathrm{SWAP}_{0,2}).
  4. Step 4:
    (QFT†)jk=Fkj‾=e−2πijk/NN\big(\mathrm{QFT}^\dagger\big)_{jk}=\overline{F_{kj}}=\frac{e^{-2\pi ijk/N}}{\sqrt N}
    Entrywise check: since FF is symmetric, F†=FˉF^\dagger=\bar F, the same matrix with ii replaced by −i-i, consistent with negating every phase angle.

Result

circuit order: SWAP(0,2); H(2); CP(2,1,−π2); H(1); CP(2,0,−π4); CP(1,0,−π2); H(0)\boxed{\text{circuit order: }\mathrm{SWAP}(0,2);\ H(2);\ \mathrm{CP}(2,1,-\tfrac\pi2);\ H(1);\ \mathrm{CP}(2,0,-\tfrac\pi4);\ \mathrm{CP}(1,0,-\tfrac\pi2);\ H(0)}

Quick check

Which circuit implements QFT4†\mathrm{QFT}_4^\dagger on two qubits (gates listed in the order they are applied)?

Select one answer

Worked example: two qubits#

For n=2n=2, N=4N=4 and ω=e2πi/4=i\omega=e^{2\pi i/4}=i, so Fjk=ijk/2F_{jk}=i^{jk}/2. The exponent jkjk only matters modulo 44: i0=1i^0=1, i1=ii^1=i, i2=−1i^2=-1, i3=−ii^3=-i. Row jj lists i0,ij,i2j,i3ji^{0},i^{j},i^{2j},i^{3j}: QFT4=12(11111i−1−i1−11−11−i−1i).\mathrm{QFT}_4=\frac12\begin{pmatrix}1&1&1&1\\1&i&-1&-i\\1&-1&1&-1\\1&-i&-1&i\end{pmatrix}. Each column has squared norm 4⋅14=14\cdot\tfrac14=1, and distinct columns are orthogonal, e.g. columns 11 and 33: 14(1⋅1+iˉ⋅(−i)+(−1)(−1)+(−i)‾⋅i)=14(1−1+1−1)=0\tfrac14\big(1\cdot1+\bar i\cdot(-i)+(-1)(-1)+\overline{(-i)}\cdot i\big)=\tfrac14(1-1+1-1)=0.

QFT4=12(11111i−1−i1−11−11−i−1i)\mathrm{QFT}_4=\frac12\begin{pmatrix}1&1&1&1\\1&i&-1&-i\\1&-1&1&-1\\1&-i&-1&i\end{pmatrix}
QFT_44×4

Unitary verified: U†U=IU^\dagger U = I (max deviation 0.0e+0)

Derivation

QFT∣01⟩\mathrm{QFT}\lvert01\rangle from the definition

  1. Step 1:
    ∣01⟩=∣j=1⟩ ⇒ QFT∣1⟩=12∑k=03ik∣k⟩\lvert01\rangle=\lvert j=1\rangle\ \Rightarrow\ \mathrm{QFT}\lvert1\rangle=\frac12\sum_{k=0}^{3}i^{k}\lvert k\rangle
    j=0⋅2+1=1j=0\cdot2+1=1 and e2πi⋅1⋅k/4=ike^{2\pi i\cdot1\cdot k/4}=i^k.
  2. Step 2:
    =12(∣0⟩+i∣1⟩−∣2⟩−i∣3⟩)=12(∣00⟩+i∣01⟩−∣10⟩−i∣11⟩)=\frac12\big(\lvert0\rangle+i\lvert1\rangle-\lvert2\rangle-i\lvert3\rangle\big)=\frac12\big(\lvert00\rangle+i\lvert01\rangle-\lvert10\rangle-i\lvert11\rangle\big)
    i0=1i^0=1, i1=ii^1=i, i2=−1i^2=-1, i3=−ii^3=-i; then write kk in binary. This is the second column of the matrix above.
Derivation

QFT∣01⟩\mathrm{QFT}\lvert01\rangle from the product formula

  1. Step 1:
    j1=0, j2=1:0.j2=0.12=12,0.j1j2=0.012=14j_1=0,\ j_2=1:\qquad 0.j_2=0.1_2=\tfrac12,\qquad0.j_1j_2=0.01_2=\tfrac14
    Binary fractions of the bits of j=01j=01.
  2. Step 2:
    QFT∣01⟩=12(∣0⟩+e2πi/2∣1⟩)⊗(∣0⟩+e2πi/4∣1⟩)=12(∣0⟩−∣1⟩)⊗(∣0⟩+i∣1⟩)\mathrm{QFT}\lvert01\rangle=\tfrac12\big(\lvert0\rangle+e^{2\pi i/2}\lvert1\rangle\big)\otimes\big(\lvert0\rangle+e^{2\pi i/4}\lvert1\rangle\big)=\tfrac12\big(\lvert0\rangle-\lvert1\rangle\big)\otimes\big(\lvert0\rangle+i\lvert1\rangle\big)
    Product representation with n=2n=2: qubit 00 carries 0.j20.j_2 and qubit 11 carries 0.j1j20.j_1j_2; eiπ=−1e^{i\pi}=-1 and eiπ/2=ie^{i\pi/2}=i.
  3. Step 3:
    =12(∣00⟩+i∣01⟩−∣10⟩−i∣11⟩)=\tfrac12\big(\lvert00\rangle+i\lvert01\rangle-\lvert10\rangle-i\lvert11\rangle\big)
    Expand the tensor product term by term. Both computations agree.

What the QFT does not give you#

Derivation

Where the QFT helps: periodic states concentrate on multiples of N/rN/r

  1. Step 1:
    ∣ψ⟩=rN∑m=0N/r−1∣x0+mr⟩(r∣N, 0≤x0<r)\lvert\psi\rangle=\sqrt{\frac rN}\sum_{m=0}^{N/r-1}\lvert x_0+mr\rangle\qquad(r\mid N,\ 0\le x_0<r)
    A state with period rr and offset x0x_0, the kind of state produced by phase kickback in the period-finding and phase-estimation algorithms. It has N/rN/r equal amplitudes.
  2. Step 2:
    ⟨k∣QFT∣ψ⟩=rN e2πix0k/N∑m=0N/r−1e2πimrk/N\langle k\vert\mathrm{QFT}\vert\psi\rangle=\frac{\sqrt r}{N}\,e^{2\pi ix_0k/N}\sum_{m=0}^{N/r-1}e^{2\pi imrk/N}
    Apply the definition to each term and pull out the factor independent of mm.
  3. Step 3:
    ∑m=0N/r−1(e2πirk/N)m={N/r,k∈{0,Nr,2Nr,… }0,otherwise\sum_{m=0}^{N/r-1}\big(e^{2\pi irk/N}\big)^m=\begin{cases}N/r,&k\in\{0,\tfrac Nr,\tfrac{2N}r,\dots\}\\0,&\text{otherwise}\end{cases}
    The orthogonality lemma with N/rN/r in place of NN: e2πirk/N=e2πik/(N/r)e^{2\pi irk/N}=e^{2\pi ik/(N/r)} equals 11 exactly when N/rN/r divides kk, and otherwise its (N/r)(N/r)-th power is 11, so the geometric series vanishes.
  4. Step 4:
    P(k=ℓN/r)=(rN⋅Nr)2=1r,ℓ=0,…,r−1P\big(k=\ell N/r\big)=\Big(\frac{\sqrt r}{N}\cdot\frac Nr\Big)^2=\frac1r,\qquad\ell=0,\dots,r-1
    The offset only contributes a phase, so the outcome is a uniformly random multiple of N/rN/r, from which rr is recovered with a few samples. The next lesson, quantum phase estimation, uses QFT†\mathrm{QFT}^\dagger in exactly this way to read an eigenphase.

Quick check

For N=8N=8, apply QFT8\mathrm{QFT}_8 to 12(∣1⟩+∣3⟩+∣5⟩+∣7⟩)\tfrac12(\lvert1\rangle+\lvert3\rangle+\lvert5\rangle+\lvert7\rangle) (period r=2r=2). What is the probability of the outcome k=4k=4?

tolerance ±0.000001
TypeScriptrunnable
import { cexpi, circuitUnitary, matrix } from '@/lib/quantum-simulator';

// Three-qubit QFT circuit, qubit 0 = most significant bit:
// on qubit q apply H, then CP(2 pi / 2^m) controlled by qubit q + m - 1 (m = 2 .. n - q); finally SWAP(0, n - 1).
const n = 3;
const N = 2 ** n;
const ops = [];
for (let q = 0; q < n; q++) {
  ops.push({ kind: 'gate', gate: 'H', targets: [q] });
  for (let m = 2; m <= n - q; m++) {
    ops.push({ kind: 'gate', gate: 'CP', targets: [q + m - 1, q], params: [(2 * Math.PI) / 2 ** m] });
  }
}
ops.push({ kind: 'gate', gate: 'SWAP', targets: [0, n - 1] });
const U = circuitUnitary({ version: 1, numQubits: n, numClbits: 0, ops });

// Definition: F[j][k] = e^{2 pi i j k / N} / sqrt(N).
const rows = [];
for (let j = 0; j < N; j++) {
  const row = [];
  for (let k = 0; k < N; k++) {
    const w = cexpi((2 * Math.PI * j * k) / N);
    row.push([w.re / Math.sqrt(N), w.im / Math.sqrt(N)]);
  }
  rows.push(row);
}
const F = matrix.fromRows(rows);
console.log('gates: ' + ops.length);
console.log('max|U - F| < 1e-12: ' + (matrix.maxAbsDiff(U, F) < 1e-12));
console.log('U is unitary: ' + matrix.isUnitary(U, 1e-12));
const e = matrix.get(U, 1, 1);
console.log('U[1][1] = ' + e.re.toFixed(4) + ' + ' + e.im.toFixed(4) + 'i');

Expected output

gates: 7
max|U - F| < 1e-12: true
U is unitary: true
U[1][1] = 0.2500 + 0.2500i

References

  1. D. Coppersmith, "An approximate Fourier transform useful in quantum factoring", IBM Research Report RC 19642 (1994); arXiv:quant-ph/0201067.arxiv.org
  2. M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, 10th anniversary ed., Cambridge University Press (2010), §5.1 The quantum Fourier transform.doi:10.1017/CBO9780511976667
  3. J. W. Cooley and J. W. Tukey, "An algorithm for the machine calculation of complex Fourier series", Mathematics of Computation 19, 297–301 (1965).doi:10.1090/S0025-5718-1965-0178586-1

Hands-on labs

  • Build the three-qubit QFT

    1 step · auto-graded

    Build the quantum Fourier transform on three qubits, F∣j⟩=18∑k=07e2πijk/8∣k⟩F\lvert j\rangle=\tfrac1{\sqrt8}\sum_{k=0}^{7}e^{2\pi ijk/8}\lvert k\rangle, from Hadamard gates, controlled-phase gates CP(φ)=diag(1,1,1,eiφ)\mathrm{CP}(\varphi)=\mathrm{diag}(1,1,1,e^{i\varphi}) and one SWAP, using at most seven gates.

    Open lab