Hi, <br><br>I am trying to speed up a vectorized R function without much luck. Perhaps it is not possible, but if it is could someone suggest how to improve the C++/Rcpp code below? The speed is approximately equal on a 64-bit version of R, but the R code is much faster when using 32-bit R (on windows). Here's an example. <br>
<br><br> # R function<br>nllR <- function(y, psi, p, nd) {<br> cp <- (p^y) * ((1 - p)^(1 - y))<br> cp <- exp(rowSums(log(cp))) # faster than: apply(cp, 1, prod)<br> loglik <- log(cp * psi + nd * (1 - psi))<br>
-sum(loglik)<br> }<br><br><br># Rcpp function<br>nllRcpp <- cxxfunction(signature(yR="matrix", psiR="numeric", pR="matrix", <br> ndR="integer"), <br> '<br> NumericMatrix y(yR);<br>
NumericVector psi(psiR);<br> NumericMatrix p(pR);<br> IntegerVector nd(ndR);<br> <br> NumericVector yV(y);<br> int R = y.nrow();<br> int J = y.ncol();<br> NumericMatrix cp(R, J);<br> NumericVector ll(R);<br>
NumericVector cpsum(R);<br> <br> for(int i=0; i<R; i++) {<br> for(int j=0; j<J; j++) {<br> cp(i,j) = pow(p(i,j), y(i,j)) * pow(1-p(i,j), 1-y(i,j));<br> }<br> NumericVector cpi = cp(i, _);<br>
cpsum[i] = exp(sum(log(cpi)));<br> ll[i] = log(cpsum[i] * psi[i] + nd[i] * (1-psi[i]));<br> }<br> double nll = -1*sum(ll); <br> return wrap(nll);<br> ', plugin="Rcpp")<br><br>
<br> <br># Test<br>M <- 10000<br>J <- 100<br>psi <- runif(M)<br>Z <- rbinom(M, 1, psi)<br>p <- matrix(runif(M*J), M, J)<br>y <- matrix(rbinom(M*J, 1, p*Z), M, J)<br>nd <- ifelse(rowSums(y, na.rm=TRUE) == 0, 1, 0)<br>
<br><br>all.equal(nllR(y, psi, p, nd), nllRcpp(y, psi, p, nd))<br><br># 64-bit<br>> benchmark(nllR(y, psi, p, nd), nllRcpp(y, psi, p, nd))<br> test replications elapsed relative user.self<br>1 nllR(y, psi, p, nd) 100 10.58 1.046489 9.35<br>
2 nllRcpp(y, psi, p, nd) 100 10.11 1.000000 9.81<br> sys.self user.child sys.child<br>1 1.14 NA NA<br>2 0.26 NA NA<br><br><br># 32-bit<br>> benchmark(nllR(y, psi, p, nd), nllRcpp(y, psi, p, nd))<br>
test replications elapsed relative user.self<br>1 nllR(y, psi, p, nd) 100 16.75 1.000000 15.54<br>2 nllRcpp(y, psi, p, nd) 100 29.97 1.789254 29.81<br> sys.self user.child sys.child<br>
1 1.08 NA NA<br>2 0.01 NA NA<br><br><br><br>> sessionInfo()<br>R version 2.12.0 (2010-10-15)<br>Platform: x86_64-pc-mingw32/x64 (64-bit)<br><br>...<br><br>attached base packages:<br>
[1] stats graphics grDevices utils datasets methods <br>[7] base <br><br>other attached packages:<br>[1] rbenchmark_0.3 inline_0.3.8 RcppArmadillo_0.2.10<br>[4] Rcpp_0.9.0 <br><br>loaded via a namespace (and not attached):<br>
[1] tools_2.12.0<br>> <br><br><br><br>Thanks,<br>Richard<br><br><br><br><br><br><br>