[Rcpp-devel] speed fix

[Richard Chandler]

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:
[1] rbenchmark_0.3       inline_0.3.8         RcppArmadillo_0.2.10
[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>