Hi Honglang, I recommend looking into using the solve() function in Armadillo, instead of inv(). Using solve() will be considerably faster.
More details here: http://arma.sourceforge.net/docs.html#solve On Thursday, December 6, 2012, Honglang Wang <wanghonglang2...@gmail.com> wrote: > The following is a full example although I don't know whether it's minimal or not: > > library(Rcpp) > library(RcppArmadillo) > sourceCpp("betahat_mod.cpp") > > #The following is data generation. > n=200 > m=20 > p=2 > t=runif(m*n,min=0, max=1) > X1=rnorm(m*n,0,1) > X1=as.matrix(1+2*exp(t)+X1) > X2=rnorm(m*n,0,1) > X2=as.matrix(3-4*t^2+X2) > X=cbind(X1,X2) > beta1=0.5*sin(t) > beta2=beta1 > rho=0.2 > sig2=1 > eps=unlist(lapply(1:n,function(x){as.matrix(arima.sim(list(ar=rho),m,sd=sqrt(sig2*(1-rho^2))))})) > y=as.matrix(X1*beta1+X2*beta2+eps) > > #A simple kernel function. > ker=function(x,h) > { > ans=x > lo=(abs(x)<h) > ans[lo]=(3/4)*(1-(x[lo]/h)^2)/h > ans[!lo]=0 > return(ans) > } > > > h=0.3 > > #assess the time for evaluating bethat of 100 t's: > system.time((betahatt=t(apply(as.matrix(t[1:100]),1,function(x) betahat(ker,x,X,y,t,h,m)$betahat)))) > > And the .cpp file is the following: > > // [[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); > } > > Best wishes! > > Honglang Wang > > Office C402 Wells Hall > Department of Statistics and Probability > Michigan State University > 1579 I Spartan Village, East Lansing, MI 48823 > wangh...@msu.edu > > > On Wed, Dec 5, 2012 at 4:07 AM, Christian Gunning <x...@unm.edu> wrote: >> >> Can you post a minimal full example? >> -Christian >> >> On Tue, Dec 4, 2012 at 8:39 PM, Honglang Wang >> <wanghonglang2...@gmail.com> wrote: >> > Yes, the main issue for my coding is the allocation of memory. And I have >> > fixed one of the biggest memory allocation issue: 4000 by 4000 diagonal >> > matrix. And since I am not familiar with Rcpp and RcppArmadillo, I have no >> > idea how to reuse the memory. I hope I can have some materials to learn >> > this. Thanks. >> > >> > >> >> >> >> What exactly do these timings show? A single call you your function? >> >> How many calls? >> >> >> > Here I called my function for 100 times. >> > >> >> >> >> Building on Romain's point: -- a portion of your function's runtime is >> >> in memory allocation >> >> (and you have a lot of allocations here). >> >> If you're calling your function thousands or millions of times, then >> >> it might pay to closely >> >> examine your memory allocation strategies and figure out what's >> >> temporary, for example. >> >> It looks like you're already using copy_aux_mem = false in a number >> >> of places, but you're >> >> allocating a lot of objects -- of approx what size? >> >> >> >> For example, wouldn't this work just as well with one less allocation? >> >> arma::vec kk = t; >> >> arma::uvec q1 = arma::find(arma::abs(tp)<h); >> >> kk.elem(q1) = ((1-arma::pow(tp.elem(q1)/h,2))/h)*0.75; >> >> // done with q1. let's reuse it. >> >> q1 = arma::find(arma::abs(tp)>=h); >> >> // was q2 >> >> kk.elem(q1).zeros(); >> >> >> >> You could potentially allocate memory for temporary working space in >> >> R, grab it with copy_aux_mem = false, write your temp results there, >> >> and reuse these objects in subsequent function calls. It doesn't make >> >> sense to go to this trouble, though, if your core algorithm consumes >> >> the bulk of runtime. >> >> >> >> Have you looked on the armadillo notes r.e. inv? Matrix inversion has >> >> O(>n^2). You may be aided by pencil-and-paper math here. >> >> http://arma.sourceforge.net/docs.html#inv >> >> >> > Here my matrix for inverse is only 4 by 4, so I think it's ok. >> > >> >> >> >> best, >> >> Christian >> >> >> >> > 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: >> >> >> >> >> >> >> >> -- >> >> A man, a plan, a cat, a ham, a yak, a yam, a hat, a canal – Panama! >> > >> > >> >> >> >> -- >> A man, a plan, a cat, a ham, a yak, a yam, a hat, a canal – Panama! > >
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
