NUTS, non-diag mass, vs. other samplers
Bob Carpenter
gave posting priveleges.
I'll have to look up the "DE" method by Vrugt and Braak
(with names like those, I'm expecting big things).
That must be one hairy model if you can only do 100K
gradient calcs/hour!
More on the content later. I'm still trying to catch
up with our initial release issues.
And thanks for posting.
- Bob
----------
I have just become a member of stan dev list, and have been looking through messages. I have some comments/questions,
and wondered if I could have posting privileges or if there was some other suitable means of communication.
Let me explain a bit - I am developing commercial software for engineering design and Bayesian uncertainty forecasting.
This involves creating a proxy/surrogate/kriging model of a large engineering simulation package, and then applying MCMC
to determine probabilistic forecasts. I have been doing this for over ten years, but am currently upgrading my MCMC
methods. NUTS seemed to fit the bill very nicely, until I started applying it to my problem this week!
I am not in academia, but in business, but I have been working on optimisation techniques for a long time.
My likelihood function is very complex, with multiple minima and lots of curved narrow valleys. I have up to 400
parameters, but normally about 50. Derivatives are available.
I have implemented basic RW, NUTS, HMCMC, DE (Vrugt and Braak) and parallel tempering (PT), and have implemented various
combinations (e.g. DE with PT). My code is fully integrated so it is very easy to switch methods, or use different
methods at different steps, or combine methods.
The goal is a software product which is completely automatic, no user tuning of parameters, the user does not even know
MCMC is being implemented. hence the attraction of NUTS. There is a practical limit on how long a user will wait - an
hour is about the maximum. This means in practice about 100,000 function/derivative calculations.
I have converted Matlab NUTS and DE code to Java.
My current findings are that NUTS is very good on nice smooth functions with a single mode, but can fail badly on more
complex functions. In contrast, DE seems to cope well with multiple modes and complex functions.
Now, the kinds of discussions I would like to have with your colleagues are ways of improving NUTS. I am not interested
in the STAN infrastructure, just the core algorithm. They fall into several categories:
1. Simple changes.
Often the step size is small and the number of leap frogs becomes large. The step size keeps on falling because it is
difficult to find better points. The algorithm seems to easily get stuck in local minima. I know there are better points
because the espilon estimation found them, but the leap frogs cannot. Stan uses a variation on epsilon, and I am
currently trying a gamma distribution on epsilon. However, Stan seems to do this when not in burn in phase, which seems
strange as I have problems within burn in.
I am also allowing a re estimation of epsilon if the limit on the max number of depths is reached. I redo the initial
epsilon estimation. This is a bit similar to DE where the step length is sometimes set to 1.0 to allow jumps to
alternative modes.
My epsilon estimation takes the average of 10 trials, because of the complex curvature.
An obvious improvement would be to first find a good starting point - this is my default - but I really want a method
which is robust enough to handle any starting point, and explore the space fully.
2. Integration with DE.
Because DE and NUTS seems to complement each other, it would be nice to have a method which combines them. For example,
run parallel chains, one of them NUTS the others DE. For the NUTS chain, it may be possible to use a DE approach to set
the momentum, instead of doing a Gaussian. I need to check how this would work, but it would ameliorate the problem of
steps which keep hiting high mountains. The DE approach is self adapting in that the points lie along valleys so
differences are also along valleys.
3. Integration with PT.
I have tried PT with DE, and it is useful, but I fear PT with NUTS would be too time consuming. Because the NUTS chain
is much more expensive, because of the multiple leap frogs, it may be useful to have other DE chains which include PT
and this will not affect performance too much.
4. Better Hessian/mass matrix.
A diagonal mass matrix is no use to me, but recent studies using BFGS and/or Reimann manifold may be beneficial. They
may not, however, help the multi mode issue. I have tried NUTS with a simple highly correlated multivariate gaussian
distribution with an exact weight matrix, and the efficiency improves drastically.
Are these kinds of discussions suitable for the mailing list? A lot of mailing list discussions are around more basic
software issues, rather than the algorithm itself.
------------------------
essence
You can call me Nigel.
My model is a kriging/surrogate/proxy/radial basis function/Guassian process model.
In fact in the example I have been using, there are about 300 models, and the overall likelihood calculates the likelihood for each individual model.
Each individual model is based on about 600 big engineering simulation runs (which can each take 24 hours in some cases).
I use Matern correlation functions which take some time, and calculating the derivatives was fun.
So for the overall likelihood function calculation I do about 300 Cholesky solves with a 600 x 600 triangular factor. The derivative calculation takes about the same time as a function evaluation.
I do it all in heavily multi threaded Java code.
Re NUTS - I did changes to the epsilon as mentioned below, and it definitely improves things. Giving the epsilon a jolt if it gets stuck at small values help a lot. The gamma distribution may also help. Being in commerce, I don't do days and days of testing, I use my instincts a lot.
Next stage I think will be to run parallel DE chains and do some cross fertilisation a la tempering.
BTW, I find that out of the box it takes of the order of 1000 leap frog steps per iteration, so there is plenty of scope for using an approximate hessian. Using the Hessian will also make it less sensitive to the size of the momentum gaussian distribution, which is very critical.
The BFGS paper is not as well known but is also very nice work - Zhang and Sutton.
essence
On Friday, 31 August 2012 20:44:31 UTC+1, Bob Carpenter wrote:
Andrew Gelman
A
>
Ben Goodrich
4. Better Hessian/mass matrix.
A diagonal mass matrix is no use to me, but recent studies using BFGS and/or Reimann manifold may be beneficial. They
may not, however, help the multi mode issue. I have tried NUTS with a simple highly correlated multivariate gaussian
distribution with an exact weight matrix, and the efficiency improves drastically.
Another idea might be to put what we call the LKJ prior on the mass matrix (presumably with shape >> 1). In that case, I think there would need to be an additional step where the algorithm decides whether to accept or reject a proposal for the mass matrix, then draws the momentum variables from a multivariate normal with that mass matrix, then update the parameters given the momentum variables ...
But limited-memory BFGS and the Reimann manifold thing sound interesting too.
Ben
essence
It can be applied to discrete parameters, although I haven;t tried this yet - I will need the capability in the near future (i.e. next week).
The main issue is that the basic approach needs 2 x N chains, but more recent work has helped. As I sometimes have 400 parameters, this is an issue. Even with 16 parameters (which is for me on the very small side), 32 chains is quite a lot. I tend to use 8 chains, it used to be so that each chain can be on a separate CPU core, but I am doing it differently now and 8 just seems a sensible number.
A lot of the ideas come from the differential evolution method in heuristic optimisation, which I have used a lot and found very useful (much better than the 'conventional' GA operators). It works, but nobody knows why!
Andrew Gelman
Matt Hoffman
I'll try to respond to some of these questions.
To clarify, Stan adapts the step size epsilon during burn-in, but not
afterwards, since adapting MCMC parameters results in a
non-time-homogeneous Markov chain which will have the wrong stationary
distribution. (I.e., you'll get samples from a distribution besides
the target posterior.) There is an option to jiggle the step size, but
the danger of doing this with NUTS is that if you use too small of a
step size then you'll wind up with a very long (and expensive)
trajectory.
NUTS tries to find the largest step size that doesn't lead to
unacceptably inaccurate simulations. The goal is to produce results
that are both as fast as possible and unlikely to produce incorrect
results—a problem with using a step size that's too large is that you
may wind up never visiting certain parts of the parameter space even
if they have a decent amount of posterior mass.
Ultimately, it's very difficult to evaluate the performance of an MCMC
algorithm on a multimodal posterior. The way I think about it is that
(a) sampling is harder than optimization (b) finding global optima is
intractable for multimodal objective functions so (c) sampling from
multimodal posteriors is intractable. So in the typical case you're
unlikely to ever get even a single "effectively independent" sample
from a multimodal distribution, since to do that you'd basically have
to be able to show that you've given each isolated mode its due
consideration.
For many practical problems, it's nonetheless useful to get a few
samples (or even one sample) from a single reasonably high-mass mode.
Even if that's a terrible summary of the posterior, it can still be a
useful result. If that's all you need, then I'd argue that you don't
actually need to use an MCMC algorithm that satisfies detailed balance
and is guaranteed to generate samples from the true posterior. Lifting
the detailed balance restriction dramatically increases the world of
algorithmic choices you can make.
Having said all that, there's certainly still plenty of room to make
NUTS better while still satisfying detailed balance! I'm not very
familiar with DE, but it's probably worth looking into. My
understanding is that it generates proposals roughly in the direction
of the vectors between points in parallel chains, an idea that's quite
closely related to Zhang and Sutton's BFGS HMC algorithm. A limitation
of these methods is that they aren't able to do any local
adaptation—if you've got a bunch of chains in different isolated
modes, then it seems likely that you're not going to get information
from chains 1 and 2 them that's relevant to the mode that chain 3 is
in. My guess is that these sorts of approaches are probably not much
more useful than using a good fixed linear preconditioner matrix
(equivalent to using a fixed covariance matrix for the momentum
distribution). But it's nice that they satisfy detailed balance and
don't require O(N^2) computation/storage.
Riemannian Manifold HMC (using the Fisher metric of Girolami and
Calderhead) is able to do local adaptation, but requires 1. a
log-concave prior (which rules out most hierarchical models) 2. the
ability to compute the derivative of the metric, which involves
computing third derivatives of the log-prior and 3. solving an N-by-N
system of equations for each step. A version of RMHMC/RMNUTS that
didn't suffer from these restrictions but still used local information
would be great, and is something I and others are actively looking at.
Matt
>
>
Nigel/essence
I am pretty familiar with NUTS, having converted it from Matlab to Java (not trivial with the recursion!).
The jiggling of epsilon I am referring to is in lines 156-161 of nuts.hpp, and if I understand it, the jiggling only happens AFTER burn in, which is what I was wondering about. Have I understood correctly, and is there a reason for this? It seems strange - I definitely want to jiggle during burn in, otherwise I will never get out of burn in.
The papers on DREAM are well worth the effort, and Gregory has found it very useful for his work finding aliens. Look up Gregory and FUSION, he is from Toronto, and has similar background to me (I also did theoretical physics PhD at Manchester).
Bob Carpenter
On 8/31/12 5:25 PM, Nigel/essence wrote:
> Thanks - there seems to be some good discussion arising. I am based in UK, so will not attempt to respond much this evening.
> I am pretty familiar with NUTS, having converted it from Matlab to Java (not trivial with the recursion!).
>
> The jiggling of epsilon I am referring to is in lines 156-161 of nuts.hpp, and if I understand it, the jiggling
> only happens AFTER burn in, which is what I was wondering about.
> Have I understood correctly, and is there a reason for this?
> It seems strange - I definitely want to jiggle during burn in, otherwise I will never get out of burn in.
in the parameter that adjusted the step size epsilon
+/- some percentage around the adapted size.
I believe the conclusion was that we thought it might
confuse the adaptation, which was already working well.
Matt didn't think moving epsilon was going to help,
but since then, we've run into some more complicated
geometries where moving it way down every few iterations
seems like it might help.
> The papers on DREAM are well worth the effort, and Gregory has
> found it very useful for his work finding aliens. Look up Gregory and FUSION,
> he is from Toronto, and has similar background to me (I also did theoretical physics PhD at Manchester).
He seems to be one emeritus professor who's still going
gangbusters, which I think is fantastic (because that's where
I want to be some day).
I'll see if I can get through the papers. These have been
more in Matt Hoffman's and Michael Betancourt's domain
(Michael's also a theoretical physicist); with their
help in front of a blackboard, I can almost feel like I'm
getting the easier parts.
What is always amazing to me is that after the dust
settles, the algorithms are usually pretty straightforward.
(NUTS is involved, but I'm a computer scientist, so recursive
binary tree algorithms are comfort food.)
- Bob
Bob Carpenter
On 8/31/12 5:12 PM, Matt Hoffman wrote:
> ... I'm not very
> understanding is that it generates proposals roughly in the direction
> of the vectors between points in parallel chains, an idea that's quite
> closely related to Zhang and Sutton's BFGS HMC algorithm.
sampler to Stan soon.
of Goodman and Weare's sampler, which
http://msp.org/camcos/2010/5-1/p04.xhtml
It samples on the ray connecting a pair of samples on
parallel chains. It does have O(N^2) memory and update
cost where N is the number of parameters, because the current points
on each chain need to span the parameter space and all need
to be updated each iteration. So thousands
of parameters will probably be a practical max on current hardware.
On the plus side, it's rotation and scale invariant,
doesn't require gradients, is very easy to parallelize,
and generates N parallel chains.
I still don't quite understand the convergence evaluation or
the effective sample size calcs in this situation, though.
There's some discussion in their paper, but I don't see how
it deals with the fact that you have correlations across chains.
Matt and Michael B. weren't too enthusiastic when I explained how
it worked, but I can't recall their concerns. But it's so
simple to implement and new samplers are so easy to integrate
into Stan (which exposes log prob and gradient functions with
unconstraining transforms built in) that I can add it and evaluate
all the BUGS models pretty easily.
- Bob
Andrew Gelman
>
Nigel/essence
My Java code for jiggling with a gamma distribution, using apache Java code for the gamma, is:
double shape = 2.0;
double scale = 1.0;
GammaDistributionImpl aGammaDistributionImpl = new GammaDistributionImpl(shape, scale);
for (int iChain = 0; iChain < nChains; iChain++) {
try {
eps[iChain] = epsilon[iChain] * aGammaDistributionImpl.inverseCumulativeProbability(Math.random()) / shape / scale;
} catch (MathException e) {
System.out.println(e.getMessage());
}
}
Does google groups have [code] tags?
I think this is a bit 'cleaner' than asking a user to specify limits. In my tests it seems to help somewhat. Occasional small steps as well as occasional long ones. It doesn't seem to affect epsilon adaptation too much. The main problem before was that it was always rejecting trees, so epsilon kept on getting smaller, and it never got out of the local minimum.
I also implement the re finding of a reasonable epsilon. I do this if (a) the depths have hit the limit (I use a limit of 8 - does this affect ergodicity etc.?) and (b) alpha/nalpha is less than 0.01 and (c) we are still in adapting phase.
I also set mu = log(epsilon) rather than log(10*epsilon) when re finding reasonable epsilon. I also have to be careful to do this not to soon before the end of burnin, or else you will have a poor epsilon after burnin.
I also have a loop around finding a reasonable epsilon, and take the average of ten attempts. Maybe I should do more. The idea is that there will be some directions with a much longer epsilon.
Of course, a very simple approach to the hessian would be to do lots of initial attempts, similar to finding reasonable epsilon, and then construct a correlation matrix from those attempts. We have a set of vectors and a set of lengths, so it should be easy to make an approximate hessian.
I looked at Lewandowski's PhD thesis, but it will take a lot for me to fully grasp the concepts. In contrast, constructing a BFGS approximation is bread and butter for me. I do a lot of regression and model fitting, but am not a trained statistician, so copulas are unknown territory.
No doubt these tweaks could be improved, but they definitely have changed my problem from one which was getting nowhere to one which gets a reasonable solution overnight.
Goodman and Weare are unfamiliar to me, thanks for the link.
Yes, it's amazing how few lines of Matlab NUTS takes (and it explodes with all the infrastructure of Stan, which is far too much for me to have time to understand). You can see all the algorithm on the screen.
In contrast to NUTS, Fortran code such as Powell's BOBYQA (200 lines) and L-BFGS-B (4800 lines) is an order or two of magnitude more complex, albeit they have a lot of comment code.
Code for BFGS MCMC is available, in Java, but it will take me a bit of work for me to understand it. There must be a lot of duplication of functions with L-BFGS-B, which I have recently converted from Fortran to Java.
It is so easy to end up with ten different cholesky factorisation methods.....
Bob Carpenter
On 9/1/12 6:39 AM, Nigel/essence wrote:
> Jiggling is always good, particularly to get out of a rut. On the other hand, jiggling could be considered synonymous with a rut....
but we haven't experimented with it in the NUTS context
much.
> My Java code for jiggling with a gamma distribution, using apache Java code for the gamma, is:
>
> double shape = 2.0;
> double scale = 1.0;
> GammaDistributionImpl aGammaDistributionImpl = new GammaDistributionImpl(shape, scale);
>
> for (int iChain = 0; iChain < nChains; iChain++) {
> try {
> eps[iChain] = epsilon[iChain] * aGammaDistributionImpl.inverseCumulativeProbability(Math.random()) / shape / scale;
> } catch (MathException e) {
> System.out.println(e.getMessage());
> }
> }
strongly advice against anywhere in the API, either
(a) printing to System.out, or (b) just swallowing
an exception without acting. What's the state
of eps[iChain] after an exception and what does that do
to the rest of the program? And if someone embeds
the code (especially in a server), where does System.out go?
See Josh Bloch's awesome book, Effective Java, for
more advice on exceptions.
> Does google groups have [code] tags?
> I think this is a bit 'cleaner' than asking a user to specify limits. In my tests it seems to help somewhat. Occasional small steps as well as occasional long ones. It doesn't seem to affect epsilon adaptation too much. The main problem before was that it was always rejecting trees, so epsilon kept on getting smaller, and it never got out of the local minimum.
> I also implement the re finding of a reasonable epsilon. I do this if (a) the depths have hit the limit (I use a limit of 8 - does this affect ergodicity etc.?) and (b) alpha/nalpha is less than 0.01 and (c) we are still in adapting phase.
a size-one tree, you just get Langevin.
> I also set mu = log(epsilon) rather than log(10*epsilon) when re finding reasonable epsilon. I also have to be careful to do this not to soon before the end of burnin, or else you will have a poor epsilon after burnin.
>
> I also have a loop around finding a reasonable epsilon, and take the average of ten attempts. Maybe I should do more. The idea is that there will be some directions with a much longer epsilon.
>
> Of course, a very simple approach to the hessian would be to do lots of initial attempts, similar to finding reasonable epsilon, and then construct a correlation matrix from those attempts. We have a set of vectors and a set of lengths, so it should be easy to make an approximate hessian.
> I looked at Lewandowski's PhD thesis, but it will take a lot for me to fully grasp the concepts. In contrast, constructing a BFGS approximation is bread and butter for me. I do a lot of regression and model fitting, but am not a trained statistician, so copulas are unknown territory.
enough for me. Ben understands it all.
> No doubt these tweaks could be improved, but they definitely have changed my problem from one which was getting nowhere to one which gets a reasonable solution overnight.
>
> Goodman and Weare are unfamiliar to me, thanks for the link.
approach the DREAM paper called DE.
> Yes, it's amazing how few lines of Matlab NUTS takes (and it explodes with all the infrastructure of Stan, which is far too much for me to have time to understand). You can see all the algorithm on the screen.
(b) computing auto-dif, (c) vectorizing, and (d) parsing the
modeling language and running it.
> In contrast to NUTS, Fortran code such as Powell's BOBYQA (200 lines) and L-BFGS-B (4800 lines) is an order or two of magnitude more complex, albeit they have a lot of comment code.
writing the language idiomatically (as opposed to
comments about code, like API specs and algorithm descriptions).
But there's a ton of things we could do to unfold and speed
up the implementations. It's like comparing Kiss-FFT to
FFTW. Same algorithm, but FFTW is much much more complex
and about twice as fast.
> Code for BFGS MCMC is available, in Java, but it will take me a bit of work for me to understand it. There must be a lot of duplication of functions with L-BFGS-B, which I have recently converted from Fortran to Java.
> It is so easy to end up with ten different cholesky factorisation methods.....
balance.
- Bob
Matt Hoffman
you're finding that setting epsilon too low is getting you stuck! One
scenario I can imagine that happening in is:
1. You start somewhere very very far away from where the sampler
ultimately ends up.
2. So the sampler has to get from a point with extremely high
potential energy to one with much lower potential energy.
3. Doing that with accurate simulations could require many, many
resamplings of the momenta, since you can only transform a limited
amount of potential energy into kinetic energy before the simulation
gets inaccurate.
4. So NUTS drives the step size ever downward while warming up,
eventually hitting the maximum number of iterations, even though using
an inaccurate simulation would be fine for the earliest "search" phase
of warmup.
Do you think that's what's happening here? This phenomenon does seem
to argue for using some kind of more intelligent initial
optimization-based search phase to initialize the sampler. Could be a
nice application for a general-purpose black-box variational inference
algorithm like the one by Gershman et al. (2012; and I'm not just
saying that because I'm the second author!).
Matt
>
>
Nigel/essence
Your epsilon estimation could also maybe benefit from some of the search length methods used in optimisation - I have had cases where your code goes into an infinite loop, so I set a max iterations in the find_reasonable_epsilon.
This behaviour can be found in the Rosenbrock example I posted on the other thread. If you start far from the minimum, it behaves badly. I was trying a dimension of 32. It would be very easy for you to try it out with Matlab NUTS.
In my real work, I do a lot of optimising to find a good starting point for sampling. But this can be dangerous, particularly if you don;t move far from that point. It results in forecast uncertainties which are too narrow.
I also do the full MCMC many times, but that is another story. But it does mean I can use my previous MCMC results to make intelligent starting points for the next time.
It's getting rather late in the UK. Gnite.
Nigel/essence
As always, it is a balance between exploration and exploitation.
There may be some useful ideas from the particle swarm literature, for which there is a well developed theory of stability and oscillations, based on eigenvalues of the transformation matrix. Some kind of inertia or dampening may be useful during burn in.
In my simple mind, if we roll a ball down the hill, we don;t want it to go all the way up the other side, but we don't want it to stop at the bottom either. We want some suitable friction or dampening. We certainly don;t want it to go round in circles like a trick cyclist.
Andrew Gelman
Nigel/essence
In my optimisation work (all unpublished, it is commercially confidential), I have integrated methods from particle swarm (PSO), differential evolution (DE) and local quasi newton minimisation. It is fully integrated, not just switching from one to the other at certain iterations. These different methods all contribute usefully.
Similarly, I see HMCMC methods as being good for local search, whereas more heuristic methods like DREAM are better for global exploration.
I did start looking at whether ideas of quantum tunneling could be used with HMCMC, to go through barriers, but didn't get very far. Parallel tempering with HMCMC was, in my initial studies, too costly, but I have often thought of having one chain as HMCMC and other chains as parallel tempered DREAM.
For multi mode problems, there may be a need for multiple tunings of NUTS.
Bob Carpenter
On 9/3/12 8:04 AM, goodwin wrote:
> I have done a lot of work developing methods will will handle very small values on normal distributions, and very fast
> approximate schemes for cum norms. Languages like R too easily set these things to zero at extremes, and then you can't
> escape. My approximations are very happy being 1000 sd's from the mean. These things are very important for robust code.
> My actual problem has a starting point of log likelihood of minus many thousands.
> Maybe something to think about for Stan if you haven't already, but you are already generally dealing with logs.
Really low log probabilities aren't a problem.
> Simple example -
> pnorm(-1000, log.p = TRUE)
> is OK in R, but
> 1.0 - pnorm(1000, log.p = TRUE)
> returns one so you have to use the tail instead of the head.
Right. We use the usual C++ log1p() functions with
first or second-order Taylor approximations for
values of x close to zero in log1p(x) = log(1 + x).
But I'm sure we have problems in places where we
haven't thought through all the algebra and recoded
things that look like log(1 + x) as log1p(x).
> But sometimes I want the area under the cum at tails, and then I need
> pnorm(-1000) - pnorm(-1001)
> so then you need to do some tricks.
This kind of thing could be a problem for us in truncated
models where we have to compute log diffs. But then
if you're that far out on a truncation/censoring model,
underflow's probably the least of our concerns.
We've calculated log sums of exponents in the usual way
in cases we've come across or are part of built-in
operations. But we haven't been too careful with the
truncated distributions involve the cumulative distribution
functions yet.
> I don't know how robust your C libraries are to these things. When these problems occur, you don't see any obvious
> manifestation.
That depends on what is happening. If we wind up
underflowing to zero for a parameter that has to be
positive, we'll trap the domain errors, throw it, and
the sampler can report it. Or we'll let it go and do
the right thing in the technically undefined limits.
At least that's the theory. I'm sure we still have a long
way to go. Hopefully the users will help us out with bug
reports.
- Bob
Nigel/essence
Garsaglia has done some interesting work in this area 'Evaluating the Normal Distribution' in JSS.
But even that was nowhere near enough for my tails, so I used the concept and wrote my own approximation.
You can see from that paper how to get log cum norm approximations at extreme tails.