Basic quantum circuit simulation in Python

I’ve always been a proponent of the idea that one of the best ways to learn about a topic is to code up a simple example that uses that idea/concept/algorithm. In conversations I’ve had with students recently, I’ve realized there is some interest in playing with quantum computing, quantum circuits, and quantum simulation without a full background in quantum physics. I thought a good way to help bridge this gap would be to provide some pedagogical examples of how one might go about simulating a quantum circuit. This won’t focus on computing performance or approaches to simulating large systems (for information on that topic, this paper is a good reference), but rather simple ways to simulate actions on a small number of qubits.

Brief background

A classical digital computer is built upon a foundation of bits, which take on values of 0 or 1, and store information. Clearly we also need to manipulate and work with that data, and the action on bits are performed by classical logical gates, such as the “AND” gate. In analogy to classical bits, in quantum computing, we have quantum bits, which take on values $\ket{0}$ and $\ket{1}$ . Similarly, we need to manipulate the data stored in quantum bits, and the equivalent operations are called quantum gates. Quantum bits have a number of nice properties absent in their classical counterparts, such as superposition and entanglement, and quantum gates are designed to make sure they can exploit these features, a property known as unitarity. For the purposes of this simple post, the exact details of these properties will be unimportant, as a correct mathematical representation of qubits and quantum gates often allows us to have these features naturally in our numerical simulations.

One Qubit

The state of a qubit may be represented as a column vector, and conventionally

(1) $\begin{align*} \ket{0} &= \left( \begin{array}{c} 1 \\ 0 \end{array} \right) \\ \ket{1} &= \left( \begin{array}{c} 0 \\ 1\end{array} \right) \end{align*}$

In Python, we can represent such states as arrays in Numpy with the correct dimensions as:

import numpy as np
import scipy as sp
import scipy.linalg

Zero = np.array([[1.0],
                 [0.0]])
One = np.array([[0.0],
                [1.0]])

which gives us a good representation of either the $\ket{0}$ (Zero) or $\ket{1}$ (One) state. Say we wanted to create a superposition of $\ket{0}$ and $\ket{1}$ , for example, the standard $\ket{+}$ state which is defined by $\ket{+}= 1/\sqrt{2} (\ket{0} + \ket{1})$ . The factor $1/\sqrt{2}$ is called a normalization factor, and helps ensure the relation to a probability distribution. If we wanted to create this state, one way of doing this in Python might be to define a normalization function, and input the raw state into it as

NormalizeState = lambda state: state / sp.linalg.norm(state)

Plus = NormalizeState(Zero + One)

Now that we know how to create single qubit states, we might want to be able to perform interesting actions on them. Quantum gates in acting on single qubits are given by two-by-two matrices. One interesting gate is the Hadamard gate, given by

(2) $\begin{align*} \frac{1}{\sqrt{2}} \left( \begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array} \right) \end{align*}$

we can represent this gate in Python and examine its action on the state $\ket{0}$ (that is, $H \ket{0} = ?$ ) through the following:

Hadamard = 1./np.sqrt(2) * np.array([[1, 1],
                                     [1, -1]])

NewState = np.dot(Hadamard, Zero)

where “np.dot” in Python serves the role of matrix multiply for Numpy arrays. Take a look at NewState, you should find it has a striking resemblance to state we just examined a few lines ago…

More than one qubit

While one qubit is interesting, more than one is certainly more interesting. How do we generalize what we just learned about representing one qubit to many qubits? States of many qubits sometimes have a few different notations, for example $\ket{001}$ or $\ket{0} \otimes \ket{0} \otimes \ket{1}$ , which are the same state. The symbol $\otimes$ is a source of confusion to many people, as rapid switches between its implicit and explicit use are rampant in quantum computing, and the symbol is essentially absent in some fields of physics, despite the operation being ubiquitous. The symbol $\otimes$ stands for the “tensor product” between two states and it allows for a correct description of many particle or many qubit states. That is, the space occupied by 2 qubits is the tensor product of the space of two single qubits, similarly the space of $N$ qubits is the N-fold tensor product of $N$ single qubit spaces. Time spent on the Wikipedia page for tensor product is prone to causing excessive head-scratching wondering what the “universal property” is, and how it’s relevant. Fear not, most of these abstract properties are not required reading for working with a practical example of a tensor product, which is what we will do here.

A practical example of a tensor product that we will use here, is the Kronecker product, which is implemented in Numpy, Matlab, and many other code packages. We will use the fact that it faithfully takes care of the important properties of the tensor product for us to represent and act on many qubit states. So how do we build a correct representation of the state $\ket{00}$ or $\ket{++}$ using it in Python? Using the Kronecker product, it’s as easy as

ZeroZero = np.kron(Zero, Zero)
OneOne = np.kron(One, One)
PlusPlus = np.kron(Plus, Plus)

What about an entangled Cat state such as $1/\sqrt{2} (\ket{00} + \ket{11})$ you ask?

CatState = NormalizeState(ZeroZero + OneOne)

And if we want to create a 5-qubit state, such as $\ket{10101}$ ? We can simply chain together Kronecker products. A simple helper function makes this look a bit less nasty

def NKron(*args):
  """Calculate a Kronecker product over a variable number of inputs"""
  result = np.array([[1.0]])
  for op in args:
    result = np.kron(result, op)
  return result

FiveQubitState = NKron(One, Zero, One, Zero, One)

Take note of how the size of the state grows as you add qubits, and be careful not to take down Python by asking it to make a 100-Qubit state in this way! This should help to give you an inkling of where quantum computers derive their advantages.

So how do we act on multi-qubit states with our familiar quantum gates like the Hadamard gate? The key observation for our implementation here is the implicit use of Identity operators and the Kronecker product. Say $H_0$ denotes a Hadamard gate on qubit 0. How do we represent $H_0$ ? The answer is, it depends, specifically on how many total qubits there are in the system. For example, if I want to act $H_0$ on $\ket{00}$ , then implicitly I am acting $H_0 \otimes I_1 \ket{00}$ and if I am acting on $\ket{10101}$ then implicitly I am performing $H_0 \otimes I_1 \otimes I_2 \otimes I_3 \otimes I_4\ket{10101}$ , and this is important for the numerical implementation. Lets see the action $H_0 \ket{10101}$ in practice,

Id = np.eye(2)
HadamardZeroOnFive = NKron(Hadamard, Id, Id, Id, Id)
NewState = np.dot(HadamardZeroOnFive, FiveQubitState)

Luckily, that’s basically all there is to single qubit gates acting on many-qubit states. All the work is in keeping track of “Id” placeholders to make sure your operators act on the same space as the qubits.

Many-qubit gates require a bit more care. In some cases, a decomposition amenable to a collection of a few qubits is available, which makes it easy to define on an arbitrary set of qubits. In other cases, it’s easier to define the operation on a fixed set of qubits, then swap the other qubits into those slots, perform the operation, and swap back. Without covering every edge case, we will use what is perhaps the most common 2-qubit gate, the CNOT gate. The action of this gate is to enact a NOT operation on a qubit, conditional on the state of another qubit. The quantum NOT operation is given by the Pauli X operator, defined by

(3) $\begin{align*} X = \left( \begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array} \right) \end{align*}$

In order to perform this action dependent on the state of another qubit, we need a sort of quantum IF-statement. We can build this from projectors of the form $P_0 = \ket{0}\bra{0}$ and $P_1 = \ket{1}\bra{1}$ that tell us if a qubit is in state $0$ or $1$ . In our representation, the so-called “bra” $\bra{0}$ is just the conjugate transpose of the “ket” $\ket{0}$ , so these projectors can be written as an outer product

(4) $\begin{align*} P_0 = \left( \begin{array}{c} 1 \\ 0 \end{array} \right) \left( \begin{array}{cc} 1 & 0 \end{array} \right) \\ P_1 = \left( \begin{array}{c} 0 \\ 1 \end{array} \right) \left( \begin{array}{cc} 0 & 1 \end{array} \right) \\ \end{align*}$

With this, we can define a CNOT gate on qubits $0$ and $1$ as $CNOT=P_0 \otimes I + P_1 \otimes X$ . So if we wanted to perform a CNOT operation, controlled on qubit 0, acting on qubit 3 of our 5 qubit state, we could do so in Python as

P0 = np.dot(Zero, Zero.T)
P1 = np.dot(One, One.T)
X = np.array([[0,1],
              [1,0]])

CNOT03 = NKron(P0, Id, Id, Id, Id) + NKron(P1, Id, Id, X, Id)
NewState = np.dot(CNOT03, FiveQubitState)

This covers almost all the requirements for simulating the action of quantum circuit of your choosing, assuming you have the computational power to do so! So long as you can reduce the circuit you are interested in into an appropriate set of 1- and 2-qubit gates. Doing this in the generic case is still a bit of a research topic, but many algorithms can already be found in this form, so it should be reasonably easy to find a few to play with!

Measurement

Now that we’ve talked about how to simulate the action of a circuit on a set of qubits, it’s important to think about how to actually get information out! This is the topic of measurement of qubits! At a basic level, we will want to be able to ask if a particular qubit is 0 or 1, and simulate the resulting state of the system after receiving the information on that measurement! Since qubits can be in superpositions between 0 and 1 (remember the state $\ket{+}$ ), when we ask that qubit for its value, we get 0 or 1 with some probability. The probability of qubit $i$ being in a state $0$ is given by the fraction of the state matching that, which can be concisely expressed as

(5) $\begin{align*} \text{Tr}\left[\ket{\Psi}\bra{\Psi} P_0^i \right] \end{align*}$

where $\ket{\Psi}$ is some generic state of qubits and $P_0$ is the same projector we had above, acting on the $i$ ‘th qubit. If we measure a 0, then the resulting state of the system is given by

(6) $\begin{align*} P_0^i \ket{\Psi} / || P_0^i \ket{\Psi}||_2 \end{align*}$

recalling that anywhere where the number of qubits don’t match up, we are following the generic rule of inserting identity operators to make sure everything works out! In Python, we could simulate measurements on a Cat state using the above rules as follows

import numpy.random

CatState = NormalizeState(ZeroZero + OneOne)
RhoCatState = np.dot(CatState, CatState.T)

#Find probability of measuring 0 on qubit 0
Prob0 = np.trace(np.dot(NKron(P0, Id), RhoCatState))

#Simulate measurement of qubit 0
if (np.random.rand() < Prob0):
    #Measured 0 on Qubit 0
    Result = 0
    ResultState = NormalizeState(np.dot(NKron(P0, Id), CatState))
else:
    #Measured 1 on Qubit 1
    Result = 1
    ResultState = NormalizeState(np.dot(NKron(P1, Id), CatState))

print "Qubit 0 Measurement Result: {}".format(Result)
print "Post-Measurement State:"
print ResultState

If you look you will see that upon measurement the state collapses to either the state $\ket{00}$ or $\ket{11}$ with equal probability, meaning that a subsequent measurement of the second qubit will yield a deterministic result! This is a signature of entanglement between the two qubits!

Summary

So in this short post we’ve covered how to represent qubits, quantum gates and their actions, and measurement of qubits in Python by using a few simple features in Numpy/Scipy. These simple actions are the foundation for numerical simulation of quantum circuits and hopefully gives you some idea of how these simulations work in general. If you have any questions or want me to elaborate on a specific point, feel free to leave a comment and I will update the post or make a new one when I get the chance!