[Rcpp-devel] [rcpp-devel] arma::mat constructor vs Rcpp::as and interaction with .insert_col

[Dale Smith]

Hello,

I had a problem this afternoon with an arma::mat constructor and the mat::insert_col method. Since I got the example from Dirk, I thought I would point out the problem I had and a solution. See

http://gallery.rcpp.org/articles/fast-linear-model-with-armadillo/

for background. Before I get to the example to reproduce the problem and the fix, let me say I'm posting from home as I can't post code from my work email; otherwise I'm fired.

I do believe Rcpp::as should be used to marshall from NumericMatrix to arma::mat, not the arma::mat constructor that accesses raw memory, as Dirk uses. I read the Armadillo documentation this afternoon and do believe the problem with .insert_col is due to the handling of memory in the constructor

mat(aux_mem*, n_rows, n_cols, copy_aux_mem = true, strict = true) 

see

http://arma.sourceforge.net/docs.html

Note that I've run similar examples on my Windows 7 box at work and on my OS X laptop (Mountain Lion) at home to confirm this is not os-dependent.

Thanks,
Dale Smith
dtsmith at mindspring.com

mat <- cbind(rnorm(50), rnorm(50), rnorm(50))
y <- rnorm(50)

cppFunction('
List fastLm(NumericVector yr, NumericMatrix Xr) {
  int n = Xr.nrow(), k = Xr.ncol();
   
  arma::mat X(Xr.begin(), n, k, false); 
  arma::colvec y(yr.begin(), yr.size(), false);
   
  arma::colvec coef = arma::solve(X, y);
  arma::colvec resid = y - X*coef;

  arma::colvec ones = arma::ones<arma::colvec>(X.n_rows);
  X.insert_cols(1, ones);
   
  double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-k));
  arma::colvec stderrest = arma::sqrt(
       sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
   
  return List::create(Named("coefficients") = coef,
                      Named("stderr")       = stderrest);
}', depends = c("RcppArmadillo"))


cppFunction('
List fastLmMod(NumericVector yr, NumericMatrix Xr) {
  int n = Xr.nrow(), k = Xr.ncol();
   
  arma::mat X = Rcpp::as<arma::mat>(Xr);
  arma::colvec y(yr.begin(), yr.size(), false);
   
  arma::colvec ones = arma::ones<arma::colvec>(X.n_rows);
  X.insert_cols(0, ones);
   
  arma::colvec coef = arma::solve(X, y);
  arma::colvec resid = y - X*coef;

  double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-k));
  arma::colvec stderrest = arma::sqrt(
       sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
   
  return List::create(Named("coefficients") = coef,
                      Named("stderr")       = stderrest);
}', depends = c("RcppArmadillo"))

fastLm(y, mat) # error is "error: Mat::init(): mismatch between size of auxiliary memory and requested size"
fastLmMod(y, mat) # coefficients are -0.15785722, -0.05668155, 0.15409366, -0.10079147
lm(y ~ mat) # confirm fastLmMod via lm function in R

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130411/82463f3b/attachment.html>