bloch-sphere-0

Quantum Fourier Transform

In order to keep advancing in our study of quantum algorithms, we now have to visit something called the fourier transform. The fourier transform is a mathematical tool that is well used outside of quantum computing but will prove useful for us.

Specifically, we will find that quantum computers are able to implement a special quantum version of the fourier transform that will be aplicable to many problems.

The discrete fourier transform is a mathematical operation that takes, as input, a sequence of numbers x0,,xN1x_0, \ldots, x_{N-1} and outputs a new sequence of numbers y0,,yN1y_0, \ldots, y_{N-1} where:

yk=1Nj=0N1xje2πijk/Ny_k = \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}x_je^{2\pi ijk/N}

Looks scary? Let's break it down.

Firstly, let's revise our complex numbers a bit. Remember that eiθe^{i \theta} merely represents a rotation of θ\theta radians around the complex unit circle. Since we know that 2π2\pi is a full rotation this means the exponential term is just jk/Njk/N full rotations around the circle.

In fact, to make this more obvious we will define a symbol to simplify our equation:

ωN=e2πi/N\omega_N = e^{2\pi i/N}

So ωN\omega_N is a 1/N1/N-th rotation around the unit circle. ωN2\omega^{2}_N is a 2/N2/N-th rotation and so on. We can then also see that ωNN\omega^{N}_N is the same as ωN0\omega^{0}_N since doing a full rotation is the same as doing no rotation. And from there it's fairly clear that we get relations like: ωN2=ωNN+2=ωN2N+2\omega^{2}_N = \omega^{N+2}_N = \omega^{2N+2}_N. So we end up with a circular pattern very similar to the powers of ii.

Just to be clear, let's now re-write our original equation with ω\omega:

yk=1Nj=0N1xjωNjky_k = \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}x_j\omega_{N}^{jk}

So our fourier transform is starting to make some more sense. Each term yky_k in the output sequence is a sum of all the input terms xjx_j multiplied by a rotation of jk/Njk/N around the unit circle.

The rotation patterns can be seen easier if we look at a concrete example in the case when N=5N=5:

y0=1N(x0ω50+x1ω50+x2ω50+x3ω50+x4ω50)y1=1N(x0ω50+x1ω51+x2ω52+x3ω53+x4ω54)y2=1N(x0ω50+x1ω52+x2ω54+x3ω56+x4ω58)y3=1N(x0ω50+x1ω53+x2ω56+x3ω59+x4ω512)y4=1N(x0ω50+x1ω54+x2ω58+x3ω512+x4ω516)y_0 = \frac{1}{\sqrt{N}}(x_0\omega^0_5 +x_1\omega^0_5 +x_2\omega^0_5 +x_3\omega^0_5 +x_4\omega^0_5) \\[3ex] y_1 = \frac{1}{\sqrt{N}}(x_0\omega^0_5 +x_1\omega^1_5 +x_2\omega^2_5 +x_3\omega^3_5 +x_4\omega^4_5) \\[3ex] y_2 = \frac{1}{\sqrt{N}}(x_0\omega^0_5 +x_1\omega^2_5 +x_2\omega^4_5 +x_3\omega^6_5 +x_4\omega^8_5) \\[3ex] y_3 = \frac{1}{\sqrt{N}}(x_0\omega^0_5 +x_1\omega^3_5 +x_2\omega^6_5 +x_3\omega^9_5 +x_4\omega^{12}_5) \\[3ex] y_4 = \frac{1}{\sqrt{N}}(x_0\omega^0_5 +x_1\omega^4_5 +x_2\omega^8_5 +x_3\omega^{12}_5 +x_4\omega^{16}_5)

We can then simplify all of our omegas using the rules we learned earlier:

y0=1N(x0+x1+x2+x3+x4)y1=1N(x0+x1ω51+x2ω52+x3ω53+x4ω54)y2=1N(x0+x1ω52+x2ω54+x3ω51+x4ω53)y3=1N(x0+x1ω53+x2ω51+x3ω54+x4ω52)y4=1N(x0+x1ω54+x2ω53+x3ω52+x4ω51)y_0 = \frac{1}{\sqrt{N}}(x_0 +x_1 +x_2 +x_3 +x_4) \\[3ex] y_1 = \frac{1}{\sqrt{N}}(x_0 +x_1\omega^1_5 +x_2\omega^2_5 +x_3\omega^3_5 +x_4\omega^4_5) \\[3ex] y_2 = \frac{1}{\sqrt{N}}(x_0 +x_1\omega^2_5 +x_2\omega^4_5 +x_3\omega^1_5 +x_4\omega^3_5) \\[3ex] y_3 = \frac{1}{\sqrt{N}}(x_0 +x_1\omega^3_5 +x_2\omega^1_5 +x_3\omega^4_5 +x_4\omega^2_5) \\[3ex] y_4 = \frac{1}{\sqrt{N}}(x_0 +x_1\omega^4_5 +x_2\omega^3_5 +x_3\omega^2_5 +x_4\omega^1_5)

It's not immediately clear what's useful about this. Usually this idea is first presented in a mathematics course where it can be used to analyze signals. In that context, the fourier transform is used to decompose a signal into its constituent frequencies. This is useful because it allows us to understand the signal better and to filter out noise.

For our case though, we'll simply stick with the definition and see how we can apply this to our quantum algorithms. It's use will hopefully become apparent later.

We will define the quantum fourier transform as follows:

y=1Nj=0N1ωNyjj\ket{y} = \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}\omega_N^{yj}\ket{j}

This looks fairly similar, but notice that we have yjyj in the omega exponent rather than jkjk like we did earlier. This is because instead of transforming a discrete sequence of numbers we are transforming a superposition quantum state.

Note also that although yy and jj are used like normal numbers in this equation, we know that the kets are actually formed of a tensor product of 1\ket{1}s and 0\ket{0}s that are the binary representation of the number. We should therefore assume that yy and jj are formed from the same number of bits and therefore N=2nN = 2^n for some nn.

To understand how we implement this in a quantum circuit, we will use some algebra to transform this summation formula into a different tensor product form, that will more easily translate into a circuit.

To do this, our first step is to split the jj summation into it's individual bits:

y=1Nj1=01jn=01ωNyjj1,,jn\ket{y} = \frac{1}{\sqrt{N}}\sum_{j_1=0}^{1} \ldots \sum_{j_n=0}^{1}\omega_N^{yj}\ket{j_1, \ldots, j_n}

Then we'll look at a slightly different representation of our ω\omega term. We know that ωNyj=e2πiyj/N\omega_N^{yj} = e^{2\pi iyj/N}. But since N=2nN = 2^n the division by NN here is simply shifting the bits of jj to after the decimal point. i.e. if n=4n = 4 and j=0101j = 0101 in binary, then we have j/2n=j/N=0.0101j/2^n = j/N = 0.0101.

In general, this shifting the bits to the right of the decimal point can be represented as:

j/N=j/2n=l=1njl2lj/N = j/2^n = \sum_{l=1}^{n}j_l2^{-l}

So our ω\omega term becomes:

ωNyj=e2πiyl=1njl2l\omega_N^{yj} = e^{2\pi iy\sum_{l=1}^{n}j_l2^{-l}}

And our main equation becomes:

y=1Nj1=01jn=01e2πiyl=1njl2lj1,,jn\ket{y} = \frac{1}{\sqrt{N}}\sum_{j_1=0}^{1} \ldots \sum_{j_n=0}^{1}e^{2\pi iy\sum_{l=1}^{n}j_l2^{-l}}\ket{j_1, \ldots, j_n}

Since the exponential co-efficient is now represented as a sum ofjjs individual bits, we can also separate the ket into a tensor product of the individual bits:

y=1Nj1=01jn=01e2πiyj121j1e2πiyjn2njn=1Nj1=01jn=01l=1ne2πiyjl2ljl\ket{y} = \frac{1}{\sqrt{N}}\sum_{j_1=0}^{1} \ldots \sum_{j_n=0}^{1}e^{2\pi iyj_12^{-1}}\ket{j_1} \otimes \ldots \otimes e^{2\pi iyj_n2^{-n}}\ket{j_n} \\[3ex] = \frac{1}{\sqrt{N}}\sum_{j_1=0}^{1} \ldots \sum_{j_n=0}^{1} \bigotimes_{l=1}^{n}e^{2\pi iyj_l2^{-l}}\ket{j_l}

And now, since we know the tensor product stitches together the 1\ket{1}s and 0\ket{0}s, we can move the tensor product to the outside to do this instead of summing over each bit:

y=1Nl=1njl=01e2πiy2ljl\ket{y} = \frac{1}{\sqrt{N}}\bigotimes_{l=1}^{n}\sum_{j_l=0}^{1}e^{2\pi iy2^{-l}}\ket{j_l}

And now our summation is just for 2 terms, so we can remove the big summation and write it explicitly:

y=1Nl=1n[e2πiy02l0+e2πiy12l1]y=1Nl=1n[0+e2πiy2l1]\ket{y} = \frac{1}{\sqrt{N}}\bigotimes_{l=1}^{n}\Bigl[e^{2\pi iy0\cdot 2^{-l}}\ket{0} + e^{2\pi iy1 \cdot 2^{-l}}\ket{1}\Bigr] \\[3ex] \ket{y} = \frac{1}{\sqrt{N}}\bigotimes_{l=1}^{n}\Bigl[\ket{0} + e^{2\pi iy2^{-l}}\ket{1}\Bigr]

The last step will involve that y2ly2^{-l} term in the exponent. If yy has bits y1,,yny_1,\ldots,y_n then multiplying by 2l2^{-l} is shifting these bits ll places to the right. e.g. if y=01011y = 01011 and l=3l=3 we have y2l=01.011y2^{-l} = 01.011. In general this is of the form: y2l=y1ynl.ynl+1yny2^{-l} = y_1 \ldots y_{n-l}.y_{n-l+1} \ldots y_n where the decimal place ends up between the nln-lth and nl+1n-l+1th digit.

But, notice that this number is multiplying 2πi2\pi i in the exponent so any whole part of this number (i.e. any number to the left of the decimal point) represents a full rotation so can be ignored. So we end up with:

y=1Nl=1n[0+e2πi0.ynl+1yn1]\ket{y} = \frac{1}{\sqrt{N}}\bigotimes_{l=1}^{n}\Bigl[\ket{0} + e^{2\pi i0.y_{n-l+1}\ldots y_n}\ket{1}\Bigr]

To clean this up, let's drop the big \otimes symbol and replace NN with 2n2^n:

y=(0+e2πi0.yn1)(0+e2πi0.yn1yn1)(0+e2πi0.y1yn1)2n/2\ket{y} = \frac{(\ket{0} + e^{2\pi i0.y_n}\ket{1})(\ket{0} + e^{2\pi i0.y_{n-1}y_n}\ket{1}) \ldots (\ket{0} + e^{2\pi i0.y_1 \ldots y_n}\ket{1})}{2^{n/2}}

This is the product representation of the quantum fourier transform. It's useful because it better shows the behaviour on individual bits and will allow us to construct a circuit easier.

To build the circuit, we'll need to invent some new gates. The quantum gates we've studied previously have mostly been focussed with discrete jumps around the Bloch Sphere. As we can see from our equation above, the quantum fourier transform will involve rotations of the phase of the 1\ket{1} state. So let's create the gates:

Rx0=0Rx1=e2πi/2x1R_x\ket{0} = \ket{0} \\[3ex] R_x\ket{1} = e^{2\pi i / 2^x}\ket{1}

So all of these gates leave the 0\ket{0} state alone and rotate the phase of the 1\ket{1} state by a 1/2x1/2^xth of a full rotation. i.e. R1R_1 rotates 1\ket{1} by 180 degress and R2R_2 rotates 1\ket{1} by 90 degrees etc.

Let's now present the quantum circuit for the quantum fourier transform and then examine why it works:

The quantum circuit for the quantum fourier transform consisting of Hadmard gates and R gates

Notice that our RR gates are controlled gates. Just like we've seen other controlled gates previously, these gates only apply their rotation if the control qubit is in the 1\ket{1} state.

To see that this works, let's start by just following the y1\ket{y_1} qubit for now.

Applying the first Hadamard gate gives us 12(0+e2πi0.y11)\frac{1}{2}(\ket{0} + e^{2\pi i0.y_1}\ket{1}) since if y1=0y_1 = 0 then e2πi0.y1=1e^{2\pi i0.y_1} = 1 and if y1=1y_1 = 1 then e2πi0.y1=1e^{2\pi i0.y_1} = -1.

The controlled R2R_2 gate then takes us to 12(0+e2πi0.y1y21)\frac{1}{2}(\ket{0} + e^{2\pi i0.y_1y_2}\ket{1}) since the 1\ket{1} state will undergo a further 90 degree rotation if y2=1y_2 = 1.

Continuing on in this fashion, the first qubit ends in the state 12(0+e2πi0.y1yn1)\frac{1}{2}(\ket{0} + e^{2\pi i0.y_1 \ldots y_n}\ket{1}). Our total state so far is then:

(0+e2πi0.y1yn1)21/2y2,,yn\frac{(\ket{0} + e^{2\pi i0.y_1 \ldots y_n}\ket{1})}{2^{1/2}}\ket{y_2, \ldots, y_n}

Continuing in this fasion we apply one fewer rotation for each qubit, so following the same steps for the second qubit brings us to:

(0+e2πi0.y1yn1)(0+e2πi0.y2yn1)22/2y3,,yn\frac{(\ket{0} + e^{2\pi i0.y_1 \ldots y_n}\ket{1})(\ket{0} + e^{2\pi i0.y_2 \ldots y_n}\ket{1})}{2^{2/2}}\ket{y_3, \ldots, y_n}

And following the pattern for the rest of our qubits we finally end up in the state:

y=(0+e2πi0.y1yn1)(0+e2πi0.y2yn1)(0+e2πi0.yn1)2n/2\ket{y} = \frac{(\ket{0} + e^{2\pi i0.y_1 \ldots y_n}\ket{1})(\ket{0} + e^{2\pi i0.y_2 \ldots y_n}\ket{1}) \ldots (\ket{0} + e^{2\pi i0.y_n}\ket{1})}{2^{n/2}}

This is almost what we want, but notice that the qubits are actually backwards compared to our equation from earlier!

This is easily solved by using n/2n/2 SWAP\text{SWAP} gates to move the positions of the qubits.

Great! So we've found a circuit that implements the quantum fourier transform!

We apply nn gates to the first qubit, n1n-1 gates to the 2nd etc. followed by n/2n/2 SWAP\text{SWAP} gates so we end up using a total of (n+(n1)++2+1)+n/2=n(n+1)2+n/2=12n2+n(n + (n - 1) + \ldots + 2 + 1) + n/2 = \frac{n(n+1)}{2} + n/2 = \frac{1}{2}n^2 + n gates. So this is fairly efficient, no exponential scaling.

We'll begin exploring in the next articles how this quantum fourier transform can be used.

> Next Article (Order Finding)