foreach

O pacote foreach possui uma forma diferente de executar códigos repetidamente no R. O principal diferencial do pacote é ter a opção de se realizar o loop em paralelo usando os processadores da sua máquina. Desta forma, operações que levam minutos utilizando a abordagem sequencial podem ser reduzidas à segundos de computação.

Primeiramente é interessante observar a forma como o pacote realiza o loop sequencial. Para isso utiliza-se o operador %do%, veja o exemplo abaixo:

raiz <- foreach(i=1:3) %do% sqrt(i)
raiz

## [[1]]
## [1] 1
## 
## [[2]]
## [1] 1.414214
## 
## [[3]]
## [1] 1.732051

Perceba que por default a saída é uma lista, diferente do loop convencional (for) que retorna um vetor. De fato, o parâmetro .combine retorna o resultado do loop da maneira que o usuário achar mais adequado. Essa funcionalidade pode ser muito útil dependendo da aplicação.

raiz_c <- foreach(i=1:3, .combine = 'c') %do% sqrt(i)
raiz_c

## [1] 1.000000 1.414214 1.732051

raiz_rbind <- foreach(i=1:3, .combine = 'rbind') %do% sqrt(i)
raiz_rbind

##              [,1]
## result.1 1.000000
## result.2 1.414214
## result.3 1.732051

Para realizar o loop em paralelo utiliza-se o operador dopar. Além disso, é necessário criar um registro de quantos processadores serão utilizados no processamento. Esse registro pode ser feito utilizando o pacote doParallel.

registerDoParallel(detectCores())

n <- 10
unif <- runif(n = n)

system.time(
  foreach(i = 1:n, .combine = 'c') %do% {
    Sys.sleep(unif)
    i
  }
)

##    user  system elapsed 
##   0.005   0.000   0.584

system.time(
  foreach(i = 1:n, .combine = 'c') %dopar% {
    Sys.sleep(unif)
    i
  }
)

##    user  system elapsed 
##   0.016   0.058   0.086

Com esse simples exemplo pode-se observar a grande diferença no tempo de execução entre o loop sequencial e o loop paralelizado.

Para mostrar a potencialidade da função vamos criar alguns exemplos. É interessante observarmos a diferença de tempo quando estamos trabalhando com algum processo custoso computacionalmente.

Normalmente em artigos, teses e dissertações em Estatística utilizam-se simulações extensivamente. Nestas simulações por muitas vezes variam-se alguns parâmetros do modelo e avaliam-se os resultados para cada configuração. Esse processo é chamado de análise de sensibilidade.

Neste sentido, vamos supor que temos acesso à \(n\) covariáveis e vamos implementar um modelo de regressão linear para cada combinação possível dessas \(n\) covariáveis (modelo guloso de seleção de covariáveis).

library(mvtnorm)

nCov <- n                                           
n <- 100000
seq_n <- 1:nCov
cov <- data.frame(rmvnorm(n = n, mean = rep(0, nCov)))
data <- cov
data$y <- rnorm(n)
comb <- Map(combn, list(seq_n), seq_along(seq_n), simplify = FALSE)
vars <- unlist(comb, recursive = FALSE)

t1 <- system.time(
  res_do <- foreach(i = 1:(2^nCov - 1), .combine = 'rbind') %do% {
    formula <- paste("y ~ ", paste0("X", vars[[i]], collapse = " + "))
    reg <- lm(formula = formula, data = data)
    df <- data.frame(i = i, AIC = AIC(reg))
    df
  }
)

t1[3]

t2 <- system.time(
  res_dopar <- foreach(i = 1:(2^nCov - 1), .combine = 'rbind') %dopar% {
    formula <- paste("y ~ ", paste0("X", vars[[i]], collapse = " + "))
    reg <- lm(formula = formula, data = data)
    df <- data.frame(i = i, AIC = AIC(reg))
    df
  }
)

t2[3]

## elapsed 
## 200.713

## elapsed 
##  57.297

Pode-se ver que o tempo de processamento do loop em paralelo é muito menor do que o loop sequencial. O tempo computacional do código paralelizado é igual a 28.55% do tempo computacional do código sequencial.

Neste caso 100% do código é paralelizável e desta forma as leis de Amdhal e Gustafson possuem a mesma estimativa de aceleração no tempo de execução. As estimativas de aceleração de ambas as leis é 4, desta forma espera-se que de fato, o código paralelizado seja 4 vezes mais rápido.

Outro contexto em que a paralelização pode ser útil é quando temos um modelo de aprendizado de máquina com diversos parâmetros de tunagem. Em geral testa-se algumas configurações destes parâmetros e escolhe-se o modelo que minimiza alguma função objetivo. Para ilustrar esse caso vamos utilizar o modelo SVM para classificação.

Para ajustar um modelo de SVM no R vamos utilizar a função svm do pacote e1071. Os hiperparâmetros que serão variados serão: scale, kernel, degree, gamma e cost. Para melhor entendimento dos parâmetros veja ?e1071::svm.

library(e1071)

#--- Criando uma base de dados
nCov <- 10
n <- 5000
X <- rmvnorm(n = n, mean = rep(0, nCov))
betas <- rnorm(n = nCov)
Xbeta <- X%*%betas 
y <- ifelse(Xbeta < 0.5, 0, 1)
data <- data.frame(y = y, X)
#--- Separando a base em treino e teste
id <- 1:nrow(data)
testindex <- sample(id, trunc(length(id)/3))
testset <- data[testindex,]
trainset <- data[-testindex,]
#--- Criando os grids para avaliação dos parâmetros de tunagem
scale = c(TRUE, FALSE)
kernel = c("linear", "polynomial", "radial", "sigmoid")
degree = 2:4
gamma = 1/seq(nrow(trainset)-(n/10), nrow(trainset) + (n/10), length.out = 4)
cost = seq(1, 10, length.out = 4)
parameters <- expand.grid(scale = scale, 
                          kernel = kernel, 
                          degree = degree,
                          gamma = gamma,
                          cost = cost)
#--- SVM sem paralelização
t1 <- system.time(
  svm.output_1 <- foreach(i = 1:nrow(parameters), .combine = 'rbind') %do% {
  #--- Modelo
  svm.model <- svm(formula = y ~ .,                     
                   data = trainset, 
                   scale = parameters[i, 1], 
                   kernel = parameters[i, 2], 
                   degree = parameters[i, 3], 
                   gamma = parameters[i, 4], 
                   cost = parameters[i, 5],
                   type = "C-classification")
  #--- Predição
  svm.pred <- predict(svm.model, testset[,-1])
  #--- Métrica de acurácia
  svm.confusao <- table(pred = svm.pred, true = testset[,1])
  svm.acuracia <- sum(diag(svm.confusao))/sum(svm.confusao)
  #--- Retornando os resultados
  cbind(parameters[i,], acuracia = svm.acuracia)
}
)
#--- SVM com paralelização
t2 <- system.time(
  svm.output_2 <- foreach(i = 1:nrow(parameters), .combine = 'rbind') %dopar% {
  #--- Modelo
  svm.model <- svm(formula = y ~ .,
                   data = trainset, 
                   scale = parameters[i, 1], 
                   kernel = parameters[i, 2], 
                   degree = parameters[i, 3], 
                   gamma = parameters[i, 4], 
                   cost = parameters[i, 5],
                   type = "C-classification")
  #--- Predição
  svm.pred <- predict(svm.model, testset[,-1])
  #--- Métrica de acurácia
  svm.confusao <- table(pred = svm.pred, true = testset[,1])
  svm.acuracia <- sum(diag(svm.confusao))/sum(svm.confusao)
  #--- Retornando os resultados
  cbind(parameters[i,], acuracia = svm.acuracia)
}
)

t1[3]
t2[3]
all.equal(svm.output_1, svm.output_2) 
best.config <- svm.output_2[which.max(svm.output_2$acuracia), ]
best.config

## elapsed 
## 231.401

## elapsed 
## 101.834

## [1] TRUE

##   scale kernel degree        gamma cost  acuracia
## 1  TRUE linear      2 0.0003528582    1 0.9957983

Neste problema vemos que o código com a paralelização teve um tempo computacional equivalente a 44.01% do tempo computacional da versão não paralelizada. O que neste contexto reflete em um ganho de 2.1594 minutos no tempo de processamento.

Pensando nas leis de Amdahl e Gustafson (ambas resultariam em uma aceleração de 4 vezes para este caso em que o código é 100% paralelizável) a estimativa para a aceleração no tempo de execução do código é acertiva.

Abaixo temos duas figuras, uma mostrando o uso dos processadores enquanto o loop sequencial rodava e outra enquanto o loop paralelizado rodava. É possível perceber que de fato o loop sequencial utiliza apenas um processador enquanto o loop paralelizado utiliza todo o potencial computacional.

Processadores durante o loop não paralelizado

Processadores durante o loop paralelizado

Com isso, é interessante ressaltar que quanto maior o número de processadores maior o ganho ao se paralelizar um código.

Uma ideia interessante é investigar o que ocorre quando declararamos mais processadores do que o próprio número de processadores da máquina. Isto pode fazer sentido quando o processo a ser executado consome pouco processamento. Obviamente, declarando muito mais processadores do que o número existente podemos cair num problema de lentidão paralela (exemplo do garçom). No exemplo abaixo temos a comparação dos tempos de execução considerando diferentes números de processadores (lembrando que a minha máquina possui 12 processadores).

function_test <- function(
  cores, 
  X = rmvnorm(n = 2000, mean = rep(0, 50)), 
  sleep.time = 0.0001, n
) {
  #--- Registrando o número de processadores
  registerDoParallel(cores)
  #--- Loop
  foreach(i = 1:n, .combine = 'c') %dopar% {
    #--- Custo de processamento
    result <- X%*%t(X)
    #--- Custo de tempo
    Sys.sleep(sleep.time)
  }
} 
X = rmvnorm(n = 150, mean = rep(0, 150))
sleep.time <- 1.5*system.time(X%*%t(X))
n <- 100

t001 <- system.time(function_test(cores = 1, X = X, sleep.time = sleep.time, n = n))[3]
t002 <- system.time(function_test(cores = 2, X = X, sleep.time = sleep.time, n = n))[3]
t004 <- system.time(function_test(cores = 4, X = X, sleep.time = sleep.time, n = n))[3]           
t008 <- system.time(function_test(cores = 8, X = X, sleep.time = sleep.time, n = n))[3]
t012 <- system.time(function_test(cores = 12, X = X, sleep.time = sleep.time, n = n))[3]
t016 <- system.time(function_test(cores = 16, X = X, sleep.time = sleep.time, n = n))[3]
t025 <- system.time(function_test(cores = 25, X = X, sleep.time = sleep.time, n = n))[3]
t032 <- system.time(function_test(cores = 32, X = X, sleep.time = sleep.time, n = n))[3]
t050 <- system.time(function_test(cores = 50, X = X, sleep.time = sleep.time, n = n))[3]
t064 <- system.time(function_test(cores = 64, X = X, sleep.time = sleep.time, n = n))[3]
t075 <- system.time(function_test(cores = 75, X = X, sleep.time = sleep.time, n = n))[3]
t090 <- system.time(function_test(cores = 90, X = X, sleep.time = sleep.time, n = n))[3]
t100 <- system.time(function_test(cores = 100, X = X, sleep.time = sleep.time, n = n))[3]

tempos <- data.frame(tempos = c(t001, t002, t004, 
                                t008, t012, t016, 
                                t025, t032, t050, 
                                t064, t075, t090, 
                                t100),
                     processadores = c(1, 2, 4, 8,
                                       12, 16, 25, 
                                       32, 50, 64, 
                                       75, 90, 100))

tempos

##    tempos processadores
## 1   1.118             1
## 2   0.507             2
## 3   0.228             4
## 4   0.129             8
## 5   0.117            12
## 6   0.127            16
## 7   0.193            25
## 8   0.213            32
## 9   0.294            50
## 10  0.360            64
## 11  0.411            75
## 12  0.485            90
## 13  0.540           100

Não utilizei a função benchmark para garantir que em nenhum momento houvessem processos concorrentes. Vejamos graficamente o que ocorreu nesta simulação:

ggTime <- ggplot(data = tempos, aes(x = processadores, y = tempos)) +
  geom_line(col = colLine, size = sizeLine) + 
  geom_point(col = colPoints, size = sizePoint) +
  xlab("Número de processadores") +
  ylab("Tempo de processamento") +
  tema

ggplotly(ggTime)

É possível ver que o tempo de execução tende a decrescer até o momento em que consideramos 12 processadores. Depois deste ponto o tempo computacional começa a crescer (inclusive podendo ser menos vantajoso que o loop sequencial). Obviamente o interessante seria realizar um processo similar ao realizado pelo benchmark, ou seja, executar algumas vezes o mesmo processo e em seguida observar alguma medida de tendência central para melhor representar o custo computacional do processo.

Diversos outros pacotes do R se propõe a realizar a paralelização tal como o pacote snow e o pacote parallel. Porém, o pacote foreach parece atender bem às expectativas além de ser muito simples e intuitivo.

Atualmente o tipo de paralelismo mais utilizado e com melhores resultados é o paralelismo utilizando placas gráficas (GPU) que vamos discutir a seguir.