Rcpp+inlineを使ってみた

まえから気になっていたRcpp+inlineを使ってみた

単純な繰り返しループの集計計算がかなり時間がかかるので、Rcppでどれくらい改善されるか挑戦

まずは通常のコード

非常に単純で以下のxとBeta,Alphaの線形結合のquantileを求めるというもの。

view plainprint?

1. x <- rnorm(100)  
2. Beta <- Alpha <- rnorm(1000000)  
3. res1 <- sapply(1:length(x), function(i)   
4.                quantile(exp(Beta*x[i]+Alpha)/(1+exp(Beta*x[i]+Alpha)),probs=0.025))  

x <- rnorm(100)
Beta <- Alpha <- rnorm(1000000)
res1 <- sapply(1:length(x), function(i) 
               quantile(exp(Beta*x[i]+Alpha)/(1+exp(Beta*x[i]+Alpha)),probs=0.025))

次にRcpp+inline

view plainprint?

1. library(Rcpp)  
2. library(inline)  
3. src <- "NumericVector xx(x);  
4.   NumericVector beta(b);  
5.   NumericVector alpha(a);  
6.   Function quantile(\"quantile\");  
7.   int xsize=xx.size(), bsize=beta.size();  
8.   NumericVector mu(bsize);  
9.   NumericVector res(xsize);  
10.   for(unsigned long int i=0; i < xsize; i++){  
11.    for(unsigned long int j=0; j < bsize; j++){  
12.     double dummy=exp(beta[j]*xx[i]+alpha[j]);  
13.     mu[j]=0.0;  
14.     mu[j]=dummy/(1.0+dummy);  
15.    }  
16.    res[i]=0.0;  
17.    res[i]=as<double>(quantile(mu, _[\"probs\"]=0.025));  
18.   }  
19.   return res;"  
20.   
21. fun.R <- cxxfunction(signature(x="numeric",a="numeric",b="numeric"),  
22.        src,plugin="Rcpp")  
23.   
24. </double>  

library(Rcpp)
library(inline)
src <- "NumericVector xx(x);
  NumericVector beta(b);
  NumericVector alpha(a);
  Function quantile(\"quantile\");
  int xsize=xx.size(), bsize=beta.size();
  NumericVector mu(bsize);
  NumericVector res(xsize);
  for(unsigned long int i=0; i < xsize; i++){
   for(unsigned long int j=0; j < bsize; j++){
    double dummy=exp(beta[j]*xx[i]+alpha[j]);
    mu[j]=0.0;
    mu[j]=dummy/(1.0+dummy);
   }
   res[i]=0.0;
   res[i]=as(quantile(mu, _[\"probs\"]=0.025));
  }
  return res;"

fun.R <- cxxfunction(signature(x="numeric",a="numeric",b="numeric"),
       src,plugin="Rcpp")

結果

view plainprint?

1. x <- rnorm(100)  
2. Beta <- Alpha <- rnorm(1000000)  
3. system.time(res1 <- sapply(1:length(x), function(i)   
4.     quantile(exp(Beta*x[i]+Alpha)/(1+exp(Beta*x[i]+Alpha)),probs=0.025)))  
5.    ユーザ   システム       経過    
6.     19.162      1.008     20.169   
7. system.time(res2 <- fun.R(x,Beta,Alpha))  
8.    ユーザ   システム       経過    
9.      8.936      0.932      9.870   

x <- rnorm(100)
Beta <- Alpha <- rnorm(1000000)
system.time(res1 <- sapply(1:length(x), function(i) 
    quantile(exp(Beta*x[i]+Alpha)/(1+exp(Beta*x[i]+Alpha)),probs=0.025)))
   ユーザ   システム       経過  
    19.162      1.008     20.169 
system.time(res2 <- fun.R(x,Beta,Alpha))
   ユーザ   システム       経過  
     8.936      0.932      9.870 

計算結果も同じだった

view plainprint?

1. summary(res1-res2)  
2.    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.   
3.       0       0       0       0       0       0   

summary(res1-res2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       0       0       0       0 

RcppでRの関数を使うのではなく、GNU Scientific Library (GSL)のquantile(gsl_stats_quantile_from_sorted_data)も試してみたけどこちらは通常のRの1.5倍ぐらいの時間がかかった....
何がボトルネックなんだろう...