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中元素明智的矩陣乘法

回答

相關問題