élément sage multiplication de la matrice dans Rcpp

Question

J'essaie d'accélérer le code R avec Rcpp qui prend un vecteur de longueur L (psi) et une matrice de dimensions (L, L) et fait quelques opérations par élément. Existe-t-il un moyen plus efficace de faire ces opérations élémentaires avec 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; 
} 

Answer 1

Vous pouvez vous remplacer apply boucles alternatives avec vectorisées.

Le premier est équivalent à:

t(phi)*psi 

Et la seconde:

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 

Eh bien, il semble que cela ne fait pas beaucoup de différence, ce qui est très surprenant que colSums est généralement beaucoup plus rapide que apply. Je ne suis pas sûr si les données de test que j'ai utilisées sont pertinentes car la sortie est 1 nombre de multiplications de nombre inférieur à 1 dans la deuxième partie. Il vaudrait peut-être mieux travailler dans une échelle logarithmique si vous voulez noter les détails d'un si petit nombre.