# Data parallelism for independent probabilistic trials

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.

• Why not leave R unevaluate and compute your result, e.g., f[r_] = (Abs[ConjugateTranspose[MatrixExp[-I r JX] . B1] . B1]^2)[[1, 1]]? Commented Nov 9, 2021 at 18:42
• @CarlWoll Thanks for your response. For reasons which stems from how the above example translates into the actual problem it's best to leave it as is regarding evaluating R. I am mostly interested in how to use data parallelism to improve the efficiency. Commented Nov 9, 2021 at 18:57

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.

• If you are interested in integrating over the whole space of unitaries, that's a bit of an interesting problem, and requires yet a different approach, I think, to avoid computing absurdly large matrix exponentials. I have a small demonstration code; please reach out! You can find my URL in my profile. Commented Nov 9, 2021 at 22:43
• Many thanks for your answer and comment, there are a number of good suggestions. I'm still revising it and as a result I may need to append the question to more accurately reflect my actual problem (work in progress). Commented Nov 10, 2021 at 7:26
• One thing you may benefit from, that I got at implicitly, is to avoid making a dedicated notebook and then evaluating the notebook N times. You can craft Modules and apply them, instead. Commented Nov 10, 2021 at 7:31