run Stan by time -- solve stuck and variable chain time problems

[Bob Carpenter]

I've been eating our dogfood a lot recently and this inevitably
brought me back to Andrew's perspective on timing chains. And it
makes me think we should run many more shorter chains given how
quickly Stan converges.

So how about this:

1. for (chain in 1:NUM_CHAINS) asynchronously [ideally in parallel]

a. run adaptation for chain for a fixed time t_a,

b. run sampling for chain for a fixed time t_s, and

2. merge the results with our jagged R-hat and n_eff.

Is there any reason this wouldn't preserve detailed balance?

Which version of jagged R-hat/n_eff to use, the one that weights
by chain or the one that weights by number of samples? The latter
would seem to be more robust in the face of stuck chains.

Running lots of chains means that each stuck chain contributes
less to the overall total. If we run 32 chains of length 100
instead of 4 chains of length 800, we get the same number of
samples. But with 32 chains, the stuck ones should wash out
in the R-hat and n_eff calcs, assuming they're stuck somewhere
near high posterior mass.

- Bob

[Andrew Gelman]

I'm not sure about this. I've usually thought that 4 or so chains should suffice. My reasoning was that we can't really verify mixing until the chains have done about 50 tours each, at which point we have something like 20 independent draws, which is good enough for usual purposes. If I wanted more draws, I'd run longer, on the theory that then the chains would explore the space further.
On the other hand, here you're talking about problems of stuck chains, and that's another story. It does seem like we need a general way to handle this problem.
A

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

[Bob Carpenter]

It's for robustness.

Is there an argument for 4? It seems like if you're
hunting for modes on init, you'd want a lot. Or if
you're looking for tails with low transition probabilities,
you'd want to run a long time.

With only 4, you do need a lot of back and forth movement
to conclude you've converged. My hope was that we'd need
less with more chains.

I think we get good convergence in most problems after 10
or 20 iterations --- some have taken up to 100 or 200. After
that, the fast to converge ones mix really well, often
near an effective sample per draw. Others not so much, with
more like 100 draws per effective sample. If we then want
100 effective samples, that's 10K iterations, and we can do them
with 30 chains of 3 samples each or fewer longer chains.

- Bob

[Andrew Gelman]

I could see that more chains could be better, we just then need some plan to do with the slow ones.

[Bob Carpenter]

The point of having more chains is that when we have
slow ones, they're only one chain of 32, so they
don't tank the R-hat and n_eff estimates. At least
that's my hope.

- Bob

[Andrew Gelman]

Yup, this is something that arises particularly with Nuts.

[Michael Betancourt]

Firstly -- let's forget about finding modes. HMC is never going to be good at transitioning between nodes even if we can find them.

Consequently the main advantage of multiple starting points is to ensure that either our adaptation hasn't gone awry or we're not missing parts of the distribution with significantly different curvature (for example, the triangle or the Cauchy). I would argue that most of the problems we're seeing are due to bad adaptation and the multiple starting point is doing exactly what it should. In particular, I really like how the combined approach of putting more thought into the initializations and choosing better priors has proved so useful. Better models for a better tomorrow!

On the other hand, if we really do have a model with very different curvature (again, think triangle or Cauchy) then running with many small chains and ignoring those with small ESS can hide that fact and the false belief that everything is automatic could lead to bias.

Personally I think the best way forward is to improve the adaption. The first easy fix is to add a period of adaptation-free warmup (which I've been playing with on and off since coming to London but just haven't had time to run significant tests). Then, once the higher-order autodiff is ready, we can starting using a running average of the SoftAbs metric instead of estimating the global covariance.

[Bob Carpenter]

On 9/9/13 4:37 PM, Michael Betancourt wrote:
> Firstly -- let's forget about finding modes. HMC is never going to be good at transitioning between nodes even if we can find them.

We're not trying to find all the modes, it's just one of
the motivations for using R-hat. If different chains are
stuck in different modes, there's not mixing.

> Consequently the main advantage of multiple starting points is to ensure that either our adaptation hasn't gone awry or we're not missing parts of the distribution with significantly different curvature (for example, the triangle or the Cauchy). I would argue that most of the problems we're seeing are due to bad adaptation and the multiple starting point is doing exactly what it should. In particular, I really like how the combined approach of putting more thought into the initializations and choosing better priors has proved so useful. Better models for a better tomorrow!
>
> On the other hand, if we really do have a model with very different curvature (again, think triangle or Cauchy) then running with many small chains and ignoring those with small ESS can hide that fact and the false belief that everything is automatic could lead to bias.

This is what I'm worried about. We get the same problem if we don't
choose an init point and get stuck.

> Personally I think the best way forward is to improve the adaption. The first easy fix is to add a period of adaptation-free warmup (which I've been playing with on and off since coming to London but just haven't had time to run significant tests). Then, once the higher-order autodiff is ready, we can starting using a running average of the SoftAbs metric instead of estimating the global covariance.

Until then, what do we do? What I think users will do is keep
running until none of the chains are stuck. It's what I'd do.

If that introduces bias, we need to warn users not to do it.

And then to do what, instead?

- Bob

[Michael Betancourt]

>> Firstly -- let's forget about finding modes. HMC is never going to be good at transitioning between nodes even if we can find them.
>
> We're not trying to find all the modes, it's just one of
> the motivations for using R-hat. If different chains are
> stuck in different modes, there's not mixing.

Sure, but making any suggestion that we can find all modes, or even some modes, by trying lots of different initialization values is a bit of a false flag. Once you're trying to deal with mode finding you really want a global, compressive algorithm like nested sampling (perhaps tempering or annealing if you have a particularly well-behaved model) and not a local algorithm like MCMC.

>> Consequently the main advantage of multiple starting points is to ensure that either our adaptation hasn't gone awry or we're not missing parts of the distribution with significantly different curvature (for example, the triangle or the Cauchy). I would argue that most of the problems we're seeing are due to bad adaptation and the multiple starting point is doing exactly what it should. In particular, I really like how the combined approach of putting more thought into the initializations and choosing better priors has proved so useful. Better models for a better tomorrow!
>>
>> On the other hand, if we really do have a model with very different curvature (again, think triangle or Cauchy) then running with many small chains and ignoring those with small ESS can hide that fact and the false belief that everything is automatic could lead to bias.
>
> This is what I'm worried about. We get the same problem if we don't
> choose an init point and get stuck.

Right, but it's easier to diagnose if we run a few chains for a while. My concern with lots of little chains is that there'd be so much Monte Carlo noise that you won't be able to diagnose anything. Give that we're seeing lots of chains getting stuck in initial models with only four chains, it seems like a reasonable choice.

>> Personally I think the best way forward is to improve the adaption. The first easy fix is to add a period of adaptation-free warmup (which I've been playing with on and off since coming to London but just haven't had time to run significant tests). Then, once the higher-order autodiff is ready, we can starting using a running average of the SoftAbs metric instead of estimating the global covariance.
>
> Until then, what do we do? What I think users will do is keep
> running until none of the chains are stuck. It's what I'd do.
>
> If that introduces bias, we need to warn users not to do it.
>
> And then to do what, instead?

I think we've had lots of success in the past few days with the usual suggestions -- think of your model more deeply and consider either better initialization (yes, this can potentially lead to bias so it's not perfect) and better priors that clean up the tail behavior. To channel Andrew, "Folk theorem folk theorem folk theorem." ;-)

[Andrew Gelman]

Just quickly: I agree with Mike that we should forget about trying to solve muliple modes. If we get lucky and can handle a problem with multiple modes, fine, but that can't be our plan.
But we have lots of non-multimodal problems with slow mixing, including various stan-users queries recently. Even the 8 schools has problems if it's not Matt-tricked. And we do want to be able to diagnose those. And, in the short and medium term, we need some plan/recommendation of what to do when we have several chains and one or more are slow moving. Maybe some way to hold the hand of HMC/Nuts to see if it can be pulled toward an area of faster motion?
A

[Michael Betancourt]

Most of these issues seem to be problems with adaptation being too aggressive in the tails. Sharing adaptation information between chains (as we have discussed before) would help, but also require a significant amount of refactoring our code (we'd basically have to implement ensemble samplers). Using SoftAbs is another option but, as previous discussed, will have to wait for the higher-order autodiff.

The easiest stopgap would be tweaked warmup -- I'll try to get to it by the end of the week.

[Bob Carpenter]

On 9/9/13 5:40 PM, Michael Betancourt wrote:
> Most of these issues seem to be problems with adaptation being too aggressive in the tails.
> Sharing adaptation information between chains (as we have discussed before) would help,
> but also require a significant amount of refactoring our code (we'd basically have to implement
> ensemble samplers). Using SoftAbs is another option but, as previous discussed, will have
> to wait for the higher-order autodiff.
>
> The easiest stopgap would be tweaked warmup -- I'll try to get to it by the end of the week.

I think it's more important right now to get the
command-line stuff working for 2.0. I'm getting really
antsy not having a release candidate that we trust.

I would rather put off making changes to adaptation or
other parts of the algorithms until we have a stable
2.0 version passing all tests. For a start, all of these
issues are w.r.t. 1.3. One advantage of 2.0, if I understand
right, is that you use all the warmup samples for adaptation --- weren't
some getting thrown away in 1.3.

So maybe if the currently running model tests pass we can
think about changing adaptation before 2.0, but I have
a strong preference for waiting until after 2.0.

- Bob

[Michael Betancourt]

> So maybe if the currently running model tests pass we can
> think about changing adaptation before 2.0, but I have
> a strong preference for waiting until after 2.0.

Oh, totally agreed. I have a bunch of back-logged features that I'll be adding but I don't expect any of them to come before the 2.0 release.