元素矩阵乘法:R与Rcpp(如何加速此代码？)

Question

我是C++编程新手(Rcpp用于无缝集成R),我很欣赏一些关于如何加速计算的建议.

请考虑以下示例:

 testmat <- matrix(1:9, nrow=3)
 testvec <- 1:3

 testmat*testvec
   #      [,1] [,2] [,3]
   #[1,]    1    4    7
   #[2,]    4   10   16
   #[3,]    9   18   27

在这里,R循环使用testvec,松散地说,testvec"变成"一个testmat与此乘法目的相同的矩阵.然后返回Hadamard产品.我希望使用Rcpp,这就是我希望i矩阵中-th行的每个元素testmat与i向量的-th元素相乘testvec.我的基准测试告诉我,我的实现速度非常慢,我很乐意建议如何加快速度.这是我的代码:

首先,使用Eigen:

#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]

using namespace Rcpp;
using namespace Eigen;

// [[Rcpp::export]]
NumericMatrix E_matvecprod_elwise(NumericMatrix Xs, NumericVector ys){
  Map<MatrixXd> X(as<Map<MatrixXd> >(Xs));
  Map<VectorXd> y(as<Map<VectorXd> >(ys));

  int k = X.cols();
  int n = X.rows();
  MatrixXd Y(n,k) ;

  // here, I emulate R's recycling. I did not find an easier way of doing this. Any hint appreciated.      
  for(int i = 0; i < k; ++i) {
   Y.col(i) = y;
  }

  MatrixXd out = X.cwiseProduct(Y);
  return wrap(out);
}

这里我的实现使用Armadillo(调整为遵循Dirk的例子,见下面的答案):

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
arma::mat A_matvecprod_elwise(const arma::mat & X, const arma::vec & y){
  int k = X.n_cols ;
  arma::mat Y = repmat(y, 1, k) ;  // 
  arma::mat out = X % Y;  
return out;
}

使用R,Eigen或Armadillo对这些解决方案进行基准测试表明,Eigen和Armadillo都比R慢约2倍.有没有办法加快这些计算速度或至少与R一样快？是否有更优雅的方式来设置它？任何建议表示赞赏和欢迎.(我也鼓励对我一般的编程风格做一些评论Rcpp / C++.)

这里有一些可重复的基准测试:

 # for comparison, define R function:
 R_matvecprod_elwise <- function(mat, vec) mat*vec

n <- 50000
k <- 50
X <- matrix(rnorm(n*k), nrow=n)
e <- rnorm(n)

benchmark(R_matvecprod_elwise(X, e), A2_matvecprod_elwise(X, e), E_matvecprod_elwise(X,e),
          columns = c("test", "replications", "elapsed", "relative"), order = "relative", replications = 1000)

这产生了

                  test replications elapsed relative
1  R_matvecprod_elwise(X, e)         1000   10.89    1.000
2  A_matvecprod_elwise(X, e)         1000   26.87    2.467
3  E_matvecprod_elwise(X, e)         1000   27.73    2.546

正如你所看到的,我的Rcpp解决方案表现得相当悲惨.有什么方法可以做得更好吗？

Answer 1

如果你想加快计算速度,你必须要小心不要复制.这通常意味着牺牲可读性.这是一个版本,它不会复制并修改矩阵X.

// [[Rcpp::export]]
NumericMatrix Rcpp_matvecprod_elwise(NumericMatrix & X, NumericVector & y){
  unsigned int ncol = X.ncol();
  unsigned int nrow = X.nrow();
  int counter = 0;
  for (unsigned int j=0; j<ncol; j++) {
    for (unsigned int i=0; i<nrow; i++)  {
      X[counter++] *= y[i];
    }
  }
  return X;
}

这是我在机器上得到的

 > library(microbenchmark)
 > microbenchmark(R=R_matvecprod_elwise(X, e), Arma=A_matvecprod_elwise(X, e),  Rcpp=Rcpp_matvecprod_elwise(X, e))

Unit: milliseconds
 expr       min        lq    median       uq      max neval
    R  8.262845  9.386214 10.542599 11.53498 12.77650   100
 Arma 18.852685 19.872929 22.782958 26.35522 83.93213   100
 Rcpp  6.391219  6.640780  6.940111  7.32773  7.72021   100

> all.equal(R_matvecprod_elwise(X, e), Rcpp_matvecprod_elwise(X, e))
[1] TRUE

Answer 2

对于初学者,我会将Armadillo版本(界面)写为

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
arama::mat A_matvecprod_elwise(const arma::mat & X, const arma::vec & y){
  int k = X.n_cols ;
  arma::mat Y = repmat(y, 1, k) ;  // 
  arma::mat out = X % Y;  
  return out;
}

因为你正在进行额外的转换(虽然wrap()通过胶水代码添加了).这const &是概念性的(正如您通过上一个问题所了解到的,a SEXP是一个轻量级的指针对象,但是更好的风格).

你没有显示你的基准测试结果所以我不能评论矩阵大小等的影响pp.我怀疑你可能会在rcpp-devel上找到比这里更好的答案.你的选择.

编辑:如果你真的想要便宜又快速的东西,我会这样做:

// [[Rcpp::export]]
mat cheapHadamard(mat X, vec y) {
    // should row dim of X versus length of Y here
    for (unsigned int i=0; i<y.n_elem; i++) X.row(i) *= y(i);
    return X;
}

它没有分配新内存,因此速度更快,可能与R竞争.

测试输出:

R> cheapHadamard(testmat, testvec)
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    4   10   16
[3,]    9   18   27
R>