Optimize evaluation time for a Quantum Walk simulation

Question

Here's the problem. I'm trying to simulate a noisy Quantum Walk to obtain results similar to these, experimentally obtained.

I wrote my code, and I get results which are in accordance to what I expected.

The problem is, my evaluation time is too high. In the article they achieve a simulation of a 100-step QW, but I have to wait up to 5 hours just to get a 16-step result.

These are the time-consuming lines of code:

finalstate = 
  ParallelTable[
    Expand[
      Nest[TotalF, 
        state[1, 1, -1, j] state[2, 0, 1, j] - 
          state[1, 0, -1, j] state[2, 1, 1, j], 
        length]], 
    {j, 1, jmax}];

probabilities = 
  ParallelTable[
    Table[{i, k, 
      Total[
        Cases[Expand @ finalstate[[j]], 
          a___ state[1, _, i, _] state[2, _, k, _] :> 
            Abs[If[a === Null, 1, a]]^2]]}, 
        {i, -length - 10, length + 10}, 
        {k, -length - 10, length + 10}], 
     {j, 1, jmax}];

Where TotalF is defined as such:

phase = 
  Table[
    Table[{Exp[I*RandomReal[{0, ϕmax}]], Exp[I*RandomReal[{0, ϕmax}]]}, 
      {i, -length - 10, length + 10}], 
    {j, 1, jmax}];

TotalF[state[index_, coin : (0 | 1), position_, conf_]] := 
  Module[{row = position + length + 3, col = If[coin == 0, 1, 2]}, 
    1/Sqrt[2] * phase[[conf, row, col]] * 
      (state[index, 0, position + 1, conf] + 
         If[coin == 0, 1, -1] state[index, 1, position - 1, conf])];

TotalF[u_ + v_] := TotalF[u] + TotalF[v];
TotalF[u_?(FreeQ[#, state] &) v_] := u TotalF[v];
TotalF[u_ v_] := TotalF[u] TotalF[v];

ϕmax, length and jmax are constants.

Now, evaluating finalstate takes 2.5 seconds for 8 steps and 680 seconds for 16 steps. Evaluating probabilities takes 37 seconds for 8 steps and a lot more (about 5 hours, didn't try it with AbsoluteTiming) for 16 steps, so if I could at least speed up this last part (although if I want to achieve the 100-steps result I'd have to optimize the first one as well) it would be great.

I hope I made myself clear enough, I'm willing to provide extra information if needed.

Answer 1

There are two parts that are immediately suspicious to me: your phase calculation and finalState.

Your calculation of phase is rather inefficient. Table is a good construct, but it is often slow. The form you have is not even allowing it to speed things up: you are nesting them instead of using both iterators with one call. Honestly, I would not even bother with Table and use the vector properties of Exp, e.g.

phase = Exp[I RandomReal[{0, \[Phi]max}, {jmax, 2 length + 21, 2}]];

which generates a packed array.

Second, in trying to determine what a single element of finalstate gives you, I ran into speed issues with only iterations. At issue is at each iteration of Nest, TotalF must go through a progressively larger expression and it is not very fast at doing so. The simpler thing to do is to expand the expression at each iteration,

Nest[
  Expand[TotalF[#]]&,
  state[1, 1, -1, j] state[2, 0, 1, j] - state[1, 0, -1, j] state[2, 1, 1, j],
  length
]

which simplifies the job TotalF has to do. On my computer, this easily works through length == 16. Of course, with a length greater than 16, I get $RecursionLimit::reclim2 messages, so I will have to think on how to better approach this.

---

The $RecursionLimit::reclim2 messages limit messages can be eliminated by changing how TotalF handles Plus. Map over Plus instead of handling it pairwise:

TotalF[a_Plus] := TotalF /@ a

We could do something similar for Times, as well, and it lends itself to a nice simplification:

Clear[TotalF];
TotalF[state[index_,coin:(0|1),position_,conf_]]:=
Module[{row=position+length+3,col=If[coin==0,1,2]},
   1/Sqrt[2]*phase[[conf,row,col]]*(state[index,0,position+1,conf]
         +If[coin==0,1,-1] state[index,1,position-1,conf])
];

TotalF[a:_Plus|_Times]:=TotalF/@a
TotalF[a_] := a

On my machine, just mapping over Plus with a length of 50 took 23 s, but eliminating the PatternTest, (u_?(FreeQ[#, state] &) v_), dropped the time it took to 15 s. A nice improvement, especially if run in parallel.

---

probabilities will require more thought.

Answer 2

Here's a possible solution. It's working (really well) for 1 particle, but needs to be checked for two particles, I guess due to minor mistakes.

Basically, I changed approach and turned to a matrix notation.

Let's start from the one-particle system. We simply write the state as a vector, and the operators as matrices, being very careful about the indexing.

So, for a system of length $2L+1$ (that means an evolution over $L$ steps starting from position 0), the state is a tensor product of the position vector and the coin vector:

posInitialState = Table[
    Normal[
     SparseArray[
       length[[v]] + 1 -> 1, 2*length[[v]] + 1]], {v, 1,lmax}];

coinInitialState = {{1, 0}, {0, 1}, {1/Sqrt[2], I/Sqrt[2]}};

where the position space has dimension $2L+1$ and the coin space has dimension $2$. length[[v]] is simply an array which determines the length of the QW, as defined in the question.

We then define the displacement and coin operators:

displacementMatrix[v_] := 
    DiagonalMatrix[Flatten[ConstantArray[{1, 0}, 2*length[[v]]]], 2] + 
    DiagonalMatrix[Flatten[ConstantArray[{0, 1}, 2*length[[v]]]], -2];


staticCoinMatrix = 
    Table[
     KroneckerProduct[
       DiagonalMatrix[Exp[I RandomReal[{0, ϕmax[[u]]},2 length[[v]] +1]]],
       HadamardMatrix[2]], {u,1,fmax},{v,1,lmax}];

Now the evolution after $L$ steps is, of course, given by the $n-th$ power of the product of these two:

unitaryStaticEvolution[u_,v_]:=
    MatrixPower[displacementMatrix[v].staticCoinMatrix[[u,v]], length[[v]]];

Now it's straightforward: let's define the initial state in the total Hilbert space:

initialState = Table[Flatten[
   KroneckerProduct[posInitialState[[v]], coinInitialState[[0]]]], 
    {v, 1, lmax}];

and

finalstate = ParallelTable[
state = 
Table[ unitaryStaticEvolution[u, v].initialState[[v]], {j, 1, jmax}];
state2 = 
Table[Abs[state[[j, 1 ;; -2 ;; 2]]]^2 + 
  Abs[state[[j, 2 ;; -1 ;; 2]]]^2, {j, 1, jmax}];
Mean[state2[[1 ;; -1]]]
, {u, 1, fmax}, {v, 1, lmax}
];

Finally, we get our probabilities in a nice table, ready to be plotted:

probabilities = Table[
    Partition[Riffle[
         Table[i,{i,-length[[v]],length[[v]]}],
         finalstate[[u,v]]],2],
    {u,1,fmax},{v,1,lmax}];

For two particles, it should be quite the same. We could define the states and operators as tensor products of the single-particle states, and do the same thing. I'm still working on it, though, as

static2pEvo[u_,v_]:
KroneckerProduct[unitaryStaticEvolution[u,v],unitaryStaticEvolution[u,v]];

does not give me something I can work on. What I get in this case are probabilities which are the same for every length or noise configuration.

But the approach is correct, as it works for 1 particle, and gives me results that are way faster than the older ones (not having to deal with pattern matchings).

---

EDIT: ok, I got the error, there was a mistake in the coin operator's definition. Now it's corrected here as well, and it works for both 1 and 2 particles.

---

RE-EDIT: alas, the approach is correct but, if the former approach got too much time to compute, this one consumes the whole RAM (I have 16GB), as for large dimensions the matrix sizes grow up to $10^{5}\times 10^{5}$, which means a very large RAM usage...