Вложенные циклы for в R с использованием функции foreach и библиотеки doParallel

Я пытаюсь вычислить косинусное сходство между столбцами в матрице. Я могу заставить его работать, используя стандартные циклы for, но когда я пытаюсь заставить его работать параллельно, чтобы код работал быстрее, он не дает мне того же ответа. Проблема в том, что я не могу получить тот же ответ, используя цикл foreach. Я подозреваю, что использую неправильный синтаксис, потому что у меня работали одиночные циклы foreach. Я попытался сделать второй цикл обычным циклом for и использовал параметр %:% с циклом foreach, но тогда функция даже не запускается.

Пожалуйста, смотрите мой прикрепленный код ниже. Заранее благодарю за любую помощь.

## Function that calculates cosine similarity using paralel functions.

#for calculating parallel processing
library(doParallel)

## Set up cluster on 8 cores

cl = makeCluster(8)

registerDoParallel(cl)

#create an example data
x=array(data=sample(1000*100), dim=c(1000, 100))

## Cosine similarity function using sequential for loops

cosine_seq =function (x) {

  co = array(0, c(ncol(x), ncol(x)))

  for (i in 2:ncol(x)) {
    for (j in 1:(i - 1)) {

      co[i, j] = crossprod(x[, i], x[, j])/sqrt(crossprod(x[, i]) * crossprod(x[, j]))
    }
  }

  co = co + t(co)

  diag(co) = 1

  return(as.matrix(co))

}

## Cosine similarity function using parallel for loops

cosine_par =function (x) {

  co = array(0, c(ncol(x), ncol(x)))

  foreach (i=2:ncol(x)) %dopar% {

    for (j in 1:(i - 1)) {

      co[i, j] = crossprod(x[, i], x[, j])/sqrt(crossprod(x[, i]) * crossprod(x[, j]))
    }
  }

  co = co + t(co)

  diag(co) = 1

  return(as.matrix(co))

}

## Calculate cosine similarity

tm_seq=system.time(
{

  x_cosine_seq=cosine_seq(x)

})

tm_par=system.time(
{

  x_cosine_par=cosine_par(x)

})

## Test equality of cosine similarity functions

all.equal(x_cosine_seq, x_cosine_par)

#stop cluster
stopCluster(cl)

r parallel-processing parallel-foreach
person Kevin    schedule 19.01.2015    source источник

Ответы (2)


arrow_upward
3
arrow_downward

Правильное распараллеливание вложенного цикла использует %:% (прочитайте здесь).

library(foreach)
library(doParallel)
registerDoParallel(detectCores())    
cosine_par1 <- function (x) {  
  co <- foreach(i=1:ncol(x)) %:%
    foreach (j=1:ncol(x)) %dopar% {    
      co = crossprod(x[, i], x[, j])/sqrt(crossprod(x[, i]) * crossprod(x[, j]))
    }
  matrix(unlist(co), ncol=ncol(x))
}

Я рекомендую вам написать его в Rcpp, а не запускать параллельно, потому что foreach(i=2:n, .combine=cbind) не всегда будет связывать столбцы в правильном порядке. Кроме того, в приведенном выше коде я удалил только нижнее треугольное условие, но время выполнения значительно медленнее, чем время выполнения непараллельного кода.

set.seed(186)
x=array(data=sample(1000*100), dim=c(1000, 100))
cseq <- cosine_seq(x)
cpar <- cosine_par1(x)
 all.equal(cpar, cseq)
#[1] TRUE
head(cpar[,1])
#[1] 1.0000000 0.7537411 0.7420011 0.7496145 0.7551984 0.7602620
head(cseq[,1])
#[1] 1.0000000 0.7537411 0.7420011 0.7496145 0.7551984 0.7602620

Дополнение: Для этой конкретной проблемы возможна (полу)векторизация cosine_seq; cosine_vec примерно в 40-50 раз быстрее, чем cosine_seq.

cosine_vec <- function(x){
  crossprod(x) / sqrt(tcrossprod(apply(x, 2, crossprod)))
}
all.equal(cosine_vec(x), cosine_seq(x))
#[1] TRUE
library(microbenchmark)
microbenchmark(cosine_vec(x), cosine_seq(x), times=20L, unit="relative")
#Unit: relative
#          expr      min       lq     mean   median       uq      max neval
# cosine_vec(x)  1.00000  1.00000  1.00000  1.00000  1.00000  1.00000    20
# cosine_seq(x) 55.81694 52.80404 50.36549 52.17623 49.56412 42.94437    20
person Khashaa    schedule 19.01.2015
comment
Спасибо за помощь, но я все еще не могу проверить параллельное представление с помощью функции all.equal R. Когда я использую синтаксис %:%, код cosine_par фактически дает только единичную матрицу, в то время как cosine_seq дает правильную матрицу, используя сходство косинуса. Также я смотрю на функции Rcpp, но мне не совсем понятно, как это сделать. Возможно, вы могли бы прикрепить код, который делает это и показывает, что 2 выхода эквивалентны? Спасибо. - person Kevin; 19.01.2015
comment
Спасибо, функция cosine_vec намного быстрее! Я пытался сделать это еще быстрее, используя mclapply вместо применения, так как я буду работать с гораздо большими матрицами, чем этот пример 1000x100. Созданная вами функция cosine_par1 имеет тот же результат, но работает намного медленнее. Однако нет необходимости перебирать все комбинации i,j, поскольку результат симметричен (именно поэтому я использовал условие нижнего треугольника). Разве функция foreach не будет работать с моими исходными операторами условия цикла i,j for, а затем добавит нижнее треугольное условие? - person Kevin; 20.01.2015
comment
Если вам нужно вставить какие-то команды между двумя foreach, как вы это сделаете? Я получаю сообщение об ошибке типа "%:%" was passed an illegal right operand - person George Dontas; 02.02.2016

arrow_upward
0
arrow_downward

Чтобы сделать вложенный цикл в foreach и использовать параллельную реализацию, есть два способа.

  1. %:% + %dopar%
  2. %dopar% + %do%

Обратите внимание, что для (1) фактически создается один объект foreach, и вы не можете ничего добавить между ними. В противном случае вы получите сообщение об ошибке: "%:%" was passed an illegal right operand.

А для (2) вы можете вставить между ними все, что захотите. Но не забудьте добавить foreach к аргументу .package во внешнем цикле, так как внутренний foreach использует пакет foreach.

Ниже приводится изящный способ решения задачи о матрице косинусов. Обратите внимание, что для иллюстрации (2) я добавил одну дополнительную строку, и, пожалуйста, не забудьте удалить ее для расчета матрицы косинуса.

testfunc <- function (x) {
  cl<-makeCluster(4)
  registerDoParallel(cl)
  co <- foreach(i=1:ncol(x), .combine = 'rbind', .packages = c('foreach', 'stats')) %dopar% {
    k <- rnorm(3)
    foreach (j=1:ncol(x), .combine = 'c') %do% {
      crossprod(x[, i], x[, j])/sqrt(crossprod(x[, i]) * crossprod(x[, j])) + k - k
    }
  }
  stopCluster(cl)
  co
}
x <- array(data=sample(20*10), dim=c(20, 10))
testfunc(x)
person Yu Yang    schedule 01.05.2020