# Multivariate normal too slow

**URL:** <https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564>\
**Category:** Modeling\
**Tags:** performance\
**Created:** [November 21, 2017, 4:47pm UTC](https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564 "2017-11-21T16:47:18Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![annptr](https://avatars.discourse-cdn.com/v4/letter/a/a88e57/32.png) [@annptr](https://discourse.mc-stan.org/u/annptr)\
**Post date:** [November 21, 2017, 4:47pm UTC](https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564/1 "2017-11-21T16:47:18Z")

</div>

Hi there,  
I’m very new to Stan and I have been trying a model with multivariate outcomes.  
Quick details about the data: I have a set of subjects and groups and for each subject-group pair there is a measure which is 6-dimensional. Because of the way these measures are taken, I assume diagonal covariance matrix.  
It takes more than 10 hours to run 500 iterations (250 warmup 250 sampling) with 2 chains on 2 cores.  
Below is my code and any helpful look into that to speed up would be really great. Thanks!

```
data {
int<lower=1> K; //dimension size (6)
int<lower=1> J; // num subjects (183)
int<lower=1> S; // num groups (65)  
int<lower=0> N; // num observations (1000)
vector[J] x[N]; // for subjects indicators
vector[S] u[N]; // for group indicators
vector[K] y[N]; // outcome measure (numeric)
}

parameters {
matrix[K, J] beta; //subject coefficients
matrix[K, S] eta; // group coefficients
vector<lower=0>[K] sigma; // std for dimesnions
}

model {
vector[K] mu[N];
matrix[K, K] Sigma;
Sigma = diag_matrix(square(sigma)); // assuming diagonal covariance matrix

for (n in 1:N)
mu[n] = beta * x[n] + eta * u[n];
y ~ multi_normal(mu , Sigma);
}
```

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [November 21, 2017, 6:04pm UTC](https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564/2 "2017-11-21T18:04:35Z")

</div>

> [@annptr](#):
>
> Sigma = diag\_matrix(square(sigma)); // assuming diagonal covariance matrix

Never use `multi_normal` with a diagonal covariance matrix, use the univariate normal distribution on the stacked observations

---

<div class="post-metadata">

**Author:** ![annptr](https://avatars.discourse-cdn.com/v4/letter/a/a88e57/32.png) [@annptr](https://discourse.mc-stan.org/u/annptr)\
**Post date:** [November 21, 2017, 6:21pm UTC](https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564/3 "2017-11-21T18:21:52Z")

</div>

Thanks Ben!

I tried that version too (code below) and it took about 2.5 hours for 500 iterations with warning about increasing max\_treedepth.

```
data {
int<lower=1> J; // num subjects (183)
int<lower=1> S; // num groups (65)  
int<lower=0> N; // num observations (6000)
matrix[N,J] x; // for subjects indicators
matrix[N,S] u; // for group indicators
real y[N]; // outcome measure (numeric)
}

parameters {
vector[J] beta;
vector[S] eta;
real<lower=0> sigma;
}

model {
y ~ normal(x* beta + u * eta, sigma);
}

```

Is there any settings that I need to do in RStan to speed this up?

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [November 21, 2017, 6:50pm UTC](https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564/4 "2017-11-21T18:50:16Z")

</div>

You probably need to use the QR reparameterization.

---

<div class="post-metadata">

**Author:** ![annptr](https://avatars.discourse-cdn.com/v4/letter/a/a88e57/32.png) [@annptr](https://discourse.mc-stan.org/u/annptr)\
**Post date:** [November 21, 2017, 8:49pm UTC](https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564/5 "2017-11-21T20:49:14Z")

</div>

Thanks so much Ben! QR reparameterization reduced the time to 16 mins from 2.5 hours.

One question, is QR reparameterization still possible when I model the multivariate outcomes with covariance parameter (similar to my initial code in the question)? Instead of a matrix for the predictors, I’ll have a set of vectors (vector[J] x[N];)

Any pointers will be very useful.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [November 21, 2017, 9:23pm UTC](https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564/6 "2017-11-21T21:23:26Z")

</div>

You are going to need a matrix to be doing the QR decomposition. I think in your case, it would be several QR decompositions.

---

<div class="post-metadata">

**Author:** ![annptr](https://avatars.discourse-cdn.com/v4/letter/a/a88e57/32.png) [@annptr](https://discourse.mc-stan.org/u/annptr)\
**Post date:** [November 28, 2017, 3:41pm UTC](https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564/7 "2017-11-28T15:41:33Z")

</div>

Following up on this question. I tried a complicated model with a larger dataset as follows:

```
data{
int<lower=1> J; // num subjects (5,000)
int<lower=1> S; // num groups (90)
int<lower=0> N; // num observations (25,000)
int<lower=1> K; // num dimensions (6)
matrix[N,J] x; // for subjects indicators
matrix[N,S] u; // for group indicators
vector[K] y[N]; // outcome measure (numeric)
int<lower=1,upper=90> group[N]; // group for subject
}

transformed data {
  // QR reparameterization of predictors
  matrix[N, J] Qx = qr_Q(x)[, 1:J] * N;
  matrix[J, J] R = qr_R(x)[1:J,] / N;
  matrix[J, J] R_inv = inverse(R);

  matrix[N, S] Qu = qr_Q(u)[, 1:S] * N;
  matrix[S, S] R2 = qr_R(u)[1:S,] / N;
  matrix[S, S] R2_inv = inverse(R2);
}

parameters {
  matrix[K,J] beta_t; // for subject
  matrix[K,S] eta_t; // for group

  cholesky_factor_corr[K] L_Omega;
  vector<lower=0>[K] L_sigma;
  vector<lower=0,upper=1>[S] alpha; // random intercept for each subject-group
}

transformed parameters {
  matrix[K,J] beta = beta_t * R_inv ;
  matrix[K,S] eta = eta_t * R2_inv;
}

model {
vector[K] mu[N];
matrix[K, K] L_Sigma;

alpha ~ beta(2,2);
L_Sigma = diag_pre_multiply(L_sigma, L_Omega);

to_vector(beta) ~ normal(0, 5);
to_vector(eta) ~ normal(0, 5);
L_Omega ~ lkj_corr_cholesky(4);
L_sigma ~ cauchy(0, 2.5);

for (n in 1:N){
for (k in 1:K){
mu[n][k] = alpha[group[n]] * dot_product(beta_t[k], Qx[n]) + dot_product(eta_t[k], Qu[n]); // to avoid transpose
}
}
y ~ multi_normal_cholesky(mu, L_Sigma);
}

```

Based on the manual, I tried few speed up trick such as removing transpose and using do\_product.  
Initial estimation for two chains in two cores is as follows:

```
Gradient evaluation took 18.4891 seconds
1000 transitions using 10 leapfrog steps per transition would take 184891 seconds.
Adjust your expectations accordingly!

Gradient evaluation took 20.4548 seconds
1000 transitions using 10 leapfrog steps per transition would take 204548 seconds.
Adjust your expectations accordingly!

```

However, this code is running for more than four days now and it has only finished the warmpup state (250/500 iterations) and started sampling (251/500).

Any suggestions for further improving the speed would be helpful!

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 2, 2017, 1:07am UTC](https://discourse.mc-stan.org/t/multivariate-normal-too-slow/2564/8 "2017-12-02T01:07:43Z")

</div>

> [@annptr](#):
>
> Any suggestions for further improving the speed would be helpful!

The next thing to do is use matrix algebra to define `mu`. Somethign like

```no-highlight
mu[n] = alpha(group[n]) + beta_t[k] * Qx + eta_t[k] * Qu;

```

Or the other way around to vectorize on `K` if that’s bigger.

Where do your estimates fall? Those hard interval constraints can be bad news, as I just blogged about:

[http://andrewgelman.com/2017/11/28/computational-statistical-issues-uniform-interval-priors/](http://andrewgelman.com/2017/11/28/computational-statistical-issues-uniform-interval-priors/)

So the best thing to do at this point might be adding tighter priors. For instance, those Cauchy can also be problematic if there’s not enough data to determine the scale of `L_sigma`.

Thanks for sharing the detailed example and reporting back.