1

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

It works fine, but I had to fix lambda to some value. For my application, I need the values in the ith row to come from a poisson distribution with mean equal to the ith element of the lambda vector. This obviously worked fine in the foreach loop., but I have yet to find a way to make it work in an apply loop.

My questions are:

  1. Is it possible to have the apply function "know" on which row it is operating and pass a corresponding argument to a function?
  2. Is there a way to work with foreach and sparse matrices without the need of creating a standard matrix and converting it into a sparse one in the next step?
  3. 如果以上都不是，我有没有办法手动将任务分配给 R 的从属进程 - 也就是说，我可以具体告诉一个进程在第 1 列上工作，另一个在第 2 列上工作，依此类推，每个创建一个稀疏向量，并且仅在最后一步中组合这些向量。
4

1 回答 1

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

于 2014-06-03T15:20:15.670 回答