Как сделать код R cpp эффективным с помощью нескольких циклов for? - PullRequest
Новая
       7

Как сделать код R cpp эффективным с помощью нескольких циклов for?

2 голосов
/ 25 января 2020

Я пытаюсь реализовать следующий код R cpp, вызывая из R. Время вычислений очень медленное. Здесь задействовано множество циклов for. 

#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::mat qpart(
      const int& n,
      const int& p,
      const int& m,
      arma::vec& G,
      arma::vec& ftime,
      arma::vec& cause,
      arma::mat& covs,
      arma::mat& S1byS0hat,
      arma::vec& S0hat,
      arma::vec& expz){

      arma::mat q(n,p);
      q.zeros();
      for(int u=0;u<n;++u){
         arma::mat q1(1,p);
         q1.zeros();
         for(int iprime=0;iprime<n;++iprime){
            for(int i=0;i<n;++i){
               if(cause(iprime)==1 & cause(i)>1 & (ftime(i) < ftime(u)) & (ftime(u) <= ftime(iprime))){
                   q1 += (covs.row(i) - S1byS0hat.row(iprime))*G(iprime)/G(i)*expz(i)/S0hat(iprime);
               }
             }
         }
         q.row(u) = q1/(m*m);
      }
return q;
}

Ниже приведен пример в R. 

#### In R ########
n = 2000
m = 500
p=3
G = runif(n)
ftime = runif(n,0.01,5)
cause = c(rep(0,600),rep(1,1000),rep(2,400))
covs = matrix(rnorm(n*p),n,p)
S1byS0hat = matrix(rnorm(n*p),p,n)
S0hat = rnorm(n)
expz = rnorm(n)

system.time( qpart(n,p,m,G,ftime,cause,covs,t(S1byS0hat),S0hat,expz))
user  system elapsed 
   21.5     0.0    21.5

Как мы видим, время вычислений очень велико.

Тот же код реализован в R, и время вычислений очень велико. 

q = matrix(0,n,p)
for(u in 1 : n){
    q1 <- matrix(0,p,1)
  for(iprime in 1 : n){
    for(i in 1 : n){
      if(cause[iprime]==1 & cause[i]>1 & (time[i]<time[u]) & (time[u] <= time[iprime])){
          q1 = q1 + (covs[i,] - S1byS0hat[,iprime])*G[iprime]/G[i]*expz[i]/S0hat[iprime]
      }
    }

  }
    q[u,] = q1/(m*m)
}

Ниже приведена формула, которую я пытаюсь реализовать. enter image description here

1 Ответ

3 голосов
/ 26 января 2020

Некоторые условия зависят только от u и iprime, поэтому вы можете проверить их намного раньше. Вы также можете заранее вычислить некоторые вещи. Это дает: 

arma::mat qpart2(
    double m,
    arma::vec& ftime,
    arma::vec& cause,
    arma::mat& covs,
    arma::mat& S1byS0hat,
    arma::vec& G_div_S0hat,
    arma::vec& expz_div_G){

  double m2 = m * m;

  int n = covs.n_rows;
  int p = covs.n_cols;

  arma::mat q(n, p, arma::fill::zeros);

  for (int u = 0; u < n; u++) {
    double ftime_u = ftime(u);
    for (int iprime = 0; iprime < n; iprime++) {
      if (cause(iprime) == 1 && ftime_u <= ftime(iprime)) {
        for (int i = 0; i < n; i++) {
          if (cause(i) > 1 && ftime(i) < ftime_u) {
            double coef = G_div_S0hat(iprime) * expz_div_G(i);
            for (int j = 0; j < p; j++) {
              q(u, j) += (covs(i, j) - S1byS0hat(iprime, j)) * coef;
            }
          }
        }
      }
    }
    for (int j = 0; j < p; j++)  q(u, j) /= m2;
  }

  return q;
}

Использование qpart2(m, ftime, cause, covs, t(S1byS0hat), G / S0hat, expz / G) занимает 3,7 se c (против 32 se c для вашего кода).

**

Небольшие замечания:

  • Есть ли причина, по которой вы используете структуры arma вместо R cpp?
  • Вы должны получить доступ к матрицам по столбцам, а не по строкам, это должно быть немного быстрее, потому что они хранятся по столбцам.
...