Saturday, 1 February 2014

quantum walking

It’s possible to build a quantum random walk simulator in Python/NumPy with code that is very close to the mathematical definitions.  Here’s how.

First, we need to import NumPy (to do the array operations) and matplotlib (to visualise the results).
from numpy import *
from matplotlib.pyplot import *
We define the number of steps, +++N+++, that we are going to walk.  We also define the total number of different positions the walker can be in after +++N+++ steps.
N = 100      # number of random steps
P = 2*N+1    # number of positions

a quantum coin

We toss a quantum coin to decide whether to go left or right (or a superposition).  In ket notation, we can write a general quantum coin as an arbitrary superposition of two states:
$$| coin \rangle = a | 0 \rangle_c + b | 1 \rangle_c; \mbox{ where} |a|^2 + |b|^2 = 1$$
The ket notation is a convenient shorthand for the actual vectors representing the state:
$$| 0 \rangle_c  = \left( \begin{array}{c}
1  \\
0  \end{array} \right);\ | 1 \rangle_c  = \left( \begin{array}{c}
0  \\
1  \end{array} \right)$$
So we can use a NumPy array to define a coin state:
coin0 = array([1, 0])  # |0>
coin1 = array([0, 1])  # |1>

the Hademard coin operator

We will see terms like +++| 0 \rangle_c \langle 0 |+++ in what follows.  This is the outer product of the relevant vectors, resulting in a matrix.  We can use NumPy to calculate these:
C00 = outer(coin0, coin0)  # |0><0| 
C01 = outer(coin0, coin1)  # |0><1| 
C10 = outer(coin1, coin0)  # |1><0| 
C11 = outer(coin1, coin1)  # |1><1| 
Quantum operators are unitary matrices.  The coin operator, that can be used to flip a quantum coin into a superposition, is:
$$\hat{C} = \frac{1}{\sqrt{2}} \left( | 0 \rangle_c \langle 0 | + | 0 \rangle_c \langle 1 | + | 1 \rangle_c \langle 0 | - | 1 \rangle_c \langle 1 | \right)$$
C_hat = (C00 + C01 + C10 - C11)/sqrt(2.)

position on a line

In ket notation, we can write the fully general position of the walker on the line as an arbitrary superposition of the +++P+++ possible states:
$$| posn \rangle = \sum_k \alpha_k | k \rangle_p ; \mbox{ where} \sum_k | \alpha_k |^2 = 1 $$

shift (step) operator

The step operator moves left or right along the line, depending on the value of the coin:
$$\hat{S} = | 0 \rangle_c \langle 0 | \otimes \sum_k | k+1 \rangle_p \langle k | + | 1 \rangle_c \langle 1 | \otimes \sum_k | k-1 \rangle_p \langle k |  $$
We assume the line is actually on a circle, so the positions at the ends wrap around.  However, we will always make the circle big enough so that this doesn’t happen during a walk.  The tensor product +++\otimes+++ is implemented with the NumPy kron operation:
ShiftPlus = roll(eye(P), 1, axis=0)
ShiftMinus = roll(eye(P), -1, axis=0)
S_hat = kron(ShiftPlus, C00) + kron(ShiftMinus, C11)

walk operator

The walk operator combines the coin operator on the coin state, and a step operator on the combined coin and position state:
$$\hat{U} = \hat{S} \left( \hat{C} \otimes \hat{\bf I}_p \right)$$
U =, C_hat))

initial state

Let’s take the initial state of the system to be a coin in a superposition of left and right, and the walker at position 0:
$$ | \psi \rangle_0 = | coin \rangle_0 \otimes | posn \rangle_0 = \frac{1}{\sqrt{2}} \left( | 0 \rangle_c + i | 1 \rangle_c\right) \otimes | 0 \rangle_p$$
posn0 = zeros(P)
posn0[N] = 1     # array indexing starts from 0, so index N is the central posn
psi0 = kron(posn0,(coin0+coin1*1j)/sqrt(2.))

state after +++N+++ steps

Then walking +++N+++ steps is just applying the walk operator +++N+++ times:
$$ | \psi \rangle_N = \hat{U}^N | \psi \rangle_0$$
psiN = linalg.matrix_power(U, N).dot(psi0)
And we’re done! +++ | \psi \rangle_N+++ is the state of the system after +++N+++ random quantum steps.

measurement operator

We can measure the state at position +++k+++ using the measurement operator
$$\hat{M}_k = \hat{\bf I}_c \otimes | k \rangle_p \langle k |$$
We can use this to build up an array of probabilities, by taking the modulus squared of the state value at each position. (We can calculate the whole distribution in one go in simulation, but we would only get one measurement per experiment on the real quantum system.)
prob = empty(P)
for k in range(P):
    posn = zeros(P)
    posn[k] = 1     
    M_hat_k = kron( outer(posn,posn), eye(2))
    proj =
    prob[k] =

plot the distribution

We can then plot the probabilities.
fig = figure()
ax = fig.add_subplot(111)

plot(arange(P), prob)
plot(arange(P), prob, 'o')
loc = range (0, P, P / 10) #Location of ticks
xlim(0, P)
ax.set_xticklabels(range (-N, N+1, P / 10))


For +++N=100+++ we get
probability distribution for a quantum random walk with +++N=100+++, symmetric initial coin

The maximum probability occurs at +++\approx N/\sqrt{2}+++.

If instead of starting with a symmetric initial coin, we start with +++|0\rangle_c+++, we get
probability distribution for a quantum random walk with +++N=100+++, initial coin +++|0\rangle+++
In neither case do these look anything like the bell curve of a classical random walk, with its maximum probability at +++0+++.  If you are drunk and trying to stagger home, best be quantum!

What I find most impressive about the Python is how closely we can make the code follow the mathematical formalism thoughout.


  1. Thanks for this...saves me writing it myself for my Computing 2 students, I'll just direct therm here instead! -- Viv

  2. Could someone help me understand the ShiftPlus and ShiftMinus operators written in terms of roll function?

  3. Both Pyhton as well as the quantum walk problem are new to me. Could someone help me to understand the ShiftPlus and ShiftMinus operators written in terms of some mysterious roll function?

    1. "roll" is a numpy function, with definition and examples defined here:

  4. Hi Susan,
    Thank you very much for sharing this great work.
    I have question please, If I want to use the same shift but add to it the mod

  5. quick fyi, the code returns a float P unless we use int, as bellow :-)
    loc = range (0, P, P // 10) #Location of ticks
    ax.set_xticklabels(range (-N, N+1, P // 10))

  6. Hi, In the classical Drunkard's Walk, the limiting expectation, n -> infinity, E(Sn) = 0, but E(|Sn|) = sqrt(2n/PI). Sn is sum of +1's and -1's at n steps. Any individual finite walk has to be positive or negative (sometimes zero). So E(|Sn|) gives an indication of long term expectation not distinguishing between directions, thus, in some sense, of any individual finite walk. However, in the quantum formulation, where any individual finite walk does not need to be positive or negative, instead a superposition of possible states, what is E(|Sn|) ? Might it be zero, rather than sqrt(2n/PI) ? Or something else ? Same as the classical value, or different ? Presumably, an expectation value of a variable would look something like the integral of an appropriately normalised absolute squared wavefunction times that variable ? Can I infer E(|Sn|) immediately from your above plots, or does it need extra work ? And is E(Sn) obvious ? Thanks, Bill

  7. Perhaps my question might be better expressed: what is the quantum analogue of the classical result E(|Sn|) = sqrt(2n/PI) and does it equal zero ? Does your code answer this question, as is ? I guess at my current level of understanding of the formalism, I cannot see immediately, how to construct the equivalent of the classical expectation value, for either Sn or |Sn|. If you have a reasonably quick answer, I'd be interested to hear it. If not, then I should set up your code myself and understand it better. Your blog is very interesting.

    1. Rather than E(|x|), it's usually easier to look at sqrt E(x^2), or E(sigma) ie, the expected std deviation. For the classical walk, this is also O(sqrt(x)) -- for the quantum case, it is O(x).

      Kempe's arXiv paper "Quantum random walks - an introductory overview" at has lots of interesting info -- including this faster than classical spreading.