Nice vectorized multi_normal speedup

[Richard McElreath]

With 2.4 out, I ran a quick test of how much speedup I could get by changing my old non-vectorized multi_normal sampling to the new vectorized form. I get a 40% time savings, without even trying hard. This is much better than I expected.

Timings with vectorized multi_normal:

#  Elapsed Time: 96.8564 seconds (Warm-up)

#                87.8308 seconds (Sampling)

#                184.687 seconds (Total)

Timings without vectorized multi_normal:

#  Elapsed Time: 168.544 seconds (Warm-up)

#                147.83 seconds (Sampling)

#                316.375 seconds (Total)

These timings are for a Gaussian mixed model with three random effects, 5000 cases, and 500 groups. The only difference between these two models was changing this line:

vary ~ multi_normal( rep_vector(0,3) , SIGMA );

to this line:

for ( j in 1:N_id ) vary[j] ~ multi_normal( rep_vector(0,3) , SIGMA );

Really nice stuff.

[Bob Carpenter]

The savings are proportionally bigger with higher dimensional
SIGMA, too.

You get the same kind of speedup vectorizing plain old normal, too.

The main reason is that the checks for well-formedness of Sigma
only need to be done once and that the autodiff is more efficient
because the expression tree getting built up has fewer nodes, and
hence fewer virtual function calls.

We still have some headroom to speed these things up --- the only
thing we're missing is programmer time!

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

[Marco Inacio]

And since you are using multi_normal with covariance parametrization,
there should be some speedup should be coming from calculating LDLT
factor of Sigma only once (and also, if Sigma is not data, its determinant).

[Bob Carpenter]

Marco's right --- though I believe the factorization and
well-formedness check are just one call to a solver. And
that's the big cost here.

- Bob

[Ben Goodrich]

On Monday, July 28, 2014 4:02:12 PM UTC-4, Richard McElreath wrote:

These timings are for a Gaussian mixed model with three random effects, 5000 cases, and 500 groups. The only difference between these two models was changing this line:

vary ~ multi_normal( rep_vector(0,3) , SIGMA );

to this line:

for ( j in 1:N_id ) vary[j] ~ multi_normal( rep_vector(0,3) , SIGMA );

Really nice stuff.

It is nice that multi_normal is vectorized now. But it is also important to keep in mind opportunities to avoid using it (can't quite tell which might be applicable in your situation, but as a general principle the following are true). If the LHS is a parameter, you can change it to a transformed parameter where

vary <- L * z;

such that L * L' = SIGMA, where L is a Cholesky factor, and z is a matrix whose elements are distributed iid standard normal a priori i.e. include

to_vector(z) ~ normal(0,1);

in the model block.

Or if SIGMA is unknown with an inverse-Wishart prior (which I am skeptical of), you can integrate it out according to Wikipedia


and obtain a prior distribution for the unknown X as a function of the degrees of freedom (nu) and scatter matrix (Psi). In the usual case where nu and Psi are known, this is proportional to a

increment_log_prob(-0.5 * (nu + n) * log_determinant(Psi + tcrossprod(X)));

So, I really don't use multi_normal() much.

Ben 

[Richard McElreath]

There were some questions about the full model code for the test I reported, so I'm attaching the complete test script here.

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

Hi, 

According to a recent question:

https://groups.google.com/forum/#!category-topic/stan-users/general/d2yWqNMcE8s

I have had a look at this older post:

 

Or if SIGMA is unknown with an inverse-Wishart prior (which I am skeptical of), you can integrate it out according to Wikipedia


and obtain a prior distribution for the unknown X as a function of the degrees of freedom (nu) and scatter matrix (Psi). In the usual case where nu and Psi are known, this is proportional to a

increment_log_prob(-0.5 * (nu + n) * log_determinant(Psi + tcrossprod(X)));

 

I don't understand, why isn't it:

increment_log_prob(-0.5 * (nu + n) * log_determinant(Psi + tcrossprod(X)) + 0.5 * nu * log_determinant(psi));

Why we can ignore the top of the fraction?

Andre

[Ben Goodrich]

Psi --- the scatter matrix in the inverse Wishart prior --- was presumed to be known (along with nu), in which case the

0.5 * nu * log_determinant(Psi)

term is a constant and does not affect the proposals for the parameters or which get accepted.

Ben