Parallel Computing in the Small

As I watched my one-month old baby sleeping, a thought ran through my mind. I wonder what technology would look like when she’s my age. There’s a lot of things going on simultaneously in technology today. I remember doing Artificial Intelligence/Machine Learning about 15 years ago. They used to be done in a university setting. They are now the new normal. Some technologies are quite recent like Big Data and Blockchain but they have become mainstream very quickly. I wonder what it will look like when my daughter grows up.

There is a technology I’m currently watching. It’s still in it’s infancy but this technology is really fascinating. It is a technology based on a remarkably counter-intuitive theory of the universe that governs the atoms and sub-atomic particles. It is called Quantum Computing. Given the high speed at which science fiction become realities, it’s possible that my baby girl will someday be programming on a quantum computer. But for now, we can only hope.

I have not seen a Quantum Computer and based on what I read so far, it’s still an engineering challenge. It might probably be a few more years before they are mass-produced but it is something to look forward to. But the good thing is, we can already talk about quantum algorithms! So let’s sharpen our pencils and get a few sheets of paper for we are going to talk about this new way of computing.

In classical computing, a bit has a value that is either one or zero. In quantum computing, a bit can now have a value that is a superposition of 1 and 0. We call this a QBit. The values 1 and 0 are now symbolized as \left| 1\right\rangle and \left| 0\right\rangle respectively. These are unit vectors in what is known as a Hilbert space. They are also known as kets, derived from the word bracket which we’ll talk more of later. An example of a vector space that you are familiar with is the color wheel. Each color is actually a combination of red, blue and green. For example, the color yellow is a combination of red and green. It’s also the same as with the qubit. A general QBit vector is of the form

\displaystyle \left| \Psi \right\rangle = a\left| 0\right\rangle + b\left| 1\right\rangle

such that \mid a\mid^2 + \mid b\mid^2 = 1. The numbers a and b are complex numbers and the symbol \mid a\mid^2 is defined as

\displaystyle \mid a\mid^2 = a^*a

where a^* is the complex conjugate of a. A QBit vector represents the state of the QBit. We shall use the term vector and state interchangeably moving forward.

However, you won’t be able to see them in superposition state. If you want to know the state of a qubit, you will have to measure it. But when you measure it, it will collapse to one of either \left| 0\right\rangle or \left| 0\right\rangle but you will have a clue as to what that state will be. The coefficients a and b will give you the probability of the result of the measurement: it will be \mid a\mid^2 in state \left| a\right\rangle and \mid b\mid^2 in state \left| b\right\rangle.

If the state of a QBit is a superposition of basis states \left| 0\right\rangle and \left| 1\right\rangle, does this mean our computation is also random? Does it mean we get different answers every time ?

It turns out that we can get a definite answer! Let’s illustrate this by a simple quantum computation.

Suppose a I have a number between 0 and 1 million. What is my number? You can guess my number and I will only tell you if you got it correctly or not. I won’t tell you if my number is higher or lower than your guess. How many tries do you think you need in order to guess my number? I’m sure our tries will be many. In the worst case, you’ll need about 1 million tries.

For a quantum computer, you only need one try to guess it and that’s what we’re going to see. However, bear with me since we have to build up our arsenal of tools first before we can even talk about it. Let’s start with how a QBit looks like as a matrix.

The basis QBits \left| 0\right\rangle and \left| 1\right\rangle can be written as 1×2 matrices:

\displaystyle \left| 0 \right\rangle =  \begin{bmatrix}  1\\  0  \end{bmatrix}

\displaystyle \left| 1 \right\rangle =  \begin{bmatrix}  0\\  1  \end{bmatrix}

Using this representation, we can write a general QBit as:

\displaystyle \mid\Psi\rangle = a\mid 0\rangle + b\mid 1\rangle = a  \begin{bmatrix}  1\\  0  \end{bmatrix}  + b  \begin{bmatrix}  0\\  1  \end{bmatrix}  = \begin{bmatrix}  a\\  b  \end{bmatrix}

So now we know that every QBit is a linear combination of ket vectors \left| 0\right\rangle and \left| 1\right\rangle. These ket vectors have a corresponding “bra” vectors \langle 0 \mid and \langle 1\mid. In matrix representation, the bra basis vectors are defined as:

\displaystyle \langle 0\mid  = \begin{bmatrix}  1 & 0  \end{bmatrix}

and

\displaystyle \langle 1\mid =  \begin{bmatrix}  0 & 1  \end{bmatrix}

For a general QBit \mid\Psi\rangle = a\mid 0\rangle + b\mid 1\rangle, the corresponding bra vector is

\displaystyle  \langle \Psi \mid = a^*\langle 0\mid + b^*\langle 1\mid =  \begin{bmatrix}  a^* & b^*  \end{bmatrix}

There is an operation between bra and ket vectors of the basis states called inner product which is defined as:

\begin{array}{rl}  \displaystyle \langle 0\mid 0 \rangle &= 1\\  \displaystyle \langle 1\mid 1 \rangle &= 1\\  \displaystyle \langle 0\mid 1 \rangle &= \langle 1\mid 0\rangle = 0  \end{array}

Using the above rules, we can define the inner product of any two QBits \mid \Psi\rangle = a\mid 0\rangle + b\mid 1\rangle and \mid\Phi\rangle = c\mid 0\rangle + d\mid 1\rangle as

\begin{array}{rl}  \displaystyle \langle \Psi\mid\Phi\rangle &= \big(a^*\langle 0\mid + b^*\langle 1\mid\big)\big(c\mid 0\rangle + d\mid 1\rangle\big)\\  &= \displaystyle a^*c \langle 0\mid 0\rangle + a^*\langle 0\mid 1\rangle + b^*c\langle 1\mid 0\rangle + b^*d\langle 1\mid 1\rangle\\  &= a^*c + b^*d  \end{array}

Using the above definition of an inner product, the norm of a vector is defined as

\displaystyle \mid\left|\Psi\right\rangle\mid = \displaystyle \sqrt{\langle \Psi \mid \Psi \rangle} =\displaystyle \sqrt{a^*a + b^*b}

Geometrically, the norm is the length of the ket \mid \Psi\rangle.

Quantum computations are done using quantum operators. They are also called quantum gates. Operators are linear mappings that take a state to another state. That is, an operator \mathbf{A} acts on a\left| 0 \right\rangle + b\left| 1\right\rangle to give another vector c\left| 0 \right\rangle + d\left| 1\right\rangle. We write this as

\displaystyle \mathbf{A}(a\left| 0\right\rangle + b\left| 1\right\rangle) = c\left| 0\right\rangle + d\left| 1\right\rangle

Since the basis vectors \left| 0\right\rangle and \left| 1\right\rangle are also vectors, operators can also operate on them. An example of a very important operator and it’s action on the basis vectors is the Hadamard operator:

\begin{array}{rl}  \displaystyle \mathbf{H} \left| 0 \right\rangle &= \displaystyle \frac{1}{\sqrt{2}}(\left| 0 \right\rangle + \left| 1\right\rangle)\\  \displaystyle \mathbf{H} \left| 1 \right\rangle &= \displaystyle \frac{1}{\sqrt{2}}(\left| 0 \right\rangle - \left| 1\right\rangle)  \end{array}

By knowing the action of an operator on the basis vectors, we can determine the matrix representation of an operator. For the Hadamard operator, the matrix representations are:

\displaystyle \mathbf{H} \left| 0 \right\rangle = \displaystyle \frac{1}{\sqrt{2}}(\left| 0 \right\rangle + \left| 1\right\rangle) = \frac{1}{\sqrt{2}}  \begin{bmatrix}  1\\  1  \end{bmatrix}

and

\displaystyle \mathbf{H} \left| 1 \right\rangle = \displaystyle \frac{1}{\sqrt{2}}(\left| 0 \right\rangle - \left| 1\right\rangle) = \frac{1}{\sqrt{2}}  \begin{bmatrix}  1\\  -1  \end{bmatrix}

Therefore the matrix representation of the Hadamard operator is

\displaystyle \mathbf{H}=  \frac{1}{\sqrt{2}}  \begin{bmatrix}  1 & 1 \\  1 & -1  \end{bmatrix}

The identity operator is defined as the operator that maps a vector into itself and is symbolized as \mathbf{I}:

\displaystyle \mathbf{I}\left|\Psi\right\rangle = \left|\Psi\right\rangle

The NOT operator is defined as

\begin{array}{rl}  \mathbf{X}\left|0\right\rangle &= \left|1\right\rangle\\  \mathbf{X}\left|1\right\rangle &= \left|0\right\rangle  \end{array}

What it does is transform the \left|0\right\rangle to \left|1\right\rangle and vice versa, just like the classical bits.

Operators are linear. We can use linearity to derive the resulting vector when an operator acts on a general QBit:

\displaystyle \mathbf{A}(a\left| 0\right\rangle + b\left| 1\right\rangle) = a\mathbf{A} \left| 0\right\rangle + b\mathbf{A} \left| 1\right\rangle)

If \mathbf{A} is the Hadamard operator, that is, \mathbf{A} = \mathbf{H}, then the above operation becomes:

\displaystyle \mathbf{H}(a\left| 0\right\rangle + b\left| 1\right\rangle) = a\mathbf{H} \left| 0\right\rangle + b\mathbf{H} \left| 1\right\rangle)

Since we know that:

\begin{array}{rl}  \displaystyle \mathbf{H} \left| 0 \right\rangle &= \displaystyle \frac{1}{\sqrt{2}}(\left| 0 \right\rangle + \left| 1\right\rangle)\\  \displaystyle \mathbf{H} \left| 1 \right\rangle &= \displaystyle \frac{1}{\sqrt{2}}(\left| 0 \right\rangle - \left| 1\right\rangle)  \end{array}

We can then write this as:

\begin{array}{rl}  \displaystyle \mathbf{H}(a\left| 0\right\rangle + b\left| 1\right\rangle) &= \displaystyle a\mathbf{H} \left| 0\right\rangle + b\mathbf{H} \left| 1\right\rangle)\\  &= \displaystyle a\frac{1}{\sqrt{2}}\big(\left| 0\right\rangle + \left| 1\right\rangle\big) + b\frac{1}{\sqrt{2}}\big(\left| 0\right\rangle - \left| 1\right\rangle\big)\\  &= \displaystyle \frac{a+b}{\sqrt{2}} \left| 0\right\rangle + \frac{a-b}{\sqrt{2}} \left| 1\right\rangle  \end{array}

Let see what the operators look like when they operate on bra vectors. Let’s start with the Hadamard operator acting on a general QBit above.

\begin{array}{rl}  \displaystyle \mathbf{H} \mid 0 \rangle &= \displaystyle \frac{1}{\sqrt{2}}(\mid 0 \rangle + \mid 1\rangle) \rightarrow \langle 0\mid \mathbf{H^\dagger} = \displaystyle \frac{1}{\sqrt{2}}(\langle 0 \mid + \langle 1\mid)\\  \displaystyle \mathbf{H} \mid 1 \rangle &= \displaystyle \frac{1}{\sqrt{2}}(\mid 0 \rangle - \mid 1\rangle) \rightarrow \langle 1\mid \mathbf{H^\dagger} = \displaystyle \frac{1}{\sqrt{2}}(\langle 0 \mid - \langle 1\mid)  \end{array}

where we symbolize \mathbf{H^\dagger} as the operator corresponding to the Hadamard operator acting on bras.

The matrix representation of the Hadamard operator acting on bras is therefore:

\mathbf{H^\dagger} = \frac{1}{\sqrt{2}}  \begin{bmatrix}  1 & 1\\  1 & -1  \end{bmatrix} = \mathbf{H}

In general, the matrix representation of \mathbf{H^\dagger} is the conjugate transpose of the matrix of \mathbf{H}. In simple terms, this means take the transpose of \mathbf{H} and apply the complex conjugates to each element. It’s amazing that the Hadamard operator is equal to its conjugate transpose. The types of operators that behave that way are called Hermitian operators.

Let \mathbf{A} be an operator and \left| \Psi\right\rangle a vector and let \Phi = \mathbf{A}\left|\Psi\right\rangle, then the norm of \Phi is

\displaystyle \langle\Phi\mid\Phi\rangle = \langle \Psi\mid \mathbf{A^\dagger}\mathbf{A}\mid \Psi\rangle = \langle\Psi\mid\Psi\rangle

This means that

\displaystyle \mathbf{A^\dagger}\mathbf{A} = \mathbf{I}

The kind of operator that satisfies the above relationship is called Unitary.

Now that we have gained knowledge of some of the tools, we can now describe how a quantum computation is set up. We need QBits for the input and QBits for the output. The simplest setup is a one QBit input and one QBit output. Let \left| x\right\rangle and \left| y\right\rangle be the input and output QBits,respectively. The general state of this quantum system is a superposition of the following states:

\left| 0 \right\rangle \otimes \left| 0 \right\rangle, \left| 0 \right\rangle \otimes \left| 1 \right\rangle, \left| 1 \right\rangle \otimes \left| 0 \right\rangle, \left| 1 \right\rangle \otimes \left| 1 \right\rangle

which could also be written as

\left|00\right\rangle, \left|01\right\rangle, \left|10\right\rangle, \left|11\right\rangle,

or, if we take them as binary representation of numbers, we can write them concisely as

\left|0\right\rangle, \left|1\right\rangle,\left|2\right\rangle,\left|3\right\rangle,

Using the above basis vectors, we can write a general 2-QBit system as:

\displaystyle \left|\Psi\right\rangle = c_1\left|0\right\rangle + c_1\left|1\right\rangle + c_1\left|2\right\rangle + c_1\left|3\right\rangle = \sum_{i=0}^{2^n-1} c_i\left|i\right\rangle

where the numbers c_i are complex and n=2. The summation is up to 2^n-1 since the maximum number n QBits can represent is 2^n-1. The summation on the right also gives us the general form of a QBit for an n-QBit system.

As we saw earlier, the quantum operators (or gates) are the ones used to compute from a given QBit inputs. Let \mathbf{U}_f represent a unitary operator implementing the function f. The function f is a function that takes the numbers from the domain \{0,1\} to \{0,1\}, that is,

\displaystyle f: \{0,1\}\rightarrow \{0,1\}

For a two-QBit system, the form of quantum computation is

\displaystyle \mathbf{U}_f\left|x\right\rangle \left|y\right\rangle = \left|x\right\rangle \left|x\oplus f(x)\right\rangle

where the symbol \oplus represent bitwise OR operation. The table below will give you the results of the bitwise OR operation for every combination of values in \{0,1\}:

\begin{array}{cc}  0 \oplus 0 &= 0\\  0 \oplus 1 &= 1\\  1 \oplus 0 &= 1\\  1 \oplus 1 &= 0  \end{array}

A simple but important example of a unitary operator acting on two QBits is the \mathbf{C_{\text{NOT}}} gate, defined as:

\displaystyle \mathbf{C_{\text{NOT}}}\left|x\right\rangle\left|y\right\rangle = \left|x\right\rangle\left|y \oplus x\right\rangle

Here is the complete action of the \mathbf{C_{\text{NOT}}}:

\begin{tabular}{c|c|c}  \text{Input} & \text{Output} & \text{Result}\\  0 & 0 & 0\\  0 & 1 & 1\\  1 & 0 & 1\\  1 & 1 & 0  \end{tabular}

If you study the table above, you’ll notice that if the input QBit is in the state \left|0\right\rangle, the output QBit stays the same. However, if the input QBit is in \left|1\right\rangle, the output QBit is flipped. The input QBit in this case controls the flipping of the second QBit. The input QBit is said to be the control bit. That’s why it’s called a \mathbf{C_{\text{NOT}}} which stands for Controlled-NOT gate.

Let’s now go back to our original problem. There is a positive number a that is less than 2^n . We want to find this number using quantum computation. We are given a black box operator \mathbf{U}_f whose action on a bra is

\displaystyle \mathbf{U}_f \left|x\right\rangle\left|y\right\rangle = \left|x\right\rangle\left| y \oplus f(x)\right\rangle

where

\displaystyle f(x) = x_0 a_0 \oplus x_1 a_1 \oplus \cdots x_{n-1}a_{n-1} = x \cdot a

where a is the unknown number we are going to guess and x_ia_i are the i-th bits of x and a.

So how many times do we have to apply \mathbf{U}_f in order to guess the number a?

To solve this, we need n input QBits and 1 output QBit. We will initialize our QBits to \left|0\right\rangle. So the initial configuration is

\displaystyle \underbrace{\left|000...0\right\rangle}_{\text{n number of zeroes}}\otimes\left|1\right\rangle

Notice how we separate the input and output qubits.

Next we apply the Hadamard operator to each QBit:

\displaystyle \mathbf{H}^{\otimes n}\otimes \mathbf{H} \left|000...0\right\rangle\left|1\right\rangle = \mathbf{H}^{\otimes n} \left|000...0\right\rangle\otimes \mathbf{H}\left|1\right\rangle

We have seen the Hadamard operator earlier. For the sake of convenience we’ll repeat it here:

\begin{array}{rl}  \displaystyle \mathbf{H} \left| 0 \right\rangle &= \displaystyle \frac{1}{\sqrt{2}}(\left| 0 \right\rangle + \left| 1 \right\rangle)\\  \displaystyle \mathbf{H} \left| 1 \right\rangle &= \displaystyle \frac{1}{\sqrt{2}}(\left| 0 \right\rangle - \left| 1\right\rangle)  \end{array}

We can generalize this by letting x be either 0 or 1:

\begin{array}{rl}  \displaystyle \mathbf{H}\left| x\right\rangle &= \frac{1}{\sqrt{2}}\big(\left|0\right\rangle + (-1)^x \left|1\right\rangle  \end{array}

We can use this to compute the action on two QBits and generalize:

\begin{array}{rl}  \displaystyle \mathbf{H}^{\otimes 2}\left| x\right\rangle &= \mathbf{H}^{\otimes 2}\left| x_0x_1\right\rangle = \mathbf{H}\left| x_0\right\rangle \otimes \mathbf{H}\left| x_1\right\rangle\\  &= \displaystyle \frac{1}{\sqrt{2}}\big(\left|0\right\rangle + (-1)^{x_0} \left|1\right\rangle \otimes \frac{1}{\sqrt{2}}\big(\left|0\right\rangle + (-1)^{x_1} \left|1\right\rangle\\  &= \displaystyle \big(\frac{1}{\sqrt{2}}\big)^2 \big( \left|00\right\rangle + (-1)^{x_0} \left|01\right\rangle + (-1)^{x_1}\left|10\right\rangle + (-1)^{x_1 + x_0}\left|11\right\rangle\\  &= \displaystyle \big(\frac{1}{\sqrt{2}}\big)^2 \big( (-1)^{0\cdot x_1 + 0\cdot x_0} \left|00\right\rangle + (-1)^{0\cdot x_1 + 1\cdot x_0} \left|01\right\rangle + (-1)^{1\cdot x_1 + 0\cdot x_0}\left|10\right\rangle + (-1)^{1\cdot x_1 + 1\cdot x_0}\left|11\right\rangle\\  &= \displaystyle \big(\frac{1}{\sqrt{2}}\big)^2 \big( (-1)^{0\cdot x_1 + 0\cdot x_0} \left|0\right\rangle + (-1)^{0\cdot x_1 + 1\cdot x_0} \left|1\right\rangle + (-1)^{1\cdot x_1 + 0\cdot x_0}\left|2\right\rangle + (-1)^{1\cdot x_1 + 1\cdot x_0}\left|3\right\rangle\\  &= \displaystyle \big(\frac{1}{\sqrt{2}}\big)^2 \sum_{y=0}^{2^2-1} (-1)^{x\cdot y}\left|y\right\rangle  \end{array}

In general,

\displaystyle \mathbf{H}^{\otimes n}\left| x\right\rangle = \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1)^{x\cdot y}\left|y\right\rangle

Using the above formula, the action of \mathbf{H}^{\otimes n} on \left|0\right\rangle is

\displaystyle \mathbf{H}^{\otimes n} \left|0\right\rangle = \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} \left|y\right\rangle

The next step is to apply our operator \mathbf{U}_{f} which was defined above:

\begin{array}{rl}  \displaystyle \mathbf{U}_f\mathbf{H}^{\otimes n}\left| 0\right\rangle \otimes \frac{1}{\sqrt{2}} \big( \left|0\right\rangle - \left|1\right\rangle\big) &= \displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} \mathbf{U}_f\left|y\right\rangle \otimes \frac{1}{\sqrt{2}} \big( \left|0\right\rangle - \left|1\right\rangle\big)\\  &= \displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} \frac{1}{\sqrt{2}} \big( \mathbf{U}_f\left|y\right\rangle \otimes  \left|0\right\rangle - \mathbf{U}_f\left|y\right\rangle \otimes \left|1\right\rangle\big)\\  &= \displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} \frac{1}{\sqrt{2}} \big( \left|y\right\rangle \otimes  \left|0\oplus f(y)\right\rangle - \left|y\right\rangle \otimes \left|1\oplus f(y)\right\rangle\big)\\  \end{array}

Let’s break this down a bit. Since the value of f(y) is either 0 or 1, we have 2 cases:

When f(y) = 0,

\displaystyle \left|y\right\rangle \otimes  \left|0\oplus f(y)\right\rangle = \left|y\right\rangle \otimes  \left|0\right\rangle
\displaystyle \left|y\right\rangle \otimes  \left|1\oplus f(y)\right\rangle = \left|y\right\rangle \otimes  \left|1\right\rangle

which means

\displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} \frac{1}{\sqrt{2}} \big( \left|y\right\rangle \otimes  \left|0\oplus f(y)\right\rangle - \left|y\right\rangle \otimes \left|1\oplus f(y)\right\rangle\big) = \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} \left|y\right\rangle \otimes \frac{1}{\sqrt{2}} \big(  \left|0\right\rangle - \left|1\right\rangle\big)

When f(y) = 1,

\displaystyle \left|y\right\rangle \otimes  \left|0\oplus f(y)\right\rangle = \left|y\right\rangle \otimes  \left|1\right\rangle
\displaystyle \left|y\right\rangle \otimes  \left|1\oplus f(y)\right\rangle = \left|y\right\rangle \otimes  \left|0\right\rangle

which means

\displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} \frac{1}{\sqrt{2}} \big( \left|y\right\rangle \otimes  \left|0\oplus f(y)\right\rangle - \left|y\right\rangle \otimes \left|1\oplus f(y)\right\rangle\big) = \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1) \left|y\right\rangle \otimes \frac{1}{\sqrt{2}} \big(  \left|0\right\rangle - \left|1\right\rangle\big)

Combining these 2 we have

\begin{array}{rl}  \displaystyle \mathbf{U}_f\mathbf{H}^{\otimes n}\left| 0\right\rangle \otimes \frac{1}{\sqrt{2}} \big( \left|0\right\rangle - \left|1\right\rangle\big)  &= \displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} \frac{1}{\sqrt{2}} \big( \left|y\right\rangle \otimes  \left|0\oplus f(y)\right\rangle - \left|y\right\rangle \otimes \left|1\oplus f(y)\right\rangle\big)\\  &= \displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1)^{f(y)} \left|y\right\rangle \otimes \frac{1}{\sqrt{2}} \big(  \left|0\right\rangle - \left|1\right\rangle\big)  \end{array}

Finally, getting the Hadamard action on this gives:

\begin{array}{rl}  \displaystyle \mathbf{H}^{\otimes n}\otimes\mathbf{H}\Big[\mathbf{U}_f\mathbf{H}^{\otimes n}\left| 0\right\rangle \otimes \frac{1}{\sqrt{2}} \big( \left|0\right\rangle - \left|1\right\rangle\big)\Big] &=\mathbf{H}^{\otimes n} \mathbf{U}_f\mathbf{H}^{\otimes n}\left| 0\right\rangle \otimes \mathbf{H}\Big[\frac{1}{\sqrt{2}} \big( \left|0\right\rangle - \left|1\right\rangle\big)\Big]\\  &= \underbrace{\mathbf{H}^{\otimes n} \Big[\displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1)^{f(y)} \left|y\right\rangle \Big]}_{A} \otimes \underbrace{\mathbf{H}\Big[\frac{1}{\sqrt{2}} \big(  \left|0\right\rangle - \left|1\right\rangle\big)\Big]}_{B}  \end{array}

Expanding underbrace B gives us:

\begin{array}{rl}  \displaystyle \mathbf{H}\Big[\frac{1}{\sqrt{2}} \big(  \left|0\right\rangle - \left|1\right\rangle\big)\Big] &= \displaystyle \frac{1}{\sqrt{2}}\big(\mathbf{H}\left|0\right\rangle - \mathbf{H}\left|1\right\rangle\big)\\  &= \displaystyle \frac{1}{\sqrt{2}}\Big[\frac{1}{\sqrt{2}} \big( \left|0\right\rangle + \left|1\right\rangle \big) - \frac{1}{\sqrt{2}}\big( \left|0\right\rangle - \left|1\right\rangle\big) \Big]\\  &= \displaystyle \frac{1}{2}\Big[ 2\left|1\right\rangle \Big]\\  &= \left|1\right\rangle  \end{array}

Expanding underbrace A above:

\begin{array}{rl}  \mathbf{H}^{\otimes n} \Big[\displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1)^{f(y)} \left|y\right\rangle \Big] &= \displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1)^{f(y)}\mathbf{H}^{\otimes n} \left|y\right\rangle  \end{array}

We have found earlier that

\displaystyle \mathbf{H}^{\otimes n}\left| y\right\rangle = \frac{1}{2^{n/2}}\sum_{x=0}^{2^n-1} (-1)^{x\cdot y}\left|x\right\rangle

Substituting this to the above and noting that

\displaystyle f(y) = y\cdot a = \sum_{i=0}^{n-1} y_ia_i \mod 2,

we get:

\begin{array}{rl}  \mathbf{H}^{\otimes n} \Big[\displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1)^{f(y)} \left|y\right\rangle \Big] &= \displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1)^{f(y)}\mathbf{H}^{\otimes n} \left|y\right\rangle \\  &= \displaystyle \frac{1}{2^{n}}\sum_{y=0}^{2^n-1} \sum_{x=0}^{2^n-1} (-1)^{f(y)}(-1)^{y\cdot x}\left|x\right\rangle\\  &= \displaystyle \frac{1}{2^{n}}\sum_{y=0}^{2^n-1} \sum_{x=0}^{2^n-1} (-1)^{y\cdot a}(-1)^{y\cdot x}\left|x\right\rangle\\  &= \displaystyle \frac{1}{2^{n}}\sum_{x=0}^{2^n-1} \underbrace{\sum_{y=0}^{2^n-1} (-1)^{(a+x)\cdot y}}\left|x\right\rangle\\  \end{array}

where we have interchanged the order of the summation. To simplify the term with the underbrace, notice that

\begin{array}{rl}  \displaystyle \Big[1+ (-1)^{z_0}\Big]\Big[1+(-1)^{z_1}\Big] &= 1 + (-1)^{z_0} + (-1)^{z_1} + (-1)^{z_0+z_1}\\  &= \displaystyle (-1)^{0\cdot z_0 + 0\cdot z_1} + (-1)^{1\cdot z_0 + 0\cdot z_1} + (-1)^{0\cdot z_0 + 1\cdot z_1} + (-1)^{1\cdot z_0+ 1\cdot z_1}\\  &= \displaystyle(-1)^{0\cdot z} + (-1)^{1\cdot z} + (-1)^{2\cdot z} + (-1)^{3\cdot z}\\  &= \displaystyle \sum_{y=0}^{2^2-1} (-1)^{y\cdot z}  \end{array}

We can therefore write the factor in underbrace as

\begin{array}{rl}  \mathbf{H}^{\otimes n} \Big[\displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1)^{f(y)} \left|y\right\rangle \Big] &= \displaystyle \frac{1}{2^{n}}\sum_{x=0}^{2^n-1} \underbrace{\sum_{y=0}^{2^n-1} (-1)^{(a+x)\cdot y}}\left|x\right\rangle\\  &= \displaystyle \frac{1}{2^{n}}\sum_{x=0}^{2^n-1} \prod_{i=0}^{n-1} \Big[ \underbrace{1 +(-1)^{(a_i+x_i)}}\Big] \left|x\right\rangle\\  \end{array}

Notice that

\displaystyle 1 +(-1)^{(a_i+x_i)} = \begin{cases}  0 & x\ne a\\  2 & x=a  \end{cases}

that is, the product \prod_{i=0}^{n-1} \Big[1 + (-1)^{(a_i+x_i)}\Big] = 0 unless x=a. Therefore

\begin{array}{rl}  \mathbf{H}^{\otimes n} \Big[\displaystyle \frac{1}{2^{n/2}}\sum_{y=0}^{2^n-1} (-1)^{f(y)} \left|y\right\rangle \Big]  &= \displaystyle \frac{1}{2^{n}}\sum_{x=0}^{2^n-1} \prod_{i=0}^{n-1} \Big[ \underbrace{1 +(-1)^{(a_i+x_i)}}\Big] \left|x\right\rangle\\  &= \displaystyle \frac{1}{2^{n}} 2^n \left|a\right\rangle\\  &= \left|a\right\rangle  \end{array}

Therefore, the quantum computation gives the answer in one application of \mathbf{U}_f:

\displaystyle \mathbf{H}^{\otimes n+1}\mathbf{U}_f \mathbf{H}^{\otimes n+1} \left|0\right\rangle = \left| a\right\rangle \left| 1\right\rangle

So what just happened? We began with a superposition of states and finished with the correct answer in one application of the black box operator \mathbf{U}_f. That’s the power of Quantum Computing!

Advertisements

Prototyping Parallel Programs on OpenShift

In the previous posts, we were introduced to parallel computing and saw how it can speed up computation and by what factor (depending on the algorithm). We also introduced Chapel, a language for parallel computation, developed by Cray. Using Chapel, we first developed a very simple parallel algorithm to compute the sum of a list of numbers. Then we learned how to implement parallel matrix multiplication. In this post, we will talk about how to run our parallel programs using a cluster of machines. The source code referenced here can be found on github.

To run our parallel program, we will need to set up a bare-metal server or a Virtual Machine (VM), install Chapel and clone it to N instances. A shared filesystem will be mounted on a pre-defined directory on each instance. This filesystem will allow all instances to load data from the same location. Any output from each instance will also be written to this filesystem (although we must ensure they write to different files or directories to avoid corruption of data).

Here is how our setup will look like:

We will also need to setup the following:

  • SSH server – this will run as an ordinary user using port 2222. The Chapel platform will connect to port 2222 on each machine.
  • Configure passwordless SSH login.

We will need to export the following environment variables to enable parallel computation across different machines. The purpose of those variables can be found here. The variable GASNET_SSH_SERVERS is a space-separated list of IP addresses or hostnames of our nodes. In our setup, we will use the following:

export CHPL_COMM=gasnet
export CHPL_LAUNCHER=amudprun
export GASNET_SSH_SERVERS="node00 node01 node10 node11"

where node00, node01, node10, node11 are hostnames of our machines/nodes.

How to Run a Parallel Program

We will run our parallel matrix_multiplication on our cluster. To launch the program, we will login to one of the nodes and change directory to /home/chapel containing the mat_mul binary and issue the following command:

chpl mat_mul -nl 4

Running in OpenShift

If having to provision a bare-metal or VM environment for parallel computing is costly, we can turn to containers for a much cheaper and faster way to provision. To make our lives easier, we will use OpenShift Origin platform to run our parallel computing environment.

In the succeeding sections, we assume that we have a working installation of OpenShift. See the instructions here on how to set up OpenShift on VirtualBox or here to install a single-node OpenShift in AWS.

Containerizing the “Node”

We will have to build a container image of the bare metal or VM we mentioned above. The format we will use is Docker. First, we clone the project from github:

git clone https://github.com/corpbob/chapel-openshift-demo.git

Change directory to chapel-openshift-demo and build the image:

docker built -t chapel .

Create the chapel project on OpenShift. We use the admin user in this example:

oc login -u admin
oc new-project chapel

We need to know the ip address of the docker-registry service in OpenShift. We execute the following command:

[root@openshift chapel-openshift-demo]# oc describe svc docker-registry -n default
Name:			docker-registry
Namespace:		default
Labels:			docker-registry=default
Selector:		docker-registry=default
Type:			ClusterIP
IP:			172.30.1.1
Port:			5000-tcp	5000/TCP
Endpoints:		172.17.0.7:5000
Session Affinity:	ClientIP
No events.

Line 7 (highlighted) of the output gives us the IP address of the docker-registry service. We will use this to push the Chapel image we just created to the registry. First we login to the docker registry:

docker login -u admin -p $(oc whoami -t) 172.30.1.1:5000

Tag the chapel image using the following format

<registry_ip>:<port>/<project_name>/<image_name>:<version>

For our example, we use the following command to tag the image:

docker tag chapel 172.30.1.1:5000/chapel/chapel:v0.1

We can now push the image to the internal docker registry:

docker push 172.30.1.1:5000/chapel/chapel

Importing the Chapel Template

In the previous section, we built the container image of Chapel and pushed it to the private docker registry. In this section, we will import a template that will do the following:

  • Set up the Chapel containers inside OpenShift
  • Create a file that will dynamically generate the variable GASNET_SSH_SERVERS containing the IP addresses of the Chapel pods that will be used in the parallel computation.

The name of the template is chapel.yml.

Import the template using the command

oc create -f chapel.yml

We need to give the default service account the view role so that it can read the IP addresses of the pods associated with chapel. To do this, execute the command:

oc policy add-role-to-user view system:serviceaccount:chapel:default 

After this we can now create the chapel application:

oc new-app chapel

This will automatically trigger a deployment of 4 chapel instances with a shared volume mounted at /home/chapel/mnt.

The “hook-post” pod is a separate instance of the chapel image that will execute the following commands

echo "export GASNET_MASTERIP=\$MY_NODE_IP" > /home/chapel/mnt/exports && \
echo "export GASNET_SSH_OPTIONS=\"-p 2222\"" >> /home/chapel/mnt/exports && \
for pod in `oc get pods -l app=chapel|grep chapel|awk '{print $1}'`; \
do \
  oc describe pod $pod |grep ^IP:|awk '{print $2}'; \
done| \
awk 'BEGIN { x="" } \
      {x = x$1" "} \
     END {print "export GASNET_SSH_SERVERS=\""x"\""}' >> \
/home/chapel/mnt/exports

The output of the above command is a file named exports and looks like the below:

Running the Sample Program

We now go to the web console-> Applications -> Pods. Select any of the pods and click Terminal. In the /home/chapel directory, there is a file named run-test.sh. This file contains the following commands:

export GASNET_SPAWNFN=S                                                                                                                          
source /home/chapel/mnt/exports                                                                                                                  
                                                                                                                                                 
./hello6-taskpar-dist -nl $* 

The commands above executes the pre-compiled chapel binary hello6-taskpar-dist which was compiled when we built the container image earlier. Executing this file gives us:

# the parameter 4 tells chapel to use 4 pods to execute the command.
sh-4.2$ ./run-test.sh 4                                                                                                                               
Warning: Permanently added '[172.17.0.16]:2222' (ECDSA) to the list of known hosts.                                                              
Warning: Permanently added '[172.17.0.18]:2222' (ECDSA) to the list of known hosts.                                                              
Warning: Permanently added '[172.17.0.15]:2222' (ECDSA) to the list of known hosts.                                                              
Warning: Permanently added '[172.17.0.17]:2222' (ECDSA) to the list of known hosts.                                                              
Hello, world! (from locale 0 of 4 named chapel-1-fgg4b)                                                                                          
Hello, world! (from locale 2 of 4 named chapel-1-pt3v9)                                                                                          
Hello, world! (from locale 1 of 4 named chapel-1-rg668)                                                                                          
Hello, world! (from locale 3 of 4 named chapel-1-rw2rc)                                                                                          

Running the Parallel Matrix Multiplication

Copy the file mat_mul.chpl to the pod and compile.

chpl mat_mul.chpl -o mnt/mat_mul

The command above will place the resulting binary inside the directory /home/chapel/mnt. This will be accessible from all pods.

Finally, execute the parallel matrix multiplication:

./run.sh mnt/mat_mul 4

Conclusion

I have not tried this in a real production environment or even on a bare-metal installation of OpenShift. This seems to be a promising use of OpenShift but I still have to find out.

Implementing Parallel Algorithms Part 2

In the previous post, we implemented a very simple parallel program to add a set of numbers. In this post, we will implement parallel matrix multiplication.

We have shown a parallel algorithm to multiply 2 big matrices using message passing. The algorithm involved block sub-matrices to be passed from node to node and multiplied within a node until the answer is found.

There will be 2 square matrices A and B. In our example, the dimension of both A and B is 4×4. We will distribute the matrices evenly to 4 nodes.

A = \left[  \begin{array}{cc|cc}  a_{00} & a_{01} & a_{02} & a_{03}\\   a_{10} & a_{11} & a_{12} & a_{13}\\ \hline  a_{20} & a_{21} & a_{22} & a_{23}\\  a_{30} & a_{31} & a_{32} & a_{33}  \end{array}  \right],    B = \left[  \begin{array}{cc|cc}  b_{00} & b_{01} & b_{02} & b_{03}\\   b_{10} & b_{11} & b_{12} & b_{13}\\ \hline  b_{20} & b_{21} & b_{22} & b_{23}\\  b_{30} & b_{31} & b_{32} & b_{33}  \end{array}  \right]

In this example, node 00 will have the following matrices:

A_{00}=\begin{bmatrix}  a_{00} & a_{01}\\  a_{10} & a_{11}  \end{bmatrix},  B_{00}=  \begin{bmatrix}  b_{00} & b_{01}\\  b_{10} & b_{11}  \end{bmatrix}

Let’s simulate each node loading entries of sub-matrices assigned to it.

const n = 4;
var vec = 1..n;
var blockSize = 2;

var A: [vec, vec] real;
var B: [vec, vec] real;
var C: [vec, vec] real;

coforall loc in Locales {
  on loc {
    var i = loc.id/2;
    var j = loc.id%2;
    var istart = i*blockSize;
    var iend = istart + blockSize;
    var jstart = j*blockSize;
    var jend = jstart + blockSize;

    for (r,s) in {istart + 1..iend, jstart + 1..jend} {
      B(r,s) = r+s;
      A(r,s) = 2*r + s;
    }
  }
}

Global Address Space

Each node has limited memory physically exclusive to itself. In order for node A to have access to the contents of the memory of another node B, node B should pass the data to node A. Fortunately, Chapel can use a library called GASNet that allows each node to have a global view of the memory of all nodes participating in the computation.

In the code above, each node loads its own data. However, the GASNet library allows each node to access the matrix elements loaded by the other nodes. Consequently, we are able to reference the sub-matrix held by each node without doing fancy message passing. The algorithm is then a straightforward implementation of

\displaystyle \mathbf{C}_{ij}=\sum_{k=0}^{2} \mathbf{A}_{ik} \mathbf{B}_{kj}

where \mathbf{A_{ij}}, \mathbf{B_{ij}} and \mathbf{C_{ij}} are submatrices of \mathbf{A}, \mathbf{B}, and \mathbf{C}, respectively.

Below is the straightforward implementation of parallel block multiplication:

coforall loc in Locales {
  on loc {
    var i = loc.id/2;
    var j= loc.id%2;
    var istart = i*blockSize;
    var iend = istart + blockSize;
    var jstart = j*blockSize;
    var jend = jstart + blockSize;
    var r = { istart + 1..iend, jstart + 1..jend };
    ref W = C[r].reindex( { 1..2,1..2 });
 
    coforall k in 0..1 {
      var U=get_block_matrix(A[vec,vec],i,k,blockSize);
      var V=get_block_matrix(B[vec,vec],k,j,blockSize);
      var P = mat_mul(U,V);
      coforall (s,t) in { 1..2,1..2 } {       
        W(s,t) += P(s,t);
      }
    }
  }
}

The procedure get_block_matrix will return the sub-matrix given the (i,j)th index and the block size.

proc get_block_matrix(A: [?D], i:int, j:int , blockSize:int) {
  var r = { i*blockSize+1 .. i*blockSize 
            +  blockSize, j*blockSize 
            + 1 .. j*blockSize + blockSize };
  return A[r];
}

The procedure mat_mul will return the matrix product of two sub-matrices:

proc mat_mul(A: [?D1], B: [?D2]) {
  var D3 = { 1..2, 1..2 };
  var C: [D3] real;
  var AA = A.reindex({1..2,1..2});
  var BB = B.reindex({1..2,1..2});

  for row in 1..2 {
    for col in 1..2 {
      var sum:real = 0;
      for k in 1..2 {
         sum += AA(row,k) * BB(k,col);
      }
      C(row,col) = sum;
    }
  }
  return C;
}

writeln(C[vec,vec]);

To run this code, we need to set the following environment variables:

source $CHPL_HOME/util/setchplenv.bash 

export CHPL_COMM=gasnet
export CHPL_LAUNCHER=amudprun

export GASNET_SSH_SERVERS="127.0.0.1 127.0.0.1 127.0.0.1 127.0.0.1"

Compiling and running this program gives the output:

#Compile
chpl mat_mul.chpl -o mat_mul

# Run using 4 nodes  
./mat_mul -nl 4

A=
3.0 4.0 5.0 6.0
5.0 6.0 7.0 8.0
7.0 8.0 9.0 10.0
9.0 10.0 11.0 12.0
B=
2.0 3.0 4.0 5.0
3.0 4.0 5.0 6.0
4.0 5.0 6.0 7.0
5.0 6.0 7.0 8.0
C=
68.0 86.0 104.0 122.0
96.0 122.0 148.0 174.0
124.0 158.0 192.0 226.0
152.0 194.0 236.0 278.0

Implementing Parallel Algorithms Part 1

Now we know that parallel algorithms allow us to make our programs run faster. So how do we implement them?

I have used mpich before, but that was more than a decade ago. Recently, I found myself looking for new ways of doing parallel programming. I discovered a very nice parallel programming language called Chapel. This is what we’ll use to implement parallel algorithms.

Algorithm 1: Parallel sum of consecutive numbers from 1 to N

To get the sum of numbers from 1 to N, where N is some integer is easy. There is a formula to do that:

\displaystyle \sum_{i=1}^N i = \frac{N(N+1)}{2}

However for the sake of illustration, we are going to compute the sum of 1 to N using a cluster of machines. Here is the Chapel code to accomplish it inside the file add_parallel.chpl.

config var N:int = 1000;
var D:domain(1) = {0..numLocales -1};
var s:[D] int;
var sum:int= 0;
var bs = N/numLocales;
coforall loc in Locales {
  on loc {
      var i = loc.id;
      var start = i*bs + 1;
      var end = start + bs -1;
      var _sum:int = 0;
      for j in start .. end {
        _sum += j;
      }
    writeln("i= " + i + ", start= "+ start + ", end=" + end + ", sum = " + _sum);
    s[i] = _sum;
  }
}

sum = + reduce s;
writeln("sum: " + sum);

This program is compiled using the command:

chpl add_parallel.chpl -o add_parallel

where add_parallel.chpl is the filename of the program and -o add_parallel specifies the filename of the binary produced after compilation.

One line 1, we have defined the default value of N to be 1000. This can be overridden on the command line by specifying the --x parameter. The number of machines we are going to use is also specified on the command line using the -nl parameter. A sample invocation of this program is the following:

./add_parallel -nl 3 --N=120

The above command means that we want to run the add_parallel binary using 3 machines with the value of N=120.

Executing the above command will give us:

./add_parallel -nl 3 --N=120
i= 0, start= 1, end=40, sum = 820
i= 1, start= 41, end=80, sum = 2420
i= 2, start= 81, end=120, sum = 4020
sum: 7260

How the program works

The program will partition N into 3 blocks. This is specified on line 5 where we divided N by numLocales to get the block size. The numLocales will contains the value of the parameter -nl which in this example is 3.

The code on line 6 tells chapel to execute a parallel for-loop executing the code inside on loc block on each Locale. A Locale has an id starting from 0. The Locale will determine it’s id and compute the starting and ending number to sum and store this value in the variable _sum.

coforall loc in Locales {
  on loc {
      var i = loc.id;
      var start = i*bs + 1;
      var end = start + bs -1;
      var _sum:int = 0;
      for j in start .. end {
        _sum += j;
      }
    writeln("i= " + i + ", start= "+ start + ", end=" + end + ", sum = " + _sum);
    s[i] = _sum;
  }
}

This _sum is stored in the array s. We take the sum of entries of the s array using the reduce keyword specifying + as the reduction operator. Finally we print the total sum across the machines.

sum = + reduce s;
writeln("sum: " + sum);

Conclusion

We have seen that that we can implement parallel programs using Chapel programming language. In part 2, we will show how to do Parallel Matrix Multiplication using Chapel.

A Very Brief Introduction to Parallel Computing

Imagine you were given a hundred 3-digit numbers to add, how much time would it take you to get the answer? If it would take you 30 seconds to add ten numbers (using a calculator), then it would take you about 300 seconds (or 5 minutes) to add 100 numbers.

Now imagine there are a hundred people and each person has a number. How would a hundred people compute the sum of all numbers? Seems like a recipe for disaster. However, we can do this:

1. Group the people by two, if there is an extra person with no group, this person can join the nearest group (to make a group of 3 people).
2. Each group will add the numbers that they have to get the sum S.
3. Each group will nominate a representative that will carry this new number S. The remaining members can sit down.
4. Repeat step 1 until there is only one person remaining.
5. The last person remaining will have the sum of the 100 numbers.

At the beginning, we have 100 people in groups of two. It takes 3 seconds for them to add their respective numbers. In the next iteration, only 50 people remain that will then add their numbers. Continuing in this way, the remaining number of people will be halved until in the 7th iteration we get our answers. So if each iteration can execute in 3 seconds, then it will take 21 seconds for a hundred people to compute the sum of 100 numbers!

I mentioned that in the 7th iteration, we are able to get the answer. There is a formula to get the number of iterations and it is

\displaystyle \lceil\log_2(n)\rceil

where n is the number of people to start with and the symbol \lceil\rceil is the ceiling function which rounds off the result of the function \log_2(n) to the next highest integer. If n=100, the number of iterations is

\displaystyle \lceil\log_2 100\rceil = 7

Having a hundred people to compute the sum of 100 numbers might be practically infeasible. We might not have the space to accommodate them. However, if we only have say 10 people, we can still have a faster computation.

Map Reduce

Given a hundred numbers, we can group the numbers by 10 and distribute it to 10 people. Each person will then add the numbers they have and combine the sum with the rest of the participants to get the sum of 100 numbers in 1/10th of the time.

Grouping the numbers into smaller subsets and distributing them to each person is called mapping. Combining the sum computed by each person to the total sum is called reduction. This is called map-reduce and the example given is a simple example of map reduce.

A More Complex Example

Let \mathbf A and \mathbf B be two matrices. The product of matrices \mathbf A and \mathbf B is the matrix \mathbf C

\displaystyle  \begin{bmatrix}  c_{00} & c_{01} & \cdots & c_{0n}\\  c_{10} & c_{11} & \cdots & c_{1n}\\  \cdots & \cdots & \cdots & \cdots\\  c_{n0} & \cdots & \cdots & c_{nn}  \end{bmatrix} =  \begin{bmatrix}  a_{00} & a_{01} & \cdots & a_{0n}\\  a_{10} & a_{11} & \cdots & a_{1n}\\  \cdots & \cdots & \cdots & \cdots\\  a_{n0} & \cdots & \cdots & a_{nn}  \end{bmatrix}  \begin{bmatrix}  b_{00} & b_{01} & \cdots & b_{0n}\\  b_{10} & b_{11} & \cdots & b_{1n}\\  \cdots & \cdots & \cdots & \cdots\\  b_{n0} & \cdots & \cdots & b_{nn}  \end{bmatrix}

such that

\displaystyle c_{ij} = \sum_{k=0}^{n-1} a_{ik}b_{kj}

where c_{ij} is the entry of the matrix C on row i and column j.

Here is a sequential algorithm to compute the matrix C:

//a and b are nxn matrices
//c[i][j] is initialized to 0 for all i,j
for(var i=0;i<n;i++){
  for(var j=0;j<n;j++){
    for(var k=0;k<n;k++){
      c[i][j] += a[i][k] * b[k][j]
    }
  }
}

There are 3 loops in the above algorithm, the innermost loop will execute n times to compute the sum of element-wise product of row i and column j. In the diagram below, the inner loop will do the following: multiply each element inside the of the box of matrix A and the elements inside the box of matrix B and take the sum. Since there are n such products, the number of addition operations is n.

Then you have to do this n times for each column of matrix B and n times for each row of matrix A, as shown below, for a total of n^3 operations.

So if you to multiply a matrix with n=100, you will need to execute 100^3 (or 1 million) operations!

Parallel Computing can help us here.

Parallel Matrix Multiplication

The good thing about matrix multiplication is that we can multiply by blocks. For the sake of simplicity, suppose we have 2 square matrices of dimension 100. We can divide the matrices into small square sub-matrices of dimension 25:

A=  \left[\begin{array}{c|c|c|c}  \mathbf A_{00} & \mathbf A_{01} & \mathbf A_{02} & \mathbf A_{03}\\  \hline  \mathbf A_{10} & \mathbf A_{11} & \mathbf A_{12} & \mathbf A_{13}\\  \hline  \mathbf A_{20} & \mathbf A_{21} & \mathbf A_{22} & \mathbf A_{23}\\  \hline  \mathbf A_{30} & \mathbf A_{31} & \mathbf A_{32} & \mathbf A_{33}  \end{array}  \right]  ,  B=  \left[\begin{array}{c|c|c|c}  \mathbf B_{00} & \mathbf B_{01} & \mathbf B_{02} & \mathbf B_{03}\\  \hline  \mathbf B_{10} & \mathbf B_{11} & \mathbf B_{12} & \mathbf B_{13}\\  \hline  \mathbf B_{20} & \mathbf B_{21} & \mathbf B_{22} & \mathbf B_{23}\\  \hline  \mathbf B_{30} & \mathbf B_{31} & \mathbf B_{32} & \mathbf B_{33}  \end{array}  \right]

where each \mathbf A_{ij} and \mathbf B_{ij} are square matrices of dimension 25.

We can then compute the resulting sub-matrix of C using the usual formula:

\displaystyle \mathbf C_{ij} = \sum_{k=0}^{p-1} \mathbf A_{ik}\mathbf B_{jk}

where p=100/25=4.

For example, to compute \mathbf C_{00} we have

\mathbf C_{00} = \mathbf A_{00} \mathbf B_{00} + \mathbf A_{01} \mathbf B_{10} + \mathbf A_{02} \mathbf B_{20} + \mathbf A_{03} \mathbf B_{30}

We then use this mechanism to distribute our sub-matrices to different computers (or in modern parlance “compute nodes”). For this example, we need 4×4=16 compute nodes. Each compute node will contain 2 sub-matrices, one for A and one for B. For easy visualization, we can make a drawing of our compute nodes arranged in a square of side 4 and labelled as shown below:

Now that we have evenly distributed our matrix data to each node, the next question is how do we compute? Notice that not one node contains all the data. Each node has a very limited subset of the data.

The Trick

We can get a clue by looking at the end result of the computation for Node 00. At the end of the computation, Node 00 should have the following result:

\mathbf C_{00} = \mathbf A_{00} \mathbf B_{00} + \mathbf A_{01} \mathbf B_{10} + \mathbf A_{02} \mathbf B_{20} + \mathbf A_{03} \mathbf B_{30}

Looking at the above formula, Node 00 only has the following data

\mathbf A_{00} \text{ and } \mathbf B_{00}

The rest of the sub-matrices are not in its memory. Therefore, Node 00 can only compute

\mathbf A_{00} \times \mathbf B_{00}

We are missing the following products:

\mathbf A_{01}\mathbf B_{10}, \mathbf A_{02}\mathbf B_{20}, \text{ and } \mathbf A_{03}\mathbf B_{30}

The matrix \mathbf A_{01} is with Node 01 and the matrix \mathbf B_{10} is with Node 10. So if Node 01 can send matrix \mathbf A_{01} to Node 00 and Node 10 can send matrix \mathbf B_{10} to Node 00, we can then get the product \mathbf A_{01}\mathbf B_{10}. In fact, if we do a slight rearrangement like the below, we can use the following algorithm to compute matrix \mathbf C:

1. Each node will send the current A sub-matrix to the node on the left and receive a new A sub-matrix from the node on the right. If there is no node on the left, it will send the sub-matrix to the last node on its row.
2. Each node will send the B sub-matrix to the node on top and receive a new B sub-matrix from the node below it. If there is no node on top, it will send the sub-matrix to the last node on its column.
3. Multiply the new sub-matrices and add the result to the current value of the sub-matrix of C that it keeps in memory.
4. Repeat until the number of iterations equal to N, where N^2 is the number of nodes.

The figure below describes this algorithm for Node 00.

Doing this for all nodes, we can visualize the parallel algorithm at work on all nodes as shown in the animation below:

This algorithm is one of my favorite algorithms since it looks like the pumping of blood in to the heart.

How fast is this algorithm?

Given a 2 square matrices of dimension N, the number of sequential computations is O(N^3). Using parallel matrix multiplication above, we divide the matrices into sub-matrices of dimension n. The resulting block matrix is of dimension N/m=p. Each node will now compute the product in parallel using O(n^3) operations per sub-matrix multiplication. Since there are p multiplications per node, the number of operations per node is O(pn^3) The ratio of these two quantities will give us how fast the parallel algorithm is

\displaystyle \frac{N^3}{pn^3} = \frac{1}{p}\frac{N^3}{n^3} = \frac{1}{p}\cdot p^3 = p^2

For our particular example, the theoretical speedup is p^2 = 4^2 = 16, that is, our parallel algorithm can compute the product 16 times faster than the sequential algorithm. If we increase p, we can increase the speedup. However, we can only increase p up to a certain point as the network will then become the bottleneck and will make things slower than the sequential algorithm.

In the next 3 posts, we will see how to program simple parallel algorithms and parallel matrix multiplication. We will also show how you can use OpenShift to prototype your parallel programs.

Have You Seen A Half-Balloon?

Balloons come in all shapes and sizes but the shape the comes to mind when I hear the word balloon is that of a spherical or round shape. Given a spherical balloon, create an imaginary partition that divides the balloon into left and right hemispheres. This effectively places the air molecules into left and right. The number of molecules in the left is more or less the same as the number of molecules on the right. This is what we see in real life.

As we pump air into the balloon, a given air molecule can randomly go left or right of the partition. If all air molecules for some reason go left (or right), then we have a situation called a “half-balloon”.

So why aren’t we seeing half-balloons? It’s because that situation is highly improbable. For illustration purposes, imagine there are 6 air molecules. Each molecule has two equally likely choices, be on the left or on the right. Now imagine each air molecule has chosen it’s “side”. A possible combination or configuration (a term we will use moving forward) will be 3 molecules choose Left and 3 choose Right:

L  L  L  R  R  R

Where L stands for Left and R stands for Right.

In fact, the number of such configurations is equal to

\displaystyle \underbrace{2\times 2\times \ldots \times 2}_{\text{6 times}} = 2^6 = 64

We enumerate below the list of possible air molecule configurations:


# This is the code in R
# x=c("L","R")
# expand.grid(m1=x,m2=x,m3=x,m4=x,m5=x,m6=x)
# d1=expand.grid(m1=x,m2=x,m3=x,m4=x,m5=x,m6=x)
   m1 m2 m3 m4 m5 m6
1   L  L  L  L  L  L
2   R  L  L  L  L  L
3   L  R  L  L  L  L
4   R  R  L  L  L  L
5   L  L  R  L  L  L
6   R  L  R  L  L  L
7   L  R  R  L  L  L
8   R  R  R  L  L  L
9   L  L  L  R  L  L
10  R  L  L  R  L  L
11  L  R  L  R  L  L
12  R  R  L  R  L  L
13  L  L  R  R  L  L
14  R  L  R  R  L  L
15  L  R  R  R  L  L
16  R  R  R  R  L  L
17  L  L  L  L  R  L
18  R  L  L  L  R  L
19  L  R  L  L  R  L
20  R  R  L  L  R  L
21  L  L  R  L  R  L
22  R  L  R  L  R  L
23  L  R  R  L  R  L
24  R  R  R  L  R  L
25  L  L  L  R  R  L
26  R  L  L  R  R  L
27  L  R  L  R  R  L
28  R  R  L  R  R  L
29  L  L  R  R  R  L
30  R  L  R  R  R  L
31  L  R  R  R  R  L
32  R  R  R  R  R  L
33  L  L  L  L  L  R
34  R  L  L  L  L  R
35  L  R  L  L  L  R
36  R  R  L  L  L  R
37  L  L  R  L  L  R
38  R  L  R  L  L  R
39  L  R  R  L  L  R
40  R  R  R  L  L  R
41  L  L  L  R  L  R
42  R  L  L  R  L  R
43  L  R  L  R  L  R
44  R  R  L  R  L  R
45  L  L  R  R  L  R
46  R  L  R  R  L  R
47  L  R  R  R  L  R
48  R  R  R  R  L  R
49  L  L  L  L  R  R
50  R  L  L  L  R  R
51  L  R  L  L  R  R
52  R  R  L  L  R  R
53  L  L  R  L  R  R
54  R  L  R  L  R  R
55  L  R  R  L  R  R
56  R  R  R  L  R  R
57  L  L  L  R  R  R
58  R  L  L  R  R  R
59  L  R  L  R  R  R
60  R  R  L  R  R  R
61  L  L  R  R  R  R
62  R  L  R  R  R  R
63  L  R  R  R  R  R
64  R  R  R  R  R  R

Looking at the data above, we can see that there are 6 configurations that contain only 1 “L”:

# This is the code in R
# cc=c()
# for(i in 1:64){
#  tmp=d1[i,]
#  cc=c(cc,sum(tmp == "L"))
# }
# d1$count = cc
# d1[d1$count == 1, ]
   m1 m2 m3 m4 m5 m6 count
32  R  R  R  R  R  L     1
48  R  R  R  R  L  R     1
56  R  R  R  L  R  R     1
60  R  R  L  R  R  R     1
62  R  L  R  R  R  R     1
63  L  R  R  R  R  R     1

This means that the probability of getting a configuration having one molecule on the left and 5 molecules on the right is:

\displaystyle \text{Probability of 1 molecule Left and 5 molecules Right} = \frac{6}{2^6} = 0.093750

If we continue summarizing our data so that we get the number of combinations containing 2, 3, 4, 5, and 6 “L”s, we can generate what is known as a probability distribution of air molecules on the left hemisphere:

# cc=c()
# for(i in 0:6){
#  cc=c(cc,length(attributes(d1[d1$count == i, ])$row.names))
# }
# data.frame(X=0:6,count=cc,probability=cc/2^6)
  X count probability
1 0     1    0.015625
2 1     6    0.093750
3 2    15    0.234375
4 3    20    0.312500
5 4    15    0.234375
6 5     6    0.093750
7 6     1    0.015625

The graph of this probability distribution is shown below:

# Here is the code that generated the plot above
barplot(cc/2^6,names.arg=0:6,main="Probability Distribution \nNumber of Air Molecules on the Left Hemisphere")

The graph above tells us that the most probable configuration of balloon molecules is that half of them are on the right and the other half on the left as shown by the high probability of the value X = 3. It also tells us the probability of all molecules choosing the Left side equal to 0.015625. When the number of air molecules is very large, this probability will turn out to be extremely small, as we will show below.

The Mathematics

At this point, manually counting the number of rows will not be practical anymore for a large number of molecules. We need to use some mathematical formula to compute the combinations. We don’t really care which molecules chose left or right. We just care about the number of molecules on the left. Given N molecules, there are 2^N possible configurations of air molecules. Of these configurations, there are

\displaystyle {N \choose m } = \frac{N!}{m! (N-m)!}

combinations that have m molecules on the left. Therefore, the probability of a configuration having m molecules on the left is

P(m) = \displaystyle \frac{\displaystyle {N\choose m}}{\displaystyle 2^N}

This is a probability density function since

\displaystyle \sum_{m=0}^N P(m) = \displaystyle \sum_{m=0}^N \frac{\displaystyle {N\choose m}}{\displaystyle 2^N} = 1

To show this, we will use the Binomial Theorem

\displaystyle \sum_{m=0}^N {N \choose m} x^m = (1+x)^N

If we let x=1, the Binomial Theorem gives us

\displaystyle \sum_{m=0}^N {N \choose m} = (1+1)^N = 2^N

Therefore

\begin{array}{rl}  \displaystyle \sum_{m=0}^N P(m) &= \displaystyle \sum_{m=0}^N \frac{\displaystyle {N\choose m}}{\displaystyle 2^N}\\  &= \displaystyle \frac{1}{2^N} \sum_{m=0}^N {N\choose m}\\  &= \displaystyle \frac{1}{2^N} 2^N\\  &= 1  \end{array}

A Mole of Air

One mole of air contains 6.022 \times 10^{23} . This means the probability of that all molecules are on the left side of the balloon is

\displaystyle \frac{1}{\displaystyle 2^{6.022 \times 10^{23}}} < 1/1024^{10^{22}} < 1/(10^3)^{10^{22}}=10^{-30,000,000,000,000,000,000,000}

This is a very small number such that it contains 30 million trillion zeros to the right of the decimal point. When you write this zeros on a piece of paper with a thickness of 0.05 millimeters, you would need to stack it up to a height 74 times the distance of earth to pluto in kilometers!

An Interview Question: Using Integer Programming

We can solve the Interview Question using a mathematical technique called Integer Programming. Let d_1, d_2, \ldots, d_N be the variables representing diskette 1, diskette 2, diskette 3, etc. The values of the d_k variables can only be 0 or 1. A 0 means the diskette is not used while a 1 means that it is used.

Each file is saved to a certain diskette. We want to know to what diskette d_i a given file f_j is assigned. To represent this, we assign the variable a_{ij} a value of 1 if file f_j is assigned to diskette d_i.

We will normalize the file sizes so that if s_i is the size of f_i, the s_i \le 1. We do this by simply dividing all file sizes by the size of the diskette. For a given diskette d_i, the following constraint should be satisfied:

d_i - s_1a_{i1} - s_2a_{i2} - \ldots - s_N a_{iN} \ge 0

for diskette i = 1, 2, \ldots, N and s_i are the normalized file sizes of file f_i for i=1,2,\ldots,N.

Since each file f_j can only be assigned to one diskette, we have the following constraint:

a_{1j} + a_{2j} + \ldots + a_{Nj} = 1

where a_{1j} is the variable representing the “file f_j is in diskette d_1“, etc.

Finally, we have to constrain the value of d_i to be either 0 or 1, that is,

d_i \le 1

for all i=1,2,\ldots,N.

Integer Programming Formulation

Given the above information, we can formulate the Integer Programming problem as

Minimize:

d_1 + d_2 + d_3 + \ldots + d_N

subject to

\begin{array}{rl}  d_1 - s_1a_{11} - s_2a_{12} - s_3a_{13} - \ldots - s_Na_{1N} &\ge 0\\  d_2 - s_1a_{21} - s_2a_{22} - s_3a_{23} - \ldots - s_Na_{2N} &\ge 0\\  :\\  d_N - s_1a_{N1} - s_2a_{N2} - s_3a_{N3} - \ldots - s_Na_{NN} &\ge 0\\  a_{11} + a_{21} + a_{31} + \ldots + a_{N1} &= 1\\  a_{12} + a_{22} + a_{32} + \ldots + a_{N2} &= 1\\  :\\  a_{1N} + a_{2N} + a_{3N} + \ldots + a_{NN} &= 1\\  d_1 &\le 1\\  d_2 &\le 1\\  :\\  d_n &\le 1  \end{array}

Solving the Problem

We will use R to solve this Integer Programming Formulation. Please see code below:

library("lpSolve")
NUMFILES=4

# Generate random file sizes between 1 and 10
FileSizes=ceiling(10*runif(NUMFILES))
x = -1*FileSizes/10
l=length(x)

# Each files can be in any of the diskettes. Suppose there are N files,
# to determine if a file j is in diskette i, the value of variable x_ij will
# 1 if file j is in diskette i, and 0 otherwise.
# Here we construct the coefficients of variables x_ij which are the
# sizes of the files (normalized to 1)
zz=c()
for(i in 1:(l-1)){
  zz=c(zz,x,rep(0,l*l))
}
zz=c(zz,x)

# Construct the coefficients of the indicator variables representing the 
# diskettes d_i
zzmatrix=matrix(zz,ncol=l*l,byrow=T)
CoefficientsOfDiskettes=c();
for(i in 1:l){
 ttt=rep(0,l)
 ttt[i] = 1
 CoefficientsOfDiskettes= c(CoefficientsOfDiskettes,ttt,zzmatrix[i,])
}

# Construct the coefficients of x_ij for constant j. These variables 
# satisfy the equation \sum_{i=1}^N x_{ij}
SumOfFileAcrossDiskettes=c()
for(i in 1:l){
  ttt=rep(0,l)
  ttt[i]=1
  SumOfFileAcrossDiskettes=c(SumOfFileAcrossDiskettes,rep(ttt,l))
}

# Prepend Coefficients of variables d_i. The value of these coefficients is 0. 
SumOfFileAcrossDiskettesMatrix=matrix(SumOfFileAcrossDiskettes,ncol=l*l,byrow=T)
PrependCoefficientsOfDiskettes=c()
for(i in 1:l){
 PrependCoefficientsOfDiskettes=c(PrependCoefficientsOfDiskettes,c(rep(0,l),SumOfFileAcrossDiskettesMatrix[i,]))
}

# Construct coefficients of d_i to construct constraint d_i <= 1
DisketteConstraints=c()
for(i in 1:l){
 ttt=rep(0,l)
 ttt[i]=1
 DisketteConstraints=c(DisketteConstraints,ttt,rep(0,l*l))
}

# Construct matrix input of lpSolve
const.mat=matrix(c(CoefficientsOfDiskettes,PrependCoefficientsOfDiskettes,DisketteConstraints),ncol=l*(l+1),byrow=T)

print("Matrix Coefficients:")
print(const.mat)

# Construct inequalities/equalities
const.dir=c(rep(">=",l),rep("=",l),rep("<=",l))

# Construct Right-Hand side
const.rhs=c(rep(0,l),rep(1,l),rep(1,l))

# Construct Objective Function
objective.in=c(rep(1,l),rep(0,l*l))

# Invoke lpSolve
mylp=lp(direction="min",objective.in=objective.in,const.mat=const.mat,const.dir=const.dir,const.rhs=const.rhs,all.int=T)

# Print Results
print(paste("Number of Diskettes: ", sum(mylp$solution[1:l])))
tz=matrix(mylp$solution,ncol=l,byrow=T)
print("File Sizes: ")
print(FileSizes)
for(i in 2:(l+1)){
   files = which(tz[i,] == 1)
   if(length(files) > 0){
      print(paste("Files in diskette ", i-1))
      print(files)
   }
}

Most of the code above is setting up the matrix of coefficients. The line 70 then calls on lpSolve to compute the optimal values of the variables

Program Output

Running this code we get the output

[1] "Matrix Coefficients:"
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
 [1,]    1    0    0    0   -1 -0.2 -0.1 -0.1    0   0.0   0.0   0.0     0   0.0   0.0   0.0     0   0.0   0.0   0.0
 [2,]    0    1    0    0    0  0.0  0.0  0.0   -1  -0.2  -0.1  -0.1     0   0.0   0.0   0.0     0   0.0   0.0   0.0
 [3,]    0    0    1    0    0  0.0  0.0  0.0    0   0.0   0.0   0.0    -1  -0.2  -0.1  -0.1     0   0.0   0.0   0.0
 [4,]    0    0    0    1    0  0.0  0.0  0.0    0   0.0   0.0   0.0     0   0.0   0.0   0.0    -1  -0.2  -0.1  -0.1
 [5,]    0    0    0    0    1  0.0  0.0  0.0    1   0.0   0.0   0.0     1   0.0   0.0   0.0     1   0.0   0.0   0.0
 [6,]    0    0    0    0    0  1.0  0.0  0.0    0   1.0   0.0   0.0     0   1.0   0.0   0.0     0   1.0   0.0   0.0
 [7,]    0    0    0    0    0  0.0  1.0  0.0    0   0.0   1.0   0.0     0   0.0   1.0   0.0     0   0.0   1.0   0.0
 [8,]    0    0    0    0    0  0.0  0.0  1.0    0   0.0   0.0   1.0     0   0.0   0.0   1.0     0   0.0   0.0   1.0
 [9,]    1    0    0    0    0  0.0  0.0  0.0    0   0.0   0.0   0.0     0   0.0   0.0   0.0     0   0.0   0.0   0.0
[10,]    0    1    0    0    0  0.0  0.0  0.0    0   0.0   0.0   0.0     0   0.0   0.0   0.0     0   0.0   0.0   0.0
[11,]    0    0    1    0    0  0.0  0.0  0.0    0   0.0   0.0   0.0     0   0.0   0.0   0.0     0   0.0   0.0   0.0
[12,]    0    0    0    1    0  0.0  0.0  0.0    0   0.0   0.0   0.0     0   0.0   0.0   0.0     0   0.0   0.0   0.0
[1] "Number of Diskettes:  2"
[1] "File Sizes: "
[1] 10  2  1  1
[1] "Files in diskette  1"
[1] 2 3 4
[1] "Files in diskette  2"
[1] 1

Interpreting the Result

Lines 2-14 of the output gives you the matrix of coefficients. Line 15 prints the number of diskettes needed to store the files. Line 17 prints the randomly generated file sizes from 1 to 10. Finally lines 18-21 prints which diskettes contain which files.

The space complexity of this solution is quite substantial. Given N files, we need to specify N^2 + N variables by 3\times N equations for a total of (N^2 + N)\times 3N memory space for coefficients.