# [Rcpp-devel] speed fix

Richard Chandler richard.chandlers at gmail.com
Mon Dec 27 13:57:41 CET 2010

```Hi,

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.

# R function
nllR <- function(y, psi, p, nd) {
cp <- (p^y) * ((1 - p)^(1 - y))
cp <- exp(rowSums(log(cp))) # faster than: apply(cp, 1, prod)
loglik <- log(cp * psi + nd * (1 - psi))
-sum(loglik)
}

# Rcpp function
nllRcpp <- cxxfunction(signature(yR="matrix", psiR="numeric", pR="matrix",
ndR="integer"),
'
NumericMatrix y(yR);
NumericVector psi(psiR);
NumericMatrix p(pR);
IntegerVector nd(ndR);

NumericVector yV(y);
int R = y.nrow();
int J = y.ncol();
NumericMatrix cp(R, J);
NumericVector ll(R);
NumericVector cpsum(R);

for(int i=0; i<R; i++) {
for(int j=0; j<J; j++) {
cp(i,j) = pow(p(i,j), y(i,j)) * pow(1-p(i,j), 1-y(i,j));
}
NumericVector cpi = cp(i, _);
cpsum[i] = exp(sum(log(cpi)));
ll[i] = log(cpsum[i] * psi[i] + nd[i] * (1-psi[i]));
}
double nll = -1*sum(ll);
return wrap(nll);
', plugin="Rcpp")

# Test
M <- 10000
J <- 100
psi <- runif(M)
Z <- rbinom(M, 1, psi)
p <- matrix(runif(M*J), M, J)
y <- matrix(rbinom(M*J, 1, p*Z), M, J)
nd <- ifelse(rowSums(y, na.rm=TRUE) == 0, 1, 0)

all.equal(nllR(y, psi, p, nd), nllRcpp(y, psi, p, nd))

# 64-bit
> benchmark(nllR(y, psi, p, nd), nllRcpp(y, psi, p, nd))
test replications elapsed relative user.self
1    nllR(y, psi, p, nd)          100   10.58 1.046489      9.35
2 nllRcpp(y, psi, p, nd)          100   10.11 1.000000      9.81
sys.self user.child sys.child
1     1.14         NA        NA
2     0.26         NA        NA

# 32-bit
> benchmark(nllR(y, psi, p, nd), nllRcpp(y, psi, p, nd))
test replications elapsed relative user.self
1    nllR(y, psi, p, nd)          100   16.75 1.000000     15.54
2 nllRcpp(y, psi, p, nd)          100   29.97 1.789254     29.81
sys.self user.child sys.child
1     1.08         NA        NA
2     0.01         NA        NA

> sessionInfo()
R version 2.12.0 (2010-10-15)
Platform: x86_64-pc-mingw32/x64 (64-bit)

...

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods
[7] base

other attached packages:
[4] Rcpp_0.9.0

loaded via a namespace (and not attached):
[1] tools_2.12.0
>

Thanks,
Richard
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20101227/70057632/attachment.htm>
```