CS 294 Fourier Transform, Period Finding and Factoring in BQP

2/11/09

Spring 2009 Lecture 4

1 Recap: Simon’s Algorithm

Recall that in the Simon’s problem, we are given a function f : Z_{2}^{n} → Z_{2}^{n} (i.e. from n-bit strings to n-bit strings), with the promise that there is a non-zero string s ∈ Z_{2}^{n}\ {0} such that

for all x ≠ y, f(x) = f(y) if and only if x ⊕ y = s .

The challenge is to determine s. As we saw last time, the problem can be solved with the following circuit.

[thin,scale=1.5] [] (0.25,-0.25) – ++(5.5,0); [] (0.25,-0.5) – ++(5.5,0); [] (0.25,-0.75) – ++(5.5,0); [] (0.25,-1.25) – ++(5.5,0); [] (0.25,-1.5) – ++(5.5,0); [] (5.75,0) – ++(0,-1) arc (-90:90:18pt and 14pt) ++(0.1,-0.25) arc (90:0:10pt) ++(-0.4,0) [-¿] – ++(0.4,0.4); [left] at (0,-0.5) |0⟩; [left] at (0,-1.375) |0⟩;

[right] at (6,-1.375) |f(x)⟩; [rectangle,draw,fill=white,minimum height=40pt,minimum width=40pt,below] at (1.5,0) H; [rectangle,draw,fill=white,minimum height=70pt,minimum width=40pt,below] at (3,0) U_{f};

[rectangle,draw,fill=white,minimum height=40pt,minimum width=40pt,below] at (4.5,0) H;

Figure 1: Circuit for Simon’s Algorithm

The above circuit corresponds to the following sequence of transformations.

|0⟩|0⟩ −→^{H^{⊗n}} \frac{1}{\sqrt{2^{n}}} ∑_{x} |x⟩|0⟩

−→^{f} \frac{1}{\sqrt{2^{n}}} ∑_{x} |x⟩|f(x)⟩

−→^{measure} \frac{1}{\sqrt{2}}(|x_{0}⟩ + |x_{0}⊕ s⟩) ⊗ |a⟩

(measuring the 2nd register, we observe a ∈ Z_{2}^{n} such that a = f(x_{0}) = f(x_{0}⊕ s))

−→^{H^{⊗n}} \frac{1}{\sqrt{2^{n}}} ∑_{y} α_{y}|y⟩|a⟩

for some numbers α_{y}.

Recall that the Hadamard transform of a general state |x⟩ is

H^{⊗n} |x⟩ = \frac{1}{\sqrt{2^{n}}} ∑_{y} (−1)^{x·y}|y⟩,

so

α_{y} = \frac{1}{\sqrt{2}}((−1)^{x_{0}·y}+ (−1)^{(x_{0}⊕s)·y}).

There are now two cases. For each y, if s · y = 1, then α_{y} = 0, whereas if s · y = 0, then α_{y} = (−1)^{x_{0}·y}\sqrt{2}.

When we observe the first register, we get a uniformly random y such that s · y = s_{1}y_{1}+ ··· + s_{n}y_{n} = 0. We repeat to collect more and more equations, and recover s from n linearly independent equations.

There is another way to view the final Hadamard transform. When we have |x_{0}⟩ + |x_{0} ⊕ s⟩, measuring immediately would destroy the state. That’s why we transform it to another basis (the Hadamard basis) before measuring.

2 Fourier Transform on Z_{M} for Integer M

Let f : Z_{M} → C be a complex-valued function on Z_{M}. Its Fourier transform f̂: Z_{M} → C is given by

f̂(t) = \frac{1}{\sqrt{M}} ∑_{x∈Z_{M}}f(x)ω^{xt}

where ω = exp(2πi/M) is a primitive M-th root of unity. If we write f as the vector

f⃗ = ( f(0) ; f(1) ; ⋮ ; f(M − 1) ) ∈ C^{M},

and similarly write f̂as f̂ ∈ C^{M}, then the vectors f⃗ and f̂ are related by a change of basis f̂ = F_{M}f⃗, where the matrix F_{M} takes the form

F_{M} = \frac{1}{\sqrt{M}} ( \begin{matrix}1 & 1 & 1 & 1 & ··· & 1 \\ 1 & ω & ω^{2} & ω^{3} & ··· & ω^{M−1} \\ 1 & ω^{2} & ω^{4} & ω^{6} & ··· & ω^{2M−2} \\ 1 & ω^{3} & ω^{6} & ω^{9} & ··· & ω^{3M−3} \\ ⋮ & ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 1 & ω^{M−1} & ω^{2M−2} & ω^{3M−3} & ··· & ω^{(M−1)(M−1)} \end{matrix} ) ,

that is, (i,j)-th entry of F_{M} is ω^{ij} (if we ignore the normalization factor 1/\sqrt{M}).

3 Classical Fast Fourier Transform

Straightforward multiplication of the vector f⃗ by F_{M} would take Ω(M^{2}) steps because multiplication of f⃗ by each row requires M multiplications. Exploiting the symmetry of F_{M}, it is possible to perform Fourier transform in O(M logM) operations when M is a power of two, i.e. M = 2^{m}. This algorithm is known as fast Fourier transform (FFT).

The idea is to rewrite the Fourier coefficients f̂(j) as

f̂(j) = ∑_{i=0}^{M−1} ω^{ij}f(i)  (where for simplicity we ignore the normalization factor 1/\sqrt{M})

= ∑_{i even}ω^{ij}f(i) + ∑_{i odd}ω^{ij}f(i) (splitting into odd and even terms)

= ∑_{i′=0}^{M/2−1}(ω^{2i′})^{j}f(2i′) + ω^{j} ∑_{i′=0}^{M/2−1}(ω^{2i′})^{j}f(2i′+ 1) (write even i as 2i′, odd as 2i′+ 1)

= ( F_{M/2}\overrightarrow{f_{even}})(j) + ω^{j} ( F_{M/2}\overrightarrow{f_{odd}})(j).

[figure: identity; _{-----------------}; v_{0}; { ; |; F T_{N/2}; ×1; +; v_{1}; { ; |; F T_{N/2}; ×w^{j}; −; multiplication by w^{j}; _{----------------}]

Figure 2: A circuit for classical fast Fourier transform

The above idea is summarized in the diagram below.

This representation gives a recursive algorithm for computing the Fourier transform in time T (M) = 2T (M/2)+ O(M) = O(M logM).

4 Quantum Fourier Transform

We continue to assume M = 2^{m}. Suppose a quantum state |f⟩ on m qubits is given as |f⟩ = ∑ _{x=0}^{M−1} f(x)|x⟩. Quantum Fourier transform (QFT) is the operation that maps |f⟩ to |f̂⟩, where |f̂⟩ = ∑ _{x=0}^{M−1} f̂(x)|x⟩ (and f̂(x) are the Fourier coefficients of f).

As we shall see, QFT can be implemented by circuit of size O(log^{2}M). However, this does not constitute an exponential speed-up over the classical algorithm because the result of quantum Fourier transform is a superposition of states which can be observed, and any measurement can extract at most m = logM bits of information.

We now describe a circuit that implements quantum Fourier transform.

Step 1: QF T_{M/2} on the first m − 1 qubits

Similar to the classical fast Fourier transform, we will split ∑_{x}f(x)|x⟩ into odd and even terms. Hence

∑_{x} f(x)|x⟩ = ∑_{i=0}^{M/2−1} f(2i)|i⟩|0⟩ + ∑_{i=0}^{M/2−1} f(2i + 1)|i⟩|1⟩.

where |2i⟩ is written as |i⟩|0⟩ because appending a zero to a binary number is the same as doubling the number (e.g. if i is 101100 in binary, then 2i is 1011000). In a quantum circuit for uantum Fourier transform, we will first apply QF T_{M/2} on the first register (i.e. the first m − 1 qubits), obtaining

∑_{i=0}^{M/2−1} α_{i}|i⟩|0⟩ + ∑_{i=0}^{M/2−1} β_{i}|i⟩|1⟩

[figure: 1st bit; R_{1}; most significant bits; { ; |; 2nd bit; QF T_{M/2}; R_{2}; least significant bit; •; •; H]

Figure 3: Circuit for quantum Fourier transform

for certain amplitudes α_{i} and β_{i}.

Step 2: Controlled phase shifts

Next, for each of the first m-1 qubits k (1 ≤ k ≤ m), if both the k-th qubit and the last qubit are 1, then we need to multiply the phase by ω^{2^{m−k}}, and otherwise leave the phase unchanged. Thus, we apply the following transformations R_{k}:

R_{k}| \overbrace{1}^{k-th qubit} \overbrace{1}^{last qubit}⟩ = ω^{2^{m−k}} |11⟩
R_{k}|01⟩  = |01⟩
R_{k}|10⟩  = |10⟩
R_{k}|00⟩  = |00⟩

Hence, R_{k} is just a controlled phase shift (with angle 2π/2^{k}). After the controlled phase shift, we get the state

∑_{i=0}^{M/2−1} α_{i}|i⟩|0⟩ + ∑_{i=0}^{M/2−1} ω^{i}β_{i}|i⟩|1⟩.

Step 3: Hadamard gate

Finally, we apply a Hadamard gate to the last qubit, and end up with the state

\frac{1}{\sqrt{2}} [ ∑_{i=0}^{M/2−1} α_{i}|i⟩(|0⟩ + |1⟩) + ∑_{i=0}^{M/2−1} ω^{i}β_{i}|i⟩(|0⟩ − |1⟩)]

= \frac{1}{\sqrt{2}} [ ∑_{i} (α_{i}+ ω^{i}β_{i}) \underbrace{|i⟩|0⟩}_{ith output}+∑_{i} (α_{i}− ω^{i}β_{i}) \underbrace{|i⟩|1⟩}_{(i + M/2)th output} ] .

Putting together, in the circuity above the quantum Fourier transform on m − 1 qubits corresponds to two Fourier transforms on m − 1 bits in the figure ??. The controlled phase shifts correspond to multiplications by ω^{j} in classical circuit. Finally, the Hadamard gate at the very end corresponds to the summation.

The number of gates T (M) satisfies the recurrence relation T (M) = T (M/2) + logM. Thus T (M) = O(log^{2}M).

5 Period Finding

Period finding is the problem in which we are given a function f : Z_{M} → C, with the promise that f is periodic with period r, i.e.

there is a r such that for all x ≠ y, f(x) = f(y) if and only if x ≡ y mod r.

The challenge is the find the period r.

This problem can be solved efficiently using the following circuit.

[thin,scale=1.5] [] (0.25,-0.25) – ++(5.5,0); [] (0.25,-0.5) – ++(5.5,0); [] (0.25,-0.75) – ++(5.5,0); [] (0.25,-1.25) – ++(5.5,0); [] (0.25,-1.5) – ++(5.5,0); [] (5.75,0) – ++(0,-1) arc (-90:90:18pt and 14pt) ++(0.1,-0.25) arc (90:0:10pt) ++(-0.4,0) [-¿] – ++(0.4,0.4); [left] at (0,-0.5) |0⟩; [left] at (0,-1.375) |0⟩;

[right] at (6,-1.375) |f(x)⟩; [rectangle,draw,fill=white,minimum height=40pt,minimum width=40pt,below] at (1.5,0) QF T_{M}; [rectangle,draw,fill=white,minimum height=70pt,minimum width=40pt,below] at (3,0) U_{f}; [rectangle,draw,fill=white,minimum height=40pt,minimum width=40pt,below] at (4.5,0) QF T_{M};

Figure 4: Circuit for period finding

The above circuit corresponds to the following sequence of transformations.

|0⟩|0⟩ −→^{QF T_{M}}\frac{1}{\sqrt{M}} ∑_{x∈Z_{M}} |x⟩|0⟩

−→^{f} \frac{1}{\sqrt{M}} ∑_{x} |x⟩|f(x)⟩

−→^{measure 2nd register} \frac{1}{\sqrt{M/r}} ∑_{k=0}^{M/r−1} |ℓ + kr⟩|f(ℓ)⟩

(Here we assume r divides M to simplify the analysis. We will remove this restriction later.)

−→^{QF T_{M}}\sqrt{\frac{r}{M}} \frac{1}{\sqrt{M}} ∑_{y} α_{y}|y⟩,

where α_{y} = ∑ _{k=0}^{M/r−1} ω^{(ℓ+kn)y} = ω^{ℓy}∑_{k}ω^{kry}.

There are two cases for y:

1. Case 1: y is a multiple of \frac{M}{r}.

In this case, then ω^{kry}= e^{2πiry/M} = 1. So α_{y} = \frac{\sqrt{r}}{M} \frac{M}{r} = \frac{1}{\sqrt{r}}.

Note that there are r multiples of M/r. The sum of the magnitudes squared for these values of y is 1. This implies that for any other y, α_{y} = 0.

2. Case 2: y is not a multiple of \frac{M}{r}.

We already showed that α_{y} must be 0 from the previous case. But we can also give an intuition for why this is the case. Note that ω^{ry},ω^{2ry},... are evenly spaced vectors of unit length around the origin. Being the sum of these complex numbers, α_{y} is 0.

In other words, if we measure the output from the second quantum Fourier transform, we get a uniformly random multiple of M/r. If we repeat the whole process t times, getting t random multiples y_{1},...,y_{t} of M/r, the greatest common divisor of the y_{j}’s is likely to be M/r. Since we know M, we can recover r from M/gcd(y_{1},...,y_{t}).

Let us compute the chance of finding the correct period with t samples. Suppose after repeating t times, we have not found the desired period M/r, but instead a multiple, say λM/r. This means that each of the t samples must be a multiple of λM/r. There are M/(λM/r) = r/λ multiples of λM/r, and since there are r multiples in total, the probability of getting a multiple of λM/r is 1/λ. Therefore,

Pr[gcdis a multiple of λM/r] = ( \frac{1}{λ}) ^{t} ≤ ( \frac{1}{2}) ^{t},

and we err with probability

Pr[gcd> M/r after t samples] ≤ M ( \frac{1}{2}) ^{t}.

So t = O(logM) measurements suffice to guarantee a solution. A more careful analysis shows that a constant number of samples is sufficient.

6 Period Finding: The General Case

For the general case where M is not a multiple of r, we will fix M = 2^{m} to be a power of two that is at least r^{2}.

The change to the above analysis is that, after measuring the 2nd register, we get

−→^{measure 2nd register} \frac{1}{\sqrt{s}} ∑_{k=0}^{s−1} |ℓ + kr⟩|f(ℓ)⟩

where s = ⌊M/r⌋ or s = ⌊M/r⌋ + 1. If we now take QFT on the first register, we get ∑_{y}α_{y}|y⟩ with

α_{y} = \frac{1}{\sqrt{sM}}ω^{ℓy}∑_{k=0}^{s−1} ω^{kry}.

Case 1 now becomes:

1. Case 1: |ry mod M| ≤ \frac{r}{2} (in this section the remainder mod M is allowed to be negative):

Intuitively, in this case the amplitudes ω^{kry}“almost line up” in the complex plane. Previously, when the period r divides M exactly, all the amplitudes for multiples of M/r “line up” at 1.

Claim: If |ry mod M| ≤ \frac{r}{2}, then |α_{y}| ≥ Cs for some constant C.

This claim implies that we have substantial probability of observing a y that falls into case 1. How many y belong to this case? When r is coprime to M (which is the usual case when we run period finding as a subroutine in factoring), the set {ry | y ∈ Z_{M}} is just Z_{M}. Put differently, as y runs through Z_{M}, the product ry also runs through Z_{M}. Hence there are about 2\frac{r}{2} = r such y. Then

Pr[Observing such a y] ≥ r · C^{2}· s^{2} \frac{1}{Ms} ≥ C^{2}· s · \frac{r}{M} ≥ const.

For the rest of the discussion, assume that we measure an y satisfying |ry mod M| ≤ r/2. How does this help us compute the period r? By assumption |ry − cM| ≤ r/2 for some integer c, and hence

| \frac{y}{M} − \frac{c}{r}| ≤ \frac{1}{2M}.

Here y and M are both known. We shall show, assuming that M > r^{2}, how to recover c/r by continued fraction.

Here is the idea: c/r is a close approximation to y/M. Is it possible to get a better rational approximation with denominator at most r? We will show it is impossible. Suppose c′/r′is a better rational approximation with denominator r′≤ r. Then

| \frac{c}{r} − \frac{c′}{r′}| = | \frac{cr′− c′r}{rr′} | ≥ \frac{1}{r^{2}}

But now it follows that

| \frac{c′}{r′}− \frac{y}{M}| ≥ \frac{1}{r^{2}} − \frac{1}{2M} > \frac{1}{2M}.

So if we compute the continued fraction expansion of y/M and look at the successive approximations to y/M, one of these must be c/r, thus yielding r. (See the section on continued fractions below.)

7 Continued Fraction

Definition 4.1 (Continued Fraction): A real number α can be approximated by an iterated fraction

CF_{n}(α) = a_{0}+ \frac{1}{a_{1}+ \frac{1}{a_{2}+ \frac{1}{...+ \frac{1}{a_{n}}}}} = \frac{P_{n}}{Q_{n}},

where a_{0},...,a_{n} (and hence P_{n} and Q_{n}) are integers.

CS 294, Spring 2009, Lecture 4 7

2 Example 1: Let us approximate π to two decimal places with a rational number. We know that

π = 3.14...
= 3 + \frac{14}{100}
= 3 + \frac{1}{\frac{100}{14}}
= 3 + \frac{1}{7 + \frac{2}{14}}
≈ 3 + \frac{1}{7}
= \frac{22}{7}

Example 2: If we decide to approximate π to four decimal places, we would have

π = 3.1415...
= 3 + \frac{1415}{10000}
= 3 + \frac{1}{\frac{10000}{1415}}
= 3 + \frac{1}{7 + \frac{95}{1415}}
= 3 + \frac{1}{7 + \frac{1}{\frac{1415}{95}}}

= 3 + \frac{1}{7 + \frac{1}{14 + \frac{85}{95}}}

≈ 3 + \frac{1}{7 + \frac{1}{14}}
= \frac{311}{99}

The following lemmas are well known facts about continued fraction that we state without proof. Lemma 4.1: CF_{n}(α) is the best rational approximation of α with denominator ≤ Q_{n}. Lemma 4.2: If α is rational then it occurs as one of the approximations CF_{n}(α).

Moreover, it is easy to see that the continued fraction is easy to compute for any rational number.

8 Shor’s Quantum Factoring Algorithm

Below we give a quantum algorithm that factors N in polylog(N) time (factoring in poly(N) time is trivial and is too slow). It turns out the problem of factoring reduces to finding a nontrivial square root.

Claim: If we can find u such that u^{2} ≡ 1 (mod N) and u ≢ ±1 (mod N), then we can factor N. (Such a number u is called a nontrivial square root of 1 (mod N).)

Proof: The condition on u is equivalent to N | u^{2} − 1 = (u + 1)(u − 1) but N ∤ u + 1 or N ∤ u − 1. So gcd(N,u + 1) and gcd(N,u − 1) are nontrivial factors of N. □

To find a non-trivial square root, we simply pick a random number x (mod N) and compute its order, where the order of x is the least positive r such that x^{r} ≡ 1 (mod N).

Claim: If N is odd, with probability at least 1/2 over a random x ∈ Z_{N}, the order r of x is even and x^{r/2} ≢ ±1 (mod N).

Example: Let N = 15. Then let’s suppose we picked x = 7. Then x = 7,x^{2} = 4,x^{3} = 13,x^{4} = 1, so x has order 4. Now, taking y = x^{r/2} = 4, notice that y − 1 = 3 and y + 1 = 5 are both factors of 15.

To compute the order, we can use a quantum circuit to compute f(a) = x^{a} (mod N). The function maps element from Z_{M} to Z_{M} with M > N^{2}. This function can be implemented efficiently with Õ(log^{2}N) gates if we do modular exponentiation with repeated squaring.

Below we give the algorithm that, given an odd integer N, outputs a nontrivial factor of N with constant probability.

[1] Pick a random x ∈ Z_{N} gcd(x,N) > 1 Output gcd(x,N) Run period finding on f(a) = x^{a} (mod N) to get the order r of x (mod N) Compute x^{r/2} (mod N) x^{r/2} ≢ ±1 (mod N) Output gcd(N,x^{r/2}± 1) Abort

Note that GCD can be computed quickly using Euclid’s algorithm.

Because of randomized nature of part of the algorithm, we may need to repeat this procedure many times. Since the algorithm succeeds with constant probability, a constant number of repititions suffice to split N into two non-trivial factors. We can then repeat the procedure on each factor until we are down to the prime factors.