Quantum Simulator 101
There’s a lot of talk surrounding quantum computers. Especially since D-Wave announced that they had produced a quantum computer1. There was a lot of controversy as to whether their machine actually counted as a quanutm computer.2 As their machine was not a general purpose quantum computer, it could only preform a few specialized tasks, but it was directly using quantum mechanical effects to produce calculations. In the years since no one has managed to create a general-purpose quantum computer. The news is usually just that a larger number was factored with Shor’s algorithm.
Discussions about quantum computers always come with comments about how they will break public-key encryption, one of the fundamental technologies employed by the internet. Despite this, a lot of the math behind quantum computing is quite simple. If you understand complex numbers, vectors and matrices then you’re almost there. The goal of this series is to show you how to use that knowledge to build a simulator of a quantum computer.
Theory is always ahead of practice. Reality is bogged down by practical concerns while theory speeds ahead unfettered. In the 1800s Lady Ada Lovelace wrote programs for a computer that has still yet to be built.3 In the 1960’s particle physicists realized that the math behind their Eightfold Way model of physics4 was well known to mathematicians by the name of group theory.5 I’ve head that math is at least 100 years ahead of physics at this point. Quanutm computers are still mostly in the theory stage which makes getting access to them difficult. Though IBM will give to access to actual quantum hardware if you sign up6 but they have issues running long quantum algorithms. Right now simulators are your best bet.
I’m going to begin with a brief section on vector spaces and basis vectors. If you’re already familiar with these then go ahead and skip it.
Vector Spaces and Basis Vectors
Two big properties of vectors is that you can add them together and multiply them by scalars. These two concepts form the basic idea behind a vector space. A vector space is a set of vectors that you can add together and multiply by scalars to get another vector in the set. This is an important property and sets of vectors that have it are said to be closed.
The xy-plane is an example of a vector space. The vectors are those that point from the origin to each point. If you add two vectors togther you’ll get a vector to another point in the plane and multiplying these vectors by scalars also gets you another vector to a point in the plane. This plane is a 2-dimensional vector space because you need two numbers to refer to a specific vector. There’s another way to think of this, with basis vectors.
You can write the vectors in the plane with an x coordinate and a y coordinate like
$$ \vec{v} = \begin{bmatrix}4\\7\end{bmatrix} $$
You can also write them like this \( \vec{v} = 4\hat{x} +7\hat{y} \), where \( \hat{x} \) and \( \hat{y} \) are unit vectors that point solely in the x and y directions. Every vector can be written in terms of \( \hat{x} \) and \( \hat{y} \) and they are called basis vectors. Since there are two basis vectors, it’s a 2-dimensional vector space. Basis vectors form the heart of any vector space.
Qubits
Classical computers operate on bits, quantum computers operate on qubits. A bit can only be in one state at any time and reading its value will always return the same result. A qubit can be in a superposition of states. When you measure the state of a qubit, you can get different values that occur with different probabilities. For example, you can have a qubit that has a 25% chance of measuring as a zero and a 75% chance of measuring as a one. However, qubits don’t have to be in a superposition. In fact, you can simulate a classical computer on a quantum computer just by using qubits that have a 100% chance to be a zero or a 100% chance to be a one. Representation
A qubit can be written as a linear combination of two basis vectors: \( \ket{0} \) and \( \ket{1} \) (\( \ket{} \) is called a ket). The notation might look scary, but they’re just vectors.
$$ \ket{0} = \begin{bmatrix}1\\0\end{bmatrix} \quad \ket{1} = \begin{bmatrix}0\\1\end{bmatrix} $$
A qubit is a linear combination of these two basis vectors
$$ \ket{N} = \alpha\ket{0} + \beta\ket{1} $$
\( \alpha \) and \( \beta \) are complex numbers and are the probability amplitudes for measuring the qubit as a zero and one. A probability is not the same as a probability amplitude. A probability is a positive real number and a probability amplitude is a complex number. It can be negative or even purely imaginary. The corresponding probabilities for \( \alpha \) and \( \beta \) are given by \( \lvert\alpha\rvert^2 \) and \( \lvert\beta\rvert^2 \). Where \( \lvert x \rvert = xx^* \) for complex numbers. Since probabilities must sum to one, we have \( \lvert\alpha\rvert^2 + \lvert\beta\rvert^2 = 1 \).
Probability amplitudes can take a a little getting used to. For example,
$$ \frac{1}{\sqrt{2}}\ket{0} + \frac{1}{\sqrt{2}}\ket{1} $$
is the one qubit state that has a 50-50 chance of a zero and one.
Basis Vectors and the Tensor Product
The state of N qubits is an 2N dimensional vector and you have 2N basis vectors. For two qubits you have
$$ \ket{00} = \begin{bmatrix}1\\0\\0\\0\end{bmatrix} \quad \ket{01} = \begin{bmatrix}0\\1\\0\\0\end{bmatrix} \quad \ket{10} = \begin{bmatrix}0\\0\\1\\0\end{bmatrix} \quad \ket{11} = \begin{bmatrix}0\\0\\0\\1\end{bmatrix} \quad $$
There’s a very handy way to generate the basis vectors for any number of qubits using an operation called the tensor product. It’s incredibly useful for qubits because
$$ \ket{10} = \ket{1} \otimes \ket{0} $$
you can keep applying the tensor product with \( \ket{0} \) and \( \ket{1} \) to get any basis vector you want.
This gives us the first hint of why simulating quantum computers on classical computers is difficult: it takes a lot of RAM. To simulate 32 qubits you need a vector with 232 entries. Each entry is a complex number and takes 8 bytes if we’re using 32 bit floats. This comes out to 32 GB of RAM. If you want to break 4096-bit RSA Shor’s algorithm requires 3 * 4096 = 12288 qubits. Storing the state of those qubits would require exponentially more space than a 64 bit address space can provide.
The Code
This post will cover the code to implement the tensor product and basis vectors for any state. This will be written in C++ using the Eigen7 library. The code can also be found in this git repo.
Let’s start with some types, a function to calculate powers of two, and a function to print binary numbers.
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/KroneckerProduct>
typedef Eigen::Matrix<float, Eigen::Dynamic, 1> Vector;
typedef Eigen::SparseMatrix<float> Matrix;
constexpr int pow2(int exp) {
return (exp == 0) ? 1 : 2 * pow2(exp - 1);
}
void printBin(int N, int BITS) {
for(int i = BITS - 1 ; i >= 0; i--) {
if((N >> i) & 0x1) {
putchar('1');
}
else {
putchar('0');
}
}
}
We’re using dense and sparse matrices here. Sparse matrices are matrices that contain mostly zeros for their entries. Because of that, they can be stored and perform operations more efficiently. Our Vector type will mostly be used to store quantum states, which tend not to be sparse. Dense matrices in eigen have three important template arguments: the type of the matrix entries, the number of rows at compile time and the number of columns at compile time. Passing Eigen::Dynamic as an argument allows for a matrix with a dynamic number of rows/columns.
We include the Eigen/KroneckerProduct header. The Kronecker product is a generalization of the tensor product that covers matrices as well as vectors. It’s an unsupported Eigen module and might be located in a separate directory from the other Eigen headers so make sure your compiler can find it.
Next we have functions to product basis vectors
//return |0>
Vector zero(void) {
Vector result(2);
result << 1.0, 0.0;
return result;
}
//return |1>
Vector one(void) {
Vector<T> result(2);
result << 0.0, 1.0;
return result;
}
//returns |N> for a space of size NQUBITS
Vector basis(int N, int NQUBITS) {
Vector basis(1);
int highBit;
int remaining = N;
basis << 1.0; //identity element for the tensor product
for(int i = NQUBITS - 1; i >= 0; i--) {
highBit = remaining >> i; //get the highest bit
remaining = ~(1 << i) & remaining; //get the remaining bits
if(highBit) {
basis = kroneckerProduct(basis, one()).eval();
}
else {
basis = kroneckerProduct(basis, zero()).eval();
}
}
return basis;
}
Lastly, we have some code to hold a quantum state.
class State {
private:
Vector state;
int NQUBITS;
public:
State(int NQUBITS, const Vector& state) : state(state), NQUBITS(NQUBITS) {
assert(state.rows() == pow2(NQUBITS));
}
~State(void) {
}
void showOutcomes(void) {
Vector probabilities = this->state.array() * this->state.conjugate().array();
printf("State | Percentage\n");
for(int i = 0; i < this->state.rows(); i++) {
printBin(i, this->NQUBITS);
printf(" | %.2f\n", (float) probabilities[i].real() * 100.0);
}
}
}
The showOutcomes function will print each state and the probability of measuring it.
Time to test our code
int main(int argc, char **argv) {
float norm = 1.0f / sqrt(2.0f);
State state(3, norm * basis(0, 3) + norm * basis(5, 3));
state.showOutcomes();
return 0;
}
That code will give you a result of
State | Percentage
000 | 50.00
001 | 0.00
010 | 0.00
011 | 0.00
100 | 0.00
101 | 50.00
110 | 0.00
111 | 0.00
which is
$$ \frac{1}{\sqrt{2}}\ket{000} + \frac{1}{\sqrt{2}}\ket{101} $$
as expected. The next post will focus on single qubit quantum gates and the hadamard gate in particular.