Rcpp中的矩阵乘法

Question

首先，我是一个新手用户，所以忘记了我的一般无知。我正在寻找一种更快的替代方法来运行R中的％*％运算符。尽管较旧的帖子提示使用RcppArmadillo，但我已经尝试了2个小时来使RcppArmadillo无法成功运行。我总是遇到产生'意外...'错误的词汇问题。我发现在RCPP下面的功能，我可以让工作：

library(Rcpp) 
func <- ' 
NumericMatrix mmult(NumericMatrix m , NumericMatrix v, bool byrow=true) 
{ 
    if(! m.nrow() == v.nrow()) stop("Non-conformable arrays") ; 
    if(! m.ncol() == v.ncol()) stop("Non-conformable arrays") ; 

    NumericMatrix out(m) ; 

    for (int i = 0; i < m.nrow(); i++) 
    { 
    for (int j = 0; j < m.ncol(); j++) 
    { 
    out(i,j)=m(i,j) * v(i,j) ; 
    } 
    } 
    return out ; 
} 
' 

这个功能，但是，进行逐元素乘法和不表现为％*％。有没有简单的方法来修改上面的代码来实现预期的结果？

编辑：

我已经想出了使用RcppEigen似乎击败％*％功能：

etest <- cxxfunction(signature(tm="NumericMatrix", 
          tm2="NumericMatrix"), 
       plugin="RcppEigen", 
       body=" 
NumericMatrix tm22(tm2); 
NumericMatrix tmm(tm); 

const Eigen::Map<Eigen::MatrixXd> ttm(as<Eigen::Map<Eigen::MatrixXd> >(tmm)); 
const Eigen::Map<Eigen::MatrixXd> ttm2(as<Eigen::Map<Eigen::MatrixXd> >(tm22)); 

Eigen::MatrixXd prod = ttm*ttm2; 
return(wrap(prod)); 
       ") 

set.seed(123) 
M1 <- matrix(sample(1e3),ncol=50) 
M2 <- matrix(sample(1e3),nrow=50) 

identical(etest(M1,M2), M1 %*% M2) 
[1] TRUE 
res <- microbenchmark(
+ etest(M1,M2), 
+ M1 %*% M2, 
+ times=10000L) 

res 

Unit: microseconds 
     expr min lq  mean median  uq max neval 
etest(M1, M2) 5.709 6.61 7.414607 6.611 7.211 49.879 10000 
    M1 %*% M2 11.718 12.32 13.505272 12.621 13.221 58.592 10000 

Answer 1

我会鼓励尝试与RcppArmadillo制定出你的问题。使用它，因为这例子也通过调用RcppArmadillo.package.skeleton()创建为简单：

// another simple example: outer product of a vector, 
// returning a matrix 
// 
// [[Rcpp::export]] 
arma::mat rcpparma_outerproduct(const arma::colvec & x) { 
    arma::mat m = x * x.t(); 
    return m; 
} 

// and the inner product returns a scalar 
// 
// [[Rcpp::export]] 
double rcpparma_innerproduct(const arma::colvec & x) { 
    double v = arma::as_scalar(x.t() * x); 
    return v; 
} 

有一个在实际例子更多的代码，但是这应该给你一个想法。

Answer 2

有充分的理由依赖现有的库/包进行标准任务。在库中的例程是

• 优化
• 彻底测试
• 的良好手段保持代码的紧凑，人类可读的，并且易于维护。

因此，我认为使用RcppArmadillo或RcppEigen应该在这里更好。但是，为了回答你的问题，下面是一个可能的RCPP代码来执行矩阵乘法：

library(Rcpp) 
cppFunction('NumericMatrix mmult(const NumericMatrix& m1, const NumericMatrix& m2){ 
if (m1.ncol() != m2.nrow()) stop ("Incompatible matrix dimensions"); 
NumericMatrix out(m1.nrow(),m2.ncol()); 
NumericVector rm1, cm2; 
for (size_t i = 0; i < m1.nrow(); ++i) { 
    rm1 = m1(i,_); 
    for (size_t j = 0; j < m2.ncol(); ++j) { 
     cm2 = m2(_,j); 
     out(i,j) = std::inner_product(rm1.begin(), rm1.end(), cm2.begin(), 0.);    
    } 
    } 
return out; 
}') 

测试一下：

A <- matrix(c(1:6),ncol=2) 
B <- matrix(c(0:7),nrow=2) 
mmult(A,B) 
#  [,1] [,2] [,3] [,4] 
#[1,] 4 14 24 34 
#[2,] 5 19 33 47 
#[3,] 6 24 42 60 
identical(mmult(A,B), A %*% B) 
#[1] TRUE 

希望这有助于。

---

作为基准测试表明，上面的RCPP代码比较慢的r内置%*%运算符。我认为，虽然我RCPP代码肯定能得到改善，这将是很难被击败的背后%*%优化的代码在性能方面：

library(microbenchmark) 
set.seed(123)  
M1 <- matrix(rnorm(1e4),ncol=100) 
M2 <- matrix(rnorm(1e4),nrow=100) 
identical(M1 %*% M2, mmult(M1,M2)) 
#[1] TRUE 
res <- microbenchmark(
      mmult(M1,M2), 
      M1 %*% M2, 
      times=1000L) 
#> res 
#Unit: microseconds 
#   expr  min  lq  mean median  uq  max neval cld 
# mmult(M1, M2) 1466.855 1484.8535 1584.9509 1494.0655 1517.5105 2699.643 1000 b 
#  M1 %*% M2 602.053 617.9685 687.6863 621.4335 633.7675 2774.954 1000 a