在R中的数组中获取最大维度答案

【问题标题】：Taking maximum over dimension in an array in R在R中的数组中获取最大维度
【发布时间】：2021-03-16 15:48:11
【问题描述】：

我目前正在使用 R 中尺寸为 5663x1000x100 的非常大的数组。我想获得 100 个最大值，这将是每个 5663x1000 矩阵的最大值。

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

到目前为止，我尝试过的两种方法包括 for 循环和应用（直觉上它不应该是最快的方法）。

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)

我认为由于数组数据结构，有一种更快的方法来运行它。我已经考虑过pmax()，但据我所知，我必须重塑我的数据（考虑到数组几乎是 4GB，我不想创建另一个对象）。此代码已经是正在并行化的代码的一部分，因此我无法进一步并行化它。

任何想法都会有很大帮助！

【问题讨论】：

"apply (从直觉上看不应该是最快的方法" 是一个强有力的陈述。是什么导致你得出这个结论？（顺便说一句：你的代码不完整，缺少一个关闭-paren.) pmax 不适合此用途。
顺便说一句，for 循环是maximas <- sapply(seq_len(dim(big_array)[3]), function(i) max(big_array[,,i]))。它比apply 更快，并且不需要预先分配maximas。 vapply(seq_len(dim(big_array)[3]), function(i) max(big_array[,,i]), numeric(1)).
每当我试图从某些东西中获得所有速度并且它乞求apply-like 逻辑时，vapply 一直是我的首选功能。它需要更多的努力，但可以获得一些好处。在这种情况下，vapply 使用小 100 倍的数据集几乎快 10%。
我看到 vapply 的速度提高了 20%，这太棒了！我开始考虑我的情况，因为矩阵对它们没有任何“好”的结构，使用 apply/vapply 是要走的路。

标签： r loops processing-efficiency memory-efficient

【解决方案1】：

为什么不直接使用Rcpp 和RcppArmadillo？试试这个

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)

输出

> 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

基准测试

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
)

结果

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

【讨论】：

感谢您的意见。我以前用过 Rcpp，看起来它可以帮助解决这个问题。