Data parallelism for independent probabilistic trials

Question

I'm considering the best way to implement the data parallelism features of Mathematica for a specific task. To that end I have a minimal code example below which contains the basic features of my actual code. The example code does the following basic task: Consider defining an initial vector of length $2^N$ (as shown only the first entry of the vector is non-zero): $$\begin{pmatrix} ~~~1&\\ ~~~0 &\\ ~~~\cdot &\\ ~~~\cdot &\\ ~~~0 \end{pmatrix}$$ The vector then undergoes a collective rotation given by the operator $\hat{J}_x := \frac{1}{2}\sum^{N}_{i}\hat{\sigma}_{i}^{x}$ (where $\hat{\sigma}_{i}^{x}$ is the Pauli matrix) about the $x$-axis where the rotation angle is probabilistically (uniformly) taken from the interval $[0,2\pi]$. Then I take the absolute value squared of the inner-product of the rotated state with the initial state. Therefore the process yields a real number at the end as a result. Since I am interested in the average of these independent probabilistic trial results this seems like a process which could benefit from data parallelism.

The main code which does the probabilistic rotation is given by (for example $N=4$):

NotebookDirectory[]
Dimension Number
Num = 4;
Identity and Pauli-X Matrix
Id = {{1, 0}, {0, 1}}; SigX = {{0, 1}, {1, 0}};

Collective X-Rotation Operator 
A = Table[Id, {i, Num - 1}]; PrependTo[A, SigX];
Pe = Permutations[A];
KP = Table[KroneckerProduct @@ Pe[[i]], {i, Length[A]}];
JX = (1/2)*Sum[KP[[i]], {i, 1, Length[A]}];

Initial Vector
B1 = ConstantArray[{0}, 2^Num - 1]; PrependTo[B1,  {1}];

Random Rotation Angle
R = RandomReal[{0, 2 \[Pi]}]

B2 = Dot[MatrixExp[-I (R) *JX], B1];
Flatten[(Abs[ConjugateTranspose[B2] . B1])^2][[1]]

and the separate code which calls the first code and stores the resultant $10$ independent trial results and calculates the mean value:

Mean[Table[NotebookEvaluate["C:the directory of the above file.nb"], 10]]  

Can anyone advise on basic parallelism schemes in Mathematica which would improve the efficiency of the outlined task? Thanks.

Answer 1

There are a few things to do to simplify your code. First, use some more Mathematica-style idioms. So,

Num = 4;
A = ConstantArray[IdentityMatrix[2], Num - 1]~Join~{PauliMatrix[1]};
JX = (1/2) Total[KroneckerProduct @@@ Permutations[A]];
B1 = {1}~Join~ConstantArray[0, 2^Num - 1];

Notice that I changed the shape of B1 from a List of Lists to just a single list, because Mathematica doesn't distinguish columns from rows. I also smoothed out some of the looping for implicit functional and array programming, for speed.

I will draw 10000 random variates from the range of R (this is the number of times you want to run your calculation, 10 in your example above),

R = RandomVariate[UniformDistribution[{0, 2 \[Pi]}], 10000];

and finally, the biggest trick for speed is using the two-input mat-vec version of MatrixExp,

overlap = Map[Re[B1\[ConjugateTranspose] . MatrixExp[-I # JX, B1]] &, R];

On my machine this evaluates almost instantaneously for Num=4. For Num=8, it took about two minutes. (If your "real example" is more involved, you can do ParallelMap.) (I added the Re because sometimes you get answers with imaginary parts that are 1e-16 or something that's coming from numerical imprecision.)

You can convince yourself that the exact answer for this example is that the overlap depends on the rotation angle via Cos[x/2]^Num. We can visualize,

Show[
    Plot[Cos[x/2]^Num, {x, 0, 2 \[Pi]}, PlotStyle -> {Black}],
    ListPlot[{R, overlap}\[Transpose],PlotStyle->{Blue}]
 ]

which gives perfect agreement

and do Mean[overlap] to get 0.379094 (in my case, yours will differ because of the RandomVariate draws). The exact answer is 1/(2 \[Pi]) Integrate[Cos[x/2]^Num, {x, 0, 2 \[Pi]}] = 3/8 = 0.375 when Num=4, but we did a monte carlo calculation with only 10000 samples so you should only expect a precision in the 1%ish range, which is right.

Further Acceleration

• If things get very big, you may benefit from making J a SparseArray, which is accepted by MatrixExp. However, if you really intend to integrate over the whole Haar measure of angles, you probably won't have much sparseness.

• You may benefit from Compile, but I am not sure MatrixExp is a compilable function.

• You may benefit from ParallelMap, depending on how many cores and licenses you have.