[Rcpp-devel] RcppArmadillo: Building with Inline
Savitsky, Terrance
savitsky at rand.org
Sun Dec 12 03:44:22 CET 2010
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
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 set-up 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.
More information about the Rcpp-devel
mailing list