# How can we speed up a matrix multiplication using ParallelTable and ParallelISum?

I do belive that Mathematica can work as fast as other packages, but it is really sensitive to the way we write the code. Here I would like to evaluate a trace resulted from a matrix multiplication at many points as follwos

KS[x_, y_,z_] := {{0, I Sin[x] + Sin[y],
13/4 - Cos[x] - Cos[y] - Cos[z] - I Sin[z], -1}, {-I Sin[x] +
Sin[y], 0, -1,
13/4 - Cos[x] - Cos[y] - Cos[z] - I Sin[z]}, {13/4 - Cos[x] -
Cos[y] - Cos[z] + I Sin[z], -1, 0, -I Sin[x] - Sin[y]}, {-1,
13/4 - Cos[x] - Cos[y] - Cos[z] + I Sin[z], I Sin[x] - Sin[y], 0}};
Oz[x_, y_,z_] :=
{{0, 0, 0, 1/2 (-Cos[z] - I Sin[z])}, {0, 0,
1/2 (Cos[z] + I Sin[z]), 0}, {0, 1/2 (Cos[z] - I Sin[z]), 0,
0}, {1/2 I (I Cos[z] + Sin[z]), 0, 0, 0}};
Ox[x_, y_,z_] := {{0, I Cos[x], Sin[x], 0}, {-I Cos[x], 0, 0,
Sin[x]}, {Sin[x], 0, 0, -I Cos[x]}, {0, Sin[x], I Cos[x], 0}};
F2[x_, y_, z_,
r_] := (Inverse[(r + I*0.01)*IdentityMatrix[4] - KS[x, y, z]]);
F1[x_, y_, z_,
r_] := (Inverse[(r - I*0.01)*IdentityMatrix[4] - KS[x, y, z]]);
Myfun[x_, y_, z_, r_] :=
Re[Tr[(Ox[x, y, z].(F2[x, y, z, r] - F1[x, y, z, r]).Oz[x, y,
z].(F2[x, y, z, r].F2[x, y, z, r] +
F1[x, y, z, r].F1[x, y, z, r]))]];


Then I used ParallelSum with Table to do the job, but it is very slow> Is there any way to speed it up?

nr=5;nc=9;
Table[ParallelSum[(nc/(2\[Pi]))^-3 1/nr Myfun[x,y,z,r],{r,-7.25,R,1/nr},{x,\[Pi]/nc,\[Pi],(2\[Pi])/nc},
{y,\[Pi]/nc,\[Pi],(2\[Pi])/nc},{z,\[Pi]/nc,\[Pi],(2\[Pi])/nc}],{R,-1,1,0.5}]//AbsoluteTiming
{9.9159,{1.64289,-1.40237,-1.40319,-1.40314,0.10968}}


I want to do it with nr=100 and nc=200.

• (1) Use nr = 5.; nc = 9.; in other words, use inexact values for nr and nc (2) let Parallelize choose the best approach by doing Parallelize[Table[Sum[.... This brings it down to 2.4s on my machine. (3) Eliminate some repeated calculations of F1 and F2 in Myfun by doing Myfun[x_, y_, z_, r_] := Block[{f1 = F1[x, y, z, r], f2 = F2[x, y, z, r]}, Re[Tr[(Ox[x, y, z].(f2 - f1).Oz[x, y, z].(f2.f2 + f1.f1))]]]; This cuts it down to just over a second. Commented Aug 14, 2020 at 21:03
• I can get it a little lower to 0.8 seconds, if I remove F1 and F2 and roll them into a single function to reduce recalculation of KS, as F1F2[x_, y_, z_, r_] := Block[{ks = KS[x, y, z], id = IdentityMatrix[4]}, {(Inverse[(r - I*0.01)*id - ks]), (Inverse[(r + I*0.01)*id -ks])}] then Myfun[x_, y_, z_, r_] := Block[{f = F1F2[x, y, z, r]}, Re[Tr[(Ox[x, y, z].(f[[2]] - f[[1]]).Oz[x, y, z].(f[[2]].f[[2]] + f[[1]].f[[1]]))]]]; then apply (1) and (2) from my previous comment. Commented Aug 14, 2020 at 21:13
• @flinty I would really appreciate if you can post it as an answer please so I can also accept it? Commented Aug 14, 2020 at 21:16
• Is 0.08 seconds on a couple year old laptop sufficiently quick?
– ciao
Commented Aug 14, 2020 at 21:16
• @ciao, how did you do it, the results I posted were performed with a workstation with 20 core and 120G as RAM-:) So, 0.08 sec would be AWESOME Commented Aug 14, 2020 at 21:18

You can compile your function to speed up the evaluation:

MyfunCompiled = Compile[{{x, _Real}, {y, _Real}, {z, _Real}, {r, _Real}},
Evaluate@Re[Tr[(Ox[x, y, z].(F2[x, y, z, r] - F1[x, y, z, r]).Oz[x, y, z].
(F2[x, y, z, r].F2[x, y, z, r] + F1[x, y, z, r].F1[x, y, z, r]))]],
CompilationTarget -> "C"];

nr = 5; nc = 9;
AbsoluteTiming@Sum[(nc/(2 \[Pi]))^-3 1/nr Myfun[x, y, z, r],
{r, -7.25, 1, 1/nr}, {x, π/nc, π, 2π/nc}, {y, π/nc, π, 2π/nc}, {z, π/nc, π, 2π/nc}]

AbsoluteTiming@Sum[(nc/(2 \[Pi]))^-3 1/nr MyfunCompiled[x, y, z, r],
{r, -7.25, 1, 1/nr}, {x, π/nc, π, 2π/nc}, {y, π/nc, π, 2π/nc}, {z, π/nc, π, 2π/nc}]

{4.8844, 0.10968}
{0.050282, 0.10968}


For further improvements, you can also compile the actual sum:

FunSum = Compile[{{R, _Real}, {nr, _Integer}, {nc, _Integer}},
Sum[(nc/(2π))^-3 1/nr MyfunCompiled[x, y, z, r],
{r, -7.25, 1, 1/nr}, {x, π/nc, π, 2π/nc}, {y, π/nc, π, 2π/nc}, {z, π/nc, π, 2π/nc}]
,CompilationTarget -> "C"];

nr = 5; nc = 9;
AbsoluteTiming@FunSum[1, nr, nc]

{0.007274, 0.10968}


Your example would then look like this (following the suggestion of flinty in comments about parallelizing over the different values of R):

nr = 5; nc = 9;
ParallelTable[FunSum[R, nr, nc], {R, -1, 1, 0.5}] // AbsoluteTiming

{0.027832, {1.64289, -1.40237, -1.40319, -1.40314, 0.10968}}


This is now fast enough that you can try going for larger values of nr and nc. For example:

AbsoluteTiming[FunSum[1, 50, 50]]
AbsoluteTiming[FunSum[1, 100, 100]]
AbsoluteTiming[FunSum[1, 100, 200]]

{5.55616, 0.452771}
{73.7241, 0.484847}
{678.142, 0.48616}


which seems to converge (at least it looks like it does).

• that is Great!, by the way, will it help if we remove F1 and F2 and roll them into a single function to reduce recalculation of KS,as suggested by "flinty" in the comments, and combine it with your approach ? Commented Aug 15, 2020 at 18:04
• It is unfortunately not directly applicable to my method, since I am basically compiling a single, long function, and all the higher-level operations like Inverse are precomputed. I did just try a version where I only precompute F1, F2, Ox and Oz, and the compiled function then assembles the Re[Tr[Dot[]]], but there is essentially no improvement in the timings. Commented Aug 15, 2020 at 19:06