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