Как сделать R код с массивом более эффективным? - PullRequest
Новая
       5

Как сделать R код с массивом более эффективным?

1 голос
/ 22 января 2020

У меня есть следующий код R, который не эффективен. Я хотел бы сделать это эффективным, используя R cpp. В частности, я не привык иметь дело с массивом в R cpp. Любая помощь будет оценена. 

myfunc <- function(n=1600,
                   m=400,
                   p = 3,
                   time = runif(n,min=0.05,max=4),
                   qi21 = rnorm(n),
                   s0c = rnorm(n),
                   zc_min_ecox_multi = array(rnorm(n*n*p),dim=c(n,n,p)),
                   qi=matrix(0,n,n),
                   qi11 = rnorm(p),
                   iIc_mat = matrix(rnorm(p*p),p,p)){

            for (j in 1:n){
              u<-time[j]
              ind<-1*(u<=time)
              locu<-which(time==u)
              qi2<- sum(qi21*ind) /s0c[locu]

              for (i in 1:n){
                qi1<-  qi11%*%iIc_mat%*%matrix(zc_min_ecox_multi[i,j,],p,1)
                qi[i,j]<- -(qi1+qi2)/m

              }
            }

}

Время вычислений составляет около 7,35 с. Мне нужно вызывать эту функцию снова и снова, может быть, 20 раз. 

system.time(myfunc())
   user  system elapsed 
   7.34    0.00    7.35

Ответы [ 2 ]

6 голосов
/ 23 января 2020

Первое, что нужно сделать, это профилировать ваш код: profvis::profvis({myfunc()}).

Что вы можете сделать, это предварительно вычислить qi11 %*% iIc_mat один раз. Вы получите (с небольшими улучшениями): 

precomp <- qi11 %*% iIc_mat

for (j in 1:n) {
  u <- time[j]
  qi2 <- sum(qi21[u <= time]) / s0c[time == u]

  for (i in 1:n) {
    qi1 <- precomp %*% zc_min_ecox_multi[i, j, ]
    qi[i, j] <- -(qi1 + qi2) / m
  }
}

, что в два раза быстрее (8 с c -> 4 с * c).

Векторизация i l oop тогда кажется простым: 

q1_all_i <- tcrossprod(precomp, zc_min_ecox_multi[, j, ])
qi[, j] <- -(q1_all_i + qi2) / m

(сейчас в 12 раз быстрее)

1 голос
/ 23 января 2020

И если вы хотите попробовать это в R cpp, вам сначала понадобится функция умножения матриц ... 

#include<Rcpp.h>
#include<numeric>
// [[Rcpp::plugins("cpp11")]]


Rcpp::NumericMatrix mult(const Rcpp::NumericMatrix& lhs,
                         const Rcpp::NumericMatrix& rhs)
{
  if (lhs.ncol() != rhs.nrow())
    Rcpp::stop ("Incompatible matrices");
  Rcpp::NumericMatrix out(lhs.nrow(),rhs.ncol());
  Rcpp::NumericVector rowvec, colvec;
  for (int i = 0; i < lhs.nrow(); ++i)
    {
      rowvec = lhs(i,Rcpp::_);
      for (int j = 0; j < rhs.ncol(); ++j)
      {
        colvec = rhs(Rcpp::_,j);
        out(i, j) = std::inner_product(rowvec.begin(), rowvec.end(),
                                      colvec.begin(), 0.);
      }
    }
  return out;
}

Затем портируйте свою функцию ... 

// [[Rcpp::export]]
Rcpp::NumericMatrix myfunc_rcpp( int n, int m, int p,
                                 const Rcpp::NumericVector& time,
                                 const Rcpp::NumericVector& qi21,
                                 const Rcpp::NumericVector& s0c,
                                 const Rcpp::NumericVector& zc_min_ecox_multi,
                                 const Rcpp::NumericMatrix& qi11,
                                 const Rcpp::NumericMatrix& iIc_mat)
{
  Rcpp::NumericMatrix qi(n, n);
  Rcpp::NumericMatrix outermat = mult(qi11, iIc_mat);

  for (int j = 0; j < n; ++j)
  {
    double qi2 = 0;
    for(int k = 0; k < n; ++k)
    {
      if(time[j] <= time[k]) qi2 += qi21[k];
    }
    qi2 /= s0c[j];
    for (int i = 0; i < n; ++i)
    {
      Rcpp::NumericMatrix tmpmat(p, 1);
      for(int z = 0; z < p; ++z)
      {
        tmpmat(z, 0) =  zc_min_ecox_multi[i + n*j + z*n*n];
      }
      Rcpp::NumericMatrix qi1 =  mult(outermat, tmpmat);
      qi(i,j) -= (qi1(0,0) + qi2)/m;
    }
  }
  return qi;
}

Тогда в R: 

my_rcpp_func <- function(n=1600,
                   m=400,
                   p = 3,
                   time = runif(n,min=0.05,max=4),
                   qi21 = rnorm(n),
                   s0c = rnorm(n),
                   zc_min_ecox_multi = array(rnorm(n*n*p),dim=c(n,n,p)),
                   qi11 = rnorm(p),
                   iIc_mat = matrix(rnorm(p*p),p,p))
{
  myfunc_rcpp(n, m, p, time, qi21, s0c, as.vector(zc_min_ecox_multi),
              matrix(qi11,1,p), iIc_mat)
}

Это, конечно, быстрее и дает те же результаты, что и ваша собственная функция, но это не быстрее, чем оптимизация в R, предложенная F Privé. Возможно, оптимизация кода C ++ могла бы привести к еще более быстрым результатам, но в конечном итоге вы умножаете 2 разумно большие матрицы вместе более чем в 2,5 миллиона раз, так что это никогда не будет таким быстрым. В конце концов, R достаточно хорошо оптимизирован для такого рода расчетов ...

...