[Rcpp-devel] How to increase the coding efficiency

[Romain Francois]

Le 04/12/12 21:27, Honglang Wang a écrit :
> Since the "ker" function defined in R is very simple and vectorized, I
> thought it was fine to call this function in C++. And I tried to code
> this function inside C++, and I found it became ever slower: from
> 47.310s to 49.536s.

Hard to say what is going on. Perhaps you are allocating too many things 
with :

     arma::vec tp = t-t0;
     arma::uvec q1 = arma::find(arma::abs(tp)<h);
     arma::uvec q2 = arma::find(arma::abs(tp)>=h);
     arma::vec kk = t;
     kk.elem(q1) = ((1-arma::pow(tp.elem(q1)/h,2))/h)*0.75;
     kk.elem(q2).zeros();


You are allocating tp, q1, q2 and kk.

Perhaps, there are ways not to do that.

arma brings you nice syntax, but sometimes this has a cost.


Sometimes, well written R code performs well and the speedup you gain 
with just translating into C++ might not be huge.


> // [[Rcpp::depends(RcppArmadillo)]]
>
> #include <RcppArmadillo.h>
>
> using namespace Rcpp;
>
> // [[Rcpp::export]]
> List betahat(double t0, NumericMatrix Xr, NumericMatrix yr,
> NumericVector tr, double h, int m) {
>    int n = Xr.nrow(), p = Xr.ncol();
>    arma::mat X(Xr.begin(), n, p, false);
>    arma::mat y(yr.begin(), n, 1, false);
>    arma::colvec t(tr.begin(), tr.size(), false);
>    arma::mat T = X;
>    T.each_col() %= (t-t0)/h;
>    arma::mat D = arma::join_rows(X,T);
>    arma::vec tp = t-t0;
>    arma::uvec q1 = arma::find(arma::abs(tp)<h);
>    arma::uvec q2 = arma::find(arma::abs(tp)>=h);
>    arma::vec kk = t;
>    kk.elem(q1) = ((1-arma::pow(tp.elem(q1)/h,2))/h)*0.75;
>    kk.elem(q2).zeros();
>    arma::mat W = (arma::diagmat(kk))/m;
>    arma::mat Inv = arma::trans(D)*W*D;
>    arma::vec betahat = arma::inv(Inv)*arma::trans(D)*W*y;
>    arma::colvec betahat0(betahat.begin(),betahat.size()/2,false);
>    return List::create(Named("betahat") = betahat0);
> }
>
>
> Best wishes!
> Honglang Wang
> Office C402 Wells Hall
> Department of Statistics and Probability
> Michigan State University
> 1579 I Spartan Village, East Lansing, MI 48823
> wangho16 at msu.edu <mailto:wangho16 at msu.edu>
>
>
>
> On Tue, Dec 4, 2012 at 2:19 PM, Romain Francois
> <romain at r-enthusiasts.com <mailto:romain at r-enthusiasts.com>> wrote:
>
>     You are calling back to R using Function. This is expensive.
>
>     What is ker ? Can you implement it in C++. This is a wild guess, but
>     that is I think where the bottleneck is.
>
>     Romain
>
>     Le 04/12/12 20:14, Honglang Wang a écrit :
>
>         Dear All,
>         I have tried out the first example by using RcppArmadillo, but I
>         am not
>         sure whether the code is efficient or not. And I did the
>         comparison of
>         the computation time.
>
>         1) R code using for loop in R: 87.22s
>         2) R code using apply: 77.86s
>         3) RcppArmadillo by using for loop in C++: 53.102s
>         4) RcppArmadillo together with apply in R: 47.310s
>
>         It is kind of not so big increase. I am wondering whether I used an
>         inefficient way for the C++ coding:
>
>
>         // [[Rcpp::depends(RcppArmadillo)__]]
>
>         #include <RcppArmadillo.h>
>
>         using namespace Rcpp;
>
>         // [[Rcpp::export]]
>         List betahat(Function ker, double t0, NumericMatrix Xr,
>         NumericMatrix
>         yr, NumericVector tr, double h, int m) {
>             int n = Xr.nrow(), p = Xr.ncol();
>             arma::mat X(Xr.begin(), n, p, false);
>             arma::mat y(yr.begin(), n, 1, false);
>             arma::colvec t(tr.begin(), tr.size(), false);
>             arma::mat T = X;
>             T.each_col() %= (t-t0)/h;
>             arma::mat D = arma::join_rows(X,T);
>             arma::vec kk =as<arma::vec>(ker(tr-t0,h));
>             arma::mat W = (arma::diagmat(kk))/m;
>             arma::mat Inv = arma::trans(D)*W*D;
>             arma::vec betahat = arma::inv(Inv)*arma::trans(D)*__W*y;
>             arma::colvec betahat0(betahat.begin(),__betahat.size()/2,false);
>             return List::create(Named("betahat") = betahat0);
>         }
>
>         Anyone could help me with how to increase the efficiency of the
>         coding?
>         Thanks.
>
>
>
>         Best wishes!
>         Honglang Wang
>         Office C402 Wells Hall
>         Department of Statistics and Probability
>         Michigan State University
>         1579 I Spartan Village, East Lansing, MI 48823
>         wangho16 at msu.edu <mailto:wangho16 at msu.edu>
>         <mailto:wangho16 at msu.edu <mailto:wangho16 at msu.edu>>
>
>
>
>     --
>     Romain Francois
>     Professional R Enthusiast
>     +33(0) 6 28 91 30 30 <tel:%2B33%280%29%206%2028%2091%2030%2030>
>
>     R Graph Gallery: http://gallery.r-enthusiasts.__com
>     <http://gallery.r-enthusiasts.com>
>     `- http://bit.ly/SweN1Z : SuperStorm Sandy
>
>     blog: http://romainfrancois.blog.__free.fr
>     <http://romainfrancois.blog.free.fr>
>     |- http://bit.ly/RE6sYH : OOP with Rcpp modules
>     `- http://bit.ly/Thw7IK : Rcpp modules more flexible
>
>     _________________________________________________
>     Rcpp-devel mailing list
>     Rcpp-devel at lists.r-forge.r-__project.org
>     <mailto:Rcpp-devel at lists.r-forge.r-project.org>
>     https://lists.r-forge.r-__project.org/cgi-bin/mailman/__listinfo/rcpp-devel
>     <https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel>
>
>


-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30

R Graph Gallery: http://gallery.r-enthusiasts.com
`- http://bit.ly/SweN1Z : SuperStorm Sandy

blog:            http://romainfrancois.blog.free.fr
|- http://bit.ly/RE6sYH : OOP with Rcpp modules
`- http://bit.ly/Thw7IK : Rcpp modules more flexible