Entendiendo Gradient Boosting Machines

Objetivos

La idea de este notebook es poder dar una intuición acerca de cómo funcionan los métodos basados en Gradient Boosting. Para ello, utilizaremos un dataset de juguete.

Gradient Descent en una carilla

Recordemos: ¿cuáles son los principales elementos de un algoritmo de Machine Learning?

• datos (\(X\) e \(y\))
• algoritmo o hipótesis (\(\hat{y}=f(X)\)) que contiene
• parámetros (por ejemplo, \(\beta_{0}, \beta_{1}, ..., \beta_{p}\) de un modelo de regresión lineal)
• métrica de costo o cost function, una medida que sirve para evaluar la performance del modelo \(\ell(y_{i} - \hat{y_{i}})\)

La función de costo debería servirnos para poder seleccionar cuál es el mejor set de parámetros. En líneas generales, queremos miniminzar la función de costo. Por eso, los problemas de machine learnin (y en buena medida, los de estadística) son problemas de optimización.

Pensemos en una regresión lineal simple

• tenemos nuestros datos \(x\) e \(y\)
• la hipótesis: \(\hat{y_{i}} = \beta_{0} + \beta_{1}X_{1}\)
• los parámetros: \((\beta_{0}, \beta_{1})\)
• función de costo: \(MSE = \frac{1}{N} \sum_{i=1}^N(y_{i} - \hat{y_{i}})^2\)

Una opción (poco recomendable) sería probar todos los valores posibles de \((\beta_{0},\beta_{1})\), digamos entre los infinitos negativos y los infitinos positivos, calcular las predicciones para cada uno, calcular \(MSE\) para cada uno y elegir el \(\beta_{1}\) que logra un valor mínimo en \(MSE\), o dicho en lenguaje matemático: \(\underset{\beta_{0},\beta_{1}}{\operatorname{argmax}}(\frac{1}{N} \sum_{i=1}^N(y_{i} - \hat{y_{i}})^2)\)

En este punto, el problema se tranforma en un problema de optimización. Dado que en general, no es computacionalmente posible probar todos los valores de un parámetro (y mucho menos de miles o millones) se recurren a otro tipo de métodos. El más utilizado en sus diversas variables es el llamado Gradient Descent o Descenso de gradiente.

Si bien, no entraremos en sus sutilezas, solamente diremos que se trata de un procedimiento iterativo que hace uso de las propiedades de las derivadas (algo así como el ratio de cambio de la función de costo al cambiar el valor del parámetro) para ir actualizando los valores de los parámetros hasta llegar al mínimo:

Calculando derivadas de la función de costo, llegamos a la siguente fórmula que permite ir actualizando los valores de \((\beta_{0}, \beta_{1})\) en un número fijo de iteraciones:

\(\beta_{0}^1 = \beta_{0}^0 - (y_{i} - \hat{y_{i}}) \times \gamma \times X\)

\(\beta_{1}^1 = \beta_{1}^0 - (y_{i} - \hat{y_{i}}) \times \gamma\)

Cuán pequeños (o grandes) son los cambios, dependerán del parámetro \(\gamma\). también llamada learning rate.

Aquí hay un tutorial sobre el tema, sumamente intuitivo.

Gradient Boosting Machines… ¿qué son?

Ahora bien, la idea general detrás de los Gradient Boosting Machines (GBM, a partir de ahora), es igual a la de los métodos de Boosting en general: entrenar modelos fuertes a partir de un ensamble de modelos débiles. Se hace, igualmente, entrenando clasificadores de manera secuencia, en dónde cada uno de los clasificadores posteriores se centra en los errores del clasificador anterior.

La principal diferencia con AdaBoost es que, en lugar de jugar con los pesos de cada instancia en cada iteración, GBM trata de fitear un modelo sobre los residuos del modelo anterior.

Así, se trata de modelos de tipo aditivo. La idea general es que la predicción final es generada por la suma de diversos modelos parciales:

\(\hat{y} = h_{1}(X) + h_{2}(X) + h_{3}(X) + ... + h_{m}(X)\)

\(\hat{y} = \sum_{m=1}^M h_{m}(X)\)

Puede pensarse de forma análoga a la descomposición de una serie temporal, donde la series es pensada como una composición aditiva de una componente de tendencia \(T_{t}\), una de ciclo \(C_{y}\), una estacional \(S_{t}\) y otra de ruido \(\mu_{t}\):

\(y_{t} = T_{t} + C_{t} + S_{t} + \mu_{t}\)

Para seguir profundizando aquí hay un tutorial de carácter sumamente intuitivo. También pueden consultar el libro de Aurelien Geron del cual tomamos y tradujimos a R este ejemplo, escrito originalmente en Python.

Un ejemplo para construir la intuición

Veamos un ejemplo usando árboles de decisión, recordando que se trata de un meta-algoritmo por lo cual podríamos usar (casi) cualquier otro método.

Generamos nuestros datos

library(tidyverse)
library(rpart)
library(caret)
set.seed(42)
X <- runif(100, 0, 1) - 0.5
y <- 3*X**2 + 0.05 * runif(100,0,4)
df <- cbind(X,y) %>% as_tibble()
rm(X,y)

En primer lugar, entrenemos un árbol de decisión sobre nuestro dataset y agreguemos las predicciones como columnas en el dataset:

tree_1 <- rpart(y~X, data=df, method='anova', control=list(cp=0.0000001, maxdepth=2))                
df <- df %>% mutate(h_1 = predict(tree_1, df),
                p_1 = h_1,
                y_1 = y - p_1)

Ahora, entrenemos un segundo árbol pero sobre los residuos del anterior, es decir, sobre \(y_{i} - \hat{y_{i}}\) o más precisamente: y - predict(tree_1, df) y agreguemos, una vez más, los resultados al dataset:

tree_2 <- rpart(y_1~X, data=df, method='anova', control=list(cp=0.0000001, maxdepth=2))                
df <- df %>% mutate(h_2 = predict(tree_2, df),
                    p_2 = p_1 + h_2,
                    y_2 = y - p_2)

Repitamos el proceso una tercera vez… y una cuarta y una quinta y una sexta…

tree_3 <- rpart(y_2~X, data=df, method='anova', control=list(cp=0.0000001, maxdepth=2))                
df <- df %>% mutate(h_3 = predict(tree_3, df),
                    p_3 = p_2 + h_3,
                    y_3 = y - p_3)
tree_4 <- rpart(y_3~X, data=df, method='anova', control=list(cp=0.0000001, maxdepth=2))                
df <- df %>% mutate(h_4 = predict(tree_4, df),
                    p_4 = p_3 + h_4,
                    y_4 = y - p_4)
tree_5 <- rpart(y_4~X, data=df, method='anova', control=list(cp=0.0000001, maxdepth=2))                
df <- df %>% mutate(h_5 = predict(tree_5, df),
                    p_5 = p_4 + h_5,
                    y_5 = y - p_5)
tree_6 <- rpart(y_5~X, data=df, method='anova', control=list(cp=0.0000001, maxdepth=2))                
df <- df %>% mutate(h_6 = predict(tree_6, df),
                    p_6 = p_5 + h_6,
                    y_6 = y - p_6)

Ahora tenemos un ensable que consiste en seis árboles de decisión. Podemos hacer una predicción sobre una instancia nueva simplemente agregando las predicciones de nuestros seis árboles.

ggpubr::ggarrange(ncol=2, nrow=3,
                  ggplot(df) + 
                          geom_point(aes(x=X, y=y), color='blue') + 
                          geom_line(aes(x=X, y=h_1), color='green'),
                  ggplot(df) + 
                          geom_point(aes(x=X, y=y), color='blue') + 
                          geom_line(aes(x=X, y=p_1), color='red'),
                  ggplot(df) + 
                          geom_point(aes(x=X, y=y_1), color='blue') + 
                          geom_line(aes(x=X, y=h_2), color='green') +
                          scale_y_continuous(limits=c(-0.4,1)),
                  ggplot(df) + 
                          geom_point(aes(x=X, y=y), color='blue') + 
                          geom_line(aes(x=X, y=p_2), color='red'),
                  ggplot(df) + 
                          geom_point(aes(x=X, y=y_2), color='blue') + 
                          geom_line(aes(x=X, y=h_3), color='green') +
                          scale_y_continuous(limits=c(-0.4,1)),
                  ggplot(df) + 
                          geom_point(aes(x=X, y=y), color='blue') + 
                          geom_line(aes(x=X, y=p_3), color='red')
)

El plot muestra cómo va evolucionando el modelo a medida que vamos agregando una mayor cantidad de iteraciones. A la izquierda tenemos el ajuste de cada modelo con su propia variable dependiente (es decir, los residuos del modelo anterior) y a la derecha, el ajuste del ensamble global sobre \(X\) e \(y\).

Puede verse cómo el ensamble se va haciendo cada vez más preciso a medida que vamos agregando modelos.

¿Cuáles son los hiperparámetros de GBM?

Vamos a ver algunos parámetros de GBM en caret y su interpretación. El método se llama “GBM” y se usa, como siempre, en la función train.

Vemos que en tuneGrid hay cuatro parámetros para evaluar:

• n.trees: es la cantidad de modelos a entrenar; es decir, la cantidad de iteraciones
• interaction.depth: la complejidad de cada árbol de decisión
• shrinkage: también llamada learning rate en otros lenguajes e implementaciones, y puede interpretarse como la contribución de cada árbol al modelo final… al setear un shrinkage bajo, por ejemplo, shrinkage=0.1 será necesario una mayor cantidad de árboles para entrenar, pero los resultados tenderán a ser mejores.
• n.minobsinnode: el tamaño mínimo para poder hacer un split en cada nodo

Entrenaremos dos GBM, uno con pocos árboles y learing rate alto y otro con las características inversas:

gbm_slow_05 <- train(y~X, method='gbm',
             data=df,
             trControl=trainControl,
             tuneGrid=data.frame(n.trees=200,
                                 interaction.depth=3,
                                 shrinkage=0.05,
                                 n.minobsinnode=2)
             )

Iter   TrainDeviance   ValidDeviance   StepSize   Improve
     1        0.0578             nan     0.0500    0.0049
     2        0.0534             nan     0.0500    0.0042
     3        0.0498             nan     0.0500    0.0043
     4        0.0464             nan     0.0500    0.0040
     5        0.0429             nan     0.0500    0.0033
     6        0.0396             nan     0.0500    0.0026
     7        0.0368             nan     0.0500    0.0031
     8        0.0339             nan     0.0500    0.0026
     9        0.0317             nan     0.0500    0.0023
    10        0.0292             nan     0.0500    0.0022
    20        0.0141             nan     0.0500    0.0011
    40        0.0048             nan     0.0500    0.0002
    60        0.0028             nan     0.0500   -0.0000
    80        0.0022             nan     0.0500   -0.0000
   100        0.0020             nan     0.0500   -0.0000
   120        0.0019             nan     0.0500   -0.0000
   140        0.0018             nan     0.0500   -0.0000
   160        0.0017             nan     0.0500   -0.0000
   180        0.0016             nan     0.0500   -0.0000
   200        0.0015             nan     0.0500   -0.0000

Agregamos las predicciones al dataset y ploteamos:

df <- df %>% mutate(y_gbm=predict(gbm, df),
                    y_gbm_n3_s01=predict(gbm_n3_s01, df),
                    y_gbm_n200_s1=predict(gbm_n200_s1, df),
                    y_gbm_slow=predict(gbm_slow, df),
                    y_gbm_slow_05=predict(gbm_slow_05, df),
                    y_gbm_slow_001=predict(gbm_slow_001, df)
                    )
ggpubr::ggarrange(ncol=2, nrow=3,
ggplot(df) + 
        geom_point(aes(x=X, y=y), color='blue') + 
        geom_line(aes(x=X,y=y_gbm)) +
        labs(title='n.trees=3, shrinkage=1'),
ggplot(df) + 
        geom_point(aes(x=X, y=y), color='blue') + 
        geom_line(aes(x=X,y=y_gbm_n3_s01)) +
        labs(title='n.trees=3, shrinkage=0.1'),
ggplot(df) + 
        geom_point(aes(x=X, y=y), color='blue') + 
        geom_line(aes(x=X,y=y_gbm_n200_s1)) +
        labs(title='n.trees=200, shrinkage=1'),
ggplot(df) + 
        geom_point(aes(x=X, y=y), color='blue') + 
        geom_line(aes(x=X,y=y_gbm_slow)) +
        labs(title='n.trees=200, shrinkage=0.1'),
ggplot(df) + 
        geom_point(aes(x=X, y=y), color='blue') + 
        geom_line(aes(x=X,y=y_gbm_slow_05)) +
        labs(title='n.trees=200, shrinkage=0.05'),
ggplot(df) + 
        geom_point(aes(x=X, y=y), color='blue') + 
        geom_line(aes(x=X,y=y_gbm_slow_001)) +
        labs(title='n.trees=200, shrinkage=0.001')
)

¿Cómo describir a cada uno? El de la izquierda parece no ser lo suficientemente flexible, mientras que el de la derecha parece generar overfitting de forma bastante clara.