Optimización bayesiana de hiperparámetros


Versión PDF: Github

Más sobre ciencia de datos: cienciadedatos.net o joaquinamatrodrigo.github.io

Introducción


La mayoría de modelos de machine learning tienen una serie de parámetros que no pueden aprenderse de los datos, sino que tienen que ser establecidos previo entrenamiento. Ha estos parámetros se les suele conocer como hiperparámetros.

Dado que no existe forma de conocer de antemano cuáles son los valores adecuados, la única forma de identificarlos es probando diferentes combinaciones y evaluarlas mediante métodos de validación por resampling (cross-validation, bootstrapping). A este proceso se le conoce como model tuning.

Dos de las estrategias más empleadas son grid search y random search. En la primera, los valores estudiados de cada hiperparámetro se distribuyen de forma uniforme dentro de un rango delimitado por el analista. En el segundo caso, los datos son aleatorios dentro de ese rango.

Aunque estas dos estrategias son totalmente válidas y generan buenos resultados, sobretodo cuando se tiene criterio para acotar el rango de búsqueda, comparten una carencia común: ninguna tiene en cuenta los resultados obtenidos hasta el momento, lo que les impide focalizar la búsqueda en las regiones de mayor interés y evitar regiones innecesarias. Una alternativa es la búsqueda de hiperparámetros con métodos de optimización bayesiana.

En términos generales, la optimización bayesiana de hiperparámetros consiste en crear un modelo probabilístico en el que el valor de la función objetivo es la métrica de validación del modelo (rmse, auc, precisión..). Con esta estrategia, se consigue que la búsqueda se vaya redirigiendo en cada iteración hacia las regiones de mayor interés. El objetivo final es reducir el número de combinaciones de hiperparámetros con las que se evalúa el modelo, eligiendo únicamente los mejores candidatos. Esto significa que, la ventaja frente a las otras estrategias mencionadas, se maximiza cuando el espacio de búsqueda es muy amplio o la evaluación del modelo es muy lenta.

Los paquetes rBayesianOptimization y ParBayesianOptimization permiten combinar fácilmente el proceso de optimización bayesiana con el entrenamiento de de cualquier modelo.

MLR3 y ParBayesianOptimization


library(ParBayesianOptimization)
library(mlr3verse)



Datos


Se emplea los datos sobre el precio de viviendas disponible en el paquete AmesHoising.

library(AmesHousing)
library(tidyverse)

datos <- make_ames()
# Se emplean únicamente las variables:
variables <- c("Sale_Price", "Bldg_Type", "Neighborhood", "Year_Built",
               "Gr_Liv_Area", "Full_Bath", "Year_Sold", "Lot_Area", "Central_Air",
               "Heating", "Latitude", "Longitude")
datos <- datos %>%
         select(one_of(variables))
task_datos <- TaskRegr$new(
                id = "task_datos",
                backend = datos,
                target = "Sale_Price"
              )

set.seed(1234)
rsmp_holdout <- rsmp("holdout", ratio = 0.8)
rsmp_holdout$instantiate(task = task_datos)

# Índices train
id_train <- rsmp_holdout$train_set(i = 1)
# Índices test
id_test <- rsmp_holdout$test_set(i = 1)
task_train <- task_datos$clone()$filter(id_train)
task_test  <- task_datos$clone()$filter(id_test)



Predictores

# Se define la variable respuesta y los predictores.
var_respuesta <- "Sale_Price"
# Para este modelo se emplean todos los predictores disponibles.
predictores   <- task_train$feature_names



Función de scoring


scoring_function <- function(num.trees, min.node.size) {
 
  learner <- lrn("regr.ranger", num.threads = 4)
  learner$param_set$values <- list(num.trees = num.trees,
                                   min.node.size = min.node.size
                               )
  resampling_cv <- rsmp("cv", folds = 5)
  metrica       <- msr("regr.rmse")
  
  resampling <- resample(
                  task         = task_train,
                  learner      = learner,
                  resampling   = resampling_cv,
                  store_models = FALSE
               )

  score <- resampling$aggregate(measures = metrica)
  
  # Si la métrica es un error, dado que el proceso es de maximización, se convierte
  # en negativo.
  score <- -score
  
  return(list(Score = score))
}



Espacio de búsqueda


limites <- list( 
              # Emplear sufijo L si son valores enteros
              num.trees = c(100L, 2000L),
              min.node.size = c(1L, 100L)
           )



Optimización


set.seed(123)
bayes_opt <- bayesOpt(
                FUN = scoring_function, 
                bounds = limites,
                initPoints = 10, 
                iters.n = 20,
                iters.k = 1,
                otherHalting = list(timeLimit = 5*60)
            )
bayes_opt
## Class: bayesOpt
## 
##                  Epochs: 18
##              Iterations: 28
##     Average FUN Seconds: 7.86
##     Highest FUN Seconds: 19.08
## Final Upper Conf. Bound: 0.08171149
##              GP Updated: TRUE
##             Stop Status: Stopped Early. See $stopStatus



Resultados


bayes_opt$scoreSummary
getBestPars(bayes_opt)
## $num.trees
## [1] 1442
## 
## $min.node.size
## [1] 1
library(gganimate)
ggplot(data = bayes_opt$scoreSummary,
       aes(x = num.trees, y = min.node.size, size = Score, color = Score)) +
  geom_point() +
  scale_color_viridis_c() +
  theme_bw()



Paralelización


library(doParallel)

cl <- makeCluster(4)
registerDoParallel(cl)
clusterExport(cl, c('task_train'))
clusterEvalQ(cl, expr = {library(mlr3verse)})

set.seed(123)
bayes_opt <- bayesOpt(
                FUN = scoring_function, 
                bounds = limites,
                initPoints = 10, 
                iters.n = 20,
                iters.k = 4,
                parallel = TRUE,
                otherHalting = list(timeLimit = 10*60)
            )

stopCluster(cl)
registerDoSEQ()



H2O y rBayesianOptimization


library(rBayesianOptimization)
library(h2o)



Datos


Se emplea los datos sobre el precio de viviendas disponible en el paquete AmesHoising.

library(AmesHousing)
library(tidyverse)

datos <- make_ames()
# Se emplean únicamente las variables:
variables <- c("Sale_Price", "Bldg_Type", "Neighborhood", "Year_Built",
               "Gr_Liv_Area", "Full_Bath", "Year_Sold", "Lot_Area", "Central_Air",
               "Heating", "Latitude", "Longitude")
datos <- datos %>%
         select(one_of(variables))



Iniciación del cluster

# Creación de un cluster local con todos los cores disponibles.
h2o.init(ip = "localhost",
         # -1 indica que se empleen todos los cores disponibles.
         nthreads = -1,
         # Máxima memoria disponible para el cluster.
         max_mem_size = "4g")

# Se eliminan los datos del cluster por si ya había sido iniciado.
h2o.removeAll()
h2o.no_progress()
datos_h2o <- as.h2o(x = datos, destination_frame = "datos_h2o")
set.seed(123)
particiones     <- h2o.splitFrame(data = datos_h2o, ratios = c(0.8), seed = 123)
datos_train_h2o <- h2o.assign(data = particiones[[1]], key = "datos_train_h2o")
datos_test_h2o  <- h2o.assign(data = particiones[[2]], key = "datos_test_h2o")



Predictores

# Se define la variable respuesta y los predictores.
var_respuesta <- "Sale_Price"
# Para este modelo se emplean todos los predictores disponibles.
predictores <- setdiff(h2o.colnames(datos_train_h2o), var_respuesta)



Función de scoring


scoring_function <- function(learn_rate, max_depth, sample_rate) {
 
  nfolds <- 3
  metrica <- "RMSE"

  modelo <- h2o.gbm(
              # Variable respuesta y predictores
              y = var_respuesta,
              x = predictores,
              # Datos de entrenamiento
              training_frame = datos_train_h2o,
              # Preprocesado
              ignore_const_cols = TRUE,
              # Hiperparámetros
              ntrees = 5000,
              learn_rate  = learn_rate,
              max_depth   = max_depth,
              sample_rate = sample_rate,
              # Estrategia de validación
              seed   = 123,
              nfolds = nfolds,
              # Parada temprana
              stopping_rounds = 5,
              stopping_metric = metrica,
              stopping_tolerance  = 0.001,
              score_tree_interval = 100
          )
  
  score <- h2o.performance(model = modelo, xval = TRUE)@metrics[[metrica]]
  
  # Si la métrica es un error, dado que el proceso es de maximización, se convierte
  # en negativo.
  score <- -score
  return(list(Score = score))
}



Espacio de búsqueda


limites <- list( 
              # Emplear sufijo L si son valores enteros
              learn_rate  = c(0.01, 0.1),
              max_depth   = c(1L, 20L), 
              sample_rate = c(0.5, 1)
           )



Optimización


set.seed(123)
bayes_opt <- BayesianOptimization(
              FUN = scoring_function, 
              bounds = limites,
              init_points = 4,
              n_iter = 4,
              acq = "ucb",
              kappa = 2.576,
              eps = 0.0,
              verbose = TRUE
            )
## elapsed = 57.26  Round = 1   learn_rate = 0.0359 max_depth = 19.0000 sample_rate = 0.7757    Value = %#5.0-1e 
## elapsed = 23.25  Round = 2   learn_rate = 0.0809 max_depth = 2.0000  sample_rate = 0.7283    Value = %#5.0-1e 
## elapsed = 23.68  Round = 3   learn_rate = 0.0468 max_depth = 11.0000 sample_rate = 0.9784    Value = %#5.0-1e 
## elapsed = 29.73  Round = 4   learn_rate = 0.0895 max_depth = 18.0000 sample_rate = 0.7267    Value = %#5.0-1e 
## elapsed = 27.83  Round = 5   learn_rate = 0.0432 max_depth = 9.0000  sample_rate = 0.9486    Value = %#5.0-1e 
## elapsed = 32.18  Round = 6   learn_rate = 0.0100 max_depth = 1.0000  sample_rate = 0.7140    Value = %#5.0-1e 
## elapsed = 34.10  Round = 7   learn_rate = 0.0698 max_depth = 19.0000 sample_rate = 0.6088    Value = %#5.0-1e 
## elapsed = 42.21  Round = 8   learn_rate = 0.0789 max_depth = 20.0000 sample_rate = 0.6695    Value = %#5.0-1e 
## 
##  Best Parameters Found: 
## Round = 2    learn_rate = 0.0809 max_depth = 2.0000  sample_rate = 0.7283    Value = %#5.0-1e



Resultados


bayes_opt$History
bayes_opt$Best_Par
##  learn_rate   max_depth sample_rate 
##  0.08094746  2.00000000  0.72830737