在R中使用多核和稀疏矩阵

Question

我正在研究一个项目，该项目需要具有较大数量的零的大矩阵。不幸的是，由于其中一些矩阵可以有超过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)

lambda都非常小，因此只有每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行的值来自泊松分布，其平均值等于向量的第I元素。这显然在foreach循环工作得很好。，但我还没有找到一种方法使它在应用程序循环中工作。

我的问题是:

1. 是否有可能让应用函数"知道"它在哪一行操作并传递相应的参数给函数?
2. 是否有一种方法可以与foreach和稀疏矩阵一起工作，而无需在下一步中创建标准矩阵并将其转换为稀疏矩阵?
3. 如果以上都没有，是否有一种方法可以让我手动将任务分配给R的从属进程-也就是说，我可以明确地告诉一个进程在列1上工作，另一个在列2上工作等等，每个进程创建一个稀疏向量，并且只在最后一步将这些组合起来。

Answer 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向量。