Vectorized function slower than its non vectorized version.

[Tammo Nikolas]

Hello, I hope this is not obvious. I am very new to STAN.

I have a user-defined function (convolution) and coded it twice, expecting the second, vectorized version to be faster. 

It turns out that it takes almost twice as long as the first version 1. Am I doing something wrong here?

I'm using pystan 2.4.0.1

Version 1 (not vectorized)

functions{
 vector convolution(int N, real tStep, vector aif, vector irf){ 
    vector<lower=0>[N] ctr; 
       for(i in 1:N){

           ctr[i] <- 0.0;

           for(j in 1:i){

               ctr[i] <- ctr[i] + aif[i - j + 1] * irf[j];

            }

         }

       return ctr*tStep;

 }

Version 2 (vectorized)

functions{
 vector convolution(int N, real tStep, vector aif, vector irf){
   vector<lower=0>[N] ctr;
   vector<lower=0>[N] summand;
   for(i in 1:N){
      summand <- rep_vector(0, N);
      for(j in 1:i){
         summand[j] <- aif[i - j + 1] * irf[j];
         }
      ctr[i] <- sum(summand);
      }
   return ctr*tStep;
 }

Cheers

Tammo

[Marco Inacio]

I suspect that the fact that you are using a summand vector of size N at each iteration is the guilt part.

Does this work:

functions{
 vector convolution(int N, real tStep, vector aif, vector irf){
   vector<lower=0>[N] ctr;
   for(i in 1:N){

      vector<lower=0>[i] summand;

      for(j in 1:i){
         summand[j] <- aif[i - j + 1] * irf[j];
      }
      ctr[i] <- sum(summand);
   }
   return ctr*tStep;
 }

?

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

[Tammo Nikolas]

Thanks! It works and it is faster than before, but it's still slower (by about 10% for the few tests I ran) than the non-vectorized Version 1.

[Marco Inacio]

You're right! And for me it's more than just 10%:

(here using RStan)

> fit1_n100 <- sampling(cmodel1, data=list(N=100), seed=10, refresh=1e5) #vectorized sum

SAMPLING FOR MODEL 'model1' NOW (CHAIN 1).

Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 0.493747 seconds (Warm-up)
#                0.402926 seconds (Sampling)
#                0.896673 seconds (Total)

> fit2_n100 <- sampling(cmodel2, data=list(N=100), seed=10, refresh=1e5) #non-vectorized sum

SAMPLING FOR MODEL 'model2' NOW (CHAIN 1).

Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 0.32559 seconds (Warm-up)
#                0.263426 seconds (Sampling)
#                0.589016 seconds (Total)


This also seems to be the case for bigger N

> fit1_n100 <- sampling(cmodel1, data=list(N=1000), seed=10)

SAMPLING FOR MODEL 'model1' NOW (CHAIN 1).

Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 47.1636 seconds (Warm-up)
#                44.21 seconds (Sampling)
#                91.3736 seconds (Total)

> fit2_n100 <- sampling(cmodel2, data=list(N=1000), seed=10)

SAMPLING FOR MODEL 'model2' NOW (CHAIN 1).

Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)
#  Elapsed Time: 31.0938 seconds (Warm-up)
#                29.0972 seconds (Sampling)
#                60.191 seconds (Total)






Source code:

require(rstan)

model1 <- "

functions {
 vector convolution(int N, real tStep, vector aif, vector irf) {
   vector<lower=0>[N] ctr;
   for(i in 1:N) {
     vector<lower=0>[i] summand;
     for(j in 1:i) {
       summand[j] <- aif[i - j + 1] * irf[j];
     }
     ctr[i] <- sum(summand);
   }
   return ctr*tStep;
 }
}

data {
  int<lower=0> N;
}
transformed data {
  vector [N] aif_;
  vector [N] irf_;
 
  aif_[1] <- 1;
  irf_[1] <- 1;
  
  for (i in 2:N) {
    aif_[i] <- 1 + aif_[i-1]*0.1 + irf_[i-1]*0.16;
    irf_[i] <- 1 + aif_[i-1]*0.06 + irf_[i-1]*0.13;
  }
}
parameters {
  vector [N] beta;
}
transformed parameters {
  vector [N] beta_t;
  beta_t <- convolution(N, 1.1, aif_, irf_);
}
model {
  beta ~ normal(0, 1);
}"

model2 <- "

functions {
 vector convolution(int N, real tStep, vector aif, vector irf) {
    vector<lower=0>[N] ctr;
       for(i in 1:N) {
           ctr[i] <- 0.0;
           for(j in 1:i) {
               ctr[i] <- ctr[i] + aif[i - j + 1] * irf[j];
            }
       }
       return ctr*tStep;
 }
}

data {
  int<lower=0> N;
}
transformed data {
  vector [N] aif_;
  vector [N] irf_;
 
  aif_[1] <- 1;
  irf_[1] <- 1;
  
  for (i in 2:N) {
    aif_[i] <- 1 + aif_[i-1]*0.1 + irf_[i-1]*0.16;
    irf_[i] <- 1 + aif_[i-1]*0.06 + irf_[i-1]*0.13;
  }
}
parameters {
  vector [N] beta;
}
transformed parameters {
  vector [N] beta_t;
  beta_t <- convolution(N, 1.1, aif_, irf_);
}
model {
  beta ~ normal(0, 1);
}"

cmodel1 <- stan_model(model_code=model1)

cmodel2 <- stan_model(model_code=model2)

fit1_n100 <- sampling(cmodel1, data=list(N=100), seed=10, refresh=1e5)

fit2_n100 <- sampling(cmodel2, data=list(N=100), seed=10, refresh=1e5)

On 14-08-08 04:35 PM, Tammo Nikolas wrote:

Thanks! It works and it is faster than before, but it's still slower (by about 10% for the few tests I ran) than the non-vectorized Version 1.

[Bob Carpenter]

On Aug 8, 2014, at 4:58 PM, Marco Inacio <marcoig...@gmail.com> wrote:



> You're right! And for me it's more than just 10%:

Thanks for investigating. Two things surprise me about the model.

1. I'd expect the sum() version to be faster for larger N.

2. I'd expect a compilation failure because we don't allow
constraints on local variable declarations because there's no
good point to validate them. I'll investigate and open an issue,
but we shouldn't have constraints on local variables.

If aif were reversed (which could be done once or the data
structures could be changed), then you cold do this
with a slicing operation in one loop, reducing the body of
the function to:

vector[N] ctr;
for (n in 1:N)
ctr[n] <- sum(tail(aif_rev,n) .* head(irf,n));
return tStep * ctr;

But that'd probably be even less efficient as you're now
generating temporary lists for head and tail. Maybe we should
think more about expression templates, which would mitigate
the copying cost. I'm sure that's what's killing the sum-based
approach --- the allocations of the vectors, that is; but here
there's no way to reduce that cost because each summand is
a different size.

- Bob

P.S. I've tried to restrict the term "vectorized" to when
a function can take a vector for an argument as well as
a scalar, and that's not what's going on here. But it's
used all over in different ways in other packages and by
other people, so it's probably a hopeless battle to restrict
its technical meaning.

> vector convolution(int N, real tStep, vector aif, vector irf) {
> vector<lower=0>[N] ctr;

> for (i in 1:N) {
> vector<lower=0>[i] summand;
> for (j in 1:i) {

[Marco Inacio]

On 14-08-10 05:14 PM, Bob Carpenter wrote:
> On Aug 8, 2014, at 4:58 PM, Marco Inacio <marcoig...@gmail.com> wrote:
>
>
>
>> You're right! And for me it's more than just 10%:
> Thanks for investigating. Two things surprise me about the model.
>
> 1. I'd expect the sum() version to be faster for larger N.
>
> 2. I'd expect a compilation failure because we don't allow
> constraints on local variable declarations because there's no
> good point to validate them. I'll investigate and open an issue,
> but we shouldn't have constraints on local variables.

I've also spotted that, but since I'm not familiar with the newly
introduced user defined functions, I thought somehow it was ok there.
I've also tested the models without that constraint on summand variable,
and the results are similar.

Also, as could be expected, if we use the function convolution inside
the generated quantities block, the results also favor the
"non-vectorized" form.

[Ben Goodrich]

On Sunday, August 10, 2014 4:15:03 PM UTC-4, Bob Carpenter wrote:

P.S.  I've tried to restrict the term "vectorized" to when
a function can take a vector for an argument as well as
a scalar, and that's not what's going on here.  But it's
used all over in different ways in other packages and by
other people, so it's probably a hopeless battle to restrict
its technical meaning.

Can we start calling such functions "vector friendly"? And some vector friendly functions are also "derivative accumulating"? And a small number of for loops in Stan are "batched" by the compiler / processor?

Ben

[Bob Carpenter]

On Aug 11, 2014, at 2:06 AM, Ben Goodrich <goodri...@gmail.com> wrote:

> On Sunday, August 10, 2014 4:15:03 PM UTC-4, Bob Carpenter wrote:
> P.S. I've tried to restrict the term "vectorized" to when
> a function can take a vector for an argument as well as
> a scalar, and that's not what's going on here. But it's
> used all over in different ways in other packages and by
> other people, so it's probably a hopeless battle to restrict
> its technical meaning.
>
> Can we start calling such functions "vector friendly"? And some vector friendly functions are also "derivative accumulating"?

OK. And there are really two different ways to be "vector friendly".

1. the way log(vector) works, is componentwise. That is, it just
repeats log on each element.

2. the way normal_log(vector,...) works, which is to add the result
for each element together.

Eigen calls the former elementwise operations and the latter reductions.
Some reductions are built-in and not just accumulations over elementwise
definitions, such as sum() and maxCoeff().

> And a small number of for loops in Stan are "batched" by the compiler / processor?

I don't know what this means.

- Bob