[Rcppdevel] RcppArmadillo: Building with Inline
Dirk Eddelbuettel
edd at debian.org
Sun Dec 12 07:32:55 CET 2010
On 11 December 2010 at 18:44, Savitsky, Terrance wrote:
 Hello, I take the same fastLm code implemented in fastLmPure in
 RcppArmadillo, but built and execute through inline. There is an
 overhead of a little over 4 seconds to execute a single iteration. By
Just to make sure we are on the same page:
 you do understand that this needs to be compiled, linked and loaded
before the first time it is run?
 you do understand that for the same very reason none of the timings we
posted and and wrote about in the post you refer were done using inline?
Dirk
 overhead, I mean that as I increase the row size of X, the marginal
 increase in runtime seems roughly comparable to the fastLmPure
 implementation in RcppArmadillo (as executed on my machine) so that the
 4 seconds seems like a fixed add. Then I write this post to make sure
 that the overhead is driven by 'plugin = RcppArmadillo' in cxxfunction
 (and not something I'm doing incorrectly). I'm in process of testing my
 Windows XP hardware/software setup before committing the work to code
 with RcppArmadillo. The code and cxxfunction use are shown, below:

 lmArma = function()
 {
 src < '
 Rcpp::NumericVector yr(ys); // creates Rcpp
 vector from SEXP
 Rcpp::NumericMatrix Xr(Xs); // creates Rcpp
 matrix from SEXP
 int n = Xr.nrow(), k = Xr.ncol();

 arma::mat X(Xr.begin(), n, k, false); // reuses memory
 and avoids extra copy
 arma::colvec y(yr.begin(), yr.size(), false);

 arma::colvec coef = arma::solve(X, y); // fit model y ~
 X
 arma::colvec res = y  X*coef; // residuals

 double s2 = std::inner_product(res.begin(), res.end(),
 res.begin(), double())/(n  k);
 // std.errors of
 coefficients
 arma::colvec std_err = arma::sqrt(s2 * arma::diagvec(
 arma::inv(arma::trans(X)*X) ));

 return Rcpp::List::create(Rcpp::Named("coefficients") = coef,
 Rcpp::Named("stderr") = std_err,
 Rcpp::Named("df") = n  k
 );
 '

 fun < cxxfunction(signature(ys="numeric", Xs="numeric"), src,
 plugin="RcppArmadillo")
 }

 checkLmArma = function(y,X)
 {
 fun = lmArma()
 res = fun(y,X)
 return(res)
 }

 set.seed(42)
 n < 10000
 k < 9
 X < cbind( rep(1,n), matrix(rnorm(n*k), ncol=k) )
 truecoef < 1:(k+1)
 y < as.numeric(X %*% truecoef + rnorm(n))
 # N < 100

 tic()
 ff = checkLmArma(y,X)
 toc()

 tic()
 gg = lm(y ~ X  1)
 toc()

 Thanks, Terrance Savitsky

 __________________________________________________________________________

 This email message is for the sole use of the intended recipient(s) and
 may contain confidential information. Any unauthorized review, use,
 disclosure or distribution is prohibited. If you are not the intended
 recipient, please contact the sender by reply email and destroy all copies
 of the original message.

 _______________________________________________
 Rcppdevel mailing list
 Rcppdevel at lists.rforge.rproject.org
 https://lists.rforge.rproject.org/cgibin/mailman/listinfo/rcppdevel

Dirk Eddelbuettel  edd at debian.org  http://dirk.eddelbuettel.com
More information about the Rcppdevel
mailing list