Rcpp 中的元素矩阵乘法

element wise matrix multiplication in Rcpp

本文关键字：元素 Rcpp      更新时间：2023-10-16

我正在尝试使用Rcpp加速一些R代码，该代码采用长度为 L (psi( 的向量和维度矩阵 (L，L( 并执行一些元素操作。有没有更有效的方法可以使用 Rcpp 执行这些元素级操作？

r： 

 UpdateLambda <- function(psi,phi){
                                      # updated full-day infection probabilites
    psi.times.phi <- apply(phi,1,function(x) x*psi)
    ## return Lambda_{i,j} = 1 - prod_{j} (1 - psi_{i,j,t} phi_{i,j})
    apply(psi.times.phi,2,function(x) 1-prod(1-x))
  }

.cpp：

#include <Rcpp.h>
#include <algorithm>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector UpdateLambdaC(NumericVector psi,
                NumericMatrix phi
                ){
  int n = psi.size();
  NumericMatrix psi_times_phi(n,n);
  NumericVector tmp(n,1.0);
  NumericVector lambda(n);
  for(int i=0; i<n;i++){
    psi_times_phi(i,_) = psi*phi(i,_);
   }
  for(int i=0; i<n;i++){
     // pi_{j} (1- lambda_{i,j,t})
    for(int j=0; j<n;j++){
      tmp[i] *= 1-psi_times_phi(i,j);
    }
     lambda[i] = 1-tmp[i];
  }
  return lambda;
}

您可以将循环apply替换为矢量化替代方案。

第一个相当于：

t(phi)*psi

第二个：

1-exp(colSums(log(1-psi.times.phi)))
#test data
phi <- matrix(runif(1e6),1e3)
psi <- runif(1e3)
#new function
UpdateLambda2 <- function(psi,phi) 1-exp(colSums(log(1-t(phi)*psi)))
#sanity check
identical(UpdateLambda(psi,phi),UpdateLambda2(psi,phi))
[1] TRUE
#timings
library(rbenchmark)
benchmark(UpdateLambda(psi,phi),UpdateLambda2(psi,phi))
                     test replications elapsed relative user.self sys.self
1  UpdateLambda(psi, phi)          100   16.05    1.041     15.06     0.93
2 UpdateLambda2(psi, phi)          100   15.42    1.000     14.19     1.19

好吧，它似乎没有太大区别，这非常令人惊讶，因为colSums通常比apply快得多。我不确定我使用的测试数据是否相关，因为输出都是第二部分中小于 1 的数字1的到期乘法次数。如果你想注意这么小的数字的细节，你可能最好在对数刻度中工作。

相关文章：