<div dir="ltr">The important question is whether the as<arma::mat> copies the contents of the matrix from the SEXP passed by R. Creating an Rcpp::NumericMatrix does not copy the contents and converting the NumericMatrix to an arma::mat with that trailing 'false' also just copies the pointer to the contents, not the contents itself. I'm not sure what the "as" method does.<div>
<br></div><div>I prefer the approach in RcppEigen where you have two different types of objects with "as" methods, the Eigen::MatrixXd object (copies the contents) and the Eigen::Map<Eigen::MatrixXd> object which only copies the pointer. I generally declare the objects created from SEXP's as</div>
<div><br></div><div>const Eigen::Map<Eigen::MatrixXd> X(Rcpp::as< ...>(Xs));</div><div><br></div><div>so the compiler can take responsibility for ensuring that I don't unintentionally modify the contents.</div>
<div><br></div><div>There are several examples in the RcppEigen vignette that Dirk and I wrote.<br><div><br></div><div><br></div></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Jan 9, 2013 at 11:58 PM, Paul Johnson <span dir="ltr"><<a href="mailto:pauljohn32@gmail.com" target="_blank">pauljohn32@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello, Rcpp land.<br>
<br>
I'm still collecting idioms and examples for the usage of Rcpp. Today<br>
I'm looking at ways that people manage the importation of R matrices.<br>
<br>
Consider fasLm.r file in the RcppArmadillo package's examples folder.<br>
It uses this this three-step process to bring in a matrix from R. This<br>
leads up to the fLmTwoCasts function:<br>
<br>
src <- '<br>
Rcpp::NumericMatrix Xr(Xs);<br>
Rcpp::NumericVector yr(ys);<br>
int n = Xr.nrow(), k = Xr.ncol();<br>
arma::mat X(Xr.begin(), n, k, false);<br>
arma::colvec y(yr.begin(), yr.size(), false);<br>
...<br>
<br>
It seems to waste some effort. It allocates Xr, then an iterator, and<br>
then writes it into X. 3 objects to get one?<br>
<br>
If I were writing this, I would take the more direct Rcpp::as route,<br>
like the varSimulation.r example does<br>
in the same folder:<br>
<br>
src <- '<br>
arma::mat X = Rcpp::as<arma::mat>(Xs);<br>
arma::colvec y = Rcpp::as<arma::colvec>(ys);<br>
int n = X.n_rows; int k = X.n_cols;<br>
...<br>
<br>
I tested this, it gives the exact same answers.<br>
<br>
I was certain it would be quicker. It creates fewer objects, it is<br>
more direct. But, as far as I can see in a simple test, the Rcpp::as<br>
way is not faster. If anything, it is slower, by a minute fraction.<br>
<br>
I really thought it would be faster, though, and lighter on memory<br>
usage. What do you think? I'm looking at as.h in the Rcpp code to try<br>
to figure it out. But, wow! It that complicated, or what?<br>
<br>
<br>
Here's my benchmark example, which I include only as evidence that I<br>
really tried before asking here.<br>
<br>
## built from RcppArmadillo examples/fastLm.r<br>
library(inline)<br>
library(rbenchmark)<br>
<br>
src <- '<br>
Rcpp::NumericMatrix Xr(Xs);<br>
Rcpp::NumericVector yr(ys);<br>
int n = Xr.nrow(), k = Xr.ncol();<br>
arma::mat X(Xr.begin(), n, k, false);<br>
arma::colvec y(yr.begin(), yr.size(), false);<br>
int df = n - k;<br>
<br>
// fit model y ~ X, extract residuals<br>
arma::colvec coef = arma::solve(X, y);<br>
arma::colvec res = y - X*coef;<br>
<br>
double s2 = std::inner_product(res.begin(), res.end(),<br>
res.begin(), 0.0)/df;<br>
// std.errors of coefficients<br>
arma::colvec sderr = arma::sqrt(s2 *<br>
arma::diagvec(arma::pinv(arma::trans(X)*X)));<br>
<br>
return Rcpp::List::create(Rcpp::Named("coefficients")=coef,<br>
Rcpp::Named("stderr") =sderr,<br>
Rcpp::Named("df") =df);<br>
'<br>
<br>
fLmTwoCasts <- cxxfunction(signature(Xs="numeric", ys="numeric"),<br>
src, plugin="RcppArmadillo")<br>
<br>
<br>
src <- '<br>
arma::mat X = Rcpp::as<arma::mat>(Xs);<br>
arma::colvec y = Rcpp::as<arma::colvec>(ys);<br>
int n = X.n_rows; int k = X.n_cols;<br>
<br>
int df = n - k;<br>
<br>
// fit model y ~ X, extract residuals<br>
arma::colvec coef = arma::solve(X, y);<br>
arma::colvec res = y - X*coef;<br>
<br>
double s2 = std::inner_product(res.begin(), res.end(),<br>
res.begin(), 0.0)/df;<br>
// std.errors of coefficients<br>
arma::colvec sderr = arma::sqrt(s2 *<br>
arma::diagvec(arma::pinv(arma::trans(X)*X)));<br>
<br>
return Rcpp::List::create(Rcpp::Named("coefficients")=coef,<br>
Rcpp::Named("stderr") =sderr,<br>
Rcpp::Named("df") =df);<br>
'<br>
<br>
pjTwoCasts <- cxxfunction(signature(Xs="numeric", ys="numeric"),<br>
src, plugin="RcppArmadillo")<br>
<br>
<br>
N <- 10000<br>
<br>
X <- cbind(1, rnorm(N), rpois(N, lambda=0.8), rnorm(N), rnorm(N),<br>
rnorm(N), rnorm(N), rnorm(N), rnorm(N), rnorm(N), rnorm(N), rnorm(N),<br>
rnorm(N), rnorm(N))<br>
y <- rnorm(N)<br>
<br>
res <- benchmark(<br>
fLmTwoCasts(X, y),<br>
pjTwoCasts(X, y),<br>
columns = c("test", "replications", "relative",<br>
"elapsed", "user.self", "sys.self"),<br>
order="relative",<br>
replications=5000)<br>
<br>
print(res)<br>
<span class="HOEnZb"><font color="#888888"><br>
--<br>
Paul E. Johnson<br>
Professor, Political Science Assoc. Director<br>
1541 Lilac Lane, Room 504 Center for Research Methods<br>
University of Kansas University of Kansas<br>
<a href="http://pj.freefaculty.org" target="_blank">http://pj.freefaculty.org</a> <a href="http://quant.ku.edu" target="_blank">http://quant.ku.edu</a><br>
_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
</font></span></blockquote></div><br></div>