Transformar loop em purrr:map()

Prezados, tudo bem?

Estou trabalhando com o seguinte código e, na minha base de dados, ele está demorando muito a rodar. Será que alguém sabe alguma forma de agilizar o processamento desse código? Uma ajuda seria conseguir transformar a função cv usando o pacote purrr mas eu não sei como fazer isso :confused:.

Alguém que possa ajudar? @clente, o mestre do purrr, será que rola?

### cross-validation function
cv <- function(data, used.function, bandwidth.grid)
{
  n <- nrow(data)
  mse <- matrix(, nrow=length(bandwidth.grid), ncol=2)
  for (b in 1:length(bandwidth.grid)){
    cv.value <- rep(0, n-1)
    for (i in 1:(n-1)){
      new.data <- data[-i,]
      funcargs <- list(reg.x=new.data[,1],reg.y=new.data[,2],x=data[i,1], h = bandwidth.grid[b])
      cv.value[i] <- do.call(used.function, funcargs)
    }
    mse[b,] <- c(bandwidth.grid[b], 1/n*sum((new.data[,2]-cv.value)^2))
  }

  ## MSE
  colnames(mse) <- c("bandwidth", "mse")
  return(mse)
}

### kernel estimator usind nadaraya-watson:
fcn1 <- function(reg.x, reg.y, x, h){
  return(ksmooth(reg.x, reg.y, x.point = x, kernel = "normal", bandwidth = h)$y)
}

attach(cars)
### CV-score for kernel estimator:
cv(cbind(speed, dist), fcn1, seq(10))

Obrigado!

Rafael,

Infelizmente o purrr não resolverá o seu problema. Na maior parte dos casos ele demora o mesmo tempo que o for e em certos casos pode demorar até mais:

microbenchmark::microbenchmark(
  `for` = for (i in 1:10000) sqrt(i),
  map = purrr::map(1:10000, sqrt)
)
#> Unit: microseconds
#>  expr      min       lq      mean    median        uq      max neval
#>   for  492.205  511.393  595.8878  534.6605  629.2885 2049.713   100
#>   map 2838.881 2929.532 3319.8938 3000.2775 3425.4065 8751.696   100

Created on 2021-10-27 by the reprex package (v2.0.1)

Se você não quiser mudar o código, acho que a melhor solução é paralelizar o loop. Antes de eu tentar implementar isso, entretanto, queria saber se você já testou uma das funções prontas do R para cross-validation. Eu dei uma pesquisada rápida e encontrei duas que parecem fazer seleção de bandwidth: extremefit::bandwidth.CV() e kedd::h.ccv(). Elas não atendem aos seus requisitos? Elas também são lentas demais?

Oi, Caio. Obrigado!

A ideia é implementar mesmo a função, como um exercício. Você tem alguma ideia de como alterar o código e ele performar mais rápido?

Eu já não sei mais o que fazer.

Eu fiz duas novas versões. Na primeira eu simplifiquei o código como eu faria em um programa meu e na segunda eu executei o loop em paralelo com o pacote furrr. Não encontrei um jeito de otimizar mais que isso… O programa parecer ser O(bn), não tem muito jeito.

# Argumentos utilizados para testar os algoritmos
data <- cbind(
  sample(cars$speed, 10000, replace = TRUE),
  sample(cars$dist, 10000, replace = TRUE)
)
bw <- seq(10)


# ORIGINAL -----------------------------------------------------------------

### cross-validation function
cv <- function(data, used.function, bandwidth.grid)
{
  n <- nrow(data)
  mse <- matrix(, nrow=length(bandwidth.grid), ncol=2)
  for (b in 1:length(bandwidth.grid)){
    cv.value <- rep(0, n-1)
    for (i in 1:(n-1)){  # *
      new.data <- data[-i,]
      funcargs <- list(reg.x=new.data[,1],reg.y=new.data[,2],x=data[i,1], h = bandwidth.grid[b])
      cv.value[i] <- do.call(used.function, funcargs)
    }
    mse[b,] <- c(bandwidth.grid[b], 1/n*sum((new.data[,2]-cv.value)^2)) # *
  }

  ## MSE
  colnames(mse) <- c("bandwidth", "mse")
  return(mse)
}

### kernel estimator usind nadaraya-watson:
fcn1 <- function(reg.x, reg.y, x, h){
  return(ksmooth(reg.x, reg.y, x.point = x, kernel = "normal", bandwidth = h)$y)
}


# MAP ----------------------------------------------------------------------

fcn2 <- function(reg_x, reg_y, x, h) {
  ksmooth(reg_x, reg_y, "normal", h, x.points = x)$y
}

make_cv_value <- function(i, fun, data, b) {
  fun(data[-i, 1], data[-i, 2], data[i, 1], b)
}

cv2 <- function(data, fun, bandwidth_grid) {
  n <- nrow(data)

  mse <- c()
  for (b in bandwidth_grid) {
    cv_value <- purrr::map_dbl(seq_len(n - 1), make_cv_value, fun, data, b)
    mse <- append(mse, 1 / n * sum((data[-n + 1, 2] - cv_value)^2))
  }

  mse <- cbind(seq_along(mse), mse)
  colnames(mse) <- c("bandwidth", "mse")
  return(mse)
}


# PARALELO -----------------------------------------------------------------

library(furrr)
plan("multicore", workers = 4)

cv3 <- function(data, fun, bandwidth_grid) {
  n <- nrow(data)

  mse <- c()
  for (b in bandwidth_grid) {
    cv_value <- future_map_dbl(seq_len(n - 1), make_cv_value, fun, data, b)
    mse <- append(mse, 1 / n * sum((data[-n + 1, 2] - cv_value)^2))
  }

  mse <- cbind(seq_along(mse), mse)
  colnames(mse) <- c("bandwidth", "mse")
  return(mse)
}


# BENCHMARKS ---------------------------------------------------------------

microbenchmark::microbenchmark(
  `for` = cv(data, fcn1, bw),
  `map` = cv2(data, fcn2, bw),
  `par` = cv3(data, fcn2, bw),
  times = 10
)
#> Unit: seconds
#>  expr       min        lq      mean    median        uq       max neval
#>   for 26.466603 26.738418 27.013609 26.996009 27.334605 27.679403    10
#>   map 23.987672 24.193279 24.411272 24.301028 24.686168 25.075567    10
#>   par  7.743104  7.812796  7.929976  7.928873  8.021716  8.136829    10

Created on 2021-10-27 by the reprex package (v2.0.1)


Só mais duas coisas:

  1. Eu marquei duas linhas da sua versão com *. Eu não sei nada sobre esse algoritmo, mas me parece que na primeira delas você quis dizer 1:n ao invés de 1:(n-1) e, na segunda, data ao invés de new.data. Eu fiz as minhas versões baseadas na sua, mas estranhei você usar data[-n + 1, 2] na hora de calcular mse; por que usar a base sem a penúltima linha?

  2. Eu sempre gosto de dar dicas de como as pessoas podem melhorar seus códigos. Tenho 4 sugestões para você baseadas nesse programa:

    1. Não use attach(). Ele é uma armadilha que polui o seu ambiente com objetos invisíveis! Prefira sempre usar o $.

    2. Evite colocar . no nome das variáveis. Essa é uma prática antiga e desencorajada hoje em dia por causa das classes S3.

    3. Não precisa usar do.call() quando você passa uma função como argumento. Ela continua funcionando do mesmo jeito (veja como eu fiz no meu código).

    4. Prefira nomes mais esclarecedores para as suas funções. É muito fácil se confundir quando as funções não indicam nada do que elas fazem.

1 Curtida

Oi, @clente. Muuuuuito obrigado!

Único problema é que aqui no meu PC demora muito a rodar esses códigos que vc fez. Era para demorar ~20s cada um?

A versão paralela deveria demorar 7 segundos como indicado pelo benchmark no final. Mas isso não funciona no RStudio! Tem que rodar como script no prompt.

Acho que estou com algum problema:

> microbenchmark::microbenchmark(
+   `for` = cv(data, fcn1, bw),
+   `map` = cv2(data, fcn2, bw),
+   `par` = cv3(data, fcn2, bw),
+   times = 1
+ )

Unit: seconds
 expr      min       lq     mean   median       uq      max neval
  for 61.79373 61.79373 61.79373 61.79373 61.79373 61.79373     1
  map 64.50819 64.50819 64.50819 64.50819 64.50819 64.50819     1
  par 62.44482 62.44482 62.44482 62.44482 62.44482 62.44482     1

Aqui fiz times=1 e demorou 61 secs! Será que to fazendo algo errado?

Você rodou no terminal?

Acho que sim.

image