如何在R中高效存储12组参数生成的1000个模拟样本并优化性能？

原代码的性能瓶颈分析 #

原代码运行缓慢主要源于以下几点：

• 全局变量赋值开销：全局数组X和Y在for循环中逐次赋值，触发频繁内存拷贝，大幅拖慢速度。
• 冗余的mvrnorm调用：由于协方差矩阵是单位矩阵，所有变量独立服从标准正态分布，完全可以用更快的rnorm替代mvrnorm，后者的矩阵运算开销远高于前者。
• 低效循环：逐次处理每个replicate的for循环未利用R的向量化运算能力，性能远低于批量计算。
• 内存浪费：replicate生成的三维数组data占用大量内存，后续拆分使用时造成冗余存储。

优化后的模拟代码 #

library(MASS)
set.seed(2000)

n = 500
p = 600
reps = 1000

params <- expand.grid(
  p1 = c(400,550,570),
  c1 = c(0,0.1),
  c2 = c(0.1,0.5)
)

# 向量化模拟函数，消除内部循环
sims_opt <- function(p1, c1, c2) {
  # 批量生成所有X矩阵和误差项：向量化操作
  X_all <- matrix(rnorm(n * p * reps), nrow = n * reps, ncol = p)
  e_all <- rnorm(n * reps)
  
  # 构造beta向量
  beta <- c(rep(c1, p1), rep(c2, p - p1))
  
  # 批量计算所有Y值：利用矩阵乘法向量化特性
  Y_all <- X_all %*% beta + e_all
  
  # 整理为要求的三维数组格式
  Y_array <- array(Y_all, dim = c(n, 1, reps))
  X_array <- array(X_all, dim = c(n, p, reps))
  
  list(Y = Y_array, X = X_array)
}

# 用lapply批量处理所有参数组，替代手动for循环
simulations_opt <- lapply(1:nrow(params), function(i) {
  pars <- params[i, ]
  sims_opt(pars$p1, pars$c1, pars$c2)
})

优化说明 #

• 向量化数据生成：直接用matrix(rnorm(...))批量生成所有重复的X和误差项，避免逐次调用mvrnorm和循环。
• 批量计算Y：利用矩阵乘法的向量化特性，一次性计算所有replicate的Y值，完全消除内部循环。
• 避免全局变量：函数内部生成并返回结果，减少内存拷贝和全局环境的副作用。
• 内存高效利用：直接构造最终的三维数组，无中间冗余大数组存储。

高效存储方案 #

12组数据的内存占用较高（每组X约2.2GB，Y约4MB），可通过以下方式优化存储：

1. 分文件存储 #

用saveRDS将每组结果单独保存，避免一次性加载所有数据到内存：

# 创建存储目录
dir.create("sim_results", showWarnings = FALSE)

# 存储每组模拟结果
for(i in 1:nrow(params)) {
  saveRDS(simulations_opt[[i]], file = sprintf("sim_results/param_group_%d.rds", i))
}

# 读取单个组数据
group1_data <- readRDS("sim_results/param_group_1.rds")

2. 高效压缩格式存储 #

使用fst包存储为高压缩比的格式，兼顾读取速度和存储空间：

library(fst)

# 存储每组的X和Y
for(i in 1:nrow(params)) {
  write_fst(as.data.frame(simulations_opt[[i]]$X), sprintf("sim_results/X_group_%d.fst", i), compress = 100)
  write_fst(as.data.frame(simulations_opt[[i]]$Y), sprintf("sim_results/Y_group_%d.fst", i), compress = 100)
}

# 读取并恢复数组格式
X_group1 <- as.matrix(read_fst("sim_results/X_group_1.fst"))
dim(X_group1) <- c(n, p, reps)

3. 按需生成数据 #

如果不需要长期存储所有数据，可在后续分析时实时生成单组数据，避免占用大量内存。

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