问Rcpp中的矩阵乘法
EN

Stack Overflow用户

提问于 2016-05-12 15:43:03

回答 3查看 7.4K关注 0票数 14

首先，我是一个新手，所以忘记我一般的无知吧。我正在寻找一个更快的替代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

c++

rcpp

回答 3

Stack Overflow用户

发布于 2016-05-13 10:48:38

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

optimized
thoroughly测试了
是一种保持代码紧凑、人类可读和易于维护的好方法。

因此，我认为在这里使用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

票数 12

Stack Overflow用户

发布于 2016-05-12 16:15:13

我鼓励你尝试用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;
}

示例中实际上有更多的代码，但这应该会让您有所了解。

票数 5

Stack Overflow用户

发布于 2020-12-21 15:33:02

也可以使用以下方法：

NumericMatrix mmult(NumericMatrix m, NumericMatrix v)
{
  Environment base("package:base");
  Function mat_Mult = base["%*%"];
  return(mat_Mult(m, v));
}