Rcpp中元素明智的矩阵乘法

我试图加速一些R代码与Rcpp需要一个长度为L（psi）的向量和一个维度矩阵（L，L）并执行一些元素操作。有没有一种更有效的方法来用Rcpp来完成这些基于元素的操作？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; 
}

来源

2013-02-14 scottyaz

它看起来像你可以避免在你的R代码中使用'apply'，并使用'colSums'和'log'ed变量来获得产品。 – James 2013-02-14 08:39:20

第一个'apply'相当于't（phi）* psi'，这应该是更快的 – James 2013-02-14 08:46:46

你认为日志，然后exp'ing和求和会比prods快得多？ – scottyaz 2013-02-14 09:10:16

可以代替你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。无论如何，如果你想记下这些小数字的细节，你可能会更好地使用对数级别的工作。

来源

2013-02-14 14:11:36 James

Rcpp中元素明智的矩阵乘法

回答

相关问题