Rcppの行列乗算

Question

まずは初心者ですので、私の一般的な無知を忘れてしまいます。 Rの％*％演算子のより高速な代替方法を探しています。以前の投稿ではRcppArmadilloの使用を提案していましたが、RcppArmadilloを成功させるために2時間試しました。私はいつも予期せぬエラーを生じる字句問題にぶつかります。

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 ; 
} 
' 

この機能は、しかし、要素ごとの乗算を実行し、％※％として動作しません：私は、私は仕事にすることができませんRcppに次の関数を発見しました。意図した結果を得るために上記のコードを変更する簡単な方法はありますか？

編集：私はRcppArmadilloであなたの問題をうまくしようとすることをお勧めします

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

：

私は％※％を打つように見えるRcppEigenを使用して関数が出ています。それを使用することも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