CODEWITHSUNDEEP
rrcpp
Hello
Hello

Reputation: 181

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

I have a following R code which is not efficient. I would like to make this efficient using Rcpp. Particularly, I am not used to dealing with array in Rcpp. Any help would be appreciated.

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

              }
            }

}

Computing time is about 7.35 secs. I need to call this function over and over again, maybe 20 times.

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

Upvotes: 1

Views: 181

Answers (2)

Allan Cameron
Allan Cameron

Reputation: 173793

And if you want to try it in Rcpp, you will first need a function to multiply the matrices...

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

Then port your function...

// [[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;
}

Then in 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)
}

This is certainly faster, and gives the same results as your own function, but it's no quicker than the in-R optimizations suggested by F Privé. Maybe optimizing the C++ code could get things even faster, but ultimately you are multiplying 2 reasonably large matrices together over 2.5 million times, so it's never going to be all that fast. R is optimized pretty well for this kind of calculation after all...

Upvotes: 1

F. Priv&#233;
F. Priv&#233;

Reputation: 11728

First thing to do would be to profile your code: profvis::profvis({myfunc()}).

What you can do is precompute qi11 %*% iIc_mat once. You get (with minor improvements):

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

that is twice as fast (8 sec -> 4 sec).

Vectorizing the i loop then seems straightforward:

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

(12 times as fast now)

Upvotes: 6

Related Questions