## Implementing the quantum Fourier transform with Qiskit

The quantum Fourier transform is a key building block of many quantum algorithms, from Shor’s factoring algorithm over matrix inversion to quantum phase estimation and simulations. Time to see how this can be implemented with Qiskit.

Recall that the quantum Fourier transform (or, depending on conventions, its inverse) is given by

$|x \rangle \mapsto \frac{1}{\sqrt{2^n}} \sum_s \eta^{xs} |s \rangle$

where $\eta$ is an 2n-th root of unity with n being the number of qubits in the register, i.e.

$\eta = \exp \frac{2\pi i}{2^n}$

How can we effectively create this state with a quantum circuit? They key to this is the observation (see the references below) that the result of the quantum Fourier transform can be written as a product state, namely as

$\sum_s \eta^{xs} |s \rangle = (|0\rangle + \eta^{2^{n-1} x}|1\rangle)(|0 \rangle + \eta^{2^{n-2}x}|1 \rangle) \cdots (|0 \rangle + \eta^x|1 \rangle)$

which you can easily verify by multiplying out the product and collecting terms.

Here we use the tensor product order that is prescribed by OpenQASM, i.e. the most significant bit is q[n-1]. This bit is therefore given by

$|0 \rangle + \eta^{2^{n-1}x} |1 \rangle$

Let us analyse this expression further. For that purpose, we decompose x into its representation as a binary number, i.e. we write

$x = \sum_{i=0}^{n-1} 2^i x_i$

with xi being the binary digits of x. If we now multiply this by 2n-1, we will get $2^{n-1} x_0$ plus a multiple of 2n. As $\eta$ is a root of unity, this multiple cancels out and we obtain that

$\eta^{2^{n-1}x} = \eta^{2^{n-1}x_0} = e^{\pi i x_0} = (-1)^{x_0}$

Thus, we can write the most significant qubit of the Fourier transform as

$|0 \rangle + (-1)^{x_0} |1 \rangle$

which is nothing but

$H |x_0\rangle$

Thus we obtain the most significant qubit of the quantum Fourier transform by simply applying a Hadamard gate to the least significant qubit of the input.

This is nice and simple, but what about the next qubit, i.e. qubit n-2? From the decomposition above, we can see that this is simply

$|0 \rangle + \eta^{2^{n-2}x} |1 \rangle$

Using exactly the same arguments as for the most significant qubit, we easily find that this is

$|0 \rangle + \eta^{2^{n-2}x_0} H |x_1 \rangle$

Thus we obtain this qubit from qubit 1 by first applying a Hadamard gate and then a conditional phase gate, i.e. a conditional rotation around the z-axis, conditioned on the value of x0. In general, qubit n-j is

$|0 \rangle + \prod_{i < j - 1} \eta^{2^{n-j+i}x_i} H |x_{j-1}\rangle$

which is a Hadamard gate followed by a sequence of conditional rotations around the z-axis, conditioned on the qubits with lower significance.

So we find that each qubit of the Fourier transform is obtained by applying a Hadamard followed by a sequence of conditional rotations. However, the order of the qubits in the output is reversed, i.e. qubit n-j is obtained by letting gates act on qubit j. Therefore, at the end of the circuit, we need to revert the order of the qubits.

In OpenQASM and Qiskit, a conditional rotation around the z-axis is called CU1, and there are swap gates that we can use to implement the final reversing of the qubits. Thus, we can use the following code to build a quantum Fourier transformation circuit acting on n qubits.

def nBitQFT(q,c,n):
circuit = QuantumCircuit(q,c)
#
#
for k in range(n):
j = n - k
circuit.h(q[j-1])
#
# there is one conditional rotation for
# each qubit with lower significance
for i in reversed(range(j-1)):
circuit.cu1(2*np.pi/2**(j-i),q[i], q[j-1])
#
# Finally we need to swap qubits
#
for i in range(n//2):
circuit.swap(q[i], q[n-i-1])
return circuit


Here is the circuit that this code produces for n=4. We can clearly see the structure – on each qubit, we first act with a Hadamard gate, followed by a sequence of conditional rotations with decreasing angle, conditioned on the less significant qubits, and finally reorder the qubits.

This is already a fairly complex circuit, and we need to find a way to test it. Let us look at the options we have. First, a quantum circuit is a unitary transformation and can be described by a matrix. In our case, it is especially easy to figure out what this matrix should be. Looking at the formula for the quantum Fourier transform, we find that the matrix describing this transformation with respect to the computational basis has the elements

$U_{ij} = \frac{1}{\sqrt{2^n}} \eta^{ij}$

The Qiskit frameworks comes with a simulator called the unitary_simulator that accepts a quantum circuit as input and returns the matrix describing that circuit. Thus, one possible test approach could be to build the circuit, run the unitary simulator on it, and to compare the resulting unitary matrix with the expected result given by the formula above. In Python, the expected result is produced by the following code

def qftMatrix(n):
qft = np.zeros([2**n,2**n], dtype=complex)
for i in range(2**n):
for j in range(2**n):
qft[i,j] = np.exp(i*j*2*1j*np.pi/(2**n))
return 1/np.sqrt(2**n)*qft


and the test can be done using the following function

def testCircuit(n):
q = QuantumRegister(n,"x")
c = ClassicalRegister(n,"c")
circuit = nBitQFT(q,c,n)

backend = Aer.get_backend('unitary_simulator')
job = execute(circuit, backend)
actual = job.result().get_unitary()
expected = qftMatrix(n)
delta = actual - expected
print("Deviation: ", round(np.linalg.norm(delta),10))
return circuit


The outcome is reassuring – we find that the matrices are the same within the usual floating point rounding differences.

After passing this test, a next reasonable validation step could be to run the algorithm on a specific input. We know that the QFT will map the state $|0 \rangle$ into an equal superposition of all elements of the computational basis. Conversely, we therefore expect that if we start with such a superposition, the QFT will map the superposition onto $|0\rangle$, at least up to a phase.

Let us try this out. Our test circuit will consist of a layer of Hadamard gates to create the equal superposition, followed by the QFT circuit, followed by a measurement. The resulting circuit for n=4 is displayed below.

It we run this circuit on the QASM simulator embedded into Qiskit, the result is as expected – for 1024 shots, we get 1024 times the output ‘0000’. So our circuit works – at least theoretically. But what about real hardware?

Let us compile and run the circuit targeting the IBM Q Experience 14 qubit device. If we dump the QASM code after compilation, we see that the overall circuit will have roughly 140 gates. This is already significant, and we expect to see some noise. To see how bad it is, I have conducted several test runs and plotted the results as a histogramm (if you want to play with this yourself, you will find the notebook on Github). Here is the output of a run with n=4.

We still see a clear peak at the expected result, but also see that the noise level is close to making the result unusable – if we did not know the result upfront, we would probably not dare to postulate anything from this output. With only three qubits, the situation becomes slightly better but is still far from being satisfactory.

Of course we could now start to optimize the circuit – remove cancelling Hadamard gates, remove the final swap gates, reorder qubits to take the coupling map into account and so on – but it becomes clear that with the current noise level, we are quickly reaching a point where even comparatively simple circuit will inflict a level of noise that is at best difficult to handle. Hopefully, this is what you expected after reading my posts on quantum error correction, but I found it instructive to see noise in action in this example.

## References

1. M. Nielsen, I. Chuang, Quantum Computation and Quantum Information, Cambridge University Press 2010
2. R. Cleve, A. Ekert, C. Macchiavello, M. Mosca, Quantum Algorithms revisited, arXiv:9708016

## Running the Deutsch-Jozsa algorithm on IBMs Q experience

In one of the previous posts, we have looked at the basics of the Qiskit package that allows us to create and run quantum algorithms in Python. In this post, we will apply this to model and execute a real quantum algorithm – the Deutsch-Jozsa algorithm.

Recall that the Deutsch-Jozsa algorithm is designed to solve the following problem. We are given a boolean function f

$f \colon \{0,1\}^n \rightarrow \{0,1\}$

in n variables. We know that the function is either constant or it is balanced (i.e. the number of times it is zero is equal to the number of times it is equal to 1). The task is to determine which of the to properties – balanced or constant – the function f has.

The first choice that we need to make is of course a suitable function f. To keep things simple, we will use a function with two variables. Maybe the most straightforward example for a balanced function in two variables is the classical XOR function.

$f(x_0, x_1) = x_0 \oplus x_1$

Thus our first task is to develop a quantum equivalent of this function, i.e. a reversible version of f.

Recall that in general, given a boolean function f, we can define a unitary transformation Uf by the following rule

$U_f(|x \rangle |y \rangle) = |x \rangle |y \oplus f(x) \rangle$

Thus Uf flips the target qubit y if f(x) is one and leaves it unchanged otherwise. In our case, combining this with the definition of the XOR function implies that we are looking for a circuit acting on three qubits – two input qubits and one target qubit – that flips the target qubit if and only if the two input qubits are not equal. From this textual description, it is obvious that this can be constructed using a sequence of two CNOT gates.

Let us now go through the individual steps of the Deutsch-Jozsa algorithm, using the optimized version published 1997 in this paper by R. Cleve, A. Ekert, C. Macchiavello and M. Mosca (this version is also described in one of my earlier posts but is does not hurt to recap some of that). The first step of the algorithm is to prepare a superposition state

$|\psi_0 \rangle = \frac{1}{2} \sum_x |x \rangle$

This is accomplished by applying a Hadamard gate to each of the two qubits that make up our primary register in which the variable x lives. We then add an ancilla qubit that is in the state $H|1 \rangle$, so that our state is now

$\frac{1}{2} \sum_x |x \rangle \otimes H |1 \rangle = \frac{1}{2\sqrt{2}} \sum_x |x , 0 \rangle - |x, 1 \rangle$

So far our circuit looks as follows.

Next we apply the oracle Uf to our state. According to the definition of the oracle, we have

$U_f |x, 0 \rangle = |x\rangle |f(x) \rangle$

and

$U_f |x, 1 \rangle = |x\rangle |f(x) \oplus 1\rangle$

Therefore we find that

$U_f \colon: |x, 0 \rangle - |x, 1 \rangle \mapsto |x\rangle |f(x)\rangle - |x\rangle |f(x) \oplus 1 \rangle$

From this formula, we can read off how Uf acts depending on the value of f(x). If f(x) is 0, the right hand side is equal to the left hand side and Uf acts trivially. If, however, f(x) is one, the right hand side is simply minus the left hand side, and Uf acts as multiplication by -1. Combining this with our previous formula, we obtain

$U_f \colon: \frac{1}{2} \sum_x |x \rangle \otimes H |1 \rangle \mapsto \frac{1}{2} \sum_x (-1)^{f(x)}|x \rangle \otimes H |1 \rangle$

Next, we apply again the Hadamard operator followed by the Pauli X gate to the third qubit (the ancilla), i.e. we uncompute the ancilla qubit. From the expression above, we can read off directly that the result left in the first two qubits will be

$|\psi \rangle = \frac{1}{2} \sum_x (-1)^{f(x)}|x \rangle$

This is the vector that we will have in the first two qubits of our quantum register when we have processed the following circuit.

Let us now calculate the overlap, i.e. the scalar product, between this vector and the initial superposition $|\psi_0 \rangle$. Clearly,

$\langle \psi_0 |\psi \rangle = \frac{1}{4} \sum_x (-1)^{f(x)}$

We find that this overlap is zero if and only if the function f is balanced. But how can we measure this scalar product? To see this, recall that the initial state $|\psi_0 \rangle$ is the result of applying the tensor product H2 of the two Hadamard operators to the first two qubits in the fiducial state. Thus we can write our scalar product as

$\langle \psi_0 |\psi \rangle = \langle 0 | H^2 |\psi \rangle = \langle 0 | H^2 \psi \rangle$

In other words, we can determine the scalar product by again applying a Hadamard operator to each of the first two qubits and measuring the overlap of the resulting state with the basis state $|00 \rangle$ which is the same thing as the probability to measure $|00 \rangle$ when we perform a measurement in the computational basis. Thus we finally obtain the following circuit

and the function f is balanced if and only if the probability to measure $|00 \rangle$ is zero.

Let us now turn this into Python code using Qiskit. I found it useful to create subroutines that act on a given circuit and build parts of the circuit which can then be tested independently from each other on a simulator before combining them. Here is the code to create the oracle for our balanced function f.

def createOracleBalanced(q,c):
circuit = QuantumCircuit(q,c)
circuit.cx(q[0], q[2])
circuit.cx(q[1], q[2])
circuit.barrier(q)
return circuit


Similarly, we can write routines that create the initial state and join the ancilla qubit.

def createInitialState(circuit):
circuit.h(q[0])
circuit.h(q[1])
circuit.barrier(q)

circuit.x(q[2])
circuit.h(q[2])
circuit.barrier(q)


Finally, we need to be able to uncompute the ancilla and to add the final measurements.

def uncomputeAncilla(circuit):
circuit.h(q[2])
circuit.x(q[2])
circuit.barrier(q)

circuit.h(q[0])
circuit.h(q[1])
circuit.barrier(q)
circuit.measure(q[0], c[0])
circuit.measure(q[1], c[1])
circuit.barrier(q)


A word on barriers. In this example code, we have added barriers after each part of the circuit. Barriers are an element of the OpenQASM specification and instruct the compiler not to combine gates placed on different sides of a barrier during optimization. In Qiskit, barriers have the additional benefit of structuring the visualization of a circuit, and this is the main reason I have included them here. In an optimized version, it would probably be safe to remove them.

After all these preparations, it is now easy to compile the full circuit. This is done by the following code snippet – note that we can apply the + operator to two circuits to tell Qiskit to concatenate the circuits.

q = QuantumRegister(3,"q")
c = ClassicalRegister(3,"c")
circuit = QuantumCircuit(q,c)
createInitialState(circuit)
circuit = circuit + (createOracleBalanced(q,c))
uncomputeAncilla(circuit)


We can then compile this code for a simulator or a backend as explained in my last post on the IBM Q experience (you can also find the full source code here). The following histogramm shows the result of the final measurement after running this on the ibmqx4 5 qubit quantum computer.

We see, as in the experiments conducted before, that the theoretically expected result (confirmed by running the circuit on a simulator) is blurred by noise – we get the value 00 in a few instances, even though the function is balanced. Even though the outcome can still be derived with a reasonable likelihood, we start to get a feeling for the issues that we might have with noise for more complex functions and circuits.

It is instructive to extract the compiled QASM from the resulting qobj and load that code into the IBM Q Experience composer. After a few beautifications, the resulting code (which you can get here) is displayed as follows in the composer.

If we compare this to the original circuit, we see that the compiler has in fact rearranged the CNOT gates that make up our oracle. The reason is that the IBMQX4 device can only realize specific CNOT gates. It can, for instance, implement a CNOT gate with q[1] as control as q[0] as target, but not vice versa. Therefore the compiler has added Hadamard gates to swap control and target qubit. We also see that the compiler has respected the barriers and not cancelled some of the double Hadamard gates. If we remove the barriers and remove all Hadamard gates that clearly cancel each other, we finally obtain a greatly simplified version of the circuit which looks as follows.

Of course this simplification hinges on the special choice of the function f. Nevertheless, it is useful – fewer gates mean fewer errors. If we compare the error rates of the optimized circuit with the new circuit, we find that while the original version had an error (i.e. an unexpected amplitude of $|00 \rangle$) of roughly 10%, a run with the optimized circuit showed an error of only 5,5%.

This completes our post. The next algorithm that we will implement is Shor’s algorithm and the quantum Fourier transform.

## Using Python to access IBMs quantum computers

In a previous post, we have looked at IBMs Q experience and the graphical composer that you can use to build simple circuits and run them on the IBM hardware. Alternatively, the quantum hardware can be addressed using an API and a Python library called Qiskit which we investigate in this post.

## Installation and setup

To be able to use Qiskit, there are some setup steps that we need to complete. First, we obviously have to install Qiskit. This is easy – we can simply use pip.

pip install qiskit


This will download the latest version of Qiskit, in my case (at the time of writing this) this was version 0.6.1. The next thing that you need is an API token to be able to access the IBM API. Assuming that you are already a registered user, you can create your token on the advanced tab of your user profile page.

In order to easily access your token from a Python script, it is useful to store the token locally on your hard drive. Qiskit uses a file qiskitrc in the ~/.qiskit folder in your home directory to store your credentials. The easiest way to create this file is the following code snippet

python -c 'from qiskit import IBMQ ; IBMQ.save_account("your_token")'


where obviously you need to replace the text inside the quotes with your token. After running this, you should find a file ~/.qiskit/qiskitrc containing your saved token (in clear text).

Once that file is in place, you can now easily load the credentials from within a Python program using the following code snippet

from qiskit import IBMQ


The method IBMQ.active_accounts() will also return a list of currently available accounts which can be useful for debugging purposes. Loading an account is only needed if we want to connect to the IBM site to use their quantum hardware or the online simulator, not if we use a local backend – more on this later.

## Circuits, gates and measurements

Let us now take a look at the basic data structures of Qiskit. A QuantumCircuit is what you expect – a collection of gates and registers. In Qiskit, a circuit operates on a QuantumRegister and optionally contains a ClassicalRegister which holds the results of a measurement. The following code snippet will create a quantum register with two qubits, a classical register with two qubits and a quantum circuit based on those registers.

from qiskit import QuantumCircuit
from qiskit import ClassicalRegister
from qiskit import QuantumRegister
q = QuantumRegister(2,"q")
c = ClassicalRegister(2,"c")
circuit = QuantumCircuit(q,c)


Next, we need to add gates to our circuit. Adding gates is done by calling the appropriate methods of the circuit object. The gate model of Qiskit is based on the OpenQASM language described in [1].

First, there are one qubit gates. The basic one qubit gates are rotations, defined as usual, for instance

$R_X(\Theta) = \exp \left( -i \frac{\Theta}{2} \sigma_X \right) = \cos \frac{\Theta}{2} - i \sigma_X \sin \frac{\Theta}{2}$

and similarly for the other Pauli matrices. Now it is well known that any rotation of the Bloch sphere can be written as a product of three rotations around y- and z-axis, i.e. in the form

$R_Z(\Phi)R_Y(\Theta)R_Z(\lambda)$

which is denoted by

$U(\Theta,\Phi,\lambda)$

in OpenQASM and Qiskit. For instance, $U(0, 0, \lambda)$ is a rotation around the z-axis and so forth. A short calculation shows that

$U(\Theta,\Phi,\lambda) = \begin{pmatrix} \exp \left(-\frac{i}{2}(\Phi + \lambda)\right) \cos \frac{\Theta}{2} & - \exp \left(-\frac{i}{2}(\Phi - \lambda)\right) \sin \frac{\Theta}{2} \\ \exp \left(\frac{i}{2}(\Phi - \lambda)\right) \sin \frac{\Theta}{2} & \exp \left(\frac{i}{2}(\Phi + \lambda)\right) \cos \frac{\Theta}{2} \end{pmatrix}$

Other gates can then be built from this family of one qubit gates. When it comes to multi-qubit gates, the only multi-qubit gate specified by OpenQASM is the CNOT gate denoted by CX, which can then again be combined with other gates to obtain gates operating on three and more qubits.

For qiskit, the available gates are specified in QASM syntax in the file qelib1.inc (see https://github.com/Qiskit/qiskit-terra/blob/master/qiskit/qasm/libs/qelib1.inc). Note that global phases are suppressed, so if you carry out the calculations, you will sometimes find a difference in the (unimportant) global phase between the result of the calculation and the results in Qiskit.

There is a couple of gates that you will often use in your circuits and that are summarized in the following table.

Gate Description
X Pauli X gate
Y Pauli Y gate
Z Pauli Z gate
S Phase gate $\text{diag}(1,i)$
T T gate $\text{diag}(1,e^{i\frac{\pi}{4}})$
CX CNOT gate

Gates take arguments that specify the qubits on which the gates operate. Individual qubits in a register can be addressed using an array-like notation. For example, to implement a circuit that applies an X gate to the first (least significant) qubit and then a controlled-NOT gate with this qubit as control qubit and the second qubit as target qubit, the following code can be used.

q = QuantumRegister(2,"q")
c = ClassicalRegister(2,"c")
circuit = QuantumCircuit(q,c)
circuit.x(q[0])
circuit.cx(q[0], q[1])


On real hardware, CNOTs can only be applied to specific combinations of qubits as specified in the coupling map of the device. However, when a circuit is prepared for execution on a specific device – a process called compilation – the compiler will deal with that by either reordering the qubits or adding additional swap operations.

Now we have gates, but to be able to run the circuit and measure the outputs, we still need measurements. These can easily been added with the measure method of a circuit, which accepts two parameters – the quantum register to measure and the corresponding classical register.

circuit.measure(q,c)


When the measurement step is reached during the processing of the circuit, the measurement is done – resulting in the projection of the state vector to the corresponding subspace – and the results of the measurements are copied into the classical register from which they can then be retrieved.

A nice property of Qiskit is its ability to visualize a quantum circuit. For that purpose, several classes called drawers are available in qiskit.tools.visualization. The circuit above, for instance, can be plotted with only one command

from qiskit.tools.visualization import matplotlib_circuit_drawer as drawer
my_style = {'cregbundle': True}
drawer(circuit, style=my_style)


and gives the nice representation

## Compiling and running circuits

Let us now actually run the circuit. To do this, we need a Qiskit backend. Qiskit offers several backends. The Aer package contains a few simulators that are installed locally and can be executed directly from a Python script or a notebook. Some of these simulators calculate the actual state vectors (the unitary simulator and the state vector simulator), but cannot deal with measurements, others – the QASM simulator – only provide statistical results but can simulate the entire circuit including measurements.

The IBMQ package can be used to connect to the devices offered by the IBM Q experience program, including an online simulator with up to 32 qubits and the actual devices. As for the composer, accessing the IBM Q experience devices does obviously require an account and available units.

In order to run a circuit, we first compile the circuit, which will create a version of the circuit that is tailored for the specific hardware, and then submit the circuit as a job.

backend = IBMQ.get_backend('ibmq_16_melbourne')
from qiskit import compile
qobj = compile(circuit, backend=backend, shots=1024)
job = backend.run(qobj)


Once the job has been submitted, we can poll its status using job.status() every few seconds. When the job has completed, we can access the results using job.result(). Every job consists of a certain number of shots, i.e. individual executions, and the method result.get_counts() will return a hash map that lists the measured outcomes along with how often that outcome was obtained. The following gist shows a basic Python script that assembles a circuit, compiles it, submits a job to the Q experience and prints the results.