Parallel optimization of Markov Chain Monte Carlo

Question

I'm trying to reduce the computation time of a MCMC simulation. Essentially I have a set of particles performing a random walk in a periodic random potential. The particles are independent and every 100 timesteps I perform some checks and modify the state of some of the particles.

This seemed to me a perfect case to be optimized by parallelizing the simulation, but whatever I try I always make things just a little bit worse.

Here is in short detail what I'm doing.

First I define my random potential function U[x] as an interpolation over a simple discrete time random walk

rw[L_] := Accumulate[RandomVariate[NormalDistribution[0, 1], L]]
myrw = rw[100];
myrw = Rescale[myrw];
ifun = Interpolation[
   Transpose[{Range[1, 200], Join[Reverse@myrw, myrw]}], 
   PeriodicInterpolation -> True];
U[x_] := ifun[x];

Here are the two functions needed to perform a single time step evolution of a particle on my potential (Metropolis-Hastings algorithm):

p[x_, prop_, T_] := Exp[-(U[prop] - U[x])/T];
MCMCEvo[x_, T_] :=
 Module[{prop},
  prop = RandomVariate[NormalDistribution[x, 0.1]];
  Return[If[RandomReal[] < p[x, prop, T], prop, x]]
  ]

Then I have a function that checks the state of the particles and with some probability modifies it and then performs 100 time steps of evolution for the whole population. It then returns the state of the particles as a function of this 100-time steps evolution.

PopulationEvolve[pop_, pdiv_, plife_] := Module[{newpop},
  (*I leave this SomeFunctionOfTheStates here for the sake of completeness, but commenting it and setting newpop=pop doesn't change anything in the timings*)
  newpop = SomeFunctionOfTheStates[#,pdiv,plife]&/@pop;
  Return[(Join[{#[[1]], #[[2]]}, {Mean[U /@ #], Last[#]} &[NestList[MCMCEvo[#, 0.1] &, If[#[[2]] == 0, RandomReal[{0, 100}], #[[4]]],100]]])& /@ newpop]
]

Then I generate 2 populations, one with 100 elements and one with 400 to perform some tests:

pop100 = Table[{1, 0, 0, RandomReal[{1, 100}]}, {i, 100}];
pop400 = Table[{1, 0, 0, RandomReal[{1, 100}]}, {i, 400}];

and try parallelized and unparallelized calculations (4-Cores i5-2.6GHz):

AbsoluteTiming[PopulationEvolve[pop400, 0.8, 1]][[1]]
AbsoluteTiming[Table[PopulationEvolve[pop100, 0.8, 1], {i, 4}]][[1]]
AbsoluteTiming[ParallelTable[PopulationEvolve[pop100, 0.8, 1], {i, 4}]][[1]]
0.952207
0.960080
1.088009

and these proportions doesn't change even if I make PopulationEvolve perform 1000 time steps instead of 100: in principle I would expect the ratios to change in favor of the parallelized version, since the ratio between data exchanged between the kernels and length of the computation changes.

(*With 1000 time steps*)
AbsoluteTiming[PopulationEvolve[pop400, 0.8, 1]][[1]]
AbsoluteTiming[Table[PopulationEvolve[pop100, 0.8, 1], {i, 4}]][[1]]
AbsoluteTiming[ParallelTable[PopulationEvolve[pop100, 0.8, 1], {i, 4}]][[1]]
9.723611
9.816607
10.676737

I'm not very experienced with parallel calculations but, since now, once I was sure that there was no interactions among the parallel calculations the only bottlenecks I've ever found were about passing too much data back and forth to the kernels, which doesn't seem to be the case in this example.

What am I missing here and what would be a good way to parallelize this?

Answer 1

Update on 2014-03-04:

This is the same problem that's described here. It is also mentioned in the documentation: the the third "possible issue" here. For more details check the these two links.

In short, the cause of the problem is that ifun will have bad performance on the subkernels and the workaround is this: DistributeDefinitions[ifun]; ParallelEvaluate[ifun = ifun;]

I'll leave the old answer below.

---

This is not a solution to the problem, nor a full answer, just tracking down the cause of the slowdown. Personally I believe it to be due to a bug.

---

This reminded me of a problem where Rule (which is supposed to be an inert head) seems to be evaluated much slower on subkernels than on the main kernel. I don't know why this happens, it's all very mysterious and slightly disturbing.

Generally, I would expect

AbsoluteTiming[expr]

and

ParallelEvaluate[AbsoluteTiming[expr], First@Kernels[]]

to take the same time to evaluate. In both cases we're measuring only the evaluation time of expr, not any time spent communicating between kernels. In both cases we're evaluating the same expression on a single kernel. If expr contains lots of Rules, the subkernel evaluation will take much longer.

If expr is the invocation of an interpolating function, it also takes much longer on a subkernel. Here's a test case:

First, make an interpolating function:

ifun = Interpolation[Accumulate@RandomVariate[NormalDistribution[], 100000]];

Launch a single subkernel and DistributeDefinitions. (Note that ParallelEvaluate doesn't automatically DistributeDefinitions.)

LaunchKernels[1]
DistributeDefinitions[ifun]

Main kernel evaluation:

AbsoluteTiming@Do[ifun[RandomReal[{1, 100000}]], {10000}]

(* ==> {0.038929, Null} *)

Subkernel evaluation

ParallelEvaluate@AbsoluteTiming@Do[ifun[RandomReal[{1, 100000}]], {10000}]

(* ==> {{4.127681, Null}} *)

The latter one is much slower!

Can anyone explain this behaviour?