Inconsistent timing for HMC runs

[Cole Monnahan]

I noticed when using Stan to run HMC that there was a lot of variation in the sampling time of my models. I ran a quick test of this, taking the 8-schools model and running it across a int_time of 1 to 20 and a 20 different seeds. As it turns out, the timings for the warmup period are very consistent and scale with int_time as expected. However, the sampling time was indeed very variable. Further, I ran NUTS with the same seeds and took the mean of the n_leapfrog__ for the warmup and sampling periods and that doesn't match the same pattern as HMC. Below is some quick simulation code.

From what I understand of these algorithms, particularly HMC, the core computational costs are the leapfrog algorithm, and even if the algorithm rejected every proposal it should take the same time as one that accepted every proposal.

This raises several questions:

1. Why would the sampling period for HMC be variable when the warmup was not (both were the same number of iterations)?
2. Shouldn't HMC be very consistent and scale directly with int_time?
3. Shouldn't the NUTS algorithm take the same time as HMC given the same number of total leapfrog steps?
4. Shouldn't the timing be independent of the tuning parameters, given total leapfrog steps? 
5. Am I missing something obvious here?

System Info:

R version 3.2.2 (2015-08-14)

Platform: x86_64-w64-mingw32/x64 (64-bit)

Running under: Windows 8 x64 (build 9200)

Rstan=2.8.2

Rtools gcc= 4.6.3

Here's code to reproduce this test, using the 8schools.stan model.

## Speed tests for Stan algorithms. HMC should scale linearly with int_time

## and NUTS should be similar for average n_leapfrog__ for the warmup and

## sampling period.

seeds <- 1:20

int_time <- seq(1,20, by=2)

Nsim <- 5000

schools_dat <- list(J=8, y=c(28,8,-3,7,-1,1,18,12), sigma=c(15,10,16,11,9,11,10,18))

fit <- stan(file='schools/schools.stan', data=schools_dat, iter=10,

                chains=1, seed=seed, algorithm='HMC', control=list(int_time=1))

speed.hmc <-

  ldply(seeds, function(seed){

          fit.hmc <- ldply(int_time, function(L){

                    fit <- stan(fit=fit, data=schools_dat, iter=Nsim,

                                chains=1, seed=seed, algorithm='HMC', control=list(int_time=L))

                    data.frame(L=L, seed=seed, rstan::get_elapsed_time(fit))

                    })})

speed.nuts <-

  ldply(seeds, function(seed){

          fit <- stan(fit=fit, data=schools_dat, iter=Nsim,

                      chains=1, seed=seed, algorithm='NUTS')

          temp <- data.frame(get_sampler_params(fit))

          data.frame(warmup.nleap=mean(temp[(1:Nsim/2), 'n_leapfrog__']),

                          sample.nleap=mean(temp[-(1:Nsim/2),'n_leapfrog__']),

                          get_elapsed_time(fit), seed=seed)})

par(mfrow=c(1,2))

plot(warmup~L, data=speed.hmc, ylab="Time", xlab="Steps", main='Warmup',

     ylim=c(0, max(speed.hmc$sample)))

points(warmup~warmup.nleap, data=speed.nuts, col='red', pch=2)

plot(sample~L, data=speed.hmc, ylab="Time", xlab="Steps", main='Sample',

     ylim=c(0, max(speed.hmc$sample)))

points(sample~sample.nleap, data=speed.nuts, col='red', pch=2)

legend('topleft', legend=c("HMC", "NUTS"), col=c('black', 'red'), pch=c(1,2))

[Michael Betancourt]

Firstly, when including figures please make sure that they are large enough to at least

read the axis labels easily enough.

Secondly, let me be a bit pedantic regarding names.  Hamiltonian Monte Carlo (HMC)

is a very generic way of generating high-performance Markov chains.  NUTS is one

_implementation_ of HMC.  As is _static HMC_, which uses a constant integration 

time and takes the final point as a proposal.  In your comparisons you are comparing

static HMC with NUTS.

Thirdly, integration time is not equal to the number of leapfrog steps.  Rather integration

time is the product of the integrator step size and total number of leapfrog steps, so

n_leapfrog = int_time / epsilon.  Also, the mean n_leapfrog in NUTS isn’t particularly

helpful due to the typical skewness of the distribution.

Finally, your expectations are almost correct _conditioned on all other tuning parameters 

being the same_.  In particular, the integrator step size and mass matrix can make huge

differences in the computational cost of a transition.  This is important because the

adaptation we do in warmup learns these other tuning parameters, and varying adaptation

leads to varying performance.  

Long story short, to do a proper comparison you’ll want to turn adaptation off to ensure

the same tuning parameters and make sure to compare the number of leapfrog steps

and not the integration time.

1. Why would the sampling period for HMC be variable when the warmup was not (both were the same number of iterations)?

Because of adaptation in warmup.  See the manual for a discussion of our windowed

adaptation procedure.  In particular, with iter=10 you have only 5 warmup iterations

which will update the mass matrix only once, after the fifth iteration.  If you use a

much more reasonable warmup then you’ll start to see variation in warmup timings,

too.

2. Shouldn't HMC be very consistent and scale directly with int_time?

Practical implementations of HMC use numerical integrators, which imply that the

algorithms should scales with the number of leapfrog steps, not integration time

directly.

3. Shouldn't the NUTS algorithm take the same time as HMC given the same number of total leapfrog steps?

Transition by transition, yes, if you ignore algorithmic overhead.  But that doesn’t

mean that they will when you average over chains set to different tunings.

4. Shouldn't the timing be independent of the tuning parameters, given total leapfrog steps? 

If you could set the total leapfrog steps then yes.  But you can’t do that in any of

our samplers right now.

[Cole Monnahan]

Hi Michael, thanks for this information. 

I'm interested in these details b/c in my field there are other software packages that use HMC/NUTS and I want to do some comparisons with Stan. I read the NUTS paper carefully and coded it and static HMC from scratch. For some reason I assumed that int_time in Stan must be the number of steps and not integration time, so that clarifies many things.  I'll respond inline below and hope you have time to answer some more of my questions.

On Thursday, December 24, 2015 at 10:25:31 AM UTC-8, Michael Betancourt wrote:

Firstly, when including figures please make sure that they are large enough to at least

read the axis labels easily enough.

I was just using the defaults for the web interface and figured a thumbnail would be fine since you could click on it and get the full version. I'll attach it full size next time.

 

Secondly, let me be a bit pedantic regarding names.  Hamiltonian Monte Carlo (HMC)

is a very generic way of generating high-performance Markov chains.  NUTS is one

_implementation_ of HMC.  As is _static HMC_, which uses a constant integration 

time and takes the final point as a proposal.  In your comparisons you are comparing

static HMC with NUTS.

Fair enough, but suspect most users won't make that distinction b/c that's how you specify them in Stan.

 

Thirdly, integration time is not equal to the number of leapfrog steps.  Rather integration

time is the product of the integrator step size and total number of leapfrog steps, so

n_leapfrog = int_time / epsilon.  Also, the mean n_leapfrog in NUTS isn’t particularly

helpful due to the typical skewness of the distribution.

It seems most other software programs with HMC (e.g., ADMB), as well as the papers describing it, you specify step size and number of steps. I think it might be worth explaining this in the adaptation section of the manual. The term int_time isn't used at all there, not really described in the help files in Rstan. I thought using mean n_leapfrog for NUTS would be reasonable b/c it is the total steps divided by the number of iterations. I reran with the n_leapfrog calculated as above and plotted that instead of int_time and it totally makes sense now (attached here).

Finally, your expectations are almost correct _conditioned on all other tuning parameters 

being the same_.  In particular, the integrator step size and mass matrix can make huge

differences in the computational cost of a transition.  This is important because the

adaptation we do in warmup learns these other tuning parameters, and varying adaptation

leads to varying performance.  

Neal (2011) describes two ways of using the mass matrix:  a Chloesky decomposition of parameter space or using it in the kinetic energy function. The first seems to be more population. Does Stan use the mass matrix to do a Cholesky decomposition? Is that why there are such computational differences depending on tuning parameters? Are the differences only in unit vs. diagonal vs. dense matrices, or can there also be differences between different dense matrices (e.g. taking advantage of sparseness)?

 

Long story short, to do a proper comparison you’ll want to turn adaptation off to ensure

the same tuning parameters and make sure to compare the number of leapfrog steps

and not the integration time.

1. Why would the sampling period for HMC be variable when the warmup was not (both were the same number of iterations)?

Because of adaptation in warmup.  See the manual for a discussion of our windowed

adaptation procedure.  In particular, with iter=10 you have only 5 warmup iterations

which will update the mass matrix only once, after the fifth iteration.  If you use a

much more reasonable warmup then you’ll start to see variation in warmup timings,

too.

 

I only used iter=10 to build the model once and then used 2500 warmup and 2500 sampling iterations for the actual analyses. If the variance in computational time is due to the adaptive mass matrix, then the last period of the warmup should be similar to the sampling period. I'll check that.

2. Shouldn't HMC be very consistent and scale directly with int_time?

Practical implementations of HMC use numerical integrators, which imply that the

algorithms should scales with the number of leapfrog steps, not integration time

directly.

3. Shouldn't the NUTS algorithm take the same time as HMC given the same number of total leapfrog steps?

Transition by transition, yes, if you ignore algorithmic overhead.  But that doesn’t

mean that they will when you average over chains set to different tunings.

4. Shouldn't the timing be independent of the tuning parameters, given total leapfrog steps? 

If you could set the total leapfrog steps then yes.  But you can’t do that in any of

our samplers right now.

I think it would be valuable for the user to be able to specify the number of steps opposed to the integration time. For instance Langevin sampling can't be done directly in Stan (if I understand correctly). I could see that being potentially useful for some people in some models. What was the reasoning behind specifying integration time instead? Also, the calculation from n_step=int_time/stepsize isn't integer.. is this value rounded up?

I'm new to Stan but really like it so far and appreciate all the development that has gone into it. I think it may be helpful to the user base to have a more thorough explanation of the adaptation procedure that is going on internally. This will be particularly true for users who come from other software platforms or have coded it themselves on specific problems. For example I couldn't find the default behavior of argument metric, or how to pass an arbitrary mass matrix.

Is it correct that this is what happens when using static HMC: (1) User gives int_time, (2) a 'reasonable' stepsize is found and from this an approximate number of leapfrog steps (L), (3) the warmup procedure begins and stepsize is updated (and hence L) using dual averaging. (4) Continue with sampling.

 Thanks,

Cole

[Michael Betancourt]

Secondly, let me be a bit pedantic regarding names.  Hamiltonian Monte Carlo (HMC)

is a very generic way of generating high-performance Markov chains.  NUTS is one

_implementation_ of HMC.  As is _static HMC_, which uses a constant integration 

time and takes the final point as a proposal.  In your comparisons you are comparing

static HMC with NUTS.

Fair enough, but suspect most users won't make that distinction b/c that's how you specify them in Stan.

It is in CmdStan, which distinguishes between algorithm (HMC) and 

engine (static or NUTS).  The conventions used in RStan (and then

PyStan) were not proposed by the entire team.  This will hopefully

be remedied in the refactor.

It seems most other software programs with HMC (e.g., ADMB), as well as the papers describing it, you specify step size and number of steps. I think it might be worth explaining this in the adaptation section of the manual. The term int_time isn't used at all there, not really described in the help files in Rstan.

And just about every implementation of HMC outside of Stan is

just a reimplementation of code in the NUTS paper or Stan itself.

Integration time is more natural, especially when trying to optimize

the various free parameters in HMC.  Trying to optimize step size

and number of leapfrog steps introduces all kinds of problems, as

evident in many, many papers.

Again, the inconsistencies between the interfaces is a known problem

and we are working to improve that.

Finally, your expectations are almost correct _conditioned on all other tuning parameters 

being the same_.  In particular, the integrator step size and mass matrix can make huge

differences in the computational cost of a transition.  This is important because the

adaptation we do in warmup learns these other tuning parameters, and varying adaptation

leads to varying performance.  

Neal (2011) describes two ways of using the mass matrix:  a Chloesky decomposition of parameter space or using it in the kinetic energy function. The first seems to be more population. Does Stan use the mass matrix to do a Cholesky decomposition? Is that why there are such computational differences depending on tuning parameters? Are the differences only in unit vs. diagonal vs. dense matrices, or can there also be differences between different dense matrices (e.g. taking advantage of sparseness)?

Those are two reparameterizations of the same problem, one is

really not different from the other in practice.  The computational

consequence is on how well-chosen the mass matrix is with

respect to the target distribution.

Stan estimates an optimal mass matrix (really Euclidean metric)

in windows during warmup, updating the mass matrix to the best

current estimate at the end of each window.

I think it would be valuable for the user to be able to specify the number of steps opposed to the integration time.

No, it would not.  The optimal value of the number of leapfrog

steps depends on the step size — in other words only the

integration time and step size can be optimized independently

of each other.  People keep making this mistake in the literature

and it causes no end of confusion that we (often I) have to clean

up.

For instance Langevin sampling can't be done directly in Stan (if I understand correctly). I could see that being potentially useful for some people in some models.

There is no reason for ever using Langevin over HMC — if there

were I can assure you that it would be in Stan.  We’re not just

implementing algorithms here, we’re also the ones actively

researching them.

What was the reasoning behind specifying integration time instead?

See above.

Also, the calculation from n_step=int_time/stepsize isn't integer.. is this value rounded up?

It rounds per the implementation of C++,

https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/static/base_static_hmc.hpp

I'm new to Stan but really like it so far and appreciate all the development that has gone into it. I think it may be helpful to the user base to have a more thorough explanation of the adaptation procedure that is going on internally. This will be particularly true for users who come from other software platforms or have coded it themselves on specific problems. For example I couldn't find the default behavior of argument metric, or how to pass an arbitrary mass matrix.

See Section 57.2 of the manual.  This is all extensively documented, and any

confusion can be resolved by looking at the code directly.

Is it correct that this is what happens when using static HMC: (1) User gives int_time, (2) a 'reasonable' stepsize is found and from this an approximate number of leapfrog steps (L), (3) the warmup procedure begins and stepsize is updated (and hence L) using dual averaging. (4) Continue with sampling.

The initial step size in (2) does not depend on the choice of integration time.

It uses the error from a single step.

[Bob Carpenter]

Thanks for taking the time to give us such detailed comments.

Did we manage to give you enough info to answer your original
list of questions about variability of run times?

If you look at the evaluation in the original NUTS paper,
you'll see that it's given in terms of epsilon * L, which
is integration time (less rounding); for each integration
time, an ideal step size was estimated based on rejection rate.
It took a lot of cluster computing to generate those figures.

As Michael indicated, we're working toward unifying the interfaces,
with the big change targeted for Stan 3. That should also allow
more control over specifying mass matrices; but we think that's
only going to be useful for stopping and restarting.

We're also continuing to work on the manual --- I think we should
be able to put in much more algorithm description before referring
people to the code (though that's the best we have now and will
always be the final arbiter of what Stan actually does in practice).

Michael --- isn't Langevin just 1 step of HMC? And if so, for a fixed
step size (discretization interval), couldn't we just fix integration time
to be one step and get Langevin?

- Bob

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <stan_speed_tests2.png>

[Michael Betancourt]

Michael --- isn't Langevin just 1 step of HMC?  

Yup.  Well, the “Metropolis-Adjusted Langevin Algorithm”

is if you really want to be pedantic.

And if so, for a fixed step size (discretization interval),

couldn't we just fix integration time to be one step and get Langevin?

If the step size is fixed then you could do int_time = epsilon.

Or, noting that the minimum number of leapfrog steps is 1,

you could also set int_time to something really small and

for any reasonable adapted step size you’d still get Langevin.

But just because you _could_ do it doesn’t mean that you should….

[Bob Carpenter]

Agreed. I pureposefully didn't doc it in the manual just because
it has a name.

The main reason we don't want to put in every algorithm we
can think of is that then there's a tendency by users to try them
all on problems. We've actually built several other named samplers
that we've not tried to include, specifically to avoid wasting our
own and our users' time.

- Bob

[Michael Betancourt]

YES.

[Krzysztof Sakrejda]

It took me a while to come around to this idea that Stan is trying to do best-in-class algorithms and while it might be nice to have pure MH for completeness it would mean a lot of maintenance and it's not worth it given that we have HMC.  A lot of people just aren't going to have the background to understand that, for example, slice sampling and MH don't make sense since we have HMC but ADVI and optimization have different enough properties that they're worth including.  

I just looked at the Stan manual and it doesn't really address the choice of algorithms in any clear section (maybe it's in the intro in story form) so maybe right before section 1.5 which launches into talking about sampling there could be a little section about choice of algorithms.  

Krzysztof

[Cole Monnahan]

Hi Michael, I'm probably being dense here but I'm still missing the point why specifying integration time is better than L (and let me know if this is the wrong forum to discuss these matters). For a given int_time there are different combinations of eps and L that get you there, and obviously we prefer smaller L computationally. The approach used by Stan (as in Alg 5 of the NUTS paper) will end with a different L after adaptation for different starting seeds/initial conditions, and some of those L's will be larger than others. Why wouldn't it be preferable to force L to be smaller (fix it) and let eps be bigger to get the same int_time? I'm not trying to say you should adapt both L and eps, but rather fix L. 

Also, could you point me toward some studies comparing efficiency of Langevin, HMC and NUTS? Outside of the NUTS paper I've only been able to find theory and derivation of the algorithms.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/FERVo5-h1L4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Cole Monnahan]

Hi Bob, yes I've got my questions answered for the most part. It's been a useful discussion in that I've learned a lot. I think it's worth noting that ADMB coupled AD and static HMC (using inverse Hessian as mass matrix) probably a decade ago and there's a community of scientists that are used to its way of doing things, who likely will migrate over to Stan. I think it's worthwhile for those users (me included, obviously) to be able to notice/understand the differences and improvements to ease that transition. 

I look forward to Stan 3 to see the changes. In the meantime I've got my original questions answered.

Thanks

Cole

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/FERVo5-h1L4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Cole Monnahan]

Hi Krzysztof,

That's a helpful way of putting it, and something I also didn't realize. From a naive user's point of view I do think that added section would be helpful. I read Neal (2011) and then Hoffman and Gelman (2014) thoroughly enough (I'm not a statistician) to code them up, but was still confused by the parameterization of static HMC in Stan. I think at least the argument int_time should be defined in the manual and what exactly is being adapted for static HMC.

That being said, I know documentation doesn't write itself and really appreciate the incredible amount of development that has gone into Stan.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/FERVo5-h1L4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

The place to put changes for the manual is either in the next
manual issue or the longer-term manual issue. This one's not
a bug, so I added it as a longer-term manual issue:

https://github.com/stan-dev/stan/issues/1238#issuecomment-167645355

The manual should be updated in a branch for the next manual issue.
I don't think that exists yet. If you want, just edit that, and if
you fix one of the comments, check its checkbox.

- Bob

[Bob Carpenter]

> On Dec 28, 2015, at 1:25 PM, Cole Monnahan <mon...@uw.edu> wrote:
>
> Hi Bob, yes I've got my questions answered for the most part. It's been a useful discussion in that I've learned a lot. I think it's worth noting that ADMB coupled AD and static HMC (using inverse Hessian as mass matrix) probably a decade ago

I didn't see any mention of HMC in that paper.
Was that the wrong link?

That paper talks about random-walk Metropolis, which is
a terrible algorithm. Also, I don't think their proposed variant
will work, because (a) you need to regularize for positive definiteness (e.g.,
see Michael's paper on Softabs), and (b) with varying curvature, the
Hessian changes, and it defeats detailed balance. There's no further description
of the algorithm or adjustments made for soundness, only a pointer to Gelman
et al.'s book in section 2.3.2.

> and there's a community of scientists that are used to its way of doing things, who likely will migrate over to Stan. I think it's worthwhile for those users (me included, obviously) to be able to notice/understand the differences and improvements to ease that transition.

That'd be great from our perspective. We can explain Stan, but
I don't know ADMB. If you want to write a Stan tutorial for ADMB
users, we're more than happy to help out and advertise on Andrew's
blog.

Andrew's pushing us to use Laplace approximations, both for
penalized/marginal maximum likelihood and for black-box
variational inference, but there's no reason to do that (and
lots of reasons not to) in cases where MCMC will work.

- Bob

[Cole Monnahan]

Hi Bob, a few thoughts inline below...

On Mon, Dec 28, 2015 at 12:34 PM, Bob Carpenter <ca...@alias-i.com> wrote:

> On Dec 28, 2015, at 1:25 PM, Cole Monnahan <mon...@uw.edu> wrote:
>
> Hi Bob, yes I've got my questions answered for the most part. It's been a useful discussion in that I've learned a lot. I think it's worth noting that ADMB coupled AD and static HMC (using inverse Hessian as mass matrix) probably a decade ago

I didn't see any mention of HMC in that paper.
Was that the wrong link?

That paper talks about random-walk Metropolis, which is
a terrible algorithm.  Also, I don't think their proposed variant
will work, because (a) you need to regularize for positive definiteness (e.g.,
see Michael's paper on Softabs), and (b) with varying curvature, the
Hessian changes, and it defeats detailed balance.  There's no further description
of the algorithm or adjustments made for soundness, only a pointer to Gelman
et al.'s book in section 2.3.2.

I was a bit too brief. ADMB was essentially written by a single person and then was recently bought and turned open source. It was/is designed for maximum likelihood estimation with MCMC capabilities tacked on (since why not?). It's amazingly robust and fast at optimization. If the Hessian at the mode is positive definite it will use the inverse for the RWM multivariate normal, or the mass matrix for the static HMC. To be clear it uses a fixed mass matrix avoiding issues above. It does indeed have static HMC it's just poorly documented, but you can see the details in the code here. I became interested in HMC because RWM was so slow for our big, tough models. Part of my dissertation research (hopefully) is to bring these new methods to our field.

 

> and there's a community of scientists that are used to its way of doing things, who likely will migrate over to Stan. I think it's worthwhile for those users (me included, obviously) to be able to notice/understand the differences and improvements to ease that transition.

That'd be great from our perspective.  We can explain Stan, but
I don't know ADMB.  If you want to write a Stan tutorial for ADMB
users, we're more than happy to help out and advertise on Andrew's
blog.

Yes I could do that, should I fork stan-dev and then do a pull request? ADMB and Stan are actually very similar so a few examples of ADMB models in Stan would go a long way. As part of my research I'll be coding lots of models in both so it would be easy to throw together a quick guide with examples.  

Andrew's pushing us to use Laplace approximations, both for
penalized/marginal maximum likelihood and for black-box
variational inference, but there's no reason to do that (and
lots of reasons not to) in cases where MCMC will work.

There's also TMB which is similar to ADMB but has Laplace approximations built in. I'm also trying to recreate basic HMC algorithms for it too. 

-Cole

[Andrew Gelman]

Stan should be just as fast for these optimization problems, no?

A

[Cole Monnahan]

Andrew: Yes, I would assume so. Caveat the following with the fact I'm a student just learning all this myself... What makes ADMB so powerful is AD, a template language, and compiled C++ -- all of which Stan also has. ADMB has some additional nice features since it was designed for optimization. I could see the fishery and ecology communities moving to Stan over time but there's a ton of legacy ADMB code so that might be a slow transition. For instance, many fishery stock assessments done by the government are using a 20,000 line ADMB model (stock synthesis; Methot and Wetzel, 2011) which realistically no one is going to reproduce in Stan anytime soon.

Our field has relied on ADMB since the 1990s because it was the best (only?)  option for robust non-linear optimization at the time, so it is really exciting to see the development of Stan and what it has to offer moving forward.

--

[Bob Carpenter]

> On Dec 31, 2015, at 2:40 AM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
>
> Stan should be just as fast for these optimization problems, no?
> A

I don't see why not. I'm assuming our autodiff and density
library is faster and more expressive, too, but we weren't able
to get ADMB's autodiff working to test it.

>> On Dec 29, 2015, at 12:44 PM, Cole Monnahan <mon...@uw.edu> wrote:
>>
>>
>> I was a bit too brief. ADMB was essentially written by a single person and then was recently bought and turned open source. It was/is designed for maximum likelihood estimation with MCMC capabilities tacked on (since why not?).

The reason not to implement Metropolis for users is that
you don't want to give them something that's going to return
the wrong answer.

This is why we're so worried about variational inference at the
moment and why it has a big "EXPERIMENTAL" tag on it. Even in
cases where it returns an answer, we don't undrestand the properties
of the answer (either the point estimate of the approximating density
or its variance).

There are models for which Metropolis will seem to mix OK,
non-robust autocorrelation estimators like the one in Coda
(a popular R package) will be used to assess effective sample
size and miss long-term corelations, and the user will be led
to believe they've converged and can estimate posterior
expectations when in fact their results will be wrong.

>> It's amazingly robust and fast at optimization.

Optimization's relatively well understood. At least for
convex problems.

>> If the Hessian at the mode is positive definite it will use the inverse for the RWM multivariate normal, or the mass matrix for the static HMC.

Is there a description of the algorithm somewhere?

The problem you run into with that approach is that the
mode often doesn't exist, or is not in the typical set, or the
posterior has varying curvature so that a static estimate
won't work well. We've found a diagonal variance to be more
robust and we have to compensate by taking small step sizes so
that an iterated first-order gradient can deal with the curvature.
But that introduces problems like stiff ODE solvers.

- Bob

[Bob Carpenter]

> On Dec 31, 2015, at 12:22 PM, Cole Monnahan <mon...@uw.edu> wrote:
>
> Andrew: Yes, I would assume so. Caveat the following with the fact I'm a student just learning all this myself... What makes ADMB so powerful is AD, a template language, and compiled C++ -- all of which Stan also has. ADMB has some additional nice features since it was designed for optimization. I could see the fishery and ecology communities moving to Stan over time but there's a ton of legacy ADMB code so that might be a slow transition.

We're already seeing people in that community use Stan.
And in the closely related population ecology fields. Always
amazing to me how intellectual communities organize.

> For instance, many fishery stock assessments done by the government are using a 20,000 line ADMB model (stock synthesis; Methot and Wetzel, 2011) which realistically no one is going to reproduce in Stan anytime soon.

How do you debug a 20K line ADMB model? We still haven't figured
that out for Stan. For something that big, we'd write the functions
in C++ anyway, and we know how to test that.

> Our field has relied on ADMB since the 1990s because it was the best (only?) option for robust non-linear optimization at the time, so it is really exciting to see the development of Stan and what it has to offer moving forward.

It's the combination of autodiff (or symbolic) gradients and
a simple language that's the key. I also don't know of anything
other than ADMB that did that. There were plenty of C++ and
autodiff packages, but nothing as far as I know putting them
together. Having a big community of users, like the fisheries and
wildlife people for ADMB is also a huge benefit for everybody.

- Bob

[Cole Monnahan]

Unfortunately not really. There's a brief description in chapter 2 in the manual: http://ftp.admb-project.org/admb-11.4/manuals/admb-11.4.1.pdf but honestly the author is the only person who really understands it. I've looked at the source code enough to get the gist but there are some pieces in there that I don't understand. I can tell you there's no adaptation of HMC and for RWM there's only a global scalar for the covariance matrix that tunes the acceptance rate.

The problem you run into with that approach is that the
mode often doesn't exist, or is not in the typical set, or the
posterior has varying curvature so that a static estimate
won't work well.  We've found a diagonal variance to be more
robust and we have to compensate by taking small step sizes so
that an iterated first-order gradient can deal with the curvature.
But that introduces problems like stiff ODE solvers.

Indeed we are familiar with these issues (which obviously affects MLE estimates as well). In these cases we simply can't do MCMC and have to change the model or abandon it. The user can specify an arbitrary covariance matrix so they can try that, but that's rare to do. We take what we can get with ADMB. This is why we're so excited by Stan.

- Bob

>
> --
> 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.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/FERVo5-h1L4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Cole Monnahan]

On Thu, Dec 31, 2015 at 9:29 AM, Bob Carpenter <ca...@alias-i.com> wrote:

> On Dec 31, 2015, at 12:22 PM, Cole Monnahan <mon...@uw.edu> wrote:
>
> Andrew: Yes, I would assume so. Caveat the following with the fact I'm a student just learning all this myself... What makes ADMB so powerful is AD, a template language, and compiled C++ -- all of which Stan also has. ADMB has some additional nice features since it was designed for optimization. I could see the fishery and ecology communities moving to Stan over time but there's a ton of legacy ADMB code so that might be a slow transition.

We're already seeing people in that community use Stan.
And in the closely related population ecology fields.  Always
amazing to me how intellectual communities organize.

Interesting indeed. The fisheries crowd will be slower to adapt than wildlife/ecology I suspect.

 

> For instance, many fishery stock assessments done by the government are using a 20,000 line ADMB model (stock synthesis; Methot and Wetzel, 2011) which realistically no one is going to reproduce in Stan anytime soon.

How do you debug a 20K line ADMB model?  We still haven't figured
that out for Stan.  For something that big, we'd write the functions
in C++ anyway, and we know how to test that.

I'm not the author (nor an expert in ADMB or C++) but I would guess it's very difficult. The traditional C++ debugging tools do work for ADMB models.  

> Our field has relied on ADMB since the 1990s because it was the best (only?)  option for robust non-linear optimization at the time, so it is really exciting to see the development of Stan and what it has to offer moving forward.

It's the combination of autodiff (or symbolic) gradients and
a simple language that's the key.  I also don't know of anything
other than ADMB that did that.  There were plenty of C++ and
autodiff packages, but nothing as far as I know putting them
together.  Having a big community of users, like the fisheries and
wildlife people for ADMB is also a huge benefit for everybody.

Yeah it's cool to see communities start to interact more. Good for science.

 

- Bob