如何使带有数组的 R 代码更高效？

Question

我有一个效率不高的以下 R 代码。我想使用 Rcpp 来提高效率。特别是，我不习惯在 Rcpp 中处理数组。任何帮助将不胜感激。

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

Answer 1

首先要做的是分析您的代码：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 秒 -> 4 秒）。

矢量化 i 循环似乎很简单：

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

（现在快 12 倍）

Answer 2

如果你想在 Rcpp 中尝试，你首先需要一个函数来乘以矩阵...

#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)
}

这当然更快，并且给出与您自己的函数相同的结果，但它并不比 F Privé 建议的 in-R 优化快。也许优化 C++ 代码可以让事情变得更快，但最终你将 2 个相当大的矩阵相乘超过 250 万次，所以它永远不会那么快。毕竟 R 对这种计算进行了很好的优化...

如何使带有数组的 R 代码更高效？

How to make R code with an array to be more efficient?

r

rcpp