Работа с несколькими ядрами и разреженными матрицами в R

Question

Я работаю над проектом, который требует больших матриц с большим количеством нулей. К сожалению, поскольку некоторые из этих матриц могут иметь более 1e10 элементов, работа со «стандартными» R-матрицами не является опцией из-за ограничений RAM. Кроме того, мне нужно работать с несколькими ядрами, так как вычисление может занять довольно много времени, а на самом деле не должно.

До сих пор я работал с пакетом 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) 

В лямбдах все довольно мало, так что только каждый 5-й элемент или так отличен от нуля, что делает его разумным магазин результаты в разреженной матрице.

К сожалению, теперь необходимо увеличить количество итераций от 1e6 до не менее 1e7, так что матрица, создаваемая петлей foreach, слишком велика для хранения в 8 ГБ ОЗУ. Теперь я хочу разделить задачи на шаги, каждая из которых имеет 1е6 итерации, и объединить их в единую разреженную матрицу.

теперь у меня есть следующие как идея:

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) 

Он отлично работает, но я должен был исправить лямбда до некоторого значения. Для моего приложения мне нужны значения в i-й строке, исходящие из распределения пуассонов со средним значением, равным i-му элементу вектора лямбда. Очевидно, он отлично работал в цикле foreach., Но мне еще предстоит найти способ заставить его работать в цикле приложения.

Мои вопросы:

1. Можно ли иметь функцию применять «знать», на котором ряд его работы и передать соответствующий аргумент функции?
2. Есть ли способ работать с foreach и разреженными матрицами без необходимости создания стандартной матрицы и преобразования ее в разреженный на следующем шаге?
3. Если не указано выше, есть ли способ вручную назначить задачи подчиненным процессам R - то есть, могу ли я конкретно сказать, что процесс работает над столбцом 1, другой - для работы с столбцом 2 и т. Д. каждый из которых создает разреженный вектор и объединяет их только на последнем шаге.

Answer 1

Я смог найти решение моей проблемы.

В моем случае я могу определить уникальный идентификатор для каждого из столбцов и задать для него параметры. Следующий код должен пояснить, что я имею в виду:

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) 

Хитрость заключается в том, чтобы применить функцию не над матрицей, а над вектором уникальных идентификаторов, которые выравнивают с индексами параметров.