More model-specific speed decreases from 1.3 to 2.1

[Jeff Goldsmith]

Hi all,

First, I want to express my appreciation for all the work that has gone into developing Stan -- it is remarkably user-friendly and flexible, and has a great community.

This post is similar in nature to that of Avi Feller on 12/30/13, but I thought it wouldn't hurt to start a new thread. I've been working on problems in function-on-scalar regression -- basically, subject-level responses are smooth functions observed on a dense grid, predictors are scalars, and coefficients are smooth functions. Coefficient functions are estimated using a B-spline basis, so estimation and inference focuses on the B-spline coefficients. I'm attaching an R file that simulates such a dataset, poses a Bayesian estimation model using RStan, and estimates the spline coefficients. This is a pretty simplified version of the model in Goldsmith & Kitago "Assessing Systematic Effects of Stroke on Motor Control Using Hierarchical Function-on-Scalar Regression", in case the attachedl code isn't particularly clear.

The issue that I'm running into is related to computation time. Direct comparisons between Stan 1.3 and 2.1 are difficult (I'm running on two different machines to keep from having to repeatedly reinstall different versions), but my best guess is that for the same effective sample sizes 2.1 is taking between 10 and 20 times longer than 1.3. Moreover, the difference seems to get worse as the number of iterations increases. For more complex models, the difference can be much more substantial.

Any thoughts about how to speed up the code in general would be appreciated, but in particular I'm wondering if there is some explanation for the surprising increase in computation time for the same effective sample sizes.

Thanks for any insights the group can provide, and happy new year!

Best,

Jeff

[Bob Carpenter]

I think this may be due to the fact that there was a bug
in the multivariate normal code that we backed out to much
slower and simpler implementations pending fixes. Marcus or
Daniel should know for sure.

Fixing them hasinvolved endless tests and design and process
discussions on our end (see the dev list and issue tracker), but
we can see the light at the end of the tunnel and it should
be faster again in 2.2.0.

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

[Marcus Brubaker]

The LDLT changes were never in 1.3.  This is most likely due to the changes in the sampler code, e.g., adaptation, acceptance rate, etc.

Looking at the model I don't think the LDLT changes would make a big difference here as the covariance matrices are pretty small.

Finally, Jeff, your model would probably be faster if you did something like the following:

transformed data {

...

cholesky_factor_cov{Kt] chol_inv_CovMat;

chol_inv_CovMat <- cholesky_decompose(inverse(CovMat));

}

model {

...

// Prior for spline coefficients for beta

for (pcur in 1:p) {

(beta[pcur])' ~ multi_normal_chol(mu_beta, beta_sig[pcur] * chol_inv_CovMat);

}

...

}

and drop the transformed parameters.

Also, drop all of the ~uniform calls, they are redundant and just slow things down.  If you want to restrict the range of these variables, change their definitions.  Or, even better, put proper priors on them as I suspect Andrew will chime in here saying any minute now.

Cheers,

Marcus

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.

[Jeff Goldsmith]

Thanks for the quick response, guys. 

With respect to using the Cholesky formulation, I'm happy to play around with it and see if it speeds things up. As a heads-up, there seems to be a typo in the compiled Stan model highlighted in red below:

• 152:         try { stan::math::check_choelsky_factor(function__,chol_inv_CovMat,"chol_inv_CovMat"); } catch (std::domain_error& e) { throw std::domain_error(std::string("Invalid value of chol_inv_CovMat: ") + std::string(e.what())); };
  

This looks like a pretty easy fix, but in case you hadn't run into the error in the past I thought I'd point it out. Also, thanks for the tip regarding removing the uniform priors -- I know the specification is redundant, but liked to keep it there to remind myself exactly what is going on. I'll take it out and see if there's a noticeable difference in speed.

With respect to where the change in speed is coming from (multivariate normal code or changes in the sampler code), I guess I'll hold out hope that bug fixes in 2.2.0 can at least help speed things up. As I mentioned, this is a pretty simplified version of the model to illustrate the issues -- in more realistic scenarios the computation time has gone from a few minutes to at least many hours.

Thanks again for your input, and for all the effort done by the Stan development team on this really impressive project. Best,

Jeff

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.

[Bob Carpenter]

There was a bit of slowdown in the underlying code in the
move from 1.3 to a more general framework in 2.0. That
should be a pretty constant slowdown of less than a factor
of 1.5.

Then there was a big slowdown (we measured a factor of 1.5 to 2 times
on the BUGS models) in the move in 2.1 of a target acceptance rate of
0.8 from 0.65. The upside is a reduction in bias in mean estimates
and better variance estimates for hierarchical models (and other models
with varying posterior scales).

Given that it's not the multivariate implementation, I'd suggest
going back to a target acceptance rate of 0.65. The command-line option
is "adapt delta=0.65" to do that. If you're using RStan, I don't know
what it is.

***BUT***, as Michael keeps stressing, you need to make sure this doesn't bias your
results. Do the 1.3 and 2.1 version get the same posterior means, std deviation
and 90% intervals? If so, you should be OK with a reduced target acceptance rate.

- Bob

On 1/4/14, 1:22 AM, Jeff Goldsmith wrote:
> Thanks for the quick response, guys.
>
> With respect to using the Cholesky formulation, I'm happy to play around with it and see if it speeds things up. As a
> heads-up, there seems to be a typo in the compiled Stan model highlighted in red below:
>

> * 152: try { stan::math::check_choelsky_factor(function__,chol_inv_CovMat,"chol_inv_CovMat"); } catch

> stan-users+...@__googlegroups.com <javascript:>.
> 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 <javascript:>.
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

[Jeff Goldsmith]

It doesn't seem to be the multivariate implementation, just as you and Marcus suspected. I tried a slightly altered model using:

(beta[pcur])' ~ normal(mu_beta, beta_sig[pcur]);

in place of the multivariate normal. For both 1.3 and 2.1 there was a bit of a slowdown compared to the multivariate approach, but the sampling in 1.3 was still around 10x faster than in 2.1.

I've tried altering the target acceptance rate using the argument 

                 control = list(adapt_delta = .65)

in R. I'm not sure if this is the correct syntax, but it is what I've been using throughout and shows up in the original script I posted. If someone with more experience has a suggestion here, it would be very much appreciated. 

Thanks for the word of warning about potential biases with lower acceptance rates. In the example posted, the results in 1.3 and 2.1 are indistinguishable; in more systematic simulations using 1.3, there hasn't been any appreciable bias across many simulated datasets. The only noticeable difference is in the computation time.

Best,

Jeff

[Bob Carpenter]

Looking more closely at your model, one issue
is that you need to constrain the parameters to their
support.

You have:

parameters {
matrix[p,Kt] beta; // matrix of fixed effect spline coefficients
real<lower=0> y_sig; // regression variance
vector<lower=0>[p] beta_sig; // tuning parameter
}

and then in the model:

y_sig ~ uniform(0,50);

// Prior for variance components controlling smoothness in beta
for(pcur in 1:p) {
beta_sig[pcur] ~ uniform(0,50);
}

You need to add upper=50 to the declaration for y_sig and beta_sig, as in:

real<lower=0,upper=50> y_sig;

Do you expect the posterior to range over this whole range of 0 to 50? If
not, tightening the interval should help sampling speed.

Your call to stan looks like this:

PFSR.fit = stan(model_name="PFSR", model_code = pfsrcode,
data=dat, iter = 10000, warmup = 4000,
control = list(adapt_delta = .65),
chains = 1, verbose = FALSE)

I'd recommend more chains of fewer iterations, which can be run in parallel.

How long does it take the model to converge? And how many effective samples
do you need?

The issue we're going to have with debugging is that the person who
undertands the samplers the best doesn't run R. I'm running the models
now to see if I can see what's going on.

- Bob



On 1/3/14, 11:41 PM, Jeff Goldsmith wrote:

[Jeff Goldsmith]

Thanks for the suggestions to improve speed in the model by constraining the range of support for some of the parameters and to run more chains in parallel. I think that will be broadly useful across a lot of models.

Since I was running comparisons on two systems and getting very different results, I found a third machine to try things out on. This was a clean install of Stan 2.1 (rather than upgrading from 1.3 to 2.1), and most of the issues have gone away -- computation times are now different by the modest 50% or so that was expected all along. My guess is that in upgrading from 1.3 to 2.1 I did something wrong enough to slow things down immensely but not so wrong that it completely broke everything. My apologies for not checking this more carefully before posting, and thank you for all the input along the way.

Best,

Jeff

[Bob Carpenter]

Ack! The problem is in RStan 2.1, not in the core
Stan samplers.

For this example, RStan is about 10 times slower than
CmdStan for version 2.1.

Time for 1000 warmup, 1000 sampling iterations with seed=1234
and target acceptance of 0.8 (new default in 2.1):

CmdStan: 40 seconds ( 26s warmup, 14s sampling)
RStan: 404 seconds (230s warmup, 175s sampling)

We'll look into what's going on and hopefully roll out a patch
for RStan ASAP (assuming it's not something that Jeff Goldsmith
and/or I did wrong in our config, in which case, we'll roll
out some better instructions).

- Bob

P.S. I'll follow up on stan-dev list with the model and
data in separate files.

> > stan-users+...@__googlegroups.com <http://googlegroups.com> <javascript:>.

> > For more options, visit https://groups.google.com/__groups/opt_out

> <https://groups.google.com/__groups/opt_out> <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 <http://googlegroups.com> <javascript:>.

> > For more options, visit https://groups.google.com/__groups/opt_out

> <https://groups.google.com/__groups/opt_out> <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 <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

[Jeff Goldsmith]

For me, it's useful to know that the mistake I made somewhere in upgrading from 1.3 to 2.1 is easy to make :-). And in general, isolating the issue to RStan rather than the core samplers is good news!

As I mentioned in another response, on a third system I installed 2.1 from scratch and didn't run into any unexpected speed increase. If I can spot a difference between the "working" RStan 2.1 configuration and the "slow" RStan 2.1 configuration I'll let you know.

Jeff

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

[Bob Carpenter]

Thanks, Jeff. You helped me realize I also had RStan 2.1
installed improperly. I carefully RTFM this time. No idea
what I did wrong the first time.

In fact, RStan 2.1 with g++ -O3 is a bit faster than CmdStan
with clang++ -O3:

RStan properly installed: 33s
CmdStan: 40s

Maybe we need to do better buffering for CmdStan --- our output
uses std::endl, which flushes every line --- we should try replacing
it with "\n" --- we probably don't want the output to be platform-dependent
anyway.

- Bob

[Jeff Goldsmith]

Well, there you go. I just went through the same process of reinstalling, discovering that everything works now, and not being sure exactly what I did wrong the first time around. It may have been in making sure the ~/.R/Makevars file had the correct form, if that's helpful.

Thanks again and best,

Jeff

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