When does NUTS fail?
403 views
Skip to first unread message
Sebastian Weber
Mar 20, 2013, 5:35:49 AM3/20/13
to stan-...@googlegroups.com
Hi!
I am wondering under what conditions does NUTS fail to adapt discretization and leapfrog integration steps between sampling points, which are crucial for HMC to work as I understood it? I.e. I never saw a diagnostic for this adaptation process. So do I just get a small effective sample if NUTS performs bad? Is it possible to diagnose problems of NUTS and if so under what circumstances do these likely occur? Is it sufficient to ensure that the support of all my variables spans always from minus to plus infinity?
I have used Stan so far as black box, but the more I use it, the more I would love to shed some light into that black box...
Best,
Sebastian
I am wondering under what conditions does NUTS fail to adapt discretization and leapfrog integration steps between sampling points, which are crucial for HMC to work as I understood it? I.e. I never saw a diagnostic for this adaptation process. So do I just get a small effective sample if NUTS performs bad? Is it possible to diagnose problems of NUTS and if so under what circumstances do these likely occur? Is it sufficient to ensure that the support of all my variables spans always from minus to plus infinity?
I have used Stan so far as black box, but the more I use it, the more I would love to shed some light into that black box...
Best,
Sebastian
Bob Carpenter
Mar 20, 2013, 12:45:31 PM3/20/13
to stan-...@googlegroups.com
We would very much like to understand some of these
issues ourselves.
The technical issues in diagnosing convergence
presented by adaptive HMC/NUTS are the same as those
presented by Metropolis/Gibbs.
If you've already studied Radford Neal's MCMC handbook
chapter on HMC, you're well on your way to understanding
what's going on. There's then the Hoffman and Gelman
paper on NUTS itself (also on arXiv, though it's about
to come out in JMLR).
The best thing to read next is Michael Betancourt's
arXiv papers on RM-HMC, which discusses the funnel example
(also discussed in the manual) in more detail. The main
pathologies for HMC/NUTS have to do with complex posterior
geometries (high curvature and varying posterior scales).
HMC is rotation invariant (can handle high corelations)
but is not scale invariant (doesn't like badly conditioned
posteriors in the matrix sense). Riemann Manifold HMC,
which we're working on now, is usually slower per iteration,
but can handle more complex posteriors.
Michael's also writing a package for diagnostics, which
should help with this, but isn't quite ready for prime
time (or so I understand --- correct me if I'm wrong here,
Michael).
What we're worried about is the case where the convergence
diagnostics like R-hat are consistent with convergence, but
the sampler hasn't actually converged. One way Michael's
been suggesting dealing with this is to run NUTS, look
at the adapted step size and average (or max) tree depth,
then run basic HMC with a step size 1/2 the size and a
number of steps equal to 2^(depth + 1), where depth is
the max tree depth in the NUTS. This should be substantially
slower, and perhaps require more iterations, but should give
the same answer. If it does, the geometry's probably OK (right,
Michael?).
If you declare variables with constraints and the model
has support in the entire constrained space, then Stan will
deal with the transform to full support on R^N (perhaps for
a different N -- for instance, K-simplexes only need K-1 unconstrained
parameters because the last value is determined by 1 minus
the sum of the others). There's a discussion of how this works
in the manual.
- 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.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>
issues ourselves.
The technical issues in diagnosing convergence
presented by adaptive HMC/NUTS are the same as those
presented by Metropolis/Gibbs.
If you've already studied Radford Neal's MCMC handbook
chapter on HMC, you're well on your way to understanding
what's going on. There's then the Hoffman and Gelman
paper on NUTS itself (also on arXiv, though it's about
to come out in JMLR).
The best thing to read next is Michael Betancourt's
arXiv papers on RM-HMC, which discusses the funnel example
(also discussed in the manual) in more detail. The main
pathologies for HMC/NUTS have to do with complex posterior
geometries (high curvature and varying posterior scales).
HMC is rotation invariant (can handle high corelations)
but is not scale invariant (doesn't like badly conditioned
posteriors in the matrix sense). Riemann Manifold HMC,
which we're working on now, is usually slower per iteration,
but can handle more complex posteriors.
Michael's also writing a package for diagnostics, which
should help with this, but isn't quite ready for prime
time (or so I understand --- correct me if I'm wrong here,
Michael).
What we're worried about is the case where the convergence
diagnostics like R-hat are consistent with convergence, but
the sampler hasn't actually converged. One way Michael's
been suggesting dealing with this is to run NUTS, look
at the adapted step size and average (or max) tree depth,
then run basic HMC with a step size 1/2 the size and a
number of steps equal to 2^(depth + 1), where depth is
the max tree depth in the NUTS. This should be substantially
slower, and perhaps require more iterations, but should give
the same answer. If it does, the geometry's probably OK (right,
Michael?).
If you declare variables with constraints and the model
has support in the entire constrained space, then Stan will
deal with the transform to full support on R^N (perhaps for
a different N -- for instance, K-simplexes only need K-1 unconstrained
parameters because the last value is determined by 1 minus
the sum of the others). There's a discussion of how this works
in the manual.
- 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.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>
Michael Betancourt
Mar 20, 2013, 12:59:02 PM3/20/13
to stan-...@googlegroups.com
First, let's distinguish between simulating Hamiltonian dynamics and NUTS, which uses those simulations to build up a binary tree from which a new point is sampled.
When can Hamiltonian dynamics fail? While there are many failure cases, the most common involve the discretization step size being too large. If the step size is too large everywhere in space then the dynamics quickly become chaotic and diverge, leading to very short trajectories. If the step size is too large only for some parts of space, then the trajectory will not explore those regions and you can get biased samples. The former is usually well controlled by adaptation algorithm, but failure would manifest as very small tree depths and high autocorrelations (small tree depths and small autocorrelations are okay). The latter is more subtle -- if you want to be safe then run an auxiliary chain (starting at one of your nominal samples) but with 1/10 or 1/100 or the nominal adaptive step size (be careful that you're not hitting the max tree depth! you might have to raise the threshold for these diagnostic chains). Any failure would manifest as inconsistencies in these two samples -- visual examination is often sufficient, but you can run a KS test if you wanted to automate (implementing these kinds of tests in Stan is on the to-do list). If you're model is hierarchical there can also be issues with limited potential variation -- see, for example, http://arxiv.org/abs/1212.4693).
So what about NUTS? Formally NUTS never really fails, but it can terminate prematurely before reaching the optimal integration time. This happens when the posterior is highly correlated and the spectral gap of the covariance is large (large spectral gap manifests has high frequency transverse oscillations which can trigger NUTS to terminate before the low frequency oscillations have a chance to move). This can be diagnosed by looking at individual trajectories or, as above, seeing small tree depths but high autocorrelations. Unfortunately there's nothing much that can be done without appealing to more advanced Hamiltonians, in particular one with a dense mass matrix (which you can conveniently use in Stan 1.2) or Riemannian HMC (coming soon).
We're working on some visual and automated diagnostics to catch these kinds of problems, but if you notice any weird behavior please don't hesitate to make a note of it on the list.
-Mike
When can Hamiltonian dynamics fail? While there are many failure cases, the most common involve the discretization step size being too large. If the step size is too large everywhere in space then the dynamics quickly become chaotic and diverge, leading to very short trajectories. If the step size is too large only for some parts of space, then the trajectory will not explore those regions and you can get biased samples. The former is usually well controlled by adaptation algorithm, but failure would manifest as very small tree depths and high autocorrelations (small tree depths and small autocorrelations are okay). The latter is more subtle -- if you want to be safe then run an auxiliary chain (starting at one of your nominal samples) but with 1/10 or 1/100 or the nominal adaptive step size (be careful that you're not hitting the max tree depth! you might have to raise the threshold for these diagnostic chains). Any failure would manifest as inconsistencies in these two samples -- visual examination is often sufficient, but you can run a KS test if you wanted to automate (implementing these kinds of tests in Stan is on the to-do list). If you're model is hierarchical there can also be issues with limited potential variation -- see, for example, http://arxiv.org/abs/1212.4693).
So what about NUTS? Formally NUTS never really fails, but it can terminate prematurely before reaching the optimal integration time. This happens when the posterior is highly correlated and the spectral gap of the covariance is large (large spectral gap manifests has high frequency transverse oscillations which can trigger NUTS to terminate before the low frequency oscillations have a chance to move). This can be diagnosed by looking at individual trajectories or, as above, seeing small tree depths but high autocorrelations. Unfortunately there's nothing much that can be done without appealing to more advanced Hamiltonians, in particular one with a dense mass matrix (which you can conveniently use in Stan 1.2) or Riemannian HMC (coming soon).
We're working on some visual and automated diagnostics to catch these kinds of problems, but if you notice any weird behavior please don't hesitate to make a note of it on the list.
-Mike
Michael Betancourt
Mar 20, 2013, 1:04:12 PM3/20/13
to stan-...@googlegroups.com
Coincident responses! Bob's response and mine are entirely consistent.
> Michael's also writing a package for diagnostics, which
> should help with this, but isn't quite ready for prime
> time (or so I understand --- correct me if I'm wrong here,
> Michael).
We're waiting on being able to output gradients which are pretty important for diagnostics. This should be pretty straightforward once I finish the refactoring. In the meantime there's a preliminary version at https://github.com/betanalpha/visual-diagnostics which requires the "graphical_diagnostics" branch of Stan which may or may not be working at the moment.
> What we're worried about is the case where the convergence
> diagnostics like R-hat are consistent with convergence, but
> the sampler hasn't actually converged. One way Michael's
> been suggesting dealing with this is to run NUTS, look
> at the adapted step size and average (or max) tree depth,
> then run basic HMC with a step size 1/2 the size and a
> number of steps equal to 2^(depth + 1), where depth is
> the max tree depth in the NUTS. This should be substantially
> slower, and perhaps require more iterations, but should give
> the same answer. If it does, the geometry's probably OK (right,
> Michael?).
Yup, although I'd recommend dropping a full order of magnitude to be safe. It will take longer but you can always run it while writing up your results.
-Mike
> Michael's also writing a package for diagnostics, which
> should help with this, but isn't quite ready for prime
> time (or so I understand --- correct me if I'm wrong here,
> Michael).
> What we're worried about is the case where the convergence
> diagnostics like R-hat are consistent with convergence, but
> the sampler hasn't actually converged. One way Michael's
> been suggesting dealing with this is to run NUTS, look
> at the adapted step size and average (or max) tree depth,
> then run basic HMC with a step size 1/2 the size and a
> number of steps equal to 2^(depth + 1), where depth is
> the max tree depth in the NUTS. This should be substantially
> slower, and perhaps require more iterations, but should give
> the same answer. If it does, the geometry's probably OK (right,
> Michael?).
-Mike
Sebastian Weber
Mar 23, 2013, 5:44:24 AM3/23/13
to stan-...@googlegroups.com
Great answers! Thanks a lot.
I am really looking forward to seeing this diagnostic tools becoming easier to access. To be honest, since I got Stan I focus much more on my models and only take a look at any convergence if my Rhat is bad, but this usually means my model is crap. In that sense Stan really helped in the spirit of Andrew to focus on modeling, which is great.
As I am trained physicist I am wondering if Stan is automatically transforming stuff into easier geometries. I.e. I remember in my physics courses to always think long and hard about possible variable transformations before solving any problem. In fact, there are approaches out there which reformulate a given mechanics problem (which is what you solve by integrating equations of motion, right?) to write a partial differential equation system which finds a transformation turning a given problem into a set of trivially solvable differentials. The trick is always to identify Noether currents, i.e. constants of the equations of motion like total energy which corresponds to probability mass. However, there are more such constants at least in physical systems, such as conservation of angular momentum and the like. All this is quite apparent if one ever tried to solve Kepler problems in euclidian coordinates, rewriting it in angular coordinates makes life easier...
I am just curios here, so does Stan take these things into account?
Best,
Sebastian
Ah, and to raise my vote here as well again: bringing ODEs into Stan would be really great.
Sebastian Weber
I am really looking forward to seeing this diagnostic tools becoming easier to access. To be honest, since I got Stan I focus much more on my models and only take a look at any convergence if my Rhat is bad, but this usually means my model is crap. In that sense Stan really helped in the spirit of Andrew to focus on modeling, which is great.
As I am trained physicist I am wondering if Stan is automatically transforming stuff into easier geometries. I.e. I remember in my physics courses to always think long and hard about possible variable transformations before solving any problem. In fact, there are approaches out there which reformulate a given mechanics problem (which is what you solve by integrating equations of motion, right?) to write a partial differential equation system which finds a transformation turning a given problem into a set of trivially solvable differentials. The trick is always to identify Noether currents, i.e. constants of the equations of motion like total energy which corresponds to probability mass. However, there are more such constants at least in physical systems, such as conservation of angular momentum and the like. All this is quite apparent if one ever tried to solve Kepler problems in euclidian coordinates, rewriting it in angular coordinates makes life easier...
I am just curios here, so does Stan take these things into account?
Best,
Sebastian
Ah, and to raise my vote here as well again: bringing ODEs into Stan would be really great.
> --
> 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.
>
>
--
> 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.
>
>
Sebastian Weber
Ben Goodrich
Mar 23, 2013, 11:19:36 AM3/23/13
to stan-...@googlegroups.com
On Saturday, March 23, 2013 5:44:24 AM UTC-4, Sebastian Weber wrote:
As I am trained physicist I am wondering if Stan is automatically transforming stuff into easier geometries. I.e. I remember in my physics courses to always think long and hard about possible variable transformations before solving any problem.
There are a lot of automatic transformations in Stan that are not obvious to the user; i.e. all parameters declared with restrictions on their support are actually transformations of primitive unbounded parameters. Your general point about geometry is an important one, but it is mostly the responsibility of the user to choose a (re)parameterization of the model that is easier to sample from. The example of the funnel in the manual is instructive, which can also be seen in RStan by doing
funnel_reparam <- stan_demo("funnel_reparam")
pairs(funnel_reparam, pars = c("x_raw[1]", "x_raw[2]", "y", "x[1]", "x[2]", "lp__"))
And there is a bunch of stuff in the Optimizing Stan Code part of the manual about reparameterizations that are usually intended to make the primitive parameter space closer to a spherical normal.
Ben
Bob Carpenter
Mar 23, 2013, 1:20:10 PM3/23/13
to stan-...@googlegroups.com
We're also working on Riemann-Manifold HMC, which does
the transformation of the geometry dynamically.
Check out Michael Betancourt's paper and refs therein:
http://arxiv.org/pdf/1212.4693.pdf
The differential geometry perspective on HMC is very
enlightening.
- Bob
the transformation of the geometry dynamically.
Check out Michael Betancourt's paper and refs therein:
http://arxiv.org/pdf/1212.4693.pdf
The differential geometry perspective on HMC is very
enlightening.
- Bob
Michael Betancourt
Mar 23, 2013, 1:36:27 PM3/23/13
to stan-...@googlegroups.com
The Hamilton-Jacobi formalism you mention provides a construction for calculating the trivially-integrable action/angle variables, but only for integrable systems. Unfortunately exactly integrable systems are rare, and the method will not prove useful for most HMC problems. About the only exactly integrable systems would be derived from mixtures of Gaussians.
I've played around with trying to identify and utilize more symmetry, whether it be through Noether currents or Killing vectors, but these kinds of approaches always end up being very application specific. As Bob notes, Riemannian HMC approximates the system as exactly integrable in a local neighborhood, which is about as best we can do in a general setting.
I've played around with trying to identify and utilize more symmetry, whether it be through Noether currents or Killing vectors, but these kinds of approaches always end up being very application specific. As Bob notes, Riemannian HMC approximates the system as exactly integrable in a local neighborhood, which is about as best we can do in a general setting.
Tim Triche, Jr.
Mar 23, 2013, 2:02:34 PM3/23/13
to stan-...@googlegroups.com
oh wow! Girolami & Calderhead's paper certainly provoked a lot of discussion, will be interesting to see what advantages it confers in practice.
thank you (plural, as in "you Stan team members and other contributors") for all of the awesome things you have been adding to Stan lately. It's hard to keep up even having seen the source material. Very interesting times.
And thanks again for the semi-supervised HMM example added for 1.2.0
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
--
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+unsubscribe@googlegroups.com.
Tim Triche, Jr.
Mar 23, 2013, 2:15:26 PM3/23/13
to stan-...@googlegroups.com
The arXiv paper with the funnel figures (providing a consistent visual cue motivating Riemannian distance metrics over Euclidean) is more straightforward than Girolami and Calderhead's original. Thank you for posting it on arXiv and tackling the practical implications. In retrospect, even the animations posted by Girolami's group are now clearer in their intent.
Per John Adams, "it's not enough to write papers, others must read them," which is a lot easier when the difficult concepts are motivated by such figures
Again thanks for the clear exposition of difficult concepts.
Bob Carpenter
Mar 23, 2013, 2:56:02 PM3/23/13
to stan-...@googlegroups.com
On 3/23/13 2:02 PM, Tim Triche, Jr. wrote:
> oh wow! Girolami & Calderhead's paper certainly provoked a lot of discussion, will be interesting to see what
> advantages it confers in practice.
going to be that if the dimensionality is modest
(less than a couple thousand parameters) and the geometry
is complex (such as in the funnel), it will help
immensely.
Elsewhere, it's going to come down to our implementation
speeds and even we don't know where that's going to fall.
> thank you (plural, as in "you Stan team members and other contributors") for all of the awesome things you have been
> adding to Stan lately. It's hard to keep up even having seen the source material. Very interesting times.
>
> And thanks again for the semi-supervised HMM example added for 1.2.0
but I haven't tried it, which is why I haven't written it
up. I shouldn't even have checked it in without testing
it. I plan to verify that it works before the next release.
- Bob
Sebastian Weber
Mar 24, 2013, 6:57:21 AM3/24/13
to stan-...@googlegroups.com
Hi!
I am more than curious to see RHMC at work! Stan already blew away everything what I know and RHMC will for most problems we deal with improve things even further! You keep on writing later this year, any idea when this is going to happen? Is there a git branch available which has RHMC already implemented in a beta like version?
Best,
Sebastian
Sebastian Weber
I am more than curious to see RHMC at work! Stan already blew away everything what I know and RHMC will for most problems we deal with improve things even further! You keep on writing later this year, any idea when this is going to happen? Is there a git branch available which has RHMC already implemented in a beta like version?
Best,
Sebastian
> --
> 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.
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
Sebastian Weber
Michael Betancourt
Mar 24, 2013, 1:24:22 PM3/24/13
to stan-...@googlegroups.com
We have quite a bit of refactoring to do before we're ready with even a beta version. If you're willing to write your own potential/gradient/hessian/third derivatives in C++, however, then you can use the code referenced in the paper, https://github.com/betanalpha/jamon.
Bob Carpenter
Mar 24, 2013, 7:55:28 PM3/24/13
to stan-...@googlegroups.com
I can't be much more specific on date because the timing
involves several factors. We need to finish the forward-mode
auto-dif and Michael Betancourt needs to do the RM-HMC
implementation. But I'm currently refactoring the math/prob
lib with Daniel Lee and Michael's busy refactoring the
HMC infrastructure.
The timing's also complicated by Michael's move to the UK
to start a postdoc with Girolami. Michael already has a standalone
implementation of his RM-HMC metrics on GitHub in his package jamon:
https://github.com/betanalpha/jamon
This list will be the first to know when the Stan integration
is done and tested!
- Bob
involves several factors. We need to finish the forward-mode
auto-dif and Michael Betancourt needs to do the RM-HMC
implementation. But I'm currently refactoring the math/prob
lib with Daniel Lee and Michael's busy refactoring the
HMC infrastructure.
The timing's also complicated by Michael's move to the UK
to start a postdoc with Girolami. Michael already has a standalone
implementation of his RM-HMC metrics on GitHub in his package jamon:
https://github.com/betanalpha/jamon
This list will be the first to know when the Stan integration
is done and tested!
- Bob
Reply all
Reply to author
Forward
0 new messages