Vectorization problem

[Tom Wallis]

Hi there, 

I am fitting hierarchical glms to binomial data. I tried vectorising one of my models to speed it up. The model compiles and runs just fine, but my "vectorised" version ("model2.txt") runs much slower than the original ("model1.txt"). It took 18 seconds for 100 samples of 4 chains for model 1, and 96 seconds for model 2. I've also attached some test data that will allow you to run the models. 

I assume there's a problem with how I've modified the model to try to use vectors and matices where possible. Where am I going wrong?

Thanks for your help. 

[Ben Goodrich]

Your code looks okay to me. I tried to run it, but it seems that gamma is omitted from testData.txt

goodrich@CYBERPOWERPC:/tmp/SU/TW$ time ./model1 --refresh=0 --init=0 --seed=12345 --data=testData.txt

Exception: variable does not exist; processing stage=data initialization; variable name=gamma; base type=double
Diagnostic information:
Dynamic exception type: std::runtime_error
std::exception::what: variable does not exist; processing stage=data initialization; variable name=gamma; base type=double


real    0m0.099s
user    0m0.083s
sys     0m0.003s
goodrich@CYBERPOWERPC:/tmp/SU/TW$ grep gammma testData.txt
goodrich@CYBERPOWERPC:/tmp/SU/TW$ grep "<-" testData.txt
N <- 7000
y <-
D <- 4
x <-
ss <-
S <- 5

Perhaps the col() function is slower than I thought (it does make a copy unfortunately). If you change the declaration to

vector[D] beta[S];

and reference beta[d] in the loop, it might go faster.

Ben

[Tom Wallis]

Woops, right you are Ben, I forgot to include gamma. I've attached a new version of the data file with gamma included.

Could this be caused by a memory problem with my use of sequential fitting rather than something to do with the stan model code? I compared the times by fitting model 1, then fitting model 2. I also fit a 3rd model that was vectorised, as model 2, but altered to use the "matt trick" (file attached). This 3rd model took nearly 12 hours to do 100 samples in 4 chains. All models were sampled in the same script, in the same R session.

I've tried two things to clarify what's going on:

1. I ran the three models again, this time removing the fit objects from (R) memory using rm(fitObject). The time taken was 46, 251 and 663 seconds. So still an accumulating performance problem, but not nearly as bad as before.

2. I changed the order of sampling the models to 3, 2, 1 instead of 1, 2, 3. The time taken for models 1, 2, and 3 under this regime was 52, 150, and 1700 seconds. These were all run with the same seed.

Should I be clearing the fit objects from R in between fits when comparing timing? Is there some way to clear stan's memory separately? In any case, I'm not seeing a speed improvement with vectorisation. Any ideas why not?

[Bob Carpenter]

On 12/13/12 11:12 AM, Tom Wallis wrote:
> Woops, right you are Ben, I forgot to include gamma. I've attached a new version of the data file with gamma included.

lambda ~ beta(2,50); # note that beta here refers to the beta CDF, not to the beta weights of the regression.

I didn't understand the comment here.

If you want to use the beta cdf you will eventually
be able to say:

lambda <- beta_cdf(outcome,2,50);

But as of now, we don't have the beta_cdf implemented.
Michael Betancourt's working on implementing all the
CDFs starting tomorrow (we needed a real math person
to work out all the crazy non-analytic derivatives
w.r.t. the distribution parameters).

> Could this be caused by a memory problem with my use of sequential fitting rather than something to do with the stan
> model code? I compared the times by fitting model 1, then fitting model 2. I also fit a 3rd model that was vectorised,
> as model 2, but altered to use the "matt trick" (file attached). This 3rd model took nearly 12 hours to do 100 samples
> in 4 chains. All models were sampled in the same script, in the same R session.

For performance comparisons, you should always try to
keep conditions equivalent. So that would mean rerunning
from a new R session without anything else going on in
the background (e.g., not watching a video or something).

> I've tried two things to clarify what's going on:
>
> 1. I ran the three models again, this time removing the fit objects from (R) memory using rm(fitObject). The time taken
> was 46, 251 and 663 seconds. So still an accumulating performance problem, but not nearly as bad as before.

The difference was 10 minutes vs. 12 hours? I wouldn't
expect a factory of 60+ unless you somehow ran out of
memory and started swapping. But then I'm not an R
expert, so maybe Jiqiang or Ben will have a better idea of
what might be going on.

> 2. I changed the order of sampling the models to 3, 2, 1 instead of 1, 2, 3. The time taken for models 1, 2, and 3 under
> this regime was 52, 150, and 1700 seconds. These were all run with the same seed.

> Should I be clearing the fit objects from R in between fits when comparing timing? Is there some way to clear stan's
> memory separately? In any case, I'm not seeing a speed improvement with vectorisation. Any ideas why not?

Always clear everything between speed tests. In fact,
I'd suggest restarting R from scratch.

Stan uses R's environment for memory, so you just need
to get rid of the fit object and anything you've attached.

I can't imagine a situation in which the vectorized
versions will be slower -- the vectorized code should
be both faster and consume less memory. As Ben said in a
previous note, you probably don't want to use col() for
optimal performance, but it should still be faster than
not vectorizing.

- Bob

[Ben Goodrich]

On Thursday, December 13, 2012 11:12:06 AM UTC-5, Tom Wallis wrote:

Should I be clearing the fit objects from R in between fits when comparing timing? Is there some way to clear stan's memory separately? In any case, I'm not seeing a speed improvement with vectorisation. Any ideas why not?

Something is definitely strange here. I am attaching the results for --test_grad for each of the three models, with inits of zero. Models 2 and 3 have the same log posterior (which makes sense) but very different gradients (possible but odd). Model 1 has a different log posterior for reasons that my eyes have not been able to catch yet (even after commenting out the uniform prior on sigma).

Ben

[Tom Wallis]

Bob, sorry about the unclear comment; that should be beta PDF, not CDF. That is, I'm trying to place a beta(2,50) prior on lambda -- giving lambda a prior mean of 0.038. 

[Tom Wallis]

Just spotted one problem in Model 3. In the model block I'm giving beta a prior based on mu and sigma, whereas I think this line should be deleted and the prior is the unit normal on e_beta; beta will get contributions from mu and sigma by the transformed parameter block. 

This can't explain the speed issues though.

[Jiqiang Guo]

I did not see the results of test gradients.  But model 1 (at least the file I am looking at) is problematic since sigma, on which a uniform prior is specified, is defined without upper constraint. 

--

Jiqiang 

--
 
 

[Ben Goodrich]

On Thursday, December 13, 2012 4:28:56 PM UTC-5, Tom Wallis wrote:

Just spotted one problem in Model 3. In the model block I'm giving beta a prior based on mu and sigma, whereas I think this line should be deleted and the prior is the unit normal on e_beta; beta will get contributions from mu and sigma by the transformed parameter block. 

This can't explain the speed issues though.

Yes, that was a problem. When using the Matt trick, always put the prior on the thing in the parameters {} block, which implies a prior on the thing in the transformed parameters {} block, but if you explicitly put a prior on the thing in the transformed parameters {} block, you will get the wrong posterior.

I have made some changes so that model1 is equivalent to model2 and the latter is slightly faster. The gain in speed should depend mostly on D, which is only 4 in this case. Originally, in model1 x was declared as

row_vector[D] x[N]; # vector of predictor values for each trial.

and accessed by

p <- gamma + (1 - gamma - lambda[ss[n]]) * Phi(x[n] * beta[ss[n]]);

whereas in model 2 (and 3) x was declared as

matrix[N,D] x; # design matrix for linear model.

and accessed by

p <- gamma + (1 - gamma - lambda[ss[n]]) * Phi(dot_product(row(x,n),row(beta,ss[n])));

Apparently the first way is faster, but I have harmonized on the second way until we sort about the bigger issue: model3 is still not equivalent to the other two for some reason that I haven't been able to spot in the .stan code or in the C++ code.

Ben

[Bob Carpenter]

On 12/13/12 4:25 PM, Tom Wallis wrote:
> Bob, sorry about the unclear comment; that should be beta PDF, not CDF. That is, I'm trying to place a beta(2,50) prior
> on lambda -- giving lambda a prior mean of 0.038.

Great example of the issues I talked about in:

http://lingpipe-blog.com/2009/10/15/the-futility-of-commenting-code/

This example only reminded me of the issue. To cite
one of my own doing that's on topic, as of Stan v1.0.3,
for the multinomial pmf stan::prob::multinomial_log(),
the doc says the outcomes are 1:K but the implementation
says 0:(K-1).

So the code will show you what Stan will do whereas the
doc is just plain wrong. As a computer scientist, my
knee-jerk reaction is to index from 0.

I'm inclined to change the code to match the doc in this case.

- Bob

> --
>
>

[Bob Carpenter]

On 12/13/12 6:10 PM, Ben Goodrich wrote:

> I have made some changes so that model1 is equivalent to model2 and the latter is slightly faster. The gain in speed
> should depend mostly on D, which is only 4 in this case. Originally, in model1 x was declared as
>
> |
>
> row_vector[D] x[N]; # vector of predictor values for each trial.
>
> and accessed by
>

> ... x[n] ...

>

> whereas in model 2 (and 3) x was declared as
>
> matrix[N,D] x; # design matrix for linear model.
>
> and accessed by
>

> ... row(x,n) ...

>
> Apparently the first way is faster,

Yes. The reason is that for arrays, Stan uses
nested std::vector objects in C++. That means
pulling out a "row" of a vector<vector<double> >
can return a reference to an existing object
rather than allocating and copying.

Although I'm the only one who finds it natural,
x[n] returns the same result as row(x,n).

Now if we didn't mangle Eigen's expression templates,
we could deal with row() and column() without copying,
too, and column() would be efficient because the data's
local due to the column-major order we use for Eigen
matrices.

- Bob

[Bob Carpenter]

The whole point of reparameterizing is to change the geometry of the
posterior in the parameters in the "parameters" block. The goal
of the transforms is to make the posterior geometry closer to
Euclidean and hence easier to traverse stably with HMC. In particular,
we want to try to remove position (i.e., parameter) dependent curvature
and scaling as much as possible.

Michael B and I are about to add auto-diff for Hessians and Hessian-vector
products, so we'll soon be able to inspect the curvature implications of
the transformations.

- Bob