在 R 中处理多核和稀疏矩阵答案

【问题标题】：Working with multiple cores and sparse matrices in R在 R 中处理多核和稀疏矩阵
【发布时间】：2014-06-03 09:23:16
【问题描述】：

我正在从事一个需要具有大量零的大型矩阵的项目。不幸的是，由于这些矩阵中的一些可能有超过 1e10 个元素，因此由于 RAM 限制，不能使用“标准”R 矩阵。此外，我需要在多个内核上工作，因为计算可能需要很长时间，而且确实不应该。

到目前为止，我一直在使用 foreach 包，然后将结果（以标准矩阵形式出现）转换为稀疏矩阵。我忍不住想，一定有更聪明的办法。

这是我迄今为止所做的一个最小示例：

cl <- makeSOCKcluster(8)
registerDoSNOW(cl)

Mat <- foreach(j=1:length(lambda), .combine='cbind') %dopar% { 
  replicate(iter, rpois(n=1, lambda[j]))
}
Mat <- Matrix(Mat, sparse=TRUE)
stopCluster(cl)

lambdas 都非常小，因此只有每 5 个左右的元素不为零，因此将结果存储在稀疏矩阵中是明智的。

不幸的是，现在必须将迭代次数从 1e6 增加到至少 1e7，这样foreach 循环生成的矩阵太大而无法存储在 8GB 的 RAM 上。我现在想要做的是将任务拆分为每个具有 1e6 次迭代的步骤，并将它们组合成一个单一的稀疏矩阵。

我现在有以下想法：

library(Matrix)
library(snow)

cl <- makeSOCKcluster(8)

iter <- 1e6

steps <- 1e5
numsteps <- iter / steps

draws <- function(x, lambda, steps){
  replicate(n=steps, rpois(n=1, lambda=lambda))
}

for(i in 1:numsteps){
  Mat <- Matrix(0, nrow=steps, ncol=96, sparse=TRUE)

  Mat <- Matrix(
    parApply(cl=cl, X=Mat, MARGIN=2, FUN=draws, lambda=0.2, steps=steps)
    , sparse = TRUE)

  if(!exists("fullmat")) fullmat <- Mat else fullmat <- rBind(fullmat, Mat)

  rm(Mat)   
}

stopCluster(cl)

它工作正常，但我必须将 lambda 修复为某个值。对于我的应用程序，我需要第 i 行中的值来自泊松分布，其均值等于 lambda 向量的第 i 个元素。这在foreach 循环中显然工作得很好，但我还没有找到一种方法让它在应用循环中工作。

我的问题是：

是否可以让 apply 函数“知道”它正在操作哪一行并将相应的参数传递给函数？
有没有一种方法可以处理 foreach 和稀疏矩阵，而无需创建标准矩阵并在下一步将其转换为稀疏矩阵？
如果以上都不是，我有没有办法手动将任务分配给 R 的从属进程 - 也就是说，我是否可以明确告诉一个进程在第 1 列上工作，另一个在第 2 列上工作，依此类推，每个都创建一个稀疏向量，并且只在最后一步中组合它们。

【问题讨论】：

标签： r parallel-processing snow parallel-foreach

【解决方案1】：

我能够找到解决问题的方法。

在我的例子中，我可以为每一列定义一个唯一的 ID，并可以通过它来处理参数。下面的代码应该能说明我的意思：

library(snow)
library(Matrix)

iter <- 1e6
steps <- 1e5

# define a unique id
SZid <- seq(from=1, to=10, by=1)

# in order to have reproducible code, generate random parameters
SZlambda <- replicate(runif(n=1, min=0, max=.5))
SZmu <- replicate(runif(n=1, min=10, max=15))
SZsigma <- replicate(runif(n=1, min=1, max=3))

cl <- makeSOCKcluster(8)
clusterExport(cl, list=c("SZlambda", "SZmu", "SZsigma"))

numsteps <- iter / steps

MCSZ <- function(SZid, steps){ # Monte Carlo Simulation 
  lambda <- SZlambda[SZid]; mu <- SZmu[SZid]; sigma <- SZsigma[SZid];

  replicate(steps, sum(rlnorm(meanlog=mu, sdlog=sigma, 
                              n = rpois(n=1, lambda))
                       ))
}

for (i in 1:numsteps){
  Mat <- Matrix(
    parSapply(cl, X=SZid, FUN=MCSZ, steps=steps), sparse=TRUE)

  if(!exists("LossSZ")) LossSZ <- Mat else LossSZ <- rBind(LossSZ, Mat)
  rm(Mat)
}

stopCluster(cl)

诀窍是不是将函数应用于矩阵，而是应用于与参数索引对齐的唯一 id 向量。

【讨论】：