Re: benchmarking NUTS vs. JAGS vs. RM-HMC
Bob Carpenter
On 11/7/11 4:52 PM, Matt Hoffman wrote:
> The email was about the bugs rats example, which I checked in (what I
> think is) a working version of. JAGS kills us on speed, but that's
> only because it's a really easy problem and JAGS can therefore get
> nearly independent samples every iteration.
I always say that benchmarking isn't easy to
make "objective".
I was expecting this sort of speed lag for simple
examples. In fact, I brought it up with Andrew when he
was asking about BUGS/JAGS vs. Stan comparisons and with
Matt when we were talking about Stan vs. RM-HMC.
We should look better on scalability in terms of
amount of data, complexity and size of parameters,
and number of samples. Judging from the lag at
the end of sampling, I think JAGS (like BUGS) buffers all
the samples in memory, but then maybe that's just an
R interface issue (which we'll probably inherit because
R's not scalable).
What I'd like to see is intrinsic model complexity on
the horizontal axis (e.g., number of predictors, data
size, etc.) versus resource usage on the vertical axis
(e.g., time or memory to get to a given ESS), with
lines for JAGS, NUTS (straight C++ implementation and
auto-dif implementation) and RM-HMC.
I think we'll look better than JAGS on more complex models
and at larger scale. And we should look better than
RM-HMC at larger-scale models.
If it's only a compile time lag, we could work on that.
If it's something deeper, I also wrote into the proposal
that we'd build Gibbs with compilation. That should run
faster than BUGS or JAGS (less compilation time).
> I just tried to write an email to stan-dev and was told I don't have
> permission. Could someone look into that?
Daniel's working on it. I hope you're at least getting these
e-mails.
- Bob
Andrew Gelman
That said, I could see how this could still take 5x as long as simple Gibbs or Metrop, so I guess that makes sense.
Bob Carpenter
libs and we're now wondering what to do for
distributions given arguments outside their
support?
Sensible options seem to be either to
1. return 0 (or -infinity on the log scale), or
2. throw an exception
If we throw an exception, we'll probably want to make
it a special one so that someone could catch it
and go on, for instance by rejecting a sample.
Some of the more serious libs, like Boost, allow
you to configure error handling at a global level to
do one or the other, but I'd rather just go with one
for simplicity.
As Daniel brought up, in our case, it's almost
always going to be the result of either (a) some
arithmetic precision error, or (b) a model
specification error. We can try to control (a) by
the same trick Matt applied to log1p(), namely
bounding returns in the expected range.
- Bob
Daniel Lee
make a mu matrix (N x T), Stan ran a lot faster. I'm not sure if the
speed difference still exists with the new "derived parameter"
section. It might be worth trying. If no one gets to it before I have
a free minute, I'll take a look.
Daniel
Bob Carpenter
0. Looking at a real example, I'm not sure putting
the derived parameters up front helps. Opinions?
1. We're not doing ourselves any favors on these
examples by using deviation as the parameter for
the normal distribution as we wind up taking sqrt()
all over the place.
It'd be interesting to see how slight changes to
non-conjugate priors would affect the comparison.
2. As Daniel pointed out, you can replace:
for (i in 1:N)
for (j in 1:T)
mu[i,j] <- alpha[i] + beta[i]*(x[j] - xbar);
and
for (i in 1:N)
for (j in 1:T)
Y[i,j] ~ normal(mu[i,j], sigma_c);
with
for (i in 1:N)
for (j in 1:T)
Y[i,j] ~ normal(alpha[i] + beta[i]*(x[j] - xbar), sigma_c);
As is, the
double mu[N,T];
declaration will get compiled to a variable
std::vector<std::vector> mu;
and the mu[i,j] accessor compiled to mu[i][j].
There's a whole lot of push_back going on to
build it all.
Matt -- is mu[i] going to return a reference
to the embedded item? If it's copying, it's
a real problem.
3. We can replace their (x[j] - xbar) with
a new defined data item x_centered[j] = x[j] - xbar.
And we could compute xbar inside the program if it's
really just the average value of the x. That's what the
"derived data" block's meant to expedite.
4. I'd like to move to the convention where our loops look like:
for (n in 1:N)
for (t in 1:T)
y[n,t] ~ ...
Also, the use of sigma_c and alpha_c and beta_c is just crazy.
At least if they're going to do this, then
alpha[i] ~ normal(alpha_c, alpha_sigma);
should be
alpha[i] ~ normal(mu_alpha, sigma_alpha);
for consistency.
- Bob
Bob Carpenter
HMC's ass!
I checked in a revised version of the Rats
model following the conventions I laid out
(and also following Andrew's desire to have
deviation be the parameter, not variance).
And also the BDA/ARM conventions of not trying to
capitalize (or boldface) random variables,
matrices, or anything other than true constants.
For Andrew's edification, I've appended the model below.
I verified Matt's time report of about 11s for
the revised model in Stan/NUTS (g++ -O3). This was after
removing the vector<vector> allocation and
inlining it into the distribution.
JAGS is about 100 times faster here,
getting 3500 samples in a tenth the time it
takes us to get around 350, and that's not even
counting model compilation time. Ouch.
Then again, we may want to be a bit easier on ourselves by
not throwing away the initial half of the samples
in a case like this that's likely to converge quickly.
At least when we're doing performance comparisons.
That might get the time closer to a factor of 50 :-)
I'll have to ask Matt about profiling when we're both in
the office -- I want to see where the time's
going now. It could be largely in auto-dif, but
I wouldn't think it'd be an overall factor of 50.
Using g++ -O3 takes about 12s -- in line with what
you're reporting. So I don't think the std::vector
accesses were the issue.
- Bob
# http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/Vol1.pdf
# Page 3: Rats
data {
int(0,) N;
int(0,) T;
double x[T];
double xbar;
double y[N,T];
}
parameters {
double alpha[N];
double beta[N];
double mu_alpha;
double mu_beta;
double(0,) sigma_y;
double(0,) sigma_alpha;
double(0,) sigma_beta;
}
derived parameters {
double(0,) sigmasq_y;
double(0,) sigmasq_alpha;
double(0,) sigmasq_beta;
sigmasq_y <- sigma_y * sigma_y;
sigmasq_alpha <- sigma_alpha * sigma_alpha;
sigmasq_beta <- sigma_beta * sigma_beta;
}
model {
mu_alpha ~ normal(0, 100);
mu_beta ~ normal(0, 100);
sigmasq_y ~ inv_gamma(0.001, 0.001);
sigmasq_alpha ~ inv_gamma(0.001, 0.001);
sigmasq_beta ~ inv_gamma(0.001, 0.001);
for (n in 1:N) {
alpha[n] ~ normal(mu_alpha, sigma_alpha);
beta[n] ~ normal(mu_beta, sigma_beta);
for (t in 1:T)
y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), sigma_y);
}
}
Andrew Gelman
Andrew Gelman
2. I agree that a factor of 50 is pretty embarrassing, although I suppose it makes sense if Nuts is doing 50 function evaluations per iteration. In other examples (e.g., involving matrices), Jags/Bugs can just grind to a halt.
3. What would happen if Stan were running simple random-walk Metropolis? I guess there's something I'm still not understanding about Nuts. Nuts is trying to optimize, so if Metrop jumps are pretty close to independent, shouldn't Nuts be just doing 1 or 2 steps per iteration? Or is it just that the gradient calculation itself is slow? If one gradient evaluation takes 10 times as long as one funciton evaluation, then any version of HMC will lose in simple problems. Which suggests that maybe Stan should have a Metropolis option for easy problems. I dunno.
Matt Hoffman
they take. The difference between waiting one second and 10 seconds is
small compared with the amount of time needed to edit a text file.
I think random-walk Metropolis would do substantially worse than
NUTS/HMC, since Langevin (i.e. HMC with 1 step) does terribly. The
"problem" is that Gibbs is actually really fast for this model. (If
only we had this problem for all models!) Since the amount of work
needed for one well implemented Gibbs step is comparable to the amount
of work needed for one likelihood evaluation for this model, it's hard
to imagine anything doing much better than Gibbs.
Matt Hoffman
wrote but couldn't send earlier tonight with some timing evaluations.
(I think maybe Bob but nobody else saw this.)
Matt
On Mon, Nov 7, 2011 at 6:39 PM, Matt Hoffman <mdho...@cs.princeton.edu> wrote:
> I looked into it a little more, and discovered that NUTS is using
> around 70 steps per iteration. I tried HMC with a few other numbers of
> steps per iteration, and got worse results:
> Sampler Time Effective Samples
> NUTS: ~11s ~350
> HMC-10: ~2s ~2
> HMC-77: ~12s ~110
> HMC-100: ~15s ~280
> HMC-200: ~30s ~314
>
> By contrast, JAGS gets about 3600 effective samples (after 10K Gibbs
> iterations with 1K burnin) in about 1s. So it's possible that we'd be
> competitive with JAGS if we weren't paying the autodif penalty. But a
> hand-coded Gibbs sampler might still beat our pants off.
>
> Basically, I think the deal is that hierarchical models make HMC need
> to take small steps. If Gibbs is able to generate effectively
> independent samples, then it's hard to see how anything could beat it.
> But of course, the point of this project is to deal with _hard_
> problems, not problems for which JAGS generates thousands of
> effectively independent samples in under a second.
>
> BTW, cutting out the mu matrix didn't seem to make a difference,
> performance-wise.
>
> Matt
Matt Hoffman
> There's a whole lot of push_back going on to
> build it all.
>
> Matt -- is mu[i] going to return a reference
> to the embedded item? If it's copying, it's
> a real problem.
I'm pretty sure it'll return a reference if it's appropriate, since
it's the same operator you'd use to say something like
mu[i][j] = 3;
> 4. I'd like to move to the convention where our loops look like:
>
> for (n in 1:N)
> for (t in 1:T)
> y[n,t] ~ ...
>
> Also, the use of sigma_c and alpha_c and beta_c is just crazy.
> At least if they're going to do this, then
>
> alpha[i] ~ normal(alpha_c, alpha_sigma);
>
> should be
>
> alpha[i] ~ normal(mu_alpha, sigma_alpha);
>
> for consistency.
Hear hear. I (and I suspect Daniel) was just copying as mindlessly
from the bugs example as possible.
Matt
Matt Hoffman
0/-inf is that at least the system will produce something (something
correct, I'm pretty sure), whereas an exception may lead to users
giving up. Maybe we could print out a warning?
I think we should definitely take some pains to prevent this situation
from coming up to the extent that's possible.
Matt
Andrew Gelman
In any case, i agree there's no need to obsess over models that are already fast. The real issue will come when we fit big hierarchical models. If these are slow with Stan, we'll have to do more thinking. And at that point it might help to have developed some intuition from the easy problems.
Finally, 1 sec vs 10 sec is no big deal, but 1 sec vs 50 sec is starting to make a difference. 50 seconds can be a long time to wait, if we're fitting models over and over again.
Matt Hoffman
and Metropolis/HMC/NUTS aren't. That could be something that's worth
looking into.
Andrew Gelman
Matt Hoffman
Matt Hoffman
took about 1/7 as long as using autodiff.
Matt
Bob Carpenter
the stamina to interleave all these comments.
0. A factor of 7 for auto-diff vs. direct seems pretty
good to me given that we're doing allocation and managing
the stack. It should get faster when we override the probability
distros in two ways. First, we'll save memory and speed by partially
evaluating results (fewer chainables on the stack holding fewer refs
with fewer steps in derivatives). Second, we won't have to do the
extra calcs for the fixed normal(0,100) mean and deviation
parameters (though that's pretty small in this case). In fact,
this reminds me that we should treat centered distros separately, as
they remove a lot of the inner arithmetic -- (x - 0) * (x - 0)
is way worse in auto-dif terms than x * x.
We could approach hand-written speed, I think, if
we did code generation rather than templated taping.
Maybe that's what Paul Hovland at DoE was talking to
Andrew about -- it's what their tools ADIFOR and ADIC
do for Fortran. It's much much harder to do for C++ than for
Fortran or even C, and requires yet another code-gen/compilation
step.
1. Matt -- that's nice analysis of Gibbs vs. Metropolis vs.
Langevin vs. HMC vs. NUTS. It should go in a textbook
somewhere!
What would be nice would be to see a setup where all
the distros are conjugate, but it's still slow w. Gibbs.
What do you mean that Gibbs is scale-invariant? That it
doesn't matter what scale the params are on because you're
doing proper conditional probs? Won't it matter in practice
if you use something like a step-out slice sampler?
And isn't the problem that the scale varies by parameter?
And isn't Riemann-manifold HMC just an attempt to automatically
calc the scale?
2. Now that I think about it, mu[i] has to be a reference,
because it can be on the LHS of an assignment.
3. I brought up the loop issue in case someone thought we
should stick closer to the way the BUGS models were coded.
I'd like to use them as examples of how to code models in
Stan. As Daniel suggested, we should supply a "translation
manual" for helping BUGS users get started. Doing it by
BUGS examples is ideal, so we'll kill two birds with one stone.
4. As to exceptions vs. 0/inf for out-of-support arguments,
you could always catch the exception inside of Stan's HMC
loop and treat it as a complete 0 probability step. It'd
probably even be more efficient than continuing even given
the cost of dispatching exceptions.
Do some of these distributions have well-defined behavior
in the limit of parameters going to 0 or inf? If so, we
could always just implement that behavior. That'd be in
keeping with IEEE math.
We can't just print something, as that's hugely time
consuming and hard to control. That's what JAGS does and it
takes ages for the megabytes of warning messages to print out.
Even if there's a log4j-like package for C++, that's
still costly compared to function evals. We could have a
DEBUG mode controlled at compile time, but that's a huge
pain for users.
There's really no way we can prevent overflows.
A user can just multiply in a loop.
5. In UI reaction time, there's a big jump from
a few milliseconds to the half-second or so when it
becomes really apparent an operation's not "instantaneous".
At that point, if you go beyond a few seconds, users will
check their e-mail, and if you go beyond a minute or two,
they'll find something else to do.
- Bob
Matt Hoffman
memory efficient, I think). It's tricky in general, but I suspect not
as hard for the average-case Stan model, which will usually just be a
sum of predefined functions of (possibly transformed) model variables.
Something to think about, especially because the memory cost of taping
poses real limitations on model size.
In this example, there's a surprising amount of overhead associated
with repeatedly computing log(sigma). It's not huge, and it'd probably
be a pain to eliminate, but we should add that sort of thing to the
list of sources of suboptimal efficiency.
BTW, on os x the "Time Profiler" in the Instruments application is a
really nice profiler. I think it's meant to replace Shark, which is
what I used to use. That's how I got alerted to the log(sigma)
business.
By "Gibbs is scale-invariant," I meant in the case where you've got
analytic conditionals for everything. It's a little less
scale-invariant when slice sampling using the doubling method, and
much less if you're using the stepping-out method. RMHMC does address
the issue of different parameters having different scales, but it goes
much further and (effectively) rotates the space so that the
parameters are locally uncorrelated. If you're dealing with a
multivariate normal target distribution, it basically tries to turn it
into a standard isotropic normal.
I think most distributions have sort of well defined behavior for
extreme parameter values—normals turn into deltas or improper uniforms
as sigma goes to 0 or inf, t's turn into normals or deltas as nu goes
to inf or 0, gammas turn into deltas as shape goes to 0 or rate goes
to inf, etc. IEEE math won't necessarily make these things happen like
they should by default. There's a little bit of a performance hit for
enforcing correct behavior, but it's probably pretty minor in c++
compared with evaluating logs and exps.
Matt
Bob Carpenter
On 11/8/11 2:22 PM, Matt Hoffman wrote:
> I think most distributions have sort of well defined behavior for
> extreme parameter values�normals turn into deltas or improper uniforms
> as sigma goes to 0 or inf, t's turn into normals or deltas as nu goes
> to inf or 0, gammas turn into deltas as shape goes to 0 or rate goes
> to inf, etc. IEEE math won't necessarily make these things happen like
> they should by default. There's a little bit of a performance hit for
> enforcing correct behavior, but it's probably pretty minor in c++
> compared with evaluating logs and exps.
We should definitely think about this for the distributions.
If the performance hit's too much, we can have two versions
of all of the functions, one that checks its params and one
that doesn't. In fact, we could make checking a template
boolean param with a default to check the ranges.
Daniel?
> Code generation would definitely make things faster (and much more
> memory efficient, I think).
We'd still have to have all the intermediate variables
and adjoints built up as variables only now they'd be
allocated on the program stack instead of our own stack.
For instance, with (a*b), we push a vari onto the stack
holding the value of (a*b) and the partial of this expression
w.r.t. the root expression (which gets instantiated in
the backward pass). We'd need both these as local variables
on the execution stack with code generation. We also
store the pointers to a and to b on the vari object. Those
would instead be linked inside the code itself, which would
almost certainly be more efficient in time (in space?) than
managing the derivative propagation stack and memory management
ourselves.
With optimization, the compiler could almost certainly optimize
our generated code to get rid of things that aren't used.
But most variables are used if they're auto-dif variables;
if they're constants, we don't put them on the stack either.
> It's tricky in general, but I suspect not
> as hard for the average-case Stan model, which will usually just be a
> sum of predefined functions of (possibly transformed) model variables.
> Something to think about, especially because the memory cost of taping
> poses real limitations on model size.
If we limit to BUGS-type models, it would be tractable to
code generate the auto-dif stuff. Where it gets tricky is
with iterative algorithms like some of the matrix solvers.
> In this example, there's a surprising amount of overhead associated
> with repeatedly computing log(sigma). It's not huge, and it'd probably
> be a pain to eliminate, but we should add that sort of thing to the
> list of sources of suboptimal efficiency.
We should be recalculating log(sigma_alpha) and log(sigma_beta)
in each loop for n in 1:N. And additionally log(sigma_y) for
each n in 1:N and each t in 1:T.
for (n in 1:N) {
alpha[n] ~ normal(mu_alpha, sigma_alpha);
beta[n] ~ normal(mu_beta, sigma_beta);
for (t in 1:T)
y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), sigma_y);
}
I just realized from your comment that we could compute these
once and reuse the computations in the loop. That, I don't know
how to automate.
> BTW, on os x the "Time Profiler" in the Instruments application is a
> really nice profiler. I think it's meant to replace Shark, which is
> what I used to use. That's how I got alerted to the log(sigma)
> business.
Thanks -- I'll check it out.
> By "Gibbs is scale-invariant," I meant in the case where you've got
> analytic conditionals for everything. ...
OK -- that's what I thought.
> It's a little less
> scale-invariant when slice sampling using the doubling method, and
> much less if you're using the stepping-out method. RMHMC
> ... basically tries to turn it into a standard isotropic normal.
That I get. In fact, that's the way I understand
covariance matrices in general and from that, it
helped immensely with the rest of linear algebra.
- Bob