Problem with low cpu utilization over mpi

[Crispin]

Hi

I'm running an analysis using uniform priors which are quick to compute, but a posterior which is expensive, so my log likelihood function looks like

def loglik(params):

    if not within_bounds(params):

        return -np.inf

    else:

        return expensive_post_likelihood_func(params)

I'm running this over MPI with a small number of walkers (16) because the burn in time is long.

Inspecting the output of top, I find that most of the time only one or two tasks are concurrently executing expensive_post_likelihood_func, rather than 16 as I would hope. I suspect this is because most proposals fail the within_bounds() prior, but the default stretch move proposer must wait for all walkers to complete the previous time step, so most cpus are idle while a couple of expensive_post_likelihood_func()s complete.

Is there a way I can modify the moves so that anything not within_bounds() is not proposed in the first place? If I were to wrap StretchMove in a class that kept calling StretchMove until the proposal passed the priors, would that work, or would it somehow break the desired properties of the proposal chain?

Or is there another way to get better cpu utilization (without increasing n_walkers or decreasing n_threads)?

thank you

Crispin

[William Heymann]

What you may want to consider is a variable transform. I usually just have emcee sample in a range of 0 to 1 for each variable and transform that to the range you want. 

I also make a slight modification to the moves. I use the de an de_snooker moves because that gave me much better convergence based on the integrated autocorrelation time and modified both of them to add the line   q[i] = q[i] % 1  to make it so that they can't generate values outside my allowed range and will instead reflect back at the boundaries. I have not done testing to see if that is better than just setting to min max instead of reflection.

This is the move set I use based on one of the papers linked to on the emcee site and found I get about 10x better convergence than the stretch move.

moves=[
                (DESnookerMove(), 0.1),
                (DEMove(), 0.9 * 0.9),
                (DEMove(gamma0=1.0), 0.9 * 0.1),

--
You received this message because you are subscribed to the Google Groups "emcee users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to emcee-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/31749a6d-13e0-4540-86ca-b9b0909ce53bn%40googlegroups.com.

[Crispin]

Thanks you, I'll check out the DE moves as you suggest.

I'm reluctant to transform priors from the range [0,1] as I have (among others) two priors a,b which must satisfy a>0, b>0, a+b<1. While I could reparameterize this as b=a*c this will cause parameter interactions between a and c. Therefore for the time being I am repeating proposals until they fit the prior, as below. As StretchMove doesn't contain any tuning I'm hoping this shouldn't affect the properties of the sample chain, though if the move is something that gets tuned I suppose it might?

class PriorRejectMoveWrapper(RedBlueMove):
    def __init__(self,move,prior,maxattempts=30, **kwargs):
        self.move = move
        self.prior = prior
        self.maxattempts = maxattempts
        super(PriorRejectMoveWrapper,self).__init__(**kwargs)
    
    def get_proposal(self,current_walkers,other_group,random):
        proposals = []
        factors = []
        for walker in current_walkers:
            for _ in range(self.maxattempts):
                proposal,factor = self.move.get_proposal(np.array([walker]),other_group,random)
                assert proposal.shape[0]==1
                assert factor.shape==(1,)
                if self.prior(proposal[0]) > -np.inf:
                    break
            proposals += [proposal[0]]
            factors += [factor[0]]
        return np.array(proposals),np.array(factors)

# ... later ...

sampler = emcee.EnsembleSampler(n_walk, n_dim, ll_func, pool=pool, moves=PriorRejectMoveWrapper(StretchMove(),prior))

Crispin

[Dan Foreman-Mackey]

This won't give you correct samples from your posterior so godspeed if you want to use this in practice! 

The suggestion about a change of variables is a good one, and I think it might help. I'm not sure I follow your argument about "parameter interactions" and you can add a Jacobian term to account for any change of variables and not change your target posterior. 

To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/f19bd879-2f54-4fb1-9937-26861d32f2ean%40googlegroups.com.

--

Dan Foreman-Mackey

Associate Research Scientist

Flatiron Institute

https://dfm.io

[Crispin]

Thanks. How about immu's suggestion above - "  q[i] = q[i] % 1  to make it so that they can't generate values outside my allowed range and will instead reflect back at the boundaries "

Will altering proposals in this way preserve posterior likelihood?

Crispin

[immu...@gmail.com]

So you are aware in my case I have tested the following 3 cases

1) Reflection q[i] = q[i] % 1

2) Clipping  q = numpy.clip(q, 0, 1)

3) No changes to proposal and instead check bounds and return -np.inf if outside bounds

For my system I don't get any difference in the posterior and the first two run within margin of error with each other but much faster than the third case. I don't know if this is true in all cases but I did test it for my case at least and it did help substantially by eliminating so many items being proposed that are not possible.

[Dan Foreman-Mackey]

Good question. You need to be a little careful with this and make sure that you check that your proposals are correctly reversible. I looked into something like this in the past, and I think the conclusion I came to was that it would be ok for the DEMove but not the StretchMove, but I wouldn't hold myself to that! I think that you'd be better off using a change of variables. Here's how I would do it for your specific constraints: https://gist.github.com/dfm/5a3daaa835ce21aa704a7b5b2d8d948f

To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/5f4d7233-8caa-4ba4-9983-4acc0dcbf7fan%40googlegroups.com.