# Parallel mc recommendation

**URL:** https://discourse.julialang.org/t/parallel-mc-recommendation/10711
**Category:** Julia at Scale
**Created:** [May 4, 2018, 5:11pm UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711 "2018-05-04T17:11:50Z")
**Posts on this page:** 9
**Page:** 1

<div class="post-metadata">

### Author: ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)
#### Post date: [May 4, 2018, 5:11pm UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711/1 "2018-05-04T17:11:50Z")

</div>

For the following code, which runs independent chains of Random Walk Metropolis, there is the obvious potential for parallelization. I am wondering what strategies people would recommend. This is related to a prior question I had about how to use `pmap` when it was being applied to something which was an in place transformation.

```julia

function Boltzmann_likelihood(x, V, beta)
    w = exp(-beta * V(x));
    return w
end

function RWM!(X0, V::Vf, beta; niterations=10^2, Dt = 1.e-1) where {Vf}

    Xp = similar(X0);

    naccept = 0;
    p0 = Boltzmann_likelihood(X0, V, beta);

    gaussian_coef = sqrt(2 * Dt/beta);

    for j = 1:niterations

        @. Xp = X0 + gaussian_coef * randn();

        pp = Boltzmann_likelihood(Xp, V, beta);
        a = min(1, pp/p0);
        #a = 1.

        if rand()<a
            naccept = naccept+1;
            @. X0 = Xp;
            p0 = pp;
        end
    end
    X0
end

nsamples = 10^2;
d= 2;
beta = 1;

function V(x)
    return (dot(x,x)-1)^2
end

X0vals = [randn(d) for i=1:nsamples];
Xbar= zeros(d);

for k = 1:nsamples
    RWM!(X0vals[k], V, beta)

    @. Xbar += X0vals[k]/nsamples;
end

println(Xbar);

```

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [May 4, 2018, 7:22pm UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711/2 "2018-05-04T19:22:00Z")

</div>

Something like (very stylized)

```julia
random_seeds = [srand() for _ in 1:5];

function do_mcmc(rng)
    # your code that does mcmc using the
    # random number generator explicitly
end

chains = pmap(do_mcmc, random_seeds)

```

That said, RWMH is the algorithm of last resort these days, I would try a variant of HMC (eg NUTS), which should be orders of magnitude faster per effective sample.

---

<div class="post-metadata">

### Author: ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)
#### Post date: [May 4, 2018, 8:12pm UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711/3 "2018-05-04T20:12:35Z")

</div>

A few things:

1. This stylization does not appear to resolve the in place action of my `RWM!` code. Is there a way around this?
2. Your suggestion just passes a different seed to each of the processes. I usually hear this approach to parallel RNG is **ok** , but not great.
3. I know that RWM isn’t great, I’m just using this as a prototype for the kinds of problems that I am interested in doing, where I might want to run multiple independent chains in parallel.

---

<div class="post-metadata">

### Author: ![Sean\_McBane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sean_mcbane/32/4580_2.png) [@Sean\_McBane](https://discourse.julialang.org/u/Sean_McBane)
#### Post date: [May 5, 2018, 12:21pm UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711/4 "2018-05-05T12:21:34Z")

</div>

As an addition to what Tamas\_Papp said, note that `Mamba.jl` implements NUTS among other algorithms, and does parallelization automatically if you have multiple Julia processes. Its plotting functionality didn’t work as advertised for me, but it’s still useful. If you can make your model fit their syntax I’d go for it.

---

<div class="post-metadata">

### Author: ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)
#### Post date: [May 5, 2018, 12:32pm UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711/5 "2018-05-05T12:32:39Z")

</div>

Tamas has [https://github.com/tpapp/DynamicHMC.jl](https://github.com/tpapp/DynamicHMC.jl), which also implements NUTS too (without need for a framework).

---

<div class="post-metadata">

### Author: ![Sean\_McBane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sean_mcbane/32/4580_2.png) [@Sean\_McBane](https://discourse.julialang.org/u/Sean_McBane)
#### Post date: [May 5, 2018, 12:44pm UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711/6 "2018-05-05T12:44:55Z")

</div>

Ah, thanks. I hadn’t come across this before but for simple problems I rather like the more straightforward approach. I wouldn’t know what to do with most of the algorithms in Mamba anyway.

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [May 7, 2018, 8:01am UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711/7 "2018-05-07T08:01:30Z")

</div>

> [@gideonsimpson](#):
>
> in place action of my RWM! code. Is there a way around this?

Presumably, you wrap the whole thing (allocating the storage for a chain, etc). That said, in my experience for anything nontrivial the log-posterior evaluation dominates in MCMC, and savings by preallocation are not worth the code complication.

> [@gideonsimpson](#):
>
> a different seed to each of the processes. I usually hear this approach to parallel RNG is ok, but not great.

OK, so use something better 😉; this is orthogonal to parallelization.

---

<div class="post-metadata">

### Author: ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)
#### Post date: [May 7, 2018, 12:00pm UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711/8 "2018-05-07T12:00:22Z")

</div>

I guess I’m thinking about problems where my random variable is a discretization of a function on some domain and my evaluation of the likelihood involves solving a PDE. To me, that seemed like a problem where I would really want to preallocate. But you’re saying that will turn out not to be the problem?

This goes beyond my simple RWM example. Suppose I am doing a study of different initial conditions to a time dependent 2D/3D PDE. I generate my initial conditions, each of which may be rather large, and then evolve each of them, independently. There, it would seem that I might want to avoid allocating (again) in order to be `pmap` compatiable.

Indeed, the RNG issue is orthogonal (or complementary). I know the `RandomNumbers.jl` includes random123 ([D.&nbsp;E.&nbsp;Shaw Research: Research](https://www.deshawresearch.com/resources_random123.html)) which is supposed to be statistically safe for multithreading/mulitprocessing environments, but the authors don’t seem to have any examples of that usage. I was hoping maybe someone here had some experience with that kind of task. My own experience has been with SPRNG ([http://www.sprng.org/](http://www.sprng.org/)) for C/C++ coding.

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [May 7, 2018, 12:07pm UTC](https://discourse.julialang.org/t/parallel-mc-recommendation/10711/9 "2018-05-07T12:07:26Z")

</div>

> [@gideonsimpson](#):
>
> seemed like a problem where I would really want to preallocate

Preallocate for the evaluation of the log posterior if that makes sense for your problem — the cleanest way I found for this is making a callable `struct`.

My point was that preallocating for the _chain_ is a very minor (if measurable at all) improvement, in most cases.

> [@gideonsimpson](#):
>
> statistically safe for multithreading/mulitprocessing environments

I am aware of the theoretical issue, but this is usually the least of my concerns when running MCMC on a nontrivial problem.
