Matrix Multiplikation in Rcpp

Question

Zunächst einmal bin ich ein Anfänger Benutzer so vergessen meine allgemeine Ignoranz. Ich suche nach einer schnelleren Alternative zum% *% Operator in R. Obwohl ältere Beiträge die Verwendung von RcppArmadillo andeuten, habe ich 2 Stunden lang versucht, RcppArmadillo ohne Erfolg zu arbeiten. Ich stoße immer auf lexikalische Probleme, die "unerwartete ..." Fehler ergeben. Ich habe die folgende Funktion in RCPP gefunden, die ich tun kann, Arbeit machen:

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

Diese Funktion ist jedoch elementweise Multiplikation durchführt und nicht als% *% verhält. Gibt es eine einfache Möglichkeit, den obigen Code zu ändern, um das gewünschte Ergebnis zu erzielen?

EDIT:

Ich habe mit einer Funktion RcppEigen kommen, das% *% zu schlagen scheint:

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

ich ermutigen würde zu versuchen, Ihre Probleme mit RcppArmadillo zu erarbeiten. Die Benutzung ist so einfach, wie dieses Beispiel erstellt auch durch RcppArmadillo.package.skeleton() Aufruf:

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

Es tatsächlich mehr Code in dem Beispiel ist, aber dies sollte Ihnen eine Idee geben.

Answer 2

Es gibt gute Gründe, sich auf vorhandene Bibliotheken/Pakete für Standardaufgaben zu verlassen. Die Routinen in den Bibliotheken sind

• optimierte
• gründlich
• ein gutes Mittel getestet, um den Code kompakt zu halten, für Menschen lesbare und leicht zu pflegen.

Deshalb denke ich, dass die Verwendung von RcppArmadillo oder RcppEigen hier vorzuziehen sein sollte. Aber um Ihre Frage zu beantworten, ist ein möglicher RCPP Code, um eine Matrix-Multiplikation auszuführen:

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

Sagen sie es testen:

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 

Hoffnung, das hilft.

---

Als Benchmark-Tests zeigen, die oben RCPP Code langsamer als R-interne %*% Betreiber. Ich gehe davon aus, dass, während mein RCPP Code sicherlich verbessert werden kann, ist es schwer sein wird, den optimierten Code hinter %*% in Bezug auf Leistung zu schlagen:

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