Supporting Population Samplers (implemented DE-MCMC)

[michaelosthege]

With this PR I would like to propose some changes to support population samplers.

Changes

• Metropolis is marked as COMPATIBLE with continuous variables (I don't understand why this was not the case already)
• the bij variable of ArrayStepShared is elevated to an attribute
• a PopulationArrayStepShared is subclassed from ArrayStepShared
• DEMetropolis is a PopulationArrayStepShared step method that implements differential-evolution metropolis
• the sample() procedure is modified and two new methods sample_population and _iter_chains are added to support sampling with population samplers
• sample() is a bit more verbose and prints out the sampler assignment and method of sampling

Performance

I have included an example that compares Metropolis, NUTS and DEMetropolis on a correlated bivariate normal target density. It calculates the effective sampling rate and normalizes it to NUTS on the uncorrelated target density. This is the result I get on my laptop:

Effective sample size [0...1]

Method \ p	0.0	0.9
Metropolis	0.176	0.024
Slice	1.000	0.102
NUTS	0.985	0.268
DEMetropolis	0.106	0.129

Runtime [s]

Method \ p	0.0	0.9
Metropolis	21.0	21.5
Slice	63.1	57.9
NUTS	37.2	63.5
DEMetropolis	18.6	19.0

Normalized effective sampling rate [0...1]

Method \ p	0.0	0.9
Metropolis	0.317	0.042
Slice	0.598	0.067
NUTS	1.000	0.159
DEMetropolis	0.214	0.257

Usage

The DEMetropolis step method can be used just like any other stepper. It supports blocking and can be combined with other steppers.

ToDo

• get all tests to pass
• add tests for PopulationArraySharedStep and DEMetropolis

References

• Cajo J. F. Ter Braak, 2006, A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces
• Supporting multichain/population-samplers #2717

[twiecki]

@michaelosthege Looks very promising!

[twiecki]

add reference

[junpenglao]

Well done @michaelosthege! Could you please add it to the release note under new feature?

Also, you should add an EXPERIMENTAL_WARNING (see eg here https://github.com/pymc-devs/pymc3/blob/master/pymc3/step_methods/sgmcmc.py#L14-L15).

[michaelosthege]

Thank you :)
I will add a reference today.
And two more things:

1. In contrast to the original publication, I'm tuning epsilon (symmetric noise that is added to the jump). The idea is that it should help the chains to diverge from potentially identical start values. Mathematically there is no problem, but I can't say if it's a good idea either.
2. I have an experimental chain-multiprocessing in the pipeline. Will do some more testing and exception handling before I push the commits.

[twiecki]

What about the code we have for SMC? Is there any overlap between these?

[michaelosthege]

The problem is the same, but the solution is different. And performance is a concern, so it's worthwhile to think about it.

It looks like SMC uses a multiprocessing.Pool and a text file MultiTrace backend to share the population state and synchronize chains. Although I have to admit that I don't completely understand how this works in combination with the Pool.

I use Pipes to create a communication channel between the processes. That channel is then used to pass the population to the chain, and the update back to the main iterator.

My guess is that the Pipes should in principle be faster. And they certainly involve less code.

[twiecki]

My concern is code duplication in the pymc3 code-base.

[junpenglao]

Maybe we can refactor SMC at some point @hvasbath @aloctavodia?

[aloctavodia]

I will love to help with SMC, "at some point" could be during January :-)

[hvasbath]

Looks like a cool sampler! Thanks for implementing @michaelosthege . If this was merged in 1 1/2 years earlier, likely I wouldnt have started working on the SMC ;) .
The algorithms themselves are quite different. I am not sure if SMC can be refactored to support the parallel sampling here, which is the only thing the samplers have in common appart from the metropolis step ;). For sure the pool also supports python2.7, which is apparently not supported in the DEM here.
The chains in the SMC do not communicate with each other during the sampling as they do here. They communicate all together every n_chains * nsteps samples. For sure would be interesting to see the test and performance you did also compared to the SMC @michaelosthege . There is the SMC2gaussians notebook where you see how to start it. Anyways if anyone has a good idea to anyway factor out common things I also could help starting earliest January.

[twiecki]

Some conflicts creeped in @michaelosthege. I'm fine merging once resolved.

[michaelosthege]

Tomorrow I will fix that and do a little refactor. Should be ready to merge then 😊

I really like the way how blocking and CompoundStep work in pymc3. With respect to SMC, it might be nice to merge PopulationArrayStep with the ArrayStepSharedLLK into a universal, parallelized population step. (Even if information is not exchanged in every step.) /cc @hvasbath @aloctavodia

The failure of my implementation to multiprocess in Py2.7 is because even in 3.5 I couldn't get the steppers to pickle without using protocol=4. I don't understand why though, because obviously Pool can do it.

[michaelosthege]

Okay, so I'd say it's ready to merge.
I benchmarked the sampler performance against expensive forward passes. I used a single 1D Normal and a custom Op with a time.sleep() to artificially slow down the forward pass. (10 chains, 5k draws, running on my 4-core laptop)

Runtime [s]
sampler \ delay          0.0   0.5   1.0   2.0   3.0   4.0    5.0   10.0
Metropolis               7.6  44.5  43.1  60.8  84.9 107.5  129.9 238.3
Slice                   11.3 110.5 110.0 184.4 248.5 316.1 1326.4 752.7
DEMetropolis_sequential  3.2  69.3  69.2 110.1 164.9 206.6  250.1 473.9
DEMetropolis_parallel   19.5  42.3  44.6  62.5  73.8  70.2   83.7 139.6

Sampling rate [1/s]
sampler \ delay           0.0   0.5   1.0   2.0   3.0   4.0   5.0   10.0
Metropolis              2623.0 449.2 464.5 328.7 235.5 186.0 153.9  83.9
Slice                   1768.3 181.0 181.8 108.4  80.5  63.3  15.1  26.6
DEMetropolis_sequential 6252.8 288.7 288.9 181.7 121.3  96.8  80.0  42.2
DEMetropolis_parallel   1025.2 472.4 448.8 320.1 271.2 285.1 239.0 143.3

So even though the parallelized DEMetropolis is much worse with very fast forward passes, it outperforms the sequential version as soon as the forward pass takes more than a few ms.

[junpenglao]

Could you please also add it to release note under feature?

[michaelosthege]

@junpenglao Just added it. Thank you for reminding me!

Note: please Squash & Merge - there are some embarrassing bugs in the history 🙃

[junpenglao]

Great job @michaelosthege!