How do you get from this: $\displaystyle y_k \mapsto \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1} x_n e^{\frac{2\pi ijk}{N}}$ to this:?

Welcome to the 5th article in my series on Quantum Computing.

In the last article we spent the majority of our time setting up Qiskit in PyCharm. We encountered a couple of hurdles along the way (like getting Qiskit to generate charts), but in the end we managed to get it all working – but didn’t really do any quantum programming. That’s about to change.

The problem is, we don’t have any tools to convert mathematical operations in to equivalent quantum circuits. I know we had a brief look at how this was done for addition on classical circuits, but did you know that even addition is forbidden in quantum circuits? More on that later.

Instead, we are going to look at a different mathematical operation known as a fourier transform and see if we can get it constructed in to a quantum circuit.

This will require quite a bit of math. But if you can get through it, I guarantee you will feel a lot better about this step in the process of quantum computing. In fact, I think this step in particular, separates those who think they know how quantum circuits work and those that do.

This article will involve a range of more advanced topics but, like before, I will introduce each one very slowly. If you are not comfortable with the tensor product representation of a many-qubit state basis vector, you will have to brush up on that. We will also need to be familiar with the Hadamard gate and the Rotation gates. We will be working in the Fourier-transformed amplitude space, so a lot of the calculations will be rotations in the complex plane – get ready for a lot of complex exponentials. There will also be a large component that needs to be converted in to binary (so that we can see exactly how each input qubit maps to an output qubit). For this we will be delving in to the realm of fractional binary notation.

Lastly, we will finish this article off by theoretically building a circuit – yes, we will be going all the way from the Fourier Transform to an actual quantum circuit. In the next article we will implement the circuit in Qiskit’s Aqua in Python.

## The Quantum Fourier Transform Circuit

What is the Quantum Fourier Transform (QFT)?

First, it is a linear transformation – which means that it preserves addition and multiplication by scalars; and that it can also be used in quantum circuits (all quantum operations must be linear).

But this is nothing special, many interesting maps are linear. So what else?

Well, it’s unitary. OK, that’s particularly useful for quantum programming, it means it has a Hermitian adjoint. It also has a matrix representation. Also cool. But what does it actually do?

It maps vectors to funny-looking vectors:

$\displaystyle y_k \mapsto \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1} x_n e^{\frac{2\pi ijk}{N}}$

That looks ugly and there are a lot of letters in there to explain. First, it involves a summation over $N$ complex numbers; $N$ can be set to absolutely anything but trust me when I say that it is useful to set $N = 2^n$ where $n$ is the dimension of the vector. We often use the abbreviation:

$\omega_N := e^{\frac{2\pi i}{N}}$

which is a rotation, and so the map looks like this:

$\displaystyle y_k \mapsto \frac{1}{\sqrt{2^n}}\sum_{j=0}^{2^n-1} x_n \omega_N^{jk}$

and so it is clear that the QFT is rotating something around something.

What else can we do?

Well, we can use qubit notation for a basis state $|x\rangle$:

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^n}}\sum_{j=0}^{2^n-1} \omega_N^{jx} |k\rangle$

Note how $x$ replaces $k$ in the exponent.

Therefore, if we have a 3-qubit system, the QFT would map a basis vector like this:

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^3}}\sum_{j=0}^{2^3-1} \omega_{2^3}^{jx} |k\rangle$

and therefore, even though we only have 3 qubits, we will need an $8\times 8$-matrix to describe what is happening term-by-term.

Let’s see what the matrix representation of this 3-qubit system looks like:

For the first row, the first entry would be:

$\displaystyle \textbf{QFT}_{00} = \frac{1}{\sqrt{2^3}}\omega_{2^3}^{0\times 0} = 1$

The second entry would be:

$\displaystyle \textbf{QFT}_{01} = \frac{1}{\sqrt{2^3}}\omega_{2^3}^{0\times 1} = 1$

In fact, every entry along the first row is 1 due to the presence of the zero in the exponent.

Moving on to the second row. The first element would be:

$\displaystyle \textbf{QFT}_{10} = \frac{1}{\sqrt{2^3}}\omega_{2^3}^{1\times 0} = 1$

The second element on the second row:

$\displaystyle \textbf{QFT}_{11} = \frac{1}{\sqrt{2^3}}\omega_{2^3}^{1\times 1} = \omega^1$

The third element:

$\displaystyle \textbf{QFT}_{12} = \frac{1}{\sqrt{2^3}}\omega_{2^3}^{1\times 2} = \omega^2$

Continuing in this way, we arrive at the full matrix representation of the QFT:

$\displaystyle \textbf{QFT} = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & \omega^1 & \omega^2 & \omega^3 & \omega^4 & \omega^5 & \omega^6 & \omega^7 \\ 1 & \omega^{2} & \omega^{4} & \omega^{6} & \omega^{8} & \omega^{10} & \omega^{12} & \omega^{14} \\ 1 & \omega^{3} & \omega^{6} & \omega^{9} & \omega^{12} & \omega^{15} & \omega^{18} & \omega^{21} \\ 1 & \omega^{4} & \omega^{8} & \omega^{12} & \omega^{16} & \omega^{20} & \omega^{24} & \omega^{28} \\ 1 & \omega^{5} & \omega^{10} & \omega^{15} & \omega^{20} & \omega^{25} & \omega^{30} & \omega^{35} \\ 1 & \omega^{6} & \omega^{12} & \omega^{18} & \omega^{24} & \omega^{30} & \omega^{36} & \omega^{42} \\ 1 & \omega^{7} & \omega^{14} & \omega^{21} & \omega^{28} & \omega^{35} & \omega^{42} & \omega^{49} \\ \end{bmatrix}$

This bulky matrix can be simplified using a bit of algebra. Note that since $\omega^n$ is an $n^{th}$ root of unity then $\omega^a = \omega^b$ for any $a, b$ if $a = b+cn$ for some integer $c \in \mathbb{Z}$.

Since we are working with eighths (i.e. $2^3 = 8$) any time we see a $\omega^8$, or indeed, any multiple of $8$, like $\omega^{16}$ or $\omega^{32}$, we can replace it with unity.

Furthermore, we can reduce rotations to modulo 8, i.e. if we see a $\omega^{10}$, for example, we can write this as $\omega^{2}$ because ten is equal to 2 modulo 8.

The resulting matrix for QFT is:

$\displaystyle \textbf{QFT} = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & \omega^1 & \omega^2 & \omega^3 & \omega^4 & \omega^5 & \omega^6 & \omega^7 \\ 1 & \omega^{2} & \omega^{4} & \omega^{6} & 1 & \omega^{2} & \omega^{4} & \omega^{6} \\ 1 & \omega^{3} & \omega^{6} & \omega^{1} & \omega^{4} & \omega^{7} & \omega^{2} & \omega^{5} \\ 1 & \omega^{4} & 1 & \omega^{4} & 1 & \omega^{4} & 1 & \omega^{4} \\ 1 & \omega^{5} & \omega^{2} & \omega^{7} & \omega^{4} & \omega^{1} & \omega^{6} & \omega^{3} \\ 1 & \omega^{6} & \omega^{4} & \omega^{2} & 1 & \omega^{6} & \omega^{4} & \omega^{2} \\ 1 & \omega^{7} & \omega^{6} & \omega^{5} & \omega^{4} & \omega^{3} & \omega^{2} & \omega^{1} \\ \end{bmatrix}$

This is what you would call a concrete example of a QFT matrix. What we need is an generalised version of this matrix that holds for any number of qubits. Only then will we be able to see how to manipulate it for any scenario.

### The General Matrix Representation

We know that the outer product of two state vectors produces a matrix. We also know that with the QFT we have the input state vector $|x\rangle$ and the output state vector $|y\rangle$ so we must be able to compute the QFT matrix using the outer product of  these two, namely $|y\rangle\langle x|$.

$\displaystyle \mathbf{QFT} := \frac{1}{\sqrt{N}} \sum_{x=0}^{N-1} \sum_{y=0}^{N-1} \omega_N^{-yx} |y\rangle\langle x|$

This matrix will compute the entries of the associated $2^n \times 2^n$-square matrix, where $n$ is the number of qubits. Therefore, and as highlighted above, a QFT on a 3-qubit system will require an $8\times 8$-square matrix.

### But is it Unitary?

It is easy to show that the QFT is unitary, and, in fact, it is customary to denote it by $U_{\mathbf{QFT}}$, so we will do so from now on. We start by multiplying the QFT matrix by its own adjoint:

$\displaystyle U_{\mathbf{QFT}}U_{\mathbf{QFT}}^{\dagger} = \frac{1}{\sqrt{N}}\sum_{x=0}^{N-1} \sum_{y=0}^{N-1} \omega_N^{yx} |x\rangle\langle y| \frac{1}{\sqrt{N}}\sum_{x'=0}^{N-1} \sum_{y'=0}^{N-1} \omega_N^{-y'x'} |y'\rangle\langle x'|$

Writing it this way lets us put $|x\rangle\langle x'|$ and, crucially, $\langle y | y' \rangle$ which is nothing but $\delta_{y,y'}$ by definition of the inner product of basis states. The exponents in the rotation simply add as one expects, and the four summations can be combined in to one as we are summing over the same range. Thus:

$\displaystyle U_{\mathbf{QFT}}U_{\mathbf{QFT}}^{\dagger} = \frac{1}{N} \sum_{x,y,x',y'}^{N-1} \omega_N^{yx-y'x'} \delta_{y,y'} |x\rangle\langle x'|$

Now, the delta function $\delta_{y,y'}$, acting on the rotation $\omega_N^{yx-y’x’}$ will simply convert any $y'$ in to a $y$. Therefore, we can absorb the delta function in to the rotation, flat-lining any $y'$ it encounters (including in the summation):

$\displaystyle U_{\mathbf{QFT}}U_{\mathbf{QFT}}^{\dagger} = \frac{1}{N} \sum_{x,y,x'}^{N-1} \omega_N^{y(x-x')} |x\rangle\langle x'|$

But now the rotation $\omega_N^{y(x-x')}$ has the representation of another delta function! I.e. $\delta_{x,x'} = \omega_N^{y(x-x')}$, thus

$\displaystyle U_{\mathbf{QFT}}U_{\mathbf{QFT}}^{\dagger} = \frac{1}{N} \sum_{x,x'}^{N-1} \delta_{(x,x')} |x\rangle\langle x'|$

And, again, the delta function is going to turn any $x'$ in to $x$, and so

$\displaystyle U_{\mathbf{QFT}}U_{\mathbf{QFT}}^{\dagger} = \frac{1}{N} \sum_{x}^{N-1} |x\rangle\langle x| = \mathbf{1}$

Therefore, the QFT matrix is unitary.

## The QFT as a Quantum Circuit

Deriving a quantum circuit for the 3-qubit case is a little difficult to start with. Let’s go back to the 2-qubit case (remember, that means a $4\times 4$ matrix!). We know that the QFT performs the following mapping on basis states:

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^n}} \sum_{y=0}^{N-1} \omega_N^{-xy}|y\rangle$

Can we do anything with this expression? Well, a quantum circuit is qubit-dependent i.e. we will need to know exactly what to do with each qubit if we want to implement a QFT. Therefore, those state vectors in the above representation for a QFT will need to be expanded:

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^n}} \sum_{y_1,y_2,\dots,y_n \in \left\{0,1\right\}} \omega_N^{-x\sum_{k=1}^n 2^{n-k}y_k}|y_1,y_2,\dots,y_n\rangle$

Here we have split out the target state vector as a product of states (standard stuff):

$\displaystyle |y\rangle := |y_1,y_2,\dots,y_n\rangle = |y_1\rangle\otimes|y_2\rangle\otimes\cdots\otimes|y_n\rangle$

and therefore, the product in the exponent becomes a sum:

$\displaystyle \omega^{xy} = \omega_N^{xy_1}\times\omega_N^{xy_2}\times\cdots\times\omega_N^{xy_n} = \omega_N^{x\sum_{k=1}^n y_k}$

But because the original sum goes to $2^n - 1$, we will essentially have $2^n - k$ copies of the above sub-summation, so we get:

$\displaystyle \sum_{y=0}^{2^n-1} \omega^{xy} = \omega_N^{x\sum_{k=1}^n 2^{n-k}y_k}$

Now for each $n$ element in the exponent summation, we can pair it with exactly one of the $y_n$‘s in the target state basis vector, i.e. we could pair the $k=1$ terms like so:

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^n}} \sum_{y_1,y_2,\dots,y_n \in \left\{0,1\right\}} \omega_N^{-x 2^{n-1}y_1}|y_1\rangle$

$\vdots$

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^n}} \sum_{y_1,y_2,\dots,y_n \in \left\{0,1\right\}} \omega_N^{-x 2^{n-n}y_n}|y_n\rangle$

…and there would be $n$ of them related by a product. Therefore, we could simply write the the exponential of sum as the product of exponents with appropriately paired terms:

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^n}} \sum_{y_1,y_2,\dots,y_n \in \left\{0,1\right\}} \bigotimes_{k=1}^n \omega_N^{-x2^{n-k}y_k}|y_k\rangle$

Rearranging the sum and product

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^n}} \bigotimes_{k=1}^n \left( \sum_{y_k \in \left\{0,1\right\}} \omega_N^{-x2^{n-k}y_k}|y_k\rangle \right)$

Since $y_k$ can either be 0 or 1, we can expand the sum now without the product getting in the way:

$\displaystyle \sum_{y_k \in \left\{0,1\right\}} \omega_N^{-x2^{n-k}y_k}|y_k\rangle = \omega_N^{-x2^n\times 0}|0\rangle + \omega_N^{-x2^{n-1}\times 1}|1\rangle$

$\quad\quad\quad\quad\quad\quad\quad = \omega_N^0|0\rangle + \omega_N^{-x2^{n-1}}|1\rangle$

$\quad\quad\quad\quad\quad\quad\quad\quad\quad = 1|0\rangle + \omega_N^{-x2^{n-1}}|1\rangle$

$\quad\quad\quad\quad\quad\quad\quad = |0\rangle + \omega_N^{-x2^{n-1}}|1\rangle$

To get:

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^n}} \bigotimes_{k=1}^n \left( |0\rangle + \omega_N^{-x2^{n-k}}|1\rangle \right)$

### Expanding the Product Using Fractional Binary Notation

We would like to now expand the

The definition of the binary fraction is this:

$\displaystyle [0.x_1\cdots x_m] := \sum_{j=1}^m x_j 2^{-j}$

To see how this is used it is extremely useful to compute some examples. Consider first $[0.x_1]$. Here, $m=1$ ($m$ counts the number of terms in the expansion, since only $x_1$ exists beyond the decimal point $m=1$), and also $j$ starts from the first index of the first term, which, in this case, is also 1. Thus:

$\displaystyle [0.x_1] := \sum_{j=1}^1 x_j 2^{-j} = x_1 2^{-1} = \frac{x_1}{2}$

Pretty simple. What about $[0.x_1 x_2]$? Here, there are two terms beyond the decimal point, so $m=2$. The first index is, again, 1 so $j$ starts at 1, and so:

$\displaystyle [0.x_1x_2] := \sum_{j=1}^2 x_j 2^{-j} = x_1 2^{-1} + x_2 2^{-2} = \frac{x_1}{2}+\frac{x_2}{4}$

Again, pretty simple.

Things get a little more tricky with $[0.x_n]$. Here, there is again only one term, so $m=1$, but now $j$ must start from $n$ as it is the first index in the list. So we get:

$\displaystyle [0.x_n] := \sum_{j=n}^1 x_j 2^{-j} = x_n 2^{-n} = \frac{x_n}{2^n}$

What about $[0.x_{n-1} x_n]$? Here, two terms, so $m=2$. First term index is now $n-1$ so $j$ starts from $n-1$. Notice now that the terms must be in lexicographical order. The binary expansion is:

$\displaystyle [0.x_{n-1}x_n] := \sum_{j=n-1}^2 x_j 2^{-j} = x_{n-1} 2^{-(n-1)} + x_n 2^{-n} = \frac{x_{n-1}}{2^{n+1}}+\frac{x_n}{2^n}$

With this little tool in mind, let’s now return to this guy:

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2^n}} \bigotimes_{k=1}^n \left( |0\rangle + \omega_N^{-x2^{n-k}}|1\rangle \right)$

Let’s expand the product in our minds and just consider what the last term would look like expanded out, recalling that $\omega_N^{-x2^{n-k}} = e^{\frac{2\pi i}{2^n}\left(-x2^{n-k}\right)}|1\rangle$:

$\displaystyle |0\rangle + \exp\left(\frac{2\pi i}{2^n} x_n 2^{n-n}\right)|1\rangle = |0\rangle + \exp\left(2\pi i \frac{x_n}{2^n}\right)|1\rangle$

…and we can see that the last part of the exponential (the $\frac{x_n}{2^n}$ bit) can now be replaced with its binary form:

$\displaystyle |0\rangle + \exp\left(\frac{2\pi i}{2^n} x_n 2^{n-n}\right)|1\rangle = |0\rangle + \exp\left(2\pi i [0.x_n]\right)$

But there are $n$ more terms in the product expansion. Here is the 2nd-last term:

$\displaystyle |0\rangle + \exp\left(\frac{2\pi i}{2^n} x_{n-1} 2^{n-(n-1)}\right)\times\exp\left(\frac{2\pi i}{2^n} x_{n} 2^{n-n}\right)|1\rangle = |0\rangle + \exp\left(2\pi i \frac{x_{n-1}}{2^{n+1}}\right)\exp\left(2\pi i \frac{x_{n}}{2^{n}}\right)|1\rangle$

Taking a common factor of $2\pi i$ out, we get

$\displaystyle = |0\rangle + \exp\left(2\pi i \left(\frac{x_{n-1}}{2^{n+1}}+\frac{x_{n}}{2^{n}}\right)\right)|1\rangle$

Which, the last part in parentheses has binary form: $\frac{x_{n-1}}{2^{n+1}}+\frac{x_{n}}{2^{n}} = [0.x_{n-1}x_n]$, and so the 2nd last term looks like this:

$\displaystyle = |0\rangle + \exp\left(2\pi i [0.x_{n-1}x_n]\right)|1\rangle$

Continuing in the fashion all the way back to the first term we find a pattern of implementing longer and longer binary representations of the exponent. Eventually we get the full binary representation of the Quantum Fourier Transform:

$U_{\mathbf{QFT}}(|x_1 x_2 \cdots x_n\rangle) = \frac{1}{\sqrt{N}} \left( \left( |0\rangle + \exp\left(2\pi i [0.x_1x_2\dots x_n]\right)|1\rangle\right) \otimes \right.$

$\displaystyle \quad\quad\quad \left.\cdots \otimes \left( |0\rangle + \exp\left(2\pi i [0.x_{n-1}x_n]\right)|1\rangle\right) \otimes \left( |0\rangle + \exp\left(2\pi i [0.x_n]\right)|1\rangle\right)\right)$

This is a very useful form of the QFT for, specifically (and don’t forget!) the case where the number of qubits is some power of two, i.e. for $N=2^n$.

Why?

Because take a look at that first qubit:

$\displaystyle \frac{1}{\sqrt{N}} \left( |0\rangle + \exp\left(2\pi i [0.x_1x_2\dots x_n]\right)|1\rangle\right)$

It’s the ONLY one that depends on ALL values of all the other input qubits $\left\{x_1,x_2,\dots,x_n\right\}$. Each further qubit depends less and less on the input qubits! Also, the last single qubit has $\exp\left( -2\pi i[0.a] \right)$ which equals either $+1$ or $-1$ which is the Hadamard gate/transform.

## Combining What We Know to Build a Circuit

We already know that the Hadamard gate produces the following in a circuit of $2^n$ qubits:

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2}} \left( |0\rangle + \exp\left( -2\pi i [0.x_1]\right)|1\rangle\right)\otimes|x_2, x_3, \dots, x_n\rangle$

Note how the $x_1$ is consumed inside the exponent.

We also know that the Rotation Gate (a.k.a Phase Shift Gate) is

$\displaystyle R_k := \begin{bmatrix} 1 & 0 \\ 0 & \exp\left(\frac{-2\pi i}{2^k}\right) \end{bmatrix}$

If we now apply controlled rotation gates $R_2, R_3,\cdots, R_n$ on the Hadamard transformed state, we get

$\displaystyle |x\rangle \mapsto \frac{1}{\sqrt{2}} \left( |0\rangle + \exp\left( -2\pi i [0.x_1x_2\dots x_n]\right)|1\rangle\right)\otimes|x_2, x_3, \dots, x_n\rangle$

and we have perfectly reproduced the first term of the above, expanded QFT state.

Then we move to the next qubit, perform another Hadmard, then the appropriate (one fewer) rotation gates, to get the second qubit.

Do this all the way to the end and you have effectively (in binary) produced the logical output of a Fourier transform on $n$ qubits!

This is what it looks like on a circuit diagram:

## Summary

In this article we made the giant leap to being able to actually come up with a quantum circuit, at least in theory, starting from a mathematical equation. We had to expand and break up the Fourier transform equation in to bite size pieces, crucially by using its linearity and unitary properties. From here we needed to find a suitable ordering of the embedded product and sums and then re-expand it in binary to get at the individual qubits.

Once we got to that point all we had to do was notice that we were really just looking at a sequence of Hadamard and Rotation transformations! This is particularly easy to do when you have expressed your equation in binary notation.

In the next article we will seek to implement this theoretical circuit of QFT in to a working prototype in Qiskit’s Aqua in Python.