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