CODEWITHSUNDEEP
rloopsprocessing-efficiencymemory-efficient
David Veitch
David Veitch

Reputation: 125

Taking maximum over dimension in an array in R

I am currently working with a very large array with dimension 5663x1000x100 in R. I would like to get 100 maximum values, which will be the maximum of each individual 5663x1000 matrix.

big_array = array(data=rnorm(566300000),dim=c(5663,1000,100))

Two methods I have tried so far include a for loop and apply (which intuitively should not be the fastest methods).

maximas = rep(0,100)

# Method 1 - Runs in 17 seconds
for(i in seq(1,100)){
    maximas[i]=max(big_array[,,i])
}

# Method 2 - Runs in 36 seconds
apply(big_array,3,max)

I would think because of the array data structure there is an even faster way to run this. I have considered pmax() but from what I see I would have to reshape my data (which given the array is almost 4GB I do not want to create another object). This code is already part of code which is being parallelized so I am unable to parallelize it any further.

Any ideas would help greatly!

Upvotes: 0

Views: 429

Answers (1)

ekoam
ekoam

Reputation: 8844

Why not just do that with Rcpp and RcppArmadillo? Try this

library(Rcpp)
library(RcppArmadillo)

cppFunction('NumericVector max_slice(const arma::cube& Q) {
  int n = Q.n_slices; 
  NumericVector out(n);
  for (int i; i < n; i++) {
    out[i] = Q.slice(i).max();
  }
  return out;
}', depends = "RcppArmadillo")

str(big_array)
max_slice(big_array)

Output

> str(big_array)
 num [1:5663, 1:1000, 1:100] -0.282 -0.166 1.114 -0.447 -0.255 ...
> max_slice(big_array)
  [1] 5.167835 4.837959 5.026354 5.211833 5.054781 5.785444 4.782578 5.169154 5.427360 5.271900 5.197460 4.994804 4.977396 5.093390 5.124796 5.221609
 [17] 5.124122 4.857690 5.230277 5.217994 4.957608 5.060677 4.943275 5.382807 5.455486 5.226405 5.598238 4.942523 5.096521 5.000764 5.257607 4.843708
 [33] 4.866905 5.125437 5.662431 5.224198 5.026749 5.349403 4.987372 5.228885 5.456373 5.576859 5.166118 5.124967 4.991101 5.210636 5.057471 5.005961
 [49] 5.223063 5.182867 5.333683 5.528648 5.015871 4.837031 5.311825 4.981555 5.876951 5.145006 5.107017 5.252450 5.219044 5.310852 5.081958 5.210729
 [65] 5.439197 5.034269 5.339251 5.567369 5.117237 5.382006 5.332199 5.032523 5.622024 5.008994 5.537377 5.279285 5.175870 5.056068 5.019422 5.616507
 [81] 5.141175 4.948246 5.262170 4.961154 5.119193 4.908987 5.175458 5.328144 5.127913 5.816863 4.745966 5.507947 5.226849 5.247738 5.336941 5.134757
 [97] 4.899032 5.067129 5.615639 5.118519

Benchmark

cppFunction('NumericVector max_slice(const arma::cube& Q) {
  int n = Q.n_slices; 
  NumericVector out(n);
  for (int i; i < n; i++) {
    out[i] = Q.slice(i).max();
  }
  return out;
}', depends = "RcppArmadillo")

max_vapply <- function(x) vapply(seq_len(dim(x)[3]), function(i) max(x[,,i]), numeric(1))

microbenchmark::microbenchmark(
  max_vapply(big_array), max_slice(big_array), 
  times = 5L
)

Result

Unit: milliseconds
                  expr       min        lq      mean   median        uq       max neval cld
 max_vapply(big_array) 4735.7055 4789.6901 5159.8319 5380.784 5428.8319 5464.1480     5   b
  max_slice(big_array)  724.8582  742.0412  800.8939  747.811  833.2658  956.4935     5  a

Upvotes: 1

Related Questions