Randomness during adapting causing slow sampling for chains in parallel?

[Ryan]

What sort of randomness is there during the adapt phase? Could this be a cause for slowness in sampling between chains?

I have two very different types of models that I've encountered slowness in some chains when sampling from 3 chains in parallel. Both models have somewhere between 50k-100k sample sizes, some form of a normal likelihood, so nothing too crazy in either case, yet in both cases I find that one chain finishes sampling much faster than the other two.

My sense is that this caused by some sort of randomness at the start of the adapt phase in each chain. Is there anyway to monitor this, or possibly set some sort of initial values in hopes of having all chains sample at the faster rate?

Thanks for any insight!

[Bob Carpenter]

On 6/30/13 6:18 PM, Ryan wrote:
> What sort of randomness is there during the adapt phase?

If you don't specify initial values, those are initialized
between -2 and 2 uniformly on the unconstrained scale (constrained
variables are transformed to unconstrained scales automatically).

Then NUTS is run during adaptation, and that has several kinds
of built-in randomness, including initial momentum and the directions
in which to update and a final Metropolis reject step.

> Could this be a cause for slowness in sampling between chains?

Yes.

> I have two very different types of models that I've encountered slowness in some chains when sampling from 3 chains in
> parallel. Both models have somewhere between 50k-100k sample sizes, some form of a normal likelihood, so nothing too
> crazy in either case, yet in both cases I find that one chain finishes sampling much faster than the other two.

We've noticed that, too, and have been considering
time-based evaluation with appropriate weighting
among the chains.

Do all of the chains get to the same place eventually
in terms of the quantiles (or mean and sd).

> My sense is that this caused by some sort of randomness at the start of the adapt phase in each chain. Is there anyway
> to monitor this, or possibly set some sort of initial values in hopes of having all chains sample at the faster rate?

Yes, you can provide initial values for each chain. If
they all start in the same place, they're much more likely
to all run at roughly the same speed, but you lose the
benefit of the conservative convergence diagnostics.

- Bob

[Ryan]

Then NUTS is run during adaptation, and that has several kinds
of built-in randomness, including initial momentum and the directions
in which to update and a final Metropolis reject step.

Ah OK, I have tried using the same initial values for all chains, which does seem to help, but I think I might be getting bit by this randomness. Are there any plans to be able to control these, or at least seed these in some way so they start the same in all chains?

 

Do all of the chains get to the same place eventually 
in terms of the quantiles (or mean and sd).

For one of the models, yes. The other has a more complicated hierarchy that doesn't seem to be well identified, so that's separate issue altogether.

Thanks for the insight Bob. Much appreciated!

[Michael Betancourt]

Ah OK, I have tried using the same initial values for all chains, which does seem to help, but I think I might be getting bit by this randomness. Are there any plans to be able to control these, or at least seed these in some way so they start the same in all chains?

If the random numbers in each chain are the same (as would be the case given identical seeds) then the dynamics would be the same and you'd just get n copies of a single chain.

The big problem here is that if you start adaptation away from the typical set and your tails are heavy then there's always a chance of the adaptation going askew.  Starting all chains at the same point helps only if that point is close to the typical set (models that seem to be "well written" are often defined such that the origin is within the typical set, which is why --init=0 performs so well most of the time).  The easiest fix right now should be running Stan without any warmup, record the final value (you should see some convergence in the trace plots), and then use that to seed a full adaptation/sampling run.

We're currently working on more robust adaptation mechanisms.

[Bob Carpenter]

On 7/1/13 9:15 AM, Ryan wrote:
>
> Then NUTS is run during adaptation, and that has several kinds
> of built-in randomness, including initial momentum and the directions
> in which to update and a final Metropolis reject step.
>
>
> Ah OK, I have tried using the same initial values for all chains, which does seem to help, but I think I might be
> getting bit by this randomness. Are there any plans to be able to control these, or at least seed these in some way so
> they start the same in all chains?

The improvements we have planned are:

* introducing an init block in the Stan language itself

* allowing partial initializations, with other parameters being random

* introducing additional randomness around inits

* controlling how much extra randomness is added to the inits

All of this can be done in R before calling Stan, but the partial inits
are particularly problematic because we currently require either all or nothing
inits that are in the support of the posterior on the constrained scale.

As Michael mentioned, we've been thinking about more robust initialization
strategies (robust in the sense of getting to the first effective sample
more quickly and finding adaptation parameters that mix well from there).

These are high enough priority that they should actually happen :-)


- Bob

[Andrew Gelman]

I assume that we can switch inits to [-.5,.5] on the very next release?

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>

[Ryan]

Thanks for all of the assistance everyone. Here's what has seemed to work well for me in regards to getting all chains to finish at about the same time (while maintaining good effective sample sizes):

- For all runs, use delta=0.25 and max_treedepth=5.

- Perform an initial small sample run from random inits using default amount of warmup (when I used no warmup, all samples were rejected).

- Use the estimates from the initial run as initial values to the full run.

My experience has been that the max_treedepth=5 helps keep NUTS from going for too long during adaption, and once sampling is started I still get good effective sample sizes. Your mileage may vary, although it has seemed to really help since I have high dimensions.

Not exactly sure how much the lower delta of 0.25 helped, but it seemed to be closer to a good value versus 0.5. That said, I found myself wanting to be able to set the initial step sizes for all variables (say with a list like initial values?) when using varying step sizes for each parameter. Perhaps this is on the TODO list?

Thanks again,

Ryan

[Bob Carpenter]

On 7/20/13 4:43 PM, Ryan wrote:
> Thanks for all of the assistance everyone. Here's what has seemed to work well for me in regards to getting all chains
> to finish at about the same time (while maintaining good effective sample sizes):
>
> - For all runs, use delta=0.25 and max_treedepth=5.

That's probably enough tree depth for most problems. It only
allows 32 leapfrog steps and some problems require more.

> - Perform an initial small sample run from random inits using default amount of warmup (when I used no warmup, all
> samples were rejected).
> - Use the estimates from the initial run as initial values to the full run.

The problem with this approach is that it doesn't use diffuse inits, which
we feel are important for convergence monitoring. The danger is
that if you start all the samples from the same point, you won't find
different modes or non-identifiability, or perhaps detect if you're
exploring enough of the parameter space. Also, if the init doesn't get
to a good place, all the chains will suffer.

We've been thinking a good strategy might be to do something like let
warmup get close to a solution, then randomize a bit. The problem is
that of scale.

> My experience has been that the max_treedepth=5 helps keep NUTS from going for too long during adaption, and once
> sampling is started I still get good effective sample sizes. Your mileage may vary, although it has seemed to really
> help since I have high dimensions.
>
> Not exactly sure how much the lower delta of 0.25 helped, but it seemed to be closer to a good value versus 0.5. That
> said, I found myself wanting to be able to set the initial step sizes for all variables (say with a list like initial
> values?) when using varying step sizes for each parameter. Perhaps this is on the TODO list?

Yes, it's on the to-do list. I'm not sure if it's in the 2.0
refactor or not. Michael's rewritten most of the command structure
to make this kind of thing easier to generalize across samplers
and optimizers.

- Bob

[Ryan]

The problem with this approach is that it doesn't use diffuse inits, which
we feel are important for convergence monitoring.  The danger is
that if you start all the samples from the same point, you won't find
different modes or non-identifiability, or perhaps detect if you're
exploring enough of the parameter space.  Also, if the init doesn't get
to a good place, all the chains will suffer.

We've been thinking a good strategy might be to do something like let
warmup get close to a solution, then randomize a bit.  The problem is
that of scale.

I wonder if instead of using the initial mean for all chains I were to use a random sample from the initial run? This wouldn't be as overdisperesed like random inits, but it also wouldn't be the same for all chains. Perhaps this could help get closer to the right place for all chains. Any ideas if this kind of strategy would be almost as good as random inits?

[Bob Carpenter]

We've been discussing various adaptation strategies like this
in our meetings and still haven't figured out the right thing
to do.

A random sample from a single chain still defeats the purpose of
overdispersed starting points, which is to provide a sanity check on the
model and its computation.

For well-behaved models, it shouldn't matter where you start.

We want to concentrate on getting from any initial point to
the typical set of the distribution so we can start sampling
and properly adapting.

But clearly we're not there yet, which is why some chains
are so much slower than others. In fact, we've found the slow
chain when running four chains often dominates the compute time.
And it makes compute time highly variable on a per-chain basis.

One of the problems with the strategy of concentrating on a single
chain from the start is that you might pick the bad starting point.

Another approach we've been discussing is running each chain
for a fixed time rather than a fixed number of iterations. Then
we get an implicit downweighting of the slow-to-converge chains if
we weight by number of iterations. We're still working out the details.

- Bob

[Michael Betancourt]

I am extremely dubious of this approach -- I'm pretty sure you're just producing biased samples. Here's what I think is going on:

The main problem I'm seeing in adaptation is when the covariance is strongly position-dependent. If NUTS can't cross the entirety of the typical set in one trajectory (as is common in hierarchical models without the Matt trick) then the empirical covariance used in the adaptation won't be the global covariance and the post-update chain will slow to a crawl once it's left the region of the previous adaptation update.

What you've done here is to essentially make the problem even worse. By decreasing the maximum tree depth and lowering the target acceptance rate you're inducing random walk behavior that drastically slows the exploration of the sampler. This keeps the chain in a small neighborhood where the covariance is approximately constant and the adaptation appears to be working fine.

[Ryan]

What you've done here is to essentially make the problem even worse.  By decreasing the maximum tree depth and lowering the target acceptance rate you're inducing random walk behavior that drastically slows the exploration of the sampler.  This keeps the chain in a small neighborhood where the covariance is approximately constant and the adaptation appears to be working fine. 

This is a good point. I tried this approach in another problem and the maximum tree depth was too small, creating small effective sample sizes. I've also tried only limiting the maximum tree depth in the initial run while leaving it higher (or even the default) in the full run. This seems to work well, but again has issues with poor convergence diagnostics because of the same initial starting values.

Up to this point it seems like lowering max tree depth does the most to help get my chains to finish at around the same time, even when using random starts to help maintain convergence diagnostics. Due to the issue you mention, how low is too low for max tree depth? What makes 2^10 ideal? Shouldn't I be able to get away with 2^1 and simply have to run for a really long time (since it is closest to a random walk)? Won't effective sample sizes illustrate this to me, or is it possible for me to be fooled without using a very high max tree depth?

[Michael Betancourt]

Yeah, you could set the tree depth to zero but then you're not taking advantage of HMC. The problem here is not NUTS but rather the adaptation -- in Stan 2.0 you can run the samplers with a global step size adaptation, but until then you can try setting the step size by hand to avoid the adaptation. Incidently, I believe 2^10 was set as the default to prevent NUTS from running for huge lengths of time without being the user being made aware of potential problems.

The biggest problem is that poor adaptation can actually bias the sampler by preventing the chain from visiting regions of high curvature (or extremely different curvature from the current point), and that won't show up in ESS/Rhat calculations.

[Bob Carpenter]

On 7/21/13 2:33 PM, Ryan wrote:
>
> What you've done here is to essentially make the problem even worse. By decreasing the maximum tree depth and
> lowering the target acceptance rate you're inducing random walk behavior that drastically slows the exploration of
> the sampler. This keeps the chain in a small neighborhood where the covariance is approximately constant and the
> adaptation appears to be working fine.
>
>
> This is a good point. I tried this approach in another problem and the maximum tree depth was too small, creating small
> effective sample sizes. I've also tried only limiting the maximum tree depth in the initial run while leaving it higher
> (or even the default) in the full run. This seems to work well, but again has issues with poor convergence diagnostics
> because of the same initial starting values.
>
> Up to this point it seems like lowering max tree depth does the most to help get my chains to finish at around the same
> time, even when using random starts to help maintain convergence diagnostics. Due to the issue you mention, how low is
> too low for max tree depth? What makes 2^10 ideal?

2^10 was just a pragmatic compromise between 2^0 and 2^infinity.

> Shouldn't I be able to get away with 2^1 and simply have to run for a
> really long time (since it is closest to a random walk)?

Yes, with enough emphasis on "really".

> Won't effective sample sizes illustrate this to me, or is it
> possible for me to be fooled without using a very high max tree depth?

No. What'll happen is you can get stuck in a given region and
explore around it and you'll get something that has a high
estimated effective sample size and a low actual effective sample
size. You can see this, for example, in the Metropolis plot
in the NUTS paper --- it doesn't explore very far, but I'm guessing
the diagnostics would look good.

- Bob

[Andrew Gelman]

On Jul 21, 2013, at 5:39 PM, Bob Carpenter wrote:

>
> 2^10 was just a pragmatic compromise between 2^0 and 2^infinity.

10 is halfway between 1 and infinity on the lllllg scale.

>
> No. What'll happen is you can get stuck in a given region and
> explore around it and you'll get something that has a high
> estimated effective sample size and a low actual effective sample
> size. You can see this, for example, in the Metropolis plot
> in the NUTS paper --- it doesn't explore very far, but I'm guessing
> the diagnostics would look good.
>

We are working on diagnostics that would flag this sort of problem, but these are still under development.