[Rcpp-devel] How to increase the coding efficiency

[Honglang Wang]

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.

// [[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



On Tue, Dec 4, 2012 at 2:19 PM, Romain Francois <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>
>>
>
>
> --
> Romain Francois
> Professional R Enthusiast
> +33(0) 6 28 91 30 30
>
> 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<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>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20121204/b90d6d46/attachment.html>