在R中对多矩阵/图像执行操作的apply()更快替代方案

在R中，apply因存在循环开销，处理大数组时速度确实不理想，以下是几种更高效的替代方案：

1. 使用matrixStats包（最推荐，简单高效） #

matrixStats包提供了高度优化的向量化统计函数，底层基于C/C++实现，比基础R的apply快数倍。核心思路是将3D数组重塑为2D矩阵，把每个像素的所有切片值转为一列，再调用colMedians计算中位数：

# 安装并加载包
install.packages("matrixStats")
library(matrixStats)

# 调整数组维度顺序并重塑为2D矩阵
myImages_mat <- aperm(myImages, c(3,1,2))  # 将维度顺序改为(切片, 行, 列)
dim(myImages_mat) <- c(100, 1024*1024)     # 转为100行、1024*1024列的矩阵

# 计算每列中位数，再重塑回原图像尺寸
medianImage <- colMedians(myImages_mat)
dim(medianImage) <- c(1024, 1024)

2. 并行计算加速 #

如果机器有多核心，可结合并行计算进一步缩短耗时，比如用parallel包的parApply：

library(parallel)

# 创建并行集群（根据机器核心数调整，示例为4核心）
cl <- makeCluster(4)
# 将数据导出到集群节点
clusterExport(cl, "myImages")

# 并行计算逐像素中位数
medianImage <- parApply(cl, myImages, c(1,2), median)

# 关闭集群
stopCluster(cl)

注意：并行计算的额外开销在数据量较小时不明显，但对超大数组能有效提速，同时内存占用会相应增加。

3. 用Rcpp自定义高效函数（极致速度） #

若追求极限速度，可编写C++代码直接操作内存数据，完全规避R层面的循环开销：

install.packages("Rcpp")
library(Rcpp)

# 编写C++中位数计算函数
cppFunction('
NumericMatrix computeMedian3D(NumericVector arr, int nrow, int ncol, int nslice) {
  NumericMatrix res(nrow, ncol);
  for (int i = 0; i < nrow; i++) {
    for (int j = 0; j < ncol; j++) {
      NumericVector vals(nslice);
      for (int k = 0; k < nslice; k++) {
        vals[k] = arr[k * nrow * ncol + i * ncol + j];
      }
      std::sort(vals.begin(), vals.end());
      if (nslice % 2 == 1) {
        res(i,j) = vals[nslice/2];
      } else {
        res(i,j) = (vals[nslice/2 - 1] + vals[nslice/2])/2.0;
      }
    }
  }
  return res;
}
')

# 调用函数计算
medianImage <- computeMedian3D(myImages, 1024, 1024, 100)

内容的提问来源于stack exchange，提问作者mri