Quantum phase estimation – the quantum algorithm Swiss army knife

When you are faced with a problem in linear algebra and have absolutely no idea what to do, an eigenvalue decomposition is the one thing that you would typically try first. In the world of quantum algorithms, the situation is similar – finding the eigenvalues of a matrix is a central building block of many quantum algorithm. The most important approach to finding the eigenvalues is called the quantum phase estimation algorithm.

The quantum phase estimation algorithm

In its most basic form, the quantum phase algorithm is used to find an eigenvalue of a unitary matrix assuming that we know an eigenvector (we will see later that this is in fact not necessarily needed for useful applications). So let us assume that we are given a unitary matrix U, implemented as a quantum circuit, and a quantum state $|\psi \rangle$ which is an eigenvector of U. As U is unitary, all eigenvalues are on the complex unit circle, and therefore the eigenvalue equation is

$U |\psi \rangle = e^{2\pi i \varphi} |\psi \rangle$

for some number $\varphi$ between zero and one. The quantum phase estimation is designed to determine (or at least approximate) the number $\varphi$, i.e. the phase of the eigenvalue (hence the name).

The algorithm consists of two parts. Both parts operate on an ancillary working register and a register that we call the primary register. Let n denote the number of qubits in the primary register and m denote the number of qubits in the working register. We assume that the initial state in the primary register is the state $|\psi \rangle$. The first step of the algorithm is to prepare the system in the state

$\frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m - 1} |k \rangle U^k |\psi \rangle$

How can we prepare this state? Given a state $|k \rangle$, with binary digits ki, we can of course write Uk as

$U^k = \prod_i \big[ U^{2^i} \big]^{k_i}$

In other words, Uk is a sequence of conditional gates U2i, conditioned on the i-th bit of k. Together with the usual approach to prepare an equal superposition of all $|k \rangle$ by using Hadamard gates, this argument shows that the desired state can be prepared by the following quantum circuit.

So far we have not yet used that the initial state $|\psi \rangle$ is an eigenstate. Let us do this now. We know that a power of U will act on $|\psi \rangle$ as multiplication by the corresponding power of the eigenvalue, hence

$\frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m - 1} |k \rangle U^k |\psi \rangle = \frac{1}{\sqrt{2^m}} \big[\sum_{k=0}^{2^m - 1} e^{2\pi i \varphi k} |k \rangle \big] \otimes |\psi \rangle$

This is interesting – our state turns out to be a product state (as the multiplication by a scalar can be pulled through the tensor product). To understand the first part of the product (i.e. the resulting state in the working register) in more detail, let us assume for a moment that the phase $\varphi$ is an exact multiple of 2-m, i.e. that

$\varphi = \frac{t}{2^m}$

for some integer t. Then the state in the working register is

$\frac{1}{\sqrt{2^m}} \sum_{k=0}^{2^m - 1} e^{2\pi i \frac{tk}{2^m}} |k \rangle$

But this is simply the inverse quantum Fourier transform applied to the state $|t \rangle$! So we can obtain t and therefore $\varphi$ by applying a quantum Fourier transform to the state in the working register and doing a measurement. This is the quantum phase estimation algorithm.

In most cases, the phase will of course not be an exact multiple of t / 2m. However, with a bit more work (see for instance arXiv:quant-ph/9708016 or my own notes) one can see that the our result is still approximately true, i.e. if we run this algorithm, measure the value in the working register and divide by 2m, we get a good approximation to the actual phase $\varphi$. The precision of the estimation depends on the number of qubits m that we have in our working register. The following diagram demonstrates this for a specific example for two different values of m.

The x-axis shows the absolute deviation from the true value, the y-axis shows the probability for this deviation. The upper image corresponds to m=3. We see that there is a peak around zero, but still a substantial probability to have a significant error. The lower diagram shows the outcome for m=8. Here we see that the peak is very sharp, and so we obtain a very good approximation with high probability.

Shor’s algorithm as an instance of quantum phase estimation

Most of the time, the quantum phase estimation algorithm is used as a subroutine in a more general algorithm. A prominent example is Shor’s algorithm which – as we will now see – can be expressed as an application of the quantum phase estimation.

Recall that the quantum part of Shor’s algorithm is concerned with finding the period of a number a modulo M, where M is some large integer that we want to factor. Here we assume that a is a unit modulo M so that one of its powers is one. As a is a unit, the mapping

$U \colon |x \rangle \mapsto |ax \mod M\rangle$

is a unitary transformation U for x < M (which we can extend by the identity to a unitary transformation of the full Hilbert space). If r denotes the period of U, then clearly the fact that ar = 1 modulo M implies that Ur = 1. Thus all eigenvalues of U are r-th roots of unity. In other words, if an eigenvalue has phase $\varphi$, then $\varphi r$ is an integer. We can therefore use the phase estimation algorithm to obtain the period r.

There is one twist, however. To be able to apply the phase estimation algorithm, we need an eigenvector as an initial state, and it is far from clear how to prepare this. Fortunately, a deeper analysis (again see arXiv:quant-ph/9708016 or my own notes for all the details of this) shows that the state $|1 \rangle$ is a sum of eigenstates. Running the algorithm on this state will – by linearity – produce a state which is a superposition of states corresponds to different multiples of 2m / r, and performing the measurement will choose one of these multiples. As in Shor’s original version of the algorithm, we can then use a continued fraction expansion to obtain the unknown period r. Thus, the following algorithm can be used to determine the period r.

1. Prepare a primary register in the state $|1 \rangle$
2. Apply the circuit shown above
3. Apply a quantum Fourier transform to the working register
4. Perform a measurement of the working register and call the result s – this will be close to a multiple of 2m / r
5. Perform a continued fraction expansion to find the period r as in Shor’s original version of the algorithm

At the first glance, it looks as if we had found an alternative period finding algorithm, but if you actually draw up the circuit that realizes this algorithm, you will find – again see the references cited before – that this is not true, in fact the resulting circuit will be almost identical to the circuit for Shor’s version of the algorithm. Thus we have not found a new algorithm, but have found a different description of the same algorithm – period finding can be described as an instance of the phase estimation algorithm.

More applications

Of course factoring integers is not the only possible application of the quantum phase estimation. An obvious application is finding eigenvalues of Hamiltonians, i.e. energy eigenvalues – and in particular ground state energies – of physical systems. We have already touched on this possible application in a previous post on the quantum variational eigensolver.

Another application is based on the observation that if we start the algorithm with a superposition of eigenvectors, the measurement at the end of the algorithm will not only select an eigenvalue, but also an eigenvector. Thus we can use the algorithm to decompose a state into a superposition of eigenvectors. This can be used to solve linear equations of the form Ax = b for a given matrix A. In fact, if we can diagonalize A, then solving this equation becomes equivalent to division. Implementing this as a quantum algorithm is not as simple and straightforward as this description might suggest, as the division is a non-unitary operation, but it turns out that it can be approximated by a unitary transformation. The resulting algorithm is called the HHL-algorithm after its inventors Harrow, Hassidim and Lloyd, and has been published in this paper on the arxiv.

The HHL algorithm does not guarantee a quantum advantage for every matrix A, as it relies on a way to efficiently calculate the corresponding time evolution operators eiAt. However, if for instance A is sparse and has low condition number (ratio between eigenvalues), the algorithm provides an exponential speedup over the best classical algorithms. An important application is deep machine learning, where efficient matrix inversion plays an important role – see for instance this paper by Zhao, Pozas-Kerstjens, Rebentrost, and Witte, as well as the corresponding source code. The details of the HHL algorithm are far from being trivial, and we will leave a more in-depth discussion to a future post.

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.

Quantum error correction with stabilizer codes

In our previous discussion of quantum error correction, we have assumed that quantum gates can act on any two physical qubits. In reality, however, this is not true – only nearby qubits and interact, and our error correction needs to take the geometric arrangements of the qubits into account. The link between these geometric constraints and the theory of quantum error correction is the stabilizer formalism described first in [2].

This post will be a bit more formal, but is a necessary preparation to be able to define and understand the surface code and other topological codes. To introduce and motivate the formalism, let us for a moment come back to the simple example studied before – the three-bit code. In this code, the logical states were

$|0 \rangle_L = |000 \rangle$

and

$|1 \rangle_L = |111 \rangle$

These logical states span our two-dimensional code space and describe one logical qubit. A single bit flip error will flip one of the three involved qubits. Suppose, for instance, that a bit flip occurs on the first or second qubit. This error will modify our logical states as follows (the first two lines represent a bit flip on the first qubit, the next two lines a bit flip on the second qubit

$|0 \rangle_L \mapsto |100 \rangle$

$|1 \rangle_L \mapsto |011 \rangle$

$|0 \rangle_L \mapsto |010 \rangle$

$|1 \rangle_L \mapsto |101 \rangle$

Now consider the operator

$M_1 = \sigma_z^1 \sigma_z^2$

where the upper index on a Pauli matrix indicates the qubit on which the operator acts. The error-free logical states are both eigenstates of this operator with eigenvalue +1. The states that result out of a bit-flip error are also eigenstates of this operator, but with eigenvalue -1. Thus, a measurement of this operator will tell us whether a bit flip on the first or second qubit has occurred. Similarly, to also detect a bit flip on the third qubit, we need a measurement for a second operator, namely

$M_2 = \sigma_z^1 \sigma_z^3$

Let us now write down a few properties of these operators. First, they are hermitian and therefore correspond to measurements. Second, the logical states are eigenvectors of these operators with eigenvalues +1 and the code space is exactly the subspace of the three-qubit state space that is fixed by M1 and M2. We can think of these operators as defining linear constraints that together define the code space, as indicated below.

Moreover, if a bit flip error occurs, the resulting state will again be an eigenvalue, but for at least one of the operators the new eigenvalue will be -1. Thus, errors correspond to violated constraints. Finally, the two operators commute and can therefore be measured simultaneously.

These properties allow us to express the three bit code entirely in terms of the operators M1 and M2. The code space is the subspace which is left invariant by these operators. The syndrome measurement can be done by measuring both operators simultaneously, which is possible because they commute. Finally, all involved states – code space and states after bit flip errors – are eigenstates and therefore the measurement process does not collapse a superposition of these states and can therefore be executed without interfering with the actual quantum algorithm that we want to run.

This correspondence between sets of operators – the Mi in our case – and quantum error correction is the core of the stabilizer formalism. In a more general setting, we are looking at the Hilbert space of dimension 2n spanned by n physical qubits. Products of n Pauli matrices (where we include the identity matrix) are acting on this Hilbert space and form a group ${\mathcal G}_n$ called the Pauli group. Put differently, the Pauli group consists of those linear operators on the Hilbert space that can be written as a product

$g = \mu (\sigma_x^1)^{a_1} (\sigma_y^1)^{b_1} (\sigma_z^1)^{c_1}\cdots (\sigma_x^n)^{a_n}(\sigma_y^n)^{b_n} (\sigma_z^n)^{c_n}$

with coefficients $a_i, b_ic, c_i$ and an overall phase factor $\mu \in \{ \pm 1, \pm i \}$ (we can even assume that all the $b_i$ are equal to one as $\sigma_x \sigma_z$ is a multiple of $\sigma_y$). Any two elements of the Pauli group either commute or anti-commute. Within this group, we now consider a finite set $\{ M_i \}$ of commuting elements and the subgroup S of the Pauli group generated by this set. This group (which is abelian as it is generated by mutually commuting elements) is called the stabilizer group. To this group, we can associate the subspace that consists of all vectors that are fixed by all elements of the group, i.e. the space

$T_S = \{ |\psi \rangle \, \text{such that} \, s|\psi \rangle = |\psi \rangle \forall s \in S \}$

This space will be the code space of our code. If the group S is generated by l independent generators, the dimension of the code space can be seen to be 2n-l.

In the example of the three bit code, we had n=3 and l=2, which gives us a two-dimensional code space, corresponding to one logical qubit (if you try work out the details, you will see that in order for this to be true, we have to assume that $-1 \notin S$ – you might want to take a look at my notes or the chapters on quantum error correction in [1] for the mathematical details).

Which errors is our code able to detect? Suppose that E is an error operator that is itself in the Pauli group. We know that any two elements of the Pauli group either commute or anti-commute. Let us assume that E does in fact anti-commute with some element s of S. If $|\psi \rangle$ is a state in the code space, we can then calculate

$s E |\psi \rangle = - E s |\psi \rangle = - E |\psi \rangle$

Therefore the state $E |\psi \rangle$ is now in the -1 eigenspace of s. Thus if we measure all elements of S, the outcome -1 for one of the measurements will tell us that an error has occurred. A similar argument shows that if the product of any two errors anti-commutes with at least one element of S, then we can also correct the error based on only the measurement results of elements in S. Mathematically speaking, the set of all elements of the Pauli group that commute with all elements of S is the centralizer of S, which in this case turns out to be equal to the normalizer N(S) of S, and a set $\{E_\alpha \}$ of errors can be detected and corrected if and only if

$E_\alpha E_\beta \in S \cup (\mathcal{G} - N(S))$

for any two error operators $E_\alpha, E_\beta$. So the elements of the Pauli group that are neither in S nor in N(S) correspond to correctable errors.

What do the elements of the normalizer correspond to? Suppose that n is an element in the normalizer and therefore commutes with all elements of S. Then, for any element s in S and any element $|\psi \rangle$ of the code space, we have

$s n |\psi \rangle = n s |\psi \rangle = n |\psi \rangle$

Consequently, the element $n |\psi \rangle$ is again fixed by all elements of S and therefore in the code space TS. Thus the elements of N(S) generate automorphisms of the code space, i.e. logical operations. It is one of the strengths of the stabilizer formalism that once we have the stabilizer group S, we can not only derive the code space and the correctable errors, but can also use group theoretic considerations to describe logical operations – this will become clearer as we study surface codes in a later post.

For later reference, let us briefly summarize what we have found. Given l independent generators si of the stabilizer S, we can define a code space as the set of all vectors that are fixed by the si. This subspace will encode n – l logical qubits in n physical qubits. To detect an error, we measure all the si. Each measurement will give us minus one or plus one. Any occurrence of minus one indicates an error (this is why the combined result of all measurements is usually called the error syndrome) and will also tell us which operation we have to apply to correct the error. To implement logical quantum gates on our code space, we can apply elements in the normalizer N(S) that map the code space into itself. Thus we have expressed error detection, error correction and logical operations entirely in terms of the stabilizer group and the language of group theory.

But how do we find useful stabilizer codes? One approach is given by the check matrix formalism, which allows us to express stabilizer codes in terms matrices and metrics, i.e. in terms of linear algebra. This approach is generalizing the parity check matrix from the classical theory of error correction. A second source of stabilizers, however, comes from a more surprising direction – algebraic topology. In fact, given a surface described in terms of a cell complex, we can associate elements of the Pauli group with every cycle and every co-cycle. For certain choices of cycles and cocycles, this will give us abelian subgroups of the Pauli group which in turn create error correction codes. This is how geometrical constraints in the actual implementation enter the scene and leads to a class of codes known as surface codes and toric codes that we will start to study in my next post.

1. E. Rieffel, W. Polak, Quantum computing – a gentle introduction, MIT Press, Cambridge 2011
2. D. Gottesman, Stabilizer Codes and Quantum Error Correction, Caltech Ph.D. Thesis, available as arXiv:quant-ph/9705052
3. J. Kempe, Approaches to quantum error correction, S\’eminaire Poincar\’e 1 (2005), pp. 65–93, available as arXiv:quant-ph/0612185

Fault tolerant quantum computing

In the previous post, we have looked at the basic ideas behind quantum error correction, namely the encoding of logical states so that we are able to detect and correct errors like bit flip and phase flip errors. Unfortunately, this is not yet good enough to implement quantum computing in a fault-tolerant way.

What is still missing? First, we have only discussed errors that impact saved quantum states, i.e. random bit flips or phase flips that invalidate a given state stored in a quantum register. We have not yet, however, understood how we can protect quantum states during processing, i.e. deal with errors that do not occur on a saved qubit but that occur while we manipulate these qubits by executing quantum gates.

Second, and maybe even worse, implementing error correction involves adding additional qubits and quantum gates to our circuit. Obviously, these additional elements will again be subject to errors and require protection. So we might need an additional correction step to account for these errors. This additional step will again introduce additional errors and so forth. It is not a priori clear that we do not introduce more errors than we can correct when following this path. Fortunately, the famous threshold theorem asserts that under certain conditions, we can win this race and correct errors faster than we create them.

Error propagation and transversal quantum gates

But let us start with the first question – how can we implement quantum gates in a fault tolerant way? To illustrate the challenge, we start with a simple (though unrealistic) example. Suppose that we wanted to perform a CNOT gate on two states encoded with the three-bit code. The obvious approach – decode the two input states, apply an ordinary CNOT and encode again – does not appear to be a good choice, as uncontrolled errors might sneak in while the states are unencoded and therefore vulnerable. So let us try to implement a CNOT operation directly on the encoded states. This can be done with the following circuit.

Let us see why this circuit works. Suppose the control word is the encoded (logical) state

$|0 \rangle_L = |000 \rangle$

Then all three qubits in the control word are zero. Consequently, the three CNOT gates do nothing and the target word is not changed. If, on the other hand, the control word is

$|1 \rangle_L = |111 \rangle$

then all three qubits in the target word will be toggled, and thus the target word $|0 \rangle_L$ will be mapped to $|1\rangle_L$. Thus our circuit does in fact implement a CNOT gate on the logical level.

However, there is a problem with this circuit. To see this, suppose that there is a bit flip error in the first qubit of the control word. Then this error will make all three CNOT gates toggle their target state, and a logical $|0 \rangle_L$ on the target qubits will be flipped into a logical $|1 \rangle_L$. Thus even a single error in the input will render the result unusable. This is an example for uncontrolled error propagation: a single error in the control word input has propagated across all three qubits that make up the second input and our code – which is only able to protect against single qubit errors – fails.

Obviously, the propagation of errors can never be entirely avoided. However, what we can avoid is the uncontrolled propagation. To explain this, let us take a look at the following circuit.

This circuit also realizes a CNOT on the level of logical states. However, for this circuit, an error in one of the qubits of the block of three qubits that encode the logical control bit can only propagate into one qubit in the second block – the block of three qubits that encode the target bit. Therefore one error per input block will only propagate into one error in each output block. If we are able to fully recover from one qubit errors per block, we can detect and correct the errors even though they have propagated through the circuit (again as long as only one error occurs). Thus the trick is to somehow spread potential errors across multiple blocks so that no single block is fully corrupted. This connectivity property of a circuit is sometimes called transversality.

Transversal gates are fault-tolerant, as an error in one qubit can only propagate into one qubit of any other code block, a situation from which we can recover. For many codes, many basic gates can be implemented qubit-wise and therefore transversally on the level of logical qubits. Unfortunately, a full set of universal quantum gates cannot be constructed in such a way – typically one gate remains for which no transversal implementation can be found. Fortunately, there are other (considerably more complicated) methods available to construct fault-tolerant logical gates. In [1], Shor has described the first fully universal set of fault-tolerant logical gates for his nine qubit code. In many modern implementations, a class of codes called surface codes are used, that also allow us to implement a universal set of quantum gates in a fault tolerant way and that we will study in a later post. Similar approaches work for measurements that can have errors as well and that need to be organized in a way so that errors can be detected and corrected, and the same applies for the actual error correction circuits.

Fault tolerant quantum computation and the threshold theorem

Even if we are able to implement quantum gates in a fault tolerant way, one challenge remains. So far, we have only discussed how to protect against single qubit errors. Of course, the rationale behind this was the assumption that errors in two or more qubits are much less probable than single qubit errors. Thus, we have reduced the error probability. For real applications, our new, reduced error probability might still not be sufficient. We therefore might have to add additional error correction circuits. Now the probability of failure for a circuit is not only a function of the reliability of every single component, but also a function of the number of components. The more error correction mechanisms we add, the higher the number of components and the higher the overall fault probability. In the worst case, the increase of the fault probability by adding the additional components needed for the error correction might eat up all the benefits of the error correction and lead to an overall higher probability of failure.

To be able to better grasp the problem, let us introduce some terminology. Suppose we are given a circuit C that describes a quantum computation. We can think of the overall computation as being organized in individual time steps. At each time step, every qubit in the circuit takes part in a basic set of possible operations (or, to mention a different term that is sometimes used, locations). This can be a gate, the preparation of a qubit in a initial state, a measurement or a “wait step”, i.e. the identity operation leaving the involved qubits unchanged. At each time step, a qubit is involved in exactly one of these operations.

Given such a circuit, we can now try to build a fault tolerant circuit performing the same computation as follows. We replace each qubit in the circuit C by a block of qubits encoding a logical qubit. Similarly, we replace each operation in C by the corresponding fault tolerant equivalent, i.e. we replace each gate by the corresponding logical gate, a measurement by a fault tolerant measurement, and a preparation by a fault tolerant state preparation. In addition, we place a (fault tolerant) error correction circuit in each block which is the output of a logical gate. The result is called a simulation of the original circuit. An example is displayed below, using (roughly) the graphical language described in [2].

The upper part of this diagram shows an original circuit, in which we prepare two qubits in a known state, then apply the unitary transformation U and finally measure both qubits. The result of the replacement procedure described above is displayed in the lower part of the picture, using again the (unrealistic) example of a three qubit code. Here we have replaced each qubit by a corresponding block of three qubits holding the code word. The operation U has been replaced by its logical implementation, acting on the encoded states, and both after the preparation and after applying U, we add a circuit performing error correction.

Let us now study how the fault probability changes when we pass to a simulation. In general, this is a long and complicated argument, and we refer the reader to [3], [4] and [5] for the details. However, there is a simple example that can be used to illustrate some of the key ideas.

In fact, let us go back to the circuit above that we have used for a fault-tolerant implementation of the CNOT gate. To analyze the probability that this circuit fails, let us assume that the only error that can occur is a bit flip error in the target qubit of each physical CNOT gate and that this error occurs with probability p. Thus we assume that error correction, state preparation and measurement work in an idealized, error free way.

Let us now look at the different combinations of faults that can occur during a specific run of the simulation circuit, i.e. a different fault paths. A first fault path is given by a single bit flip error in one of the target qubits after executing the CNOT operation. However, this affects only one qubit and can therefore be corrected by the error correction circuit. Thus, for this fault path, the simulation works correctly.

There is, however, a different fault path, given by a fault in target qubit one and target qubit two, whereas the CNOT acting on target qubit three operates correctly. Assuming that the individual faults are statistically independent, the probability for this fault path is

$p \cdot p \cdot (1 - p) = p^2 - p^3$

This would introduce an error from which we could not recover, as our code is only able to correct errors in one qubit. But this is not yet the overall probability of failure for our new circuit, as there is a second fault path that we need to consider – we could as well have a failure in the first and the third instead of the first and second qubit. Obviously, the probability for this fault path is the same. And there is a third fault path with two errors, given by an error in the second and third target qubit. These three different fault paths, all with two errors, correspond to the

${{3} \choose {2}} = 3$

possibilities to choose two error locations out of the three target qubits. And finally, it might happen that faults occur in all three locations, an event which clearly has probability p3. Thus the overall probability to have at least two faults is

$3 (p^2 - p^3) + p^3 = 3 p^2 - 2 p^3 \leq 3 p^2$

Thus the overall probability that the circuit does not operate correctly does not reduce from p – for the physical CNOT operation without error protection – to p2, but to 3 p2, where the factor 3 reflects the growth of the circuit when passing from the physical circuit to the simulation, i.e. the impact of the additional gates.

Of course this is a very simplified example, but the hard work done in the original papers in the references demonstrates that the result is true in general: if we pass from a circuit to its simulation, the change of the error rate is described by a rule of the form

$p \mapsto c p^{t+1}$

where c is a constant that reflects the growth of the circuit and t is the number of errors that the code that we use can still correct. Obviously, this does not help unless the right hand side is smaller than the left hand side, i.e. unless p is smaller than a threshold given by

$p < c^{-\frac{1}{t}} = p_{crit}$

If this is true, we can even iterate the procedure, and reduce the error rate below every prescribed level by adding additional layers of simulation. Thus given physical operations with error rate less than pcrit, we can implement each quantum circuit in a fault-tolerant way with arbitrary reliability.

Unfortunately, the exact value of the threshold is not known. The initial estimate provided in [3] was in the range of 10-6. Later, this was improved by several authors, and is mostly assumed to be in the range between 10-4 and 10-2 (see [6] for an overview of some results in this direction).

As c is in the order of pcrit-t, this implies that we need as many as a few hundred to a few thousand physical gates and qubits to simulate one logical qubit. Note that in general, the overhead depends on the desired accuracy and the level of simulations needed as well as on the factors entering the constant c, i.e. the encoding, the architecture of the logical gates, as well as how close the actual physical error rate is to the threshold. Improving the threshold and thereby reducing the overhead remains one of the most active areas of current research in quantum computing and is essential for paving the road towards a large-scale fully fault tolerant universal quantum computer.

At the time of writing, a class of codes known as surface codes appear to be good candidates for the implementation of fault-tolerant quantum computers, not only because these codes have good error correction properties, but also because they are well adapted to the geometry of the actual quantum computer. In the next few posts, we will work towards understanding surface codes and other topological codes, starting with the formalism of stabilizer codes that will be the focus of the next post.

1. P. Shor, Scheme for reducing decoherence in quantum computer memory, Phys. Rev. A Vol. 52 (1995) Issue 4, pp. 2493–2496
2. D. Gottesman, An Introduction to Quantum Error Correction and Fault-Tolerant Quantum Computation, arXiv:0904.2557v1
3. D. Arahonov, M. Ben-Or, Fault Tolerant Quantum Computation with Constant Error, arXiv:quant-ph/9611025
4. E. Knill, R. Laflamme, W.H. Zurek, Threshold accuracy for quantum computation, arXiv:quant-ph/9610011,
5. P. Aliferis, D. Gottesman, J. Preskill, Quantum accuracy threshold for concatenated distance-3 codes, arXiv:quant-ph/0504218v3
6. S. J. Devitt, W. J. Munro, K. Nemoto, Quantum Error Correction for Beginners, arXiv:0905.2794v4

Basics of quantum error correction

Do usable universal quantum computers exist today? If you follow the recent press releases, you might believe that the answer is “yes”, with IBM announcing a 50 qubit quantum computer and Google promoting its Bristlecone architecture with up to 72 qubits. Unfortunately, the world is more complicated than this – time to demystify the hype a bit.

The need for error correction

The important point is that it is not just the number of qubits that matters, but also their quality. When we study quantum algorithms like Shor’s algorithm, we are working with idealized qubits that behave exactly the way an isolated two-state quantum system is supposed to behave – these idealized qubits are often called logical qubits. However, in a real world implementation, there is no fully isolated two-state quantum system. Every system interacts to some extent with the environment, an interaction that we can try to reduce to a minimum, for instance by cooling our device down to very low temperatures, but never fully avoid. A trapped ion could, for instance, interact with radiation entering our device from the environment, and suddenly is part of a larger quantum system, consisting of the ion, the photons making up the radiation and maybe even the source of the photon. This will introduce errors into our system, i.e. deviations of the behavior of the system from the idealized theoretical model.

In addition to unwanted interactions with the environment, other errors could creep into our computation. When manipulating qubits to realize gates, we might make mistakes, for instance by directing a microwave pulse with a slightly incorrect frequency at our qubit, and we can make mistakes during each measurement.

Thus the real qubits in a quantum computer – called physical qubits – are prone to errors. These errors might be small for one qubit, but they tend to propagate through the circuit and add up to a significant error that will render the result of our quantum computation unusable. Thus we need error correction, i.e. the ability to detect and correct errors during our computation.

So how would you do this? Of course, errors can also occur in classical systems, and there are well developed methods to detect and correct them. Unfortunately, these approaches typically rely on the ability to copy and measure individual bits. This is easy for a classical bit, but more complicated for a qubit, as a measurement will collapse our system into an eigenstate of the observed operator and thus interfere with quantum algorithm.

The good news is that quantum error correction is still possible. In this post, I will try to explain the basics, before we then dive into more advanced topics in the next few posts.

Encoding logical states

In order to understand the basic ideas and structures behind quantum error correction, it is useful to study a simplified example – the three qubit code. To introduce this code, let us suppose that we have access to a communication channel across which we can send individual qubits from one quantum device to another one. Suppose further that this transmission is not perfect, but is subject to a bit flip error with a probability p. Thus, with probability p, a one-qubit state $|\psi \rangle$ will be changed to $X |\psi \rangle$ during transmission, where $X$ is the usual bit flip operator, and with probability 1-p, the transmission does not change the state. The aim is to construct an encoding of a qubit such that these errors can be detected and corrected.

To achieve this, we encode every single qubit state in a three-qubit state before transmitting it. Thus we use the following encoding

$|0 \rangle \mapsto |000 \rangle$

$|1 \rangle \mapsto |111 \rangle$

It is not difficult to see that this encoding can in fact be realized by a unitary circuit – the circuit below will do the trick.

To transmit one qubit, we use this encoding to obtain a three-qubit message. We then send those three qubits through our communication channel. After the transmission, we apply a procedure known as syndrome measurement, using the following circuit.

Let us see what this circuit is doing. First, let us suppose that the original qubit was $|0 \rangle$, encoded as $|000 \rangle$, and no error did occur during the transmission. Then the three qubits at the top of the circuit will still be $|000 \rangle$. In this case, the CNOT gates act as the identity, and the overall state after passing the circuit is

$|000 \rangle |00 \rangle$

Similarly, if the original qubit is $|1 \rangle$ and no error occurred, all CNOT gates will act as inversion. Thus the ancilla qubits will be inverted twice, and we end up with the state

$|111 \rangle |00 \rangle$

The situation is a bit different if a bit flip error has affected one of the qubits during the transmission, say the first one. Suppose the original state was again $|0 \rangle$. After encoding and transmission with a bit flip on the first qubit, we will receive the state $|100 \rangle$. Therefore both ancilla qubits will be inverted, and we obtain the state

$|100 \rangle |11 \rangle$

Let us now generalize these considerations. Assume that we are encoding the state $a |0 \rangle + b |1 \rangle$, so that the encoded state will be $a |000 \rangle + b |111 \rangle$. When we go through the above exercise for all possible cases, we arrive at the following table that shows the transmission error and the resulting state after passing the syndrome measurement circuit.

Error Resulting state
No error $(a |000 \rangle + b |111 \rangle ) |00 \rangle$
Bit flip on first qubit $(a |100 \rangle + b |011 \rangle ) |11 \rangle$
Bit flip on second qubit $(a |010\rangle + b |101 \rangle ) |10 \rangle$
Bit flip on third qubit $(a |001\rangle + b |110 \rangle ) |01 \rangle$

Now let us see what happens if we measure the ancilla qubits. First, note that all the states are already eigenstates for the corresponding measurement operator. Thus measuring the ancilla qubits will not change the state of the first three qubits and it will not reveal any information on the encoded state. The second important observation is that the value of the ancilla qubits tells us the exact error that has occurred. Thus we have found a way to not only find out that an error has occurred without destroying our superposition, but also to figure out which qubit was flipped. Given that information, we can now apply a bit flip operator once more to the affected qubit to correct the error. Again, this will not reveal the values of a and b and not collapse our state, and we can therefore continue to work with the encoded quantum state, for instance by running it through the inverse of the encoding circuit to get our original state back.

More general error models

So what we have found so far is that it is possible, also in the quantum world, to detect and protect against pure bit flip errors without destroying the superposition. But there is more we can learn from that example. In fact, let us revisit our original assumption that the only thing that can go wrong is that the operator X is applied to some of the qubits and allow a more general operator. Suppose, for instance, that our error is represented by applying to at most one of our qubits the operator

$E = (1-\epsilon) 1 + \epsilon X$

for a small $\epsilon$. In contrast to our earlier assumption that the error operator is discrete, i.e. is either applied or not applied, this operator is now continuous, depending on the parameter $\epsilon$, i.e. it looks as if we had to deal with a full continuous spectrum of errors. A short calculation shows that after transmitting and applying the syndrome measurement circuit above, the state of our quantum system will now be

$(1-\epsilon)(\alpha |000 \rangle + \beta |111 \rangle)|00\rangle + \epsilon (\alpha |100 \rangle + \beta |011\rangle) |11 \rangle$

Now let us again apply a measurement of the ancilla qubits. Then, according to the laws of quantum mechanics, the system will collapse onto an eigenstate, i.e. it will – up to normalization – end up in one of the states

$(\alpha |000 \rangle + \beta |111 \rangle)|00\rangle$

and

$(\alpha |100 \rangle + \beta |011\rangle) |11 \rangle$

But these are exactly the states in which we end up if no error occurs or a single bit flip error occurs. Thus, our measurement forces the system to somehow decide whether an error occurred or not and if yes, which error occurred – we are opening the box in which Schrödingers cat is hidden. This is a very important observation often referred to as the digitization of errors – it suffices to protect against discrete errors as the syndrome measurement will collapse any superposition of different errors states.

So far we have worked with a code which is able to protect against a bit flip error. But of course this is not the only type of error that can occur. At the first glance, it looks like there is a vast universe of potential errors that we have to account for, as in theory, the error could be any unitary operator. However, using arguments similar to the discussion in the last paragraph, one can show that it suffices to protect against two types of errors: the bit flip error discussed above and the phase flip error, represented by the matrix

$Z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$

Note that the phase flip error has no classical equivalent, other than the bit flip error which can classically be interpreted as a bit flipping from one to zero or vice versa randomly. The first error correction code that was able to handle both, bit flip errors and phase flip errors (and combinations thereof, ), and therefore any possible type of error (as long as the number of errors is limited) was described in 1995 by P. Shor. This code uses nine qubits to encode one logical qubit and is somehow a repeated application of the three bit code, based on the observation that the Hadamard transform turns a bit flip error into a phase flip error and vice versa. Later, different types of codes were discovered that are also universal in the sense that they protect against any potential one qubit error, i.e. against any error that only affects one qubit.

Let us now look at the structure that these codes have in common. First, the encoding can be described as identifying a 2-dimensional subspace C, called the code space, in the larger space spanned by all qubits. In the case of the three qubit code, the code space is spanned by the states $|000 \rangle$ and $|111 \rangle$, corresponding to a logical zero and a logical one. More generally, the space C can have dimension 2n and the full Hilbert space can have dimension 2k, in which case we can encode n qubits in a larger set of k qubits (in the nine qubit code example, k = 9 and n = 1).

Hence the code space is spanned by a set of 2n basis vectors $|\psi_i \rangle$ that we call code words. In the case n = 1, it is common practice to denote the codewords by $|0 \rangle_L$ and $|1 \rangle_L$ to indicate that they represent a logical qubit encoded using k physical qubits.

In addition, there is a set of error operators, i.e. a finite set of operators $\{ E_\alpha\}$ like the phase flip or bit flip operators acting on the larger Hilbert space. These operators represent the impact of discretized noise and will generally move the code words out of the code space, i.e. they will rotate the code space onto different subspaces of the entire Hilbert space.

In order for a code to be useful, we of course need a relation between code space and errors that tells us that the errors can be detected and corrected. What are these conditions? A first condition is that no matter which errors occur, we can still tell the code words apart. This is guaranteed if the various error operators map different code words onto mutually orthogonal subspaces, in other words if

$\langle \psi_i | E_\alpha^{+} E_\beta | \psi_j \rangle = 0$

whenever $i \neq j$ and for all $\alpha, \beta$. Thus, even in the presence of errors, the different code words will never overlap.

What about the case that both code words are the same? It is tempting to ask for the condition

$\langle \psi_i | E_\alpha^{+} E_\beta | \psi_i \rangle = \delta_{\alpha \beta}$

i.e. to require that different errors map the same codeword to orthogonal subspaces. This would make recovery very easy. We could perform the measurements that correspond to projections onto these subspaces to detect the error and correct them by applying the inverse of the error operator. However, this condition is too restrictive. In the case of the nine qubit code, for example, it might very well happen that two different errors map a code word to the same state. However, the same applies for the correction, i.e. we do not have to distinguish between these two errors as the act similarly on the code space. Therefore, a more general condition is usually used which captures this case:

$\langle \psi_i | E_\alpha^{+}E_\beta | \psi_j \rangle = C_{\alpha\beta} \delta_{ij}$

for a hermitian matrix $C_{\alpha \beta}$. If the matrix has full rank, the code is called non-degenerate, otherwise – as in the case of the nine qubit code – the code is called degenerate.

Of course this encoding generates some overhead. To represent one logical qubit, we need more than one physical qubit. For the Shor code, we have to use nine physical qubits to encode one logical qubit. It is natural to ask what the minimum overhead is that we need. In 1996, Steane discovered a code that requires only seven physical qubits. In the same year, a code that requires only five qubits was presented and it was shown that this is a lower bound, i.e. there is no error correction code that requires less than five qubits.

So there is an unavoidable overhead – and the situation is even worse. To implement error correction, you need again quantum gates, which can of course also experience errors. Thus you need additional circuitry to protect the error correction against errors, which can again introduce errors and so forth. That there is a way out of this vicious circle is not obvious and the content of the famous threshold theorem that we will study in a later post – but even this way out is very hard to implement and might require thousands of physical qubits to implement one single logical qubit.

So even with a 72 qubit device, we are still far away from implementing only one logical qubit – and having a few thousand logical qubits to use Short’s algorithm to break an RSA key with a realistic key length is yet another story. So it is probably a good idea to take claims about universal supremacy of quantum computing within a few years with a grain of salt.

In this post, we have looked at some of the essential ideas behind quantum error correction, i.e. the ability to detect and correct errors. However, this is not enough to build a reliable quantum computer – after all, adding an error correction circuit introduces additional qubits that can also create new errors. In addition, we need to be able to perform calculations on our encoded states. So there is more that we need for fault-tolerant quantum computing which is the topic of the next post.

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.

Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

A few more features of the Qiskit package, including working with different simulators and visualization options as well as QASM output, are demonstrated in this script on my GitHub page. In one of the next posts, we will try to implement a real quantum algorithm, namely the Deutsch-Jozsa algorithm, and run it on IBMs device.

References

1. Andrew W. Cross, Lev S. Bishop, John A. Smolin, Jay M. Gambetta, Open Quantum Assembly Language, arXiv:1707.03429v2 [quant-ph]
2. The Qiskit tutorial on basic quantum operations
3. The Qiskit documentation

Quantum simulation

In his famous lecture Simulating Physics with computers, Nobel laureate Richard Feynman argued that non-trivial quantum systems cannot efficiently be simulated on a classical computer, but on a quantum computer – a claim which is widely considered to be one of the cornerstones in the development of quantum computing. Time to ask whether a universal quantum computer can deliver on that promise.

At the first glance, this seems to be obvious. Of course, there is one restriction – similar to classical computing, the resources of a quantum computer are finite, se we can only hope to be successful if we can at least approximately reduce the Hilbert space of the original problem to a finite dimensional space. But once this hurdle is taken, everything else looks easy. The time evolution of the system is governed by its Hamiltonian H – over a period of time t, the evolution of the system is described by the unitary operator

$U(t) = e^{iHt}$

This is a unitary operator, so given a universal set of quantum gates, we should – at least approximately – be able to implement this unitary evolution as a sequence of quantum gates. And we are able to store a state in a Hilbert space of dimension 2N in N qubits, in contrast to the exponential number of classical bits that we would need. So what is the problem?

The problem is hidden in the complexity to construct an arbitrary unitary operator from a finite set of quantum gates. Suppose we are considering a quantum system with a finite number of qubits N so that a unitary operator is described by an 2N x 2N matrix. As argued in [1], if we assemble unitary operators from a set of universal gates, we need in the order of 22N steps, as this number of parameters is required to specify an arbitrary (binary) matrix of size N. Thus we have mastered the memory bottleneck only to run into a time bottleneck, as the number of steps still grows exponentially with the dimension of the Hilbert space.

In addition, figuring out how to build up the unitary matrix from the given set of universal gates – a task typically accomplished with a classical computer – can still be extremely difficult and scale exponentially with N. So it seems that we are fighting a lost battle (see also section 4.5.4 of [2] for an argument why building general unitary matrices from an elementary set of gates is hard).

Fortunately, this is only true if we insist on simulating any conceivable unitary matrix. In practice, however, most Hamiltonians are much simpler and built from local interactions – and this turns out to be the key to efficiency.

As an example to illustrate this point, let us consider a typical system of interest – a lattice in space that has a particle of spin 1/2 attached to each of its vertices. Suppose further that any particle only interacts with adjacent particles. Then, in the easiest case of constant coupling, the Hamiltonian of this system can be written as

$H = \sum_{\text{pairs} \, i,j} S_i S_j$

with spin operators Si. The number of terms in the sum is proportional to N as we only allow pairs of adjacent spins to interact – actually it is roughly equal to 2N (each spin site pairs up with four other spins, and we need to divide by two as we double count – and there will be slight deviations from this rule depending on the boundary conditions). At every spin site, we have a two-dimensional Hilbert space describing that spin. If two sites i and j are adjacent, the operator SiSj acts on the tensor product of the Hilbert spaces attached to the involved sites, i.e. on a subspace of dimension four of the overall Hilbert space. The crucial point is that even though the dimension of the overall Hilbert space is huge, each term in the Hamiltonian only acts on a rather small subspace.

Let us now consider more general Hamiltonians of this form. Specifically, let us assume that we can write our Hamiltonian as a sum

$H = \sum_{s=1}^l H_s$

of l = 2N hermitian operators, where each Hs only acts on a subspace of dimension ms. Let m be an upper bound on the ms.

Of course, it is general not true that the exponential of H is the product of the exponentials of the Hi, unless the Hi commute. So

$e^{iHt} \neq \prod_{s=1}^l e^{iH_st}$

However, it turns out that this equality becomes approximately true for very small times t. In fact, the Trotter formula states that

$e^{iHt} = \lim_{n \rightarrow \infty} (\prod_{s=1}^l e^{iH_st/n})^n$

Thus, to simulate the development of the system over a certain time t, we could proceed as follows.

• Choose a large number n
• Initialize the quantum computer in a state corresponding to the initial state
• Apply each of the l operators eiHst/n in turn
• Repeat this n times
• Measure to determine the result

Note that in the last step, we typically do not measure all the individual qubits, but instead we are interested in the value of an observable A and measure that observable. In general, we will have to repeat the entire simulation several times to determine the expectation value of the measurement outcome.

We can visualize the evolution of the quantum state under this process similar to the process of approximating a smooth curve by moving along the coordinate axis. If, for instance, l = 2, so that H = H1 + H2, we will periodically first apply eiH1t, then eiH2t and so forth. We can think of the two operators as coordinates, and then applying these operators corresponds to moving along the coordinate axis. If we move in sufficiently short intervals, we will at least approximately move along the smooth curve.

Now let us try to understand the complexity of this process. Each individual Hs acts on a space with dimension at most m. The important part is that if we enlarge the system – for instance by adding additional sites to a spin lattice – the number m is not substantially changed. Thus the unitary operator eiHst/n can be approximated by universal quantum gates in the order of m2 steps, and the entire process takes in the order of nlm2 steps. If l grows at most polynomially with N, as in our example of the spin lattice, the complexity therefore grows at most polynomially with N as well. In general, this will be exponential faster than on a classical computer.

The upshot of the discussion is that even though general unitary transformations are hard to approximate with a universal set of gates, local Hamiltonians – and therefore the Hamiltonians of a large and very rich class of interesting quantum systems – can be efficiently simulated by a universal quantum computer.

The method presented above (sometimes called the Trotter-based quantum simulation) is quite universal, but in some cases, special problems admit more specialized solutions that are also often termed as simulations. Suppose, for instance, we wanted to know the lowest energy eigenvalue of a quantum system. We could do this with our simulation approach – we start with an arbitrary state and simulate the system for some time, hoping that the system will settle in the ground state.

Obviously, this is not very efficient, and we might want to employ different methods. We have already seen that one possible approach to this is the variational quantum eigensolver treated in a previous post. However, there is also an approach related to the universal simulation described above. Basically, the idea described in [4] is to use quantum phase estimation to find the eigenvalues of U = eiHt and hence of H. To do this, one needs to conditionally apply powers of the matrix U and this in turn can be done using the Trotter-based approach. The result will be an eigenvector, i.e. an energy eigenstate, and the corresponding eigenvalue. The eigenvector can then be employed further to evaluate other observables or to start a simulation from there.

In all approaches, the quantum system at hand has to be mapped to the model of a quantum computer, i.e. to qubits and unitary quantum gates. The efficiency and feasibility of this mapping often depends on the chosen representation (for instance first versus second quantization) of the original system, and also determines the resource requirements of the used quantum computer. Finding suitable mappings and being able to prepare good initial states efficiently are key challenges in real applications, especially as the number of available qubits remains limited, and a very active area of research – see [5] for an example.

Finally, what we have discussed so far is often called digital simulation or universal simulation, as it simulates an arbitrary system with the same, fixed computational model of qubits and quantum gates. A different approach to simulation is to build a system which is sufficiently close to the system under interest, but easier to control and which can therefore be placed in a specific state to observe the evolution. This is sometimes called analog quantum simulation.

The ideas explained in this post have been developed in the late nineties, and since then several researchers have demonstrated working prototypes on small quantum computers with a few qubits (and without error correction). We refer the reader to [3] for an overview of some recent developments and results and a collection of references for further reading.

References

1. S. Lloyd, Universal quantum simulators, Science, New Series, Vol. 273, No. 5278 (Aug. 23, 1996), pp. 1073-1078
2. M.Nielsen, I. Chuang, Quantum computation and quantum information theory, Cambridge University Press
3. J. Olsen et. al., Quantum Information and Computation for Chemistry, arXiv:1706.05413
4. D.S. Adams, S. Lloyd, A quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors, Phys.Rev.Lett.83:5162-5165,1999 or arXiv:quant-ph/9807070
5. A. Aspuru-Guzik, A.D. Dutoi, P.J. Love, M. Head-Gordon. Simulated Quantum Computation of Molecular Energies. Science 309, no. 5741 (September 9, 2005), pp. 1704–1707 or arXiv:quant-ph/0604193

Navigating downhill: the quantum variational eigensolver

In quantum mechanics, the dynamics of a system is determined by its Hamiltonian, which is a hermitian operator acting on the Hilbert space that describes the system at hand. The eigenstates and eigenvalues of the Hamiltonian then correspond to stationary states and their energies, and finding these eigenstates and the corresponding eigenvalues is the central task in computational quantum mechanics.

Unfortunately, in most cases, this is very hard. First, the Hilbert space of most physically relevant systems is infinite-dimensional. However, even if we are able to approximate the Hamiltonian by a hermitian operator acting on a finite dimensional subspace, finding the eigenvalues and eigenstates by applying classical methods from numerical linear algebra like the QR method is computationally challenging due to the high dimensional spaces involved. It is natural to ask whether a quantum computer can help.

In fact, there are several methods available for finding eigenvalues of a hermitian matrix on a quantum computer. In 1995, A. Kitaev described an algorithm which is now known as the quantum phase estimation (see [2] and [3]) which can be used to estimate eigenvalues of unitary matrices (and can be applied to hermitian matrices as well noting that if A is a hermitian matrix, U = eiAt is unitary and that there is an obvious relation between the eigenvalues of U and A). Unfortunately, this algorithm might require a large number (millions or even billions) of fault-tolerant quantum gates which is currently completely out of reach. This is the reason why a second algorithm which was first described in [1] and is designed to run on a quantum computer with a low number of qubits attracted a lot of attention. This algorithm, called the variational quantum eigensolver, might be able to deliver improvements over what is possible with classical hardware with somewhere between a few tens and one hundred qubits.

Setup and preliminaries

To explain the approach, let us assume that we are given a quantum register with n qubits and some hermitian operator H acting on the corresponding 2n-dimensional Hilbert space. We assume further that we can write our operator as a linear combination

$H = \sum_i h_i H_i$

of a small number of operators Hi which correspond to measurements that we can easily implement. An example could be a decomposition

$H = \sum_i h_i P_i$

where each Pi is a tensor product of Pauli matrices (it is not difficult to see that such a decomposition exists for every hermitian operator. First, write the operator as a linear combination of projections $|x_i \rangle \langle x_i |$ and the vectors as linear combinations of tensor products to show that each hermitian operator is a linear combination of tensor products of hermitian single-qubit operators, and then use the fact that each hermitian single qubit operator is a linear combination of Pauli matrices). A tensor product of Pauli matrices corresponds to a measurement that can be easily implemented by a combination of a measurement in the standard computational basis and a unitary transformation, see for instance this discussion of Pauli measurements.

Now if the operators Hi correspond to measurements that can efficiently be implemented, it is also possible to efficiently evaluate the expectation value

$\langle H_i \rangle = \frac{\langle \psi | H_i | \psi \rangle }{\langle \psi | \psi \rangle}$

for a given state $|\psi \rangle$, assuming that we can efficiently prepare the state N times and then conduct N independent measurements – in fact, if N is sufficiently large, the expectation value will simply be the average over the measurements (see [4], section IV, for a more detailed discussion) – here we see a trade-off between the number of times we need to prepare the state $|\psi \rangle$ and measure and the resulting precision.

Once we are able to measure the expectation value of each Hi, we can easily obtain the expectation value of H as the weighted average

$\langle H \rangle = \sum h_i \langle H_i \rangle$

A second observation that we need in order to understand the variational quantum eigensolver is the fact that the problem of finding the eigenvalues can be reduced to a variational problem. In fact, suppose that we are given a hermitian matrix A on a finite dimensional Hilbert space with eigenvalues $\lambda_1, \lambda_2, \dots$. Let us order the eigenvalues (which are of course real) such that $\lambda_1 \leq \lambda_2 \leq \dots$, and let us assume that $|\psi_1 \rangle, |\psi_2 \rangle \dots$ is an orthonormal basis of corresponding eigenvectors.

If we are now given an arbitrary non-zero vector $|\psi \rangle$, we can of course expand this vector as a linear combination

$|\psi \rangle = \sum_i a_i |\psi_i \rangle$

which immediately gives us the following expression for the expectation value

$\langle A \rangle = \frac{1}{\langle \psi | \psi \rangle} \sum_i |a_i|^2 \lambda_i$

Given the ordering of the eigenvalues, we can now estimate this number as follows.

$\langle A \rangle = \frac{1}{\langle \psi | \psi \rangle} \sum_i |a_i|^2 \lambda_i \geq \frac{1}{\langle \psi | \psi \rangle} \sum_i |a_i|^2 \lambda_1 = \lambda_1$

It is also obvious that we have equality if $|\psi \rangle$ is actually an eigenvector of A with eigenvalue $\lambda_1$. Thus finding an eigenvector for the lowest eigenvalue of A is equivalent to minizing the expectation value of A!

This is at the heart of many variational methods like the Ritz method that have always been part of the toolset of computational quantum mechanics. In most of these methods, one considers state vectors parametrized by a finite dimensional parameter set, i.e. states of the form $|\psi(\Theta) \rangle$ where the parameter $\Theta$ ranges over some subset of a finite dimensional euclidian space. One then tries to minimize the expectation value

$\langle H \rangle (\Theta) = \frac{\langle \psi(\Theta) | H | \psi(\Theta) \rangle}{\langle \psi(\Theta) | \psi(\Theta) \rangle}$

using classical methods from mathematical optimization. The state vector minimizing this expectation value is then taken as an approximation for an eigenvector of H for its lowest eigenvalue.

Many classical optimiziation approaches that one might want to employ for this task work iteratively. We start with some value of the parameter $\Theta$, then determine the expectation value, adjust the parameter, determine the next expectation value and so forth. Unfortunately, calculating the expectation value of a matrix in a high dimensional Hilbert space is computationally very hard, which makes this algorithm difficult to apply to quantum systems with more than a few particles.

The algorithm

This is the point where the quantum variational eigensolver comes into play. If the operator H can – as assumed above – be decomposed into a sum of operators for which finding the expectation value can efficiently be done, then we can efficiently determine the expectation value of H as well and can use that in combination with a classical optimization algorithm to find an approximate eigenvector for H.

More precisely, here is the outline of the variational quantum eigensolver algorithm.

1. Decompose the operator H as a linear combination $H = \sum_\alpha h_\alpha H_\alpha$ where the expectation value of each $H_\alpha$ can be efficiently determined
2. Choose an initial value for the parameter $\Theta$
3. Prepare the quantum computer in the state $|\psi(\Theta) \rangle$
4. For each $H_\alpha$, use measurements to determine the expectation value $\langle H_\alpha \rangle$
5. Calculate the expectation value of H as the weighted average $\langle H \rangle = \sum_\alpha h_\alpha \langle H_\alpha \rangle$
6. Apply a classical optimization step to determine a new value for the parameter $\Theta$
7. Start over with step 3 and repeat

So the algorithm switches forth and back between a quantum part – preparing the state and finding the expectation values – and a classical part, i.e. performing the actual optimization, and is therefore a prime example for a hybrid approach. This is visualized in the diagram below (which is a variation of the diagram presented in [1]).

Of course, to fully specify this algorithm, we have to define a method how to prepare the states efficiently and need to pick a classical optimization algorithm. In [1], the so-called coupled cluster method is chosen to prepare a state. In this method, the prepared state is given by applying a (parametrized) unitary transformation to a fixed reference state, and this unitary transformation can be efficiently implemented on a quantum computer. As an optimization algorithm, the downhill simplex method, also known as the Nelder-Mead method, is employed.

As noted in [1], the algorithm makes a certain trade-off between the number of iterations required to estimate the expectation values and perform the optimization step and the time required to run a single iteration on the quantum device. One of the main problems that todays quantum computers face, however, is to keep a quantum state stable over a longer period of time. Therefore the bottleneck is the run time or the number of operations that we spend in the quantum part of the algorithm, which this approach tries to minimize. As explained in [4], the variational principle also remains valid if we consider quantum systems with noise, which is an indication that the algorithm is comparatively resistant to noise. All this makes the quantum variational eigensolver a good candidate for a true quantum advantage on a near-term, noisy and not fully fault-tolerant quantum device.

Applications in computational chemistry

So far we have implicitly dealt with Hamiltonians H on a Hilbert space representing a quantum register, i.e. a tensor product of one qubit systems. For applications, there is therefore one additional step that we need to carry out in order to be able to apply the algorithm – we need to map the system and the Hamiltonian in question to an n-qubit quantum system. This involves in particular a projection of the (typically) infinite dimensional Hilbert space of the original problem onto a finite dimensional system.

A prime example that is detailed in [4] and in [1] is the Hamiltonian describing an atom with N electrons and their interactions with the nucleus and other electrons. Finding the eigenstates and eigenvalues of this Hamiltonian is known as the electronic structure problem.

In general, this is a difficult problem. Of course, the easiest case of this – the hydrogen atom – can be solved analytically and is treated in most beginner courses in quantum mechanics. However, already for the Helium atom, no closed solution is known and approximations have to be used.

In order to apply the variational quantum eigensolver to this problem, a mapping to a Hilbert space describing a quantum register needs to be done. As explained in [1], this can be done using the mechanism of second quantization and the Jordan-Wigner transform. In this paper, the authors have applied the method to estimate the strength of the chemical bond in a helium hydride ion HeH+, which is formed by the reaction of a proton with a helium atom. They show that in order to obtain a very good approximation, a four-dimensional Hilbert space is sufficient, so that the algorithm can be carried out on a quantum computer with only two qubits. Recently ([5]), this has been extended to a six-qubit system used to estimate the ground state energy of a BeH2 molecule. Of course these are still comparatively small systems, but these experimental results encourage the hope that the approach scales and can be used on medium scale near-term quantum devices to go beyond what is possible with classical computers.

This work is also described in a post on the IBM website. IBM has made the Python sourcecode for similar applications available as part of its QISKit quantum computation library, so that everybody with a Q Experience account can run the algorithm and play with it. QISKit is a very interesting piece of software, and I will most likely devote one or more future posts on QISKit and similar frameworks like pyQUIL.

References

1. A. Peruzzo et. al., A variational eigenvalue solver on a photonic quantum processor, Nature Communications volume 5 (2014), available at www.nature.com/articles/ncomms5213
2. A. Kitaev, Quantum measurements and the Abelian Stabilizer Problem, arXiv:quant-ph/9511026
3. R. Cleve, A. Ekert, C. Macchiavello, M. Mosca, Quantum Algorithms Revisited, arXiv:quant-ph/9708016
4. J. R. McClean, J. Romero, R. Babbush, A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, arXiv:1509.04279
5. A. Kandala, A. Mezzacapo, K. Temme,
M. Takita, M. Brink, J. M. Chow, J. M. Gambetta, Hardware-efficient Variational Quantum Eigensolver for Small Molecules and Quantum Magnets, Nature volume 549, pages 242–246 (14 September 2017), available as arXiv:1704.05018v2

Into the quantum lab – first steps with IBMs Q experience

Even though physical implementations of quantum computers make considerable progress, it is not likely that you will have one of them under your desk in the next couple of years. Fortunately, some firms like IBM and Rigetti have decided to make some of their quantum devices available only so that you can play with them. In this post, we will have a first look at IBMs cloud offering – the IBM Q experience.

IBM offers different ways to interact with their quantum computers. The easiest access is provided by the interactive composer which allows you to create a quantum circuit online and to run it on a 5 qubit quantum computer based on superconducting qubits after registering.

The screenshot below shows the editor screen. At the top, you find an overview of the available quantum devices, both realising five (physical) qubits, and some of the main characteristics like error rates and decoherence times.

In the lower area of the screen, you can compose quantum circuits graphically. You can add quantum gates that manipulate the five qubits q[0] – q[4] and measurements. After setting up your circuit, you can either simulate the circuit or you can run it on the actual device. It is also possible to export the circuit as a QASM code snippet and to save experiments and results as so-called quantum scores.

Let us now run this circuit on the actual quantum hardware. Doing this consumes units, i.e. credits. A standard run consists of 1024 repetitions and one measurement after each repetition and consumes three units (which are replenished after the execution has completed). The execution is queued and, typically after a few minutes, the results are made available online (and you receive a mail informing you about the execution). When running a circuit, the platform also checks whether that circuit was executed before (by you our someone else) and if yes offers you immediate access to the previous results without consuming any units. The screenshot below shows the result of one of these runs.

The interesting part is the diagram at the top that shows the measurement outcomes in the computational basis, i.e. for each member of the computational basis on the x-axis, the value shows the percentage of measurements having that basis element as result. In this example, we see that almost always, the result is 00000, as expected from the circuit in which we do a double inversion of the first qubit. However, in some cases – corresponding to 1,3 percent of all runs – the outcome is 00001. This is a nice example for an error that occurs whenever we switch from a simulator to a real physical device. Note that in the IBM Q experience, the qubit q[0] is the least significant qubit and the rightmost qubit in the notation of the computational basis.

Now let us try something else – we want to put our system into an equal superposition

$\sum_i |i \rangle$

which is the usual starting point for most quantum algorithms. We know that we can achieve this by applying a Hadamard gate on each qubit. The following screenshot shows the resulting circuit for three qubits and the results of a (cached) run.

As expected, we see that our measurement results are spread across eight (23) vectors of the computational basis. However, we also see that again, the result does not exactly match the theoretical prediction – the values are not exactly equally distributed but we see slight deviations as before.

Finally, let us try out a multi-qubit gate. The screenshot below shows the result of a run with two CNOT gates. The first CNOT gate receives a logical one on its control bit, the second CNOT gate a logical zero. Thus, the result should be 00011, i.e. the first (least significant) qubits are inverted. In reality, we again get this result plus some noise, represented by unexpected outcomes in the measurements.

Playing with the composer is fun, and makes it easy to create simple circuits from scratch. However, more complicated circuits soon become difficult to define in this way. For those circuits, the IBM platform offers an API that can be used to run experiments using Python scripts that can much easier be developed and debugged step by step. We will look at this option in a later post. Until then, I recommend you get an account and create some circuits yourself – happy hacking!