# Reduce\_sum performance

**URL:** <https://discourse.mc-stan.org/t/reduce-sum-performance/15330>\
**Category:** Modeling\
**Created:** [May 22, 2020, 3:43pm UTC](https://discourse.mc-stan.org/t/reduce-sum-performance/15330 "2020-05-22T15:43:08Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Louis](https://avatars.discourse-cdn.com/v4/letter/l/ecb155/32.png) [@Louis](https://discourse.mc-stan.org/u/Louis)\
**Post date:** [May 22, 2020, 3:43pm UTC](https://discourse.mc-stan.org/t/reduce-sum-performance/15330/1 "2020-05-22T15:43:08Z")

</div>

I have a large spatiotemporal hierarchical model that takes hours to run, so I am now considering the use of reduce\_sum to try to speed things up a bit. However, in my first attempt, the model that uses reduce\_sum is actually slower. To try to understand things better, I’ve tested reduce\_sum on a simple example. Basically, I generate 50,000 samples from a GEV distribution and then try to infer the GEV parameters from the simulated data. In this example, the model that uses reduce\_sum is three times faster (1 min vs 3 min). However, note that I have about 144 ‘CPUs’ available for parallel threads: 4 sockets x 18 cores per socket x 2 threads per core = 144 CPUs. I set grainsize = 1, and with this Stan seems to use 60 CPUs when on reduce\_sum. A 3x speed increase is significant but I was expecting more given that Stan was using 60 CPUs. I guess, is this expected? I just want to make sure I am doing things right before I go through the effort of introducing reduce\_sum into my larger model.

Stan code below:

```stan
functions{
    real gev_lpdf(real[] y,real mu,real sigma,real xi,int n){
        vector[n] z;
            
        for(i in 1:n){
            z[i] = (1 + xi * (y[i] - mu) / sigma)^(-1/xi);
        }
        return -sum(log(sigma) - (1 + xi) * log(z) + z); 
    }
        
    real partial_sum(real[] ys,int is,int ie,real mu,real sigma,real xi){
        int ns = ie - is + 1;
        return gev_lpdf(ys | mu,sigma,xi,ns);   
    }
}  
data{
    int<lower=0> n;
    real y[n];
}
transformed data{
    int grainsize = 1;
}
parameters{
    real mu;
    real<lower=0> sigma;
    real xi;
}
model{
    mu ~ normal(0,2);
    sigma ~ normal(0,1);
    xi ~ normal(0,1);
        
    //y ~ gev(mu,sigma,xi,n); // no threading
    target += reduce_sum(partial_sum,y,grainsize,mu,sigma,xi);
}

```

---

<div class="post-metadata">

**Author:** ![wds15](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wds15/32/908_2.png) [@wds15](https://discourse.mc-stan.org/u/wds15)\
**Post date:** [May 22, 2020, 6:55pm UTC](https://discourse.mc-stan.org/t/reduce-sum-performance/15330/2 "2020-05-22T18:55:07Z")

</div>

What a machine! First, hyperthreading usually never buys you anything since Stan needs a lot of cache. The fact that things are spread on many sockets is also not so good, I think.

You will probably get 3x speedup with much less CPUs. You should try to replace the large for loop by a single vectorized expression, if you can.

The grainsize of 1 is too small for you problem. I would try 500 or so as a start.

Also, have a read of Amdahls law on wikipedia. You may spend 1000 cores and only get 20x speedup in some problems.

---

<div class="post-metadata">

**Author:** ![Louis](https://avatars.discourse-cdn.com/v4/letter/l/ecb155/32.png) [@Louis](https://discourse.mc-stan.org/u/Louis)\
**Post date:** [May 22, 2020, 8:21pm UTC](https://discourse.mc-stan.org/t/reduce-sum-performance/15330/3 "2020-05-22T20:21:16Z")

</div>

Thanks for your reply, wds15. I see now that choosing a sensible grainsize is important. I tried several values for the grainsize but, in this case, I can’t do better than ~1 min. And thanks also for the reference to Amdahl’s law, which makes a lot of sense.

There is no vectorized version of pow() in Stan, so I am afraid that the loop is unavoidable.

---

<div class="post-metadata">

**Author:** ![nhuurre](https://avatars.discourse-cdn.com/v4/letter/n/ad7895/32.png) [@nhuurre](https://discourse.mc-stan.org/u/nhuurre)\
**Post date:** [May 22, 2020, 8:40pm UTC](https://discourse.mc-stan.org/t/reduce-sum-performance/15330/4 "2020-05-22T20:40:08Z")

</div>

There’s a vectorized `log`, that loop is equivalent to

```stan
real gev_lpdf(real[] y, real mu, real sigma, real xi, int n){
    vector[n] log_z = log1p(xi * (to_vector(y) - mu) / sigma) * (-1/xi);
    return -n * log(sigma) + (1 + xi) * sum(log_z) - sum(exp(log_z));
}

```

---

<div class="post-metadata">

**Author:** ![wds15](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wds15/32/908_2.png) [@wds15](https://discourse.mc-stan.org/u/wds15)\
**Post date:** [May 22, 2020, 9:04pm UTC](https://discourse.mc-stan.org/t/reduce-sum-performance/15330/5 "2020-05-22T21:04:26Z")

</div>

1m is already very fast…you would need a larger problem to see larger speed ups. Getting multiple h or days down to something to less than or just a few h is what I think makes sense.

The suggestion from @nhuurre sounds very good.

---

<div class="post-metadata">

**Author:** ![Louis](https://avatars.discourse-cdn.com/v4/letter/l/ecb155/32.png) [@Louis](https://discourse.mc-stan.org/u/Louis)\
**Post date:** [May 22, 2020, 9:14pm UTC](https://discourse.mc-stan.org/t/reduce-sum-performance/15330/6 "2020-05-22T21:14:55Z")

</div>

Thanks nhuurre for the suggestion. That’s interesting because I already tried the trick exp(n\*log(x)) in the bigger version of my model, as I discussed in this [thread](https://discourse.mc-stan.org/t/vectorization-of-exp-and-log/13683), and this seemed to make the code slower than looping. Surprisingly, when I try the trick in the simple model above the model runs in only 12 s! It’s amazing the large effect small changes can have, and such effect can be different depending on the model.