ParallelDo barely faster than Do

Question

I am new to parallel computing, so I am trying to get a grasp as to how this works. Many posts seem to indicate that overhead (like here) is a common pitfall, but I do not see why this would be the case here. Here are the following two codes, one using Do, and the other ParallelDo

Set-up:

lmax = 5;
expr =((1 + f0)*q*Cos[ϕ]*Sin[θ])/(Sqrt[f0]*r0) - 
      ((Cos[ϕ]*Sin[θ] + f0*Cos[ϕ]*Sin[θ])*
      (2*q*Cos[θ]^2 - q*Sin[θ]^2))/(Sqrt[f0]*r0) + 
      ((-1 + f0)*Cos[ϕ]*Sin[θ]*(q*Cos[θ]^2*Cos[ϕ]^2 - 
      2*q*Cos[ϕ]^2*Sin[θ]^2 + q*Sin[ϕ]^2))/(Sqrt[f0]*r0) + 
      (2*Sqrt[f0]*Cos[ϕ]*Sin[θ]*(q*Cos[ϕ]^2 + 
      q*Cos[θ]^2*Sin[ϕ]^2 - 2*q*Sin[θ]^2*Sin[ϕ]^2))/r0;

Do code:

result = Table[0, {l, 0, lmax}, {m, 0, l}];
(Do[result[[l + 1, m + 1]] = 
    Integrate[(expr SphericalHarmonicY[l, m, θ, ϕ] Sin[θ]), 
       {θ, 0, π}, {ϕ, 0, 2 π}], {l, 0, lmax}, {m, 0, l}]) // AbsoluteTiming
 result  // Simplify // MatrixForm

On my computer, this takes about ~50 seconds.

ParallelDo code:

result = Table[0, {l, 0, lmax}, {m, 0, l}];
SetSharedVariable[result]
(ParallelDo[result[[l + 1, m + 1]] = 
    Integrate[(expr SphericalHarmonicY[l, m, θ, ϕ] Sin[θ]), 
       {θ, 0, π}, {ϕ, 0, 2 π}], {l, 0, lmax}, {m, 0, l}]) // AbsoluteTiming
 result  // Simplify // MatrixForm

This takes ~37 seconds. I have $ProcessorCount = 4, and Mathematica sets $KernelCount=4 by default when I ask it to parallelize something. why is the gain so little? Is there a way to make the computation faster?

---

Running the modified code proposed by @ShinKim answer, the computation is in my case barely faster; about ~ 31 seconds.

Following the advice of @HenrikSchumacher and @SjoerdSmit, I attempted:

(result = ParallelTable[
Integrate[(expr SphericalHarmonicY[l, m, θ, ϕ] 
Sin[θ]), {θ, 0, π}, 
{ϕ, 0, 2 π}], {l, 0, lmax}, {m, 0, l}]) // AbsoluteTiming

This also takes ~39 seconds to complete.

Note: I always use Quit[] before trying each piece of code for a fair comparison.

Answer 1

As usual, using a better algorithm vastly beats parallelization.

Define a function that integrates trigonometric functions over the sphere: using pattern matching instead of Integrate,

sphericalintegrate[a_] := 
  2π^2 Expand[TrigToExp[a]] /. {E^(u_ θ + v_ ϕ) -> 0,
                                E^(v_ ϕ) -> 0,
                                E^(u_ θ) /; OddQ[u/I] -> -2/(π u),
                                E^(u_ θ) -> 0}

Let's compare:

computeresult[l_, m_] := Integrate[
  expr SphericalHarmonicY[l, m, θ, ϕ] Sin[θ], {θ, 0, π}, {ϕ, 0, 2π}]

fastcomputeresult[l_, m_] := sphericalintegrate[
  expr SphericalHarmonicY[l, m, θ, ϕ] Sin[θ]]

lmax = 5;
result = Table[computeresult[l, m], {l, 0, lmax}, {m, 0, l}]; //AbsoluteTiming//First
(*    30.8732    *)

fastresult = Table[fastcomputeresult[l, m], {l, 0, lmax}, {m, 0, l}]; //AbsoluteTiming//First
(*    0.922231    *)

result == fastresult // FullSimplify
(*    True    *)

Answer 2

Reduce the sizes of computation 'chunks' to distribute the workload more optimally by defining index set:

ind = Flatten[Table[{l, m}, {l, 0, lmax}, {m, 0, l}], 1];
result = Table[0, {l, 0, lmax}, {m, 0, l}];
SetSharedVariable[result];
ParallelDo[
    {l, m} = ind[[i]];
    result[[l + 1, m + 1]] = Integrate[
        (expr SphericalHarmonicY[l, m, \[Theta], \[Phi]] Sin[\[Theta]])
    , {\[Theta], 0, \[Pi]}, {\[Phi], 0, 2 \[Pi]}]
, {i, Length@ind}] // AbsoluteTiming
ClearAll[ind];

(* {27.4837, Null} *)

This makes it about 7.6 times only $\times 2$ faster. My PC is nothing fancy, but has $ProcessorCount = 4 and $KernelCount = 8.

---

UPDATE. I can confirm the 'caching' behavior in @Michael's answer. Here I improved the parallelization that gives consistent performance independent of caching. First, FullSimplifying the expr(in fact I just used expr*Sin[θ]) apparently cuts down the computation time in half. Also I further reduced the computation size by separating the integration domains w.r.t. the variables:

CloseKernels[];LaunchKernels[];
lmax = 5;
expr = ((1 + f0)*q*Cos[ϕ]*Sin[θ])/(Sqrt[f0]*r0) -
       ((Cos[ϕ]*Sin[θ] + f0*Cos[ϕ]*Sin[θ])*(2*q*Cos[θ]^2 - q*Sin[θ]^2))/(Sqrt[f0]*r0) +
       ((-1 + f0)*Cos[ϕ]*Sin[θ]*(q*Cos[θ]^2*Cos[ϕ]^2 - 2*q*Cos[ϕ]^2*Sin[θ]^2 + q*Sin[ϕ]^2))/(Sqrt[f0]*r0) +
       (2*Sqrt[f0]*Cos[ϕ]*Sin[θ]*(q*Cos[ϕ]^2 + q*Cos[θ]^2*Sin[ϕ]^2 - 2*q*Sin[θ]^2*Sin[ϕ]^2))/r0;
expr = FullSimplify[expr*Sin[θ]];


(*prep for parallelization*)
ind = RandomSample@Flatten[Table[{l, m}, {l, 0, lmax}, {m, 0, l}], 1];
comps1 = {};
comps2 = {};
SetSharedVariable[comps1, comps2];


(*parallel compute*)
ParallelDo[
    {l, m} = ind[[i]];
    val = Integrate[(expr*SphericalHarmonicY[l, m, θ, ϕ]), {θ, 0, Pi}];
    AppendTo[comps1, {l + 1, m + 1, val}];
, {i, Length@ind}] // AbsoluteTiming;
t1 = First@%;
ParallelDo[
    {l, m, func} = comps1[[i]];
    val = Integrate[func, {ϕ, 0, 2 Pi}];
    AppendTo[comps2, {l, m, val}];
, {i, Length@ind}] // AbsoluteTiming;
First@% + t1

(* 10.0961 *)


(*make result*)
result = GatherBy[Sort[comps2], First][[All, All, 3]];
result // Simplify // MatrixForm

ClearAll[ind,comps1,comps2,t1];

This makes it about 5 times faster. IMHO this is an acceptable multiplier provided that $KernelCount = 8.

Answer 3

[System: V13.0.0, Mac M1 Max, 10 cores]

You can observe a rough analysis with the menu command Evaluation > Parallel Kernel Status.... We can do a bit better by timing each integral. It will be seen that the integrals take relatively long time, and the overhead of parallelization of just about any parallel setup won't be very significant. The real overhead cost here is the initialization of the Integrate code on each kernel. (If I run @ShinKim's SharedVariable method twice, the second run takes less than a second, unless I permute ind, in which case, it runs only a second or two faster than the original run. This shows, I think, that in addition to whatever initialization there is, results or parts of them are cached.) Since the longer integrals take 3-5s, a kernel that happens to get two such integrals might take 10s to get its job done, though in fact, as luck would have it, it's more like 8-9s. If you have only 4 kernels, then some are going to get 4 integrals, maybe four long ones.

Now, the interesting part is the cacheing. I don't know what is cached (or perhaps just autoloaded and initialized). If we do the same analysis on the unparallelized integrations, we see only a couple of long integrals and the rest are faster. Thus, in parallelizing, we lose the benefit of this cacheing (if that is what it is).

Before each run:

Quit[]

LaunchKernels[];

lmax = 4; (* reduced to save time :) *)
expr = ((1 + f0)*q*Cos[ϕ]*Sin[θ])/(Sqrt[f0]*r0) - 
       ((Cos[ϕ]*Sin[θ] + f0*Cos[ϕ]*Sin[θ]) *
         (2*q*Cos[θ]^2 - q*Sin[θ]^2)) / 
         (Sqrt[f0]*r0) + 
       ((-1 + f0)*Cos[ϕ]*Sin[θ] * 
         (q*Cos[θ]^2*Cos[ϕ]^2 - 
          2*q*Cos[ϕ]^2*Sin[θ]^2 + 
          q*Sin[ϕ]^2)) / (Sqrt[f0]*r0) + 
       (2*Sqrt[f0]*Cos[ϕ]*Sin[θ] * 
         (q*Cos[ϕ]^2 + 
          q*Cos[θ]^2*Sin[ϕ]^2 - 
          2*q*Sin[θ]^2*Sin[ϕ]^2)) / r0;

With ParallelTable (no SharedVariable). Only the outer loop is parallelized, which means one kernel gets 5 integrals (and more of the longer ones):

(times = Flatten[
    ParallelTable[
     (Integrate[expr SphericalHarmonicY[l, m, θ, ϕ] Sin[θ],
       {θ, 0, π}, {ϕ, 0, 2 π}];
     $KernelID) // 
      AbsoluteTiming,
     {l, 0, lmax}, {m, 0, l}],
    1]) // AbsoluteTiming
(*
{16.955, {{1.0218, 10}, {1.05096, 9}, {4.39347, 9}, {1.08692, 
   8}, {4.31533, 8}, {3.46061, 8}, {1.0968, 7}, {4.706, 7}, {3.45951, 
   7}, {3.59418, 7}, {1.05576, 6}, {4.38848, 6}, {3.4545, 
   6}, {3.47644, 6}, {3.29451, 6}}}
*)

Merge[Apply[Rule]@*Reverse /@ times, List]
Total@*Flatten /@ %
(*  $KernelID -> times of integrals (not all kernels used)
<|10 -> {{1.0218}},
   9 -> {{1.05096, 4.39347}}, 
   8 -> {{1.08692, 4.31533, 3.46061}}, 
   7 -> {{1.0968,  4.706,   3.45951, 3.59418}}, 
   6 -> {{1.05576, 4.38848, 3.4545,  3.47644, 3.29451}}|>

<|10 -> 1.0218, 9 -> 5.44443, 8 -> 8.86287, 7 -> 12.8565, 6 -> 15.6697|>
*)

@ShinKim's method with SharedVariable; the table just returns the time and kernel of each result.

ind = Flatten[Table[{l, m}, {l, 0, lmax}, {m, 0, l}], 1];
result = Table[0, {l, 0, lmax}, {m, 0, l}];
SetSharedVariable[result];
(times = ParallelTable[
    ({l, m} = ind[[i]];
      result[[l + 1, m + 1]] = 
       Integrate[expr SphericalHarmonicY[l, m, θ, ϕ] Sin[θ],
        {θ, 0, π}, {ϕ, 0, 2 π}];
      $KernelID) // 
     AbsoluteTiming,
    {i, Length@ind}]) // AbsoluteTiming
(*
{9.12166, {{1.3988, 10}, {0.336149, 10}, {5.10557, 9}, {0.597471, 
   9}, {5.03329, 8}, {3.36767, 8}, {1.4241, 7}, {4.75394, 
   7}, {4.78472, 6}, {3.54514, 6}, {1.44249, 5}, {5.10244, 
   4}, {4.73389, 3}, {4.82558, 2}, {4.79167, 1}}}
*)

Merge[Apply[Rule]@*Reverse /@ times, List]
Total@*Flatten /@ %
(*  $KernelID -> times of integrals
<|10 -> {{1.3988,  0.33615}},  9 -> {{5.10557, 0.597471}}, 
   8 -> {{5.03329, 3.36767}},  7 -> {{1.4241,  4.75394}}, 
   6 -> {{4.78472, 3.54514}},  5 -> {{1.44249}},
   4 -> {{5.10244}},           3 -> {{4.73389}},
   2 -> {{4.82558}},           1 -> {{4.79167}}|>

<|10 -> 1.73495, 9 -> 5.70304, 8 -> 8.40097, 7 -> 6.17804, 6 -> 8.32986,
   5 -> 1.44249, 4 -> 5.10244, 3 -> 4.73389, 2 -> 4.82558, 1 -> 4.79167|>
*)

Unparallelized:

(times = Flatten[
    Table[(
     Integrate[(expr SphericalHarmonicY[l, m, θ, ϕ] Sin[θ]),
      {θ, 0, π}, {ϕ, 0, 2 π}];
     $KernelID) // 
      AbsoluteTiming,
     {l, 0, lmax}, {m, 0, l}],
    1]) // AbsoluteTiming
(*
{20.596, {{0.923452, 0}, {0.242653, 0}, {4.09701, 0}, {0.301539, 
   0}, {0.865541, 0}, {3.20008, 0}, {0.291668, 0}, {1.18621, 
   0}, {0.579081, 0}, {3.30374, 0}, {0.254854, 0}, {0.894277, 
   0}, {0.600397, 0}, {0.670284, 0}, {3.18521, 0}}}
*)

Merge[Apply[Rule]@*Reverse /@ times, List]
Total@*Flatten /@ %
(*
<|0 ->
  {{0.923452,
    0.242653, 4.09701,                  <-- head of column is longest
    0.301539, 0.865541, 3.20008,          ↙  ↓
    0.291668, 1.18621,  0.579081, 3.30374, 
    0.254854, 0.894277, 0.600397, 0.670284, 3.18521}}|>

<|0 -> 20.596|>
*)

My guess is that each new value of m causes something to be cached or code to be initialized. You lose this efficiency in parallelization. For instance, switching the iterators in the first ParallelTable above to {m, 0, lmax}, {l, m, lmax} cuts the time in half.