Understanding Model-Specific Speed Decrease between Stan 1.3 and Stan 2.x

[Avi Feller]

Dear wonderful Stan list,

I've spent the past several months using Stan in the context of causal inference models. Overall, Stan has been a life-saver, allowing me to fit vastly more complicated models than I would be able to quickly code by hand. So thanks again for making such a great product.

I'm writing because the upgrade from Stan 1.3 to Stan 2.x (both 2.0 and 2.1) proved to be a mixed blessing for the kinds of models I've been focusing on. On the one hand, effective sample sizes for the same models are marginally higher overall under Stan 2.1 than Stan 1.3. On the other, total computation time for the same number of iterations increased dramatically between v.1.3 and v.2.x, often by a factor of 3 or 4 in some of the more complicated models. While I haven't done a particularly scientific side-by-side (it's a pain to uninstall and re-install different versions), I'm confident that there's been a sizable decrease in ESS/sec when using Stan 2.x for these kinds of models.

I have two related questions.

First, is there a way to re-create the old Stan 1.3 default sampler behavior using the current version? I realize that this is not a sustainable question for future versions, but if someone's already thought about this, it would make it much easier for me to isolate the source of the slow-down and better understand the trade-offs in this case. (I've indeed tinkered with the sampler options but have yet to find a consistent speed-up)

Second, I've attached some much-simplified example Stan code, which is based on Hirano, Imbens, Rubin, and Zhou (2000; Biostatistics). Unfortunately, I can't post the original data set, but am happy to email it directly to folks if that's useful. My actual data set of interest involves many more bells-and-whistles, but the overall speed change seems to hold for this simpler model. I very much welcome suggestions on how to speed up this base model (note that I've been hard-coding the 'softmax' transform given the persistent memory leak issues there). 

Thanks, as always, for the help. Session info below.

Best,

Avi

> sessionInfo()

R version 3.0.2 (2013-09-25)

Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:

[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

[1] rstan_2.1.0   inline_0.3.13 Rcpp_0.10.5  

loaded via a namespace (and not attached):

[1] stats4_3.0.2 tools_3.0.2 

[Andrew Gelman]

Avi

I don't know if this would speed things up, but I think you could vectorize just about all your model statements (thus no longer needing all those little loops).

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

[Michael Betancourt]

Avi,

One suggestion on your model -- using softmax to compute pi will be more efficient. Or it should be if we didn't have a memory leak in the softmax code right now. Once its fixed you should definitely make the change.

The biggest difference in 2.1 that would effect performance is the increasing of the target acceptance probability to improve the robustness of our sampler, particularly for hierarchical models. You can recreate the 1.3 behavior by setting adapt delta=0.65, but be careful that you're not biasing the samples! 1.3 might be faster but also might be biased in the tails. For more information about the possible issues and some relevant diagnostics, see http://arxiv.org/abs/1312.0906.

[Avi Feller]

A,

Thanks for the response. I've messed around with a few vectorization approaches for the big 'for' loop. But these tend to be clunky due to the need to add (i.e., 'log_sum_exp') mixture terms. More importantly, doing so doesn't seem to increase the speed all that much, since 'for' loops are pretty damn fast on their own!

Of course, the 'softmax' part is still an extra loop (though N_type is typically very small; e.g., 2 or 3):

    for(k in 1:N_type){

      pi[k] <- exp(gamma_type[k] + delta[k] * x[n]);

    }

    pi <- pi/sum(pi);

Unfortunately, when I've tried to directly use 'softmax' here, I didn't see any speed-up. It's possible that this could change with the latest hot fix, but I'm not too hopeful.

Of course, I'll do whatever works. If there's an elegant vectorization that I've missed, I'd love to hear about it.

Thanks,

Avi

--
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/sxuTmhm6oCs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Michael Betancourt]

Oh, yeah! If you preprocess your data into each group beforehand then you can call the vectorized version of the distributions for significantly better performance.

[Michael Betancourt]

> More importantly, doing so doesn't seem to increase the speed all that much, since 'for' loops are pretty damn fast on their own!

Vectorization isn't trying to avoid the for loop, it has to do with our autodiff implementation. Needless to say, for the current version
of Stan creating a temporary array to hold each type of data and then using the vectorized calls will be MUCH faster. In your case
you can do even better since the groups are partitioned by data that can be done in preprocessing, either external to Stan or in
the transformed data block.

What you'll want to do here is compute the sums over each group first (using the vectorized bernoulli_logit_log and separate
arrays for each group/alpha_type) then sum over the mixture components in increment_log_prob.

(Actually, looking at the model again you may be boned. If the mixture is over each datum individually then there's nothing you can do).

>
> Of course, the 'softmax' part is still an extra loop (though N_type is typically very small; e.g., 2 or 3):
>
> for(k in 1:N_type){
> pi[k] <- exp(gamma_type[k] + delta[k] * x[n]);
> }
> pi <- pi/sum(pi);
>
> Unfortunately, when I've tried to directly use 'softmax' here, I didn't see any speed-up. It's possible that this could change with the latest hot fix, but I'm not too hopeful.

Hypothetically, softmax will make the code easier to read so there's really no reason not to use it. Of course you shouldn't use it
yet because of the bug. The sampling issue here is that the gamma_type and delta are non-identified because of the constraint
and the way you've written your model there's not way to place the more-appropriate simplex data type in here.

[Avi Feller]

Very helpful thoughts, as always. Responses inline:

On Mon, Dec 30, 2013 at 3:03 PM, Michael Betancourt <betan...@gmail.com> wrote:

> More importantly, doing so doesn't seem to increase the speed all that much, since 'for' loops are pretty damn fast on their own!

Vectorization isn't trying to avoid the for loop, it has to do with our autodiff implementation.  Needless to say, for the current version
of Stan creating a temporary array to hold each type of data and then using the vectorized calls will be MUCH faster.  In your case
you can do even better since the groups are partitioned by data that can be done in preprocessing, either external to Stan or in
the transformed data block.

What you'll want to do here is compute the sums over each group first (using the vectorized bernoulli_logit_log and separate
arrays for each group/alpha_type) then sum over the mixture components in increment_log_prob.

(Actually, looking at the model again you may be boned.  If the mixture is over each datum individually then there's nothing you can do).

Completely understood on autodiff vs. the 'for' loop. And I agree that it's straightforward to pre-process the data. But I think your last parenthetical is exactly on point---I'm not sure how to vectorize this since each person has a different mixture probability (i.e., 'pi' differs for each individual). Again, I've been staring at this same model for far too long, so it's likely that I'm missing something obvious.

For what it's worth, I basically implemented the mixture model in Section 11 of the manual. If there is a good way to vectorize this, it'd be great to add to that section of the manual. And, certainly, if there's a quick fix "under the hood" for this type of mixture modeling, I imagine other folks would benefit as well (though I have no idea what that fix might be).

>
> Of course, the 'softmax' part is still an extra loop (though N_type is typically very small; e.g., 2 or 3):
>
>     for(k in 1:N_type){
>       pi[k] <- exp(gamma_type[k] + delta[k] * x[n]);
>     }
>     pi <- pi/sum(pi);
>
> Unfortunately, when I've tried to directly use 'softmax' here, I didn't see any speed-up. It's possible that this could change with the latest hot fix, but I'm not too hopeful.

Hypothetically, softmax will make the code easier to read so there's really no reason not to use it.  Of course you shouldn't use it
yet because of the bug.  The sampling issue here is that the gamma_type and delta are non-identified because of the constraint
and the way you've written your model there's not way to place the more-appropriate simplex data type in here.

Agreed, and I'll wait until things are fixed.

I understand the non-identifiability issue here, but I don't follow your last part. Are you suggesting there's a better way to parameterize this?

[Avi Feller]

Also, changing delta to 0.65 does not appear to change the overall speed that much (I had tried this previously as well). But will investigate a bit more and report back.

[Michael Betancourt]

Hmm, that's odd. Do ESS and time change, but ESS/time is the same or is everything the same?
Also double check to see what average acceptance probabilities you get in your final chains to
ensure that the adaptation is narrowing down on the right values.

We refactored the samplers which induced a small slow down going from 1.3 to 2+
but it was less than 1.5 for all of the BUGS models and that was due to overhead and the slowdown
vanished for models that took any nontrivial amount of time to run.

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

[Bob Carpenter]

[changed subject for a slight digression]

There's also a basic problem with non-identifiability of
softmax-based parameterizations which requires one of
two solutions to identify well: (1) keep one of the
values fixed, as is typically done in a K-1-ary logistic
regression in the statistics world, or (2) provide strong
priors for whatever goes into determining the components.

In the loop under discussion:

>> for(k in 1:N_type){
>> pi[k] <- exp(gamma_type[k] + delta[k] * x[n]);
>> }
>> pi <- pi/sum(pi);

This means either setting one of the pi[k] = 0, or
putting strong priors on gamma_type[k] in particular to
avoid the translation non-identifiability (add z to all
of the pi and you get the same answer.

For solution 2, you can also simplify the whole thing and
write:

pi <- softmax(gamma_type + delta * x[n]);

- Bob

[Bob Carpenter]

I. There's nothing to be gained speedwise from doing this:

for (n in 1:N)
... set log_p_vec[n] <- ...val comp for n ...
}

increment_log_prob( sum(log_p_vec) );

instead of:

for (n in 1:N)
... increment_log_prob(...val comp for n...); ...

The increment_log_prob builds up its own internal
vector and applies a big sum at the end that includes all
of the contributions of Jacobians and sampling statements and
other increment_log_prob calls. In fact, it'll be faster
than writing out into an array because (a) no array or
allocation or copies, and (b) flatter auto-diff expression
tree because the val will be flattened into every other
component of the log_prob (such as the priors here).

II. What would speed things up is only computing
pi[k] / sum(pi) for the k that you need. You only use
one per loop and then you take the log, so that'd be:

for(k in 1:N_type)
pi[k] <- exp(gamma_type[k] + delta[k] * x[n]);

as is, and then instead of using

pi <- pi/sum(pi);

and then pulling out one of the pi and logging it, as in
log(pi[k]), use:

log(pi[k] / sum(pi))

which is the same as:

log(pi[k]) - log(sum(pi)).

You could keep going and instead of taking the
log of an exp for log(pi[k]), just use

for(k in 1:N_type)
pi[k] <- gamma_type[k] + delta[k] * x[n];

(FWIW, I don't like setting variables and the renormalizing; I
know we do it in our examples, but here, two variables pi and exp_pi
would be preferable --- I should update our doc.)

Then, for whichever component j you use, just do:

pi[j] - log(sum(exp(pi))

which is more efficiently and arithmetically stably coded as:

pi[i] - log_sum_exp(pi)

Normalizing a log-scale multinomial is the bread-and-butter of
the log_sum_exp operation.

- Bob







- Bob

On 12/30/13, 10:31 PM, Avi Feller wrote:
> Dear wonderful Stan list,
>
> I've spent the past several months using Stan in the context of causal inference models. Overall, Stan has been a
> life-saver, allowing me to fit vastly more complicated models than I would be able to quickly code by hand. So thanks
> again for making such a great product.
>
> I'm writing because the upgrade from Stan 1.3 to Stan 2.x (both 2.0 and 2.1) proved to be a mixed blessing for the kinds
> of models I've been focusing on. On the one hand, effective sample sizes for the same models are marginally higher
> overall under Stan 2.1 than Stan 1.3. On the other, total computation time for the same number of iterations increased
> dramatically between v.1.3 and v.2.x, often by a factor of 3 or 4 in some of the more complicated models. While I
> haven't done a particularly scientific side-by-side (it's a pain to uninstall and re-install different versions), I'm
> confident that there's been a sizable decrease in ESS/sec when using Stan 2.x for these kinds of models.
>
> I have two related questions.
>
> First, is there a way to re-create the old Stan 1.3 default sampler behavior using the current version? I realize that
> this is not a sustainable question for future versions, but if someone's already thought about this, it would make it
> much easier for me to isolate the source of the slow-down and better understand the trade-offs in this case. (I've
> indeed tinkered with the sampler options but have yet to find a consistent speed-up)
>
> Second, I've attached some much-simplified example Stan code, which is based on Hirano, Imbens, Rubin, and Zhou (2000;

> /Biostatistics/). Unfortunately, I can't post the original data set, but am happy to email it directly to folks if

> that's useful. My actual data set of interest involves many more bells-and-whistles, but the overall speed change seems
> to hold for this simpler model. I very much welcome suggestions on how to speed up this base model (note that I've been
> hard-coding the 'softmax' transform given the persistent memory leak issues there).
>
> Thanks, as always, for the help. Session info below.
>
> Best,
> Avi
>
>
>
> > sessionInfo()
> R version 3.0.2 (2013-09-25)
> Platform: x86_64-apple-darwin10.8.0 (64-bit)
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] rstan_2.1.0 inline_0.3.13 Rcpp_0.10.5
>
> loaded via a namespace (and not attached):
> [1] stats4_3.0.2 tools_3.0.2
>
>

[Avi Feller]

Fantastically useful comments. Responses below.

On Mon, Dec 30, 2013 at 4:35 PM, Bob Carpenter <ca...@alias-i.com> wrote:

I.  There's nothing to be gained speedwise from doing this:

  for (n in 1:N)
    ... set log_p_vec[n] <- ...val comp for n ...
  }

  increment_log_prob( sum(log_p_vec) );

instead of:

  for (n in 1:N)
    ... increment_log_prob(...val comp for n...); ...

The increment_log_prob builds up its own internal
vector and applies a big sum at the end that includes all
of the contributions of Jacobians and sampling statements and
other increment_log_prob calls.    In fact, it'll be faster
than writing out into an array because (a) no array or
allocation or copies, and (b) flatter auto-diff expression
tree because the val will be flattened into every other
component of the log_prob (such as the priors here).

Excellent advice. And an easy change to make. Thanks.

This makes sense. Do you have a sense of how the speed compares to the full 'softmax' model you reference in the other email chain? 

Speaking of which, thanks for that helpful summary of the two softmax approaches. I've been going back and forth between them for some time now and haven't noticed a difference in speed or stability either way. I now imagine that, once all the bugs are updated, the full vector form will be more efficient, right?

 

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 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/sxuTmhm6oCs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Michael Betancourt]

Really the probabilities need to be parameterized in the (n-1)-dimensional
latent space, but that's non-trivial for the model as defined here. Really
you'd need to define some product of simplices that yields another simplex.

[Bob Carpenter]

Why is that a problem? Just make gamma_type and delta
of dimension (N_type - 1) and use:

for(k in 1:(N_type - 1))

pi[k] <- exp(gamma_type[k] + delta[k] * x[n]);

pi[N_type] <- 0;

At least that's the usual multi-logit used in the GLM
world in stats. If you read the machine learning literature,
they usually use N dimensions instead of (N - 1) and don't
bother with identifiability.

- Bob

[Bob Carpenter]

Vectorization only makes things faster in Stan if it
reduces the size of the expression tree. It's a bit
subtle to tell where that is, but it's basically when
you reuse the same input for multiple outputs (as in
scalar * vector) or reuse multiple inputs for the same
output (as in the vectorized probability functions).
In contrast, exp(vector) is no faster than the loop, just
easier to write.

- Bob

On 12/31/13, 1:41 AM, Avi Feller wrote:
> Fantastically useful comments. Responses below.
>

> [1] en_US.UTF-8/en_US.UTF-8/en_US.__UTF-8/C/en_US.UTF-8/en_US.UTF-__8

>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] rstan_2.1.0 inline_0.3.13 Rcpp_0.10.5
>
> loaded via a namespace (and not attached):
> [1] stats4_3.0.2 tools_3.0.2
>
>
> --
> 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 <mailto:stan-users%2Bunsu...@googlegroups.com>.
>
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

>
>
> --
> 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/sxuTmhm6oCs/__unsubscribe
> <https://groups.google.com/d/topic/stan-users/sxuTmhm6oCs/unsubscribe>.
> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@__googlegroups.com
> <mailto:stan-users%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/__groups/opt_out <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.

[Andrew Gelman]

I'd prefer to parameterize in K dimensions, that just seems cleaner.

[Michael Betancourt]

> Why is that a problem? Just make gamma_type and delta
> of dimension (N_type - 1) and use:
>
> for(k in 1:(N_type - 1))
> pi[k] <- exp(gamma_type[k] + delta[k] * x[n]);
> pi[N_type] <- 0;
>
> At least that's the usual multi-logit used in the GLM
> world in stats.

But then pi[N_type] is always zero -- isn't that an odd modeling choice?

> If you read the machine learning literature,
> they usually use N dimensions instead of (N - 1) and don't
> bother with identifiability.

Sure, but when you only care about point prediction identifiability isn't so important!

[Bob Carpenter]

On 12/31/13, 8:34 PM, Michael Betancourt wrote:
>> Why is that a problem? Just make gamma_type and delta
>> of dimension (N_type - 1) and use:
>>
>> for(k in 1:(N_type - 1))
>> pi[k] <- exp(gamma_type[k] + delta[k] * x[n]);
>> pi[N_type] <- 0;
>>
>> At least that's the usual multi-logit used in the GLM
>> world in stats.
>
> But then pi[N_type] is always zero -- isn't that an odd modeling choice?

It's standard to identify the model. They even discuss it in
the Wikipedia entry under identifiability:

http://en.wikipedia.org/wiki/Multinomial_logistic_regression

saying:

Because only differences of vectors of regression coefficients are used,
adding an arbitrary constant to all coefficient vectors has no effect on the
model. This means that, just as in the log-linear model, only K-1 of the
coefficient vectors are identifiable, and the last one can be set to an arbitrary
value (e.g. 0).

The downside is that then the priors don't make a lot of sense
due to the asymmetry of the parameterization.

>> If you read the machine learning literature,
>> they usually use N dimensions instead of (N - 1) and don't
>> bother with identifiability.
>
> Sure, but when you only care about point prediction identifiability isn't so important!

Because only differences of vectors of regression coefficients are used, adding an arbitrary constant to all coefficient
vectors has no effect on the model. This means that, just as in the log-linear model, only K-1 of the coefficient
vectors are identifiable, and the last one can be set to an arbitrary value (e.g. 0).

I assume that's why they don't care about it. It
does make error checking problematic as well as convergence monitoring
if you look at parameters and have some randomization. Hastie et al.'s
machine learning book uses the (K-1)-parameterization, whereas Bishop's
uses the K-parameterization. This confused me so much a few years ago that
I had to work through all the math myself:

We should really be providing a (K-1)-parameterized softmax along with
a complete multi-logit model. The gradients are really simple and it
would be much much faster than using softmax the way we do now.

- Bob

[Bob Carpenter]

On 12/31/13, 6:51 PM, Andrew Gelman wrote:
> I'd prefer to parameterize in K dimensions, that just seems cleaner.

I agree that it's cleaner. But you have to have stronger priors
to avoid all the problems of (additively) non-identifiable models (or singular
models in Watanabe's terminology), such as effective sample size
and R-hat measurements. For example, your general Cauchy-prior recommendation
is a disaster in terms of sampling efficiency.

You can mitigate the convergence monitoring problem by creating generated quantities
by subtracting the K-th coefficient from all coefficients, but that's not so
easy in a GLM. And it doesn't help with adaptation.

I've played around with this a bit and the K-dimensional parameterization
requires MANY more samples to converge or much stronger priors.

- Bob

[Andrew Gelman]

I have some thoughts on this but maybe we could discuss side by side, it will be easier than email.

[Michael Betancourt]

Yeah, I think a parameterization in terms of (K - 1) deltas is the right approach,
both in terms of computational efficiency and defining interpretable priors.
It also has the right relationship to a simplex (K - 1 translation-invariant deltas
exponentiate to K - 1 scale-invariant probabilities) -- this is the intuition for which
I was looking earlier.

[andre.p...@googlemail.com]

Dear Stan list, 

>> for(k in 1:(N_type - 1)) 
>>    pi[k] <- exp(gamma_type[k] + delta[k] * x[n]); 
>> pi[N_type] <- 0;

I've been looking how they resolve this in mnrfit in matlab (which I'm familiar with):

http://www.mathworks.ch/ch/help/stats/multinomial-models-for-nominal-responses.html

When I'm not totally wrong that would 

pi[N_type] <- 1 - sum(pi[1:(N_type-1)]); 

Or maybe we could use the softmax where when we set gamma_type[k] + delta[k] * x[n]  == (equvialent for all) 1

for k = N_type?

Andre

[Bob Carpenter]

Yes, the mathworks link is for the usual K-1 parameter
approach to K-ary logistic regression. We'll probably
be adding that to Stan, soon.

- 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/d/optout.