11
$\begingroup$

Difficulty in parallelizing Packages using DistributeDefinitions

Context

I am trying to parallelize some parameter estimation using MCMC.

As a test run I expect it to find the width of a Gaussian.

If I use the built in Mathematica function FindDistributionParameters, it parallelizes nicely:

Table[x = RandomVariate[NormalDistribution[0, 2/10], 50000]; 
  FindDistributionParameters[x, NormalDistribution[0, s]], {i, 
   40}] // Timing

{0.324758, {{s -> 0.200035},…}

runs much slower than

ParallelTable[x = RandomVariate[NormalDistribution[0, 2/10], 50000]; 
  FindDistributionParameters[x, NormalDistribution[0, s]], {i, 
   40}] // Timing

{0.016444, {{s -> 0.199899}, …}

If I now use a MCMC package (following this answer)

<<MCMC.m

Unprotect[Monitor];
Monitor = # &;
Protect[Monitor];
DistributeDefinitions[MCMC];

For some reason this (non parallel) implementation is as slow

Table[x = RandomVariate[NormalDistribution[0, 2/10], 500]; 
  f[s_] = LogLikelihood[NormalDistribution[0, s], x]; 
  MCMC[f[s], {{s, 0.1, 0.4, Reals}}, 1000]["BestFitParameters"]
  , {i, 4}] // Timing

{2.26298,{{s->0.188522},{s->0.208759},{s->0.208809},{s->0.1983}}

as that (parallel) implementation:

ParallelTable[x = RandomVariate[NormalDistribution[0, 2/10], 500]; 
  f[s_] = LogLikelihood[NormalDistribution[0, s], x]; 
  MCMC[f[s], {{s, 0.1, 0.4, Reals}}, 1000]["BestFitParameters"]
  , {i, 4}] // Timing

{3.38599,{{s->0.199691},{s->0.200183},{s->0.206191},{s->0.201176}}

Question

Why does ParallelTable not scale up on this problem given that the different jobs are fully independent?

PS: note that an alternative option could involve parallelizing the function itself, following e.g. this answer

$\endgroup$
5
  • 1
    $\begingroup$ I recently found that in using packages I would need to use Needs[] as well as ParallelNeeds[] to have the definitions distributed correctly to the Kernels. $\endgroup$
    – gwr
    Commented Mar 3, 2017 at 15:58
  • $\begingroup$ Suprised the heck out of me, too... just when one thought, one knew all the tricks. :) $\endgroup$
    – gwr
    Commented Mar 3, 2017 at 16:02
  • $\begingroup$ When I do this with MCMC, I get Min::nord and FrontEndObject::notavail (albeit getting the correct result). Do you get the same messages? $\endgroup$
    – gwr
    Commented Mar 3, 2017 at 16:17
  • 1
    $\begingroup$ @gwr yes unless I do ParallelDo[Unprotect[Monitor]; Monitor = # &; Protect[Monitor]; , {4}]; because the function uses Monitor. $\endgroup$
    – chris
    Commented Mar 3, 2017 at 16:36
  • $\begingroup$ I'm not an expert but "parallel" and "random numbers" ring some bells. I think there are some peculiarities with that combination. If unsure you might want to make sure that what you do is correct... $\endgroup$ Commented Mar 4, 2017 at 20:19

1 Answer 1

10
$\begingroup$

Loading Packages with ParallelNeeds

As Szabolcs has nicely posted here distributing definitions to the Kernels by using DistributeDefinitions like so:

DistributeDefinitions[ "PackageContext`"] 

will not work. Since definitions of the default context are distributed automatically, one needs to make sure, that packages loaded via Needs also "reach" the Kernels.

This can be achieved by using ParallelNeeds. But one still needs to load the package into the global context thus in your case:

Needs[ "MCMC`"];
ParallelNeeds[ "MCMC`"];

That should do the trick.

In practice for the problem at hand:

ParallelDo[Unprotect[Monitor];Monitor = # &; Protect[Monitor];, {4}]

Then we can compute best estimate and error on estimate as follow:

res = ParallelTable[
   x = RandomVariate[NormalDistribution[0, 2/10], 5000]; 
   f[m_, s_] = LogLikelihood[NormalDistribution[m, s], x]; 
   MCMC[f[m, s], {{m, 0.1, 0.01, Reals},
      {s, 0.1, 0.01, Reals}}, 5000][{"BestFitParameters", 
     "ParameterErrors"}]
   , {i, 32*8}];

dat = Join[{m, s} /. Map[#[[1, 2]] &, res] // 
     Transpose, {m, s} /. Map[#[[2, 2]] &, res] // Transpose] // 
   Transpose;
Needs["ErrorBarPlots`"];
ErrorListPlot[{{#1, #2}, ErrorBarPlots`ErrorBar @@ {#3, #4}} & @@@ 
  dat, PlotRangePadding -> {Scaled[0.15], Automatic}]

Mathematica graphics

Indeed the Maximum likelihood method is unbiased

 dat1 = {m, s} /. Map[#[[1, 2]] &, res]; dat1 // DensityHistogram

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.