Gradiente descendente

Laboratório de Estatística e Geoinformação - LEG/UFPR

Gradiente descendente é um dos algoritmos de maior sucesso em problemas de Machine Learning. O método consiste em encontrar, de forma iterativa, os valores dos parâmetros que minimizam determinada função de interesse. Para entender esse método, vamos relembrar, inicialmente, como é a obtenção dos parâmetros de forma analítica. Considere um problema de regressão linear em que \(\mathbf{X} \in \mathrm{M}_{n \times p}(\mathrm{R})\), com \(n > p\) e \(posto(\mathbf{X}) = p\). Seja \(Y_i\), \(i=1,\dots,n\), dado por: \[ Y_i = \mathbf{\beta}^t\mathbf{x_i} + \varepsilon_i, \ \ \textrm{com}\ \mathbf{\varepsilon} \sim N(0,\sigma^2).\nonumber \]

Nosso interesse é encontrar os \(\beta's\) que minimize a função \(J\) \[\begin{eqnarray} J(\mathbf{\beta}) &=& \sum_{i=1}^{n} (y_i - \mathbf{\beta}^t\mathbf{x}_i)^2 \nonumber \\ &=& \langle \mathbf{X}\mathbf{\beta} - \mathbf{y},\mathbf{X}\mathbf{\beta} - \mathbf{y} \rangle \nonumber \end{eqnarray}\]

Nesse caso, a solução por mínimos quadrados é dada por: \[ \mathbf{\hat \beta} = \left(\mathbf{X}^t\mathbf{X}\right)^{-1}\mathbf{X}^t \mathbf{y}. \]

Teorema: Seja \(\mathbf{X} \in \mathrm{M}_{n \times p}(\mathrm{R})\), com \(n > p\) e \(posto(\mathbf{X}) = p\). Definimos \(J: \mathrm{R}^{p} \rightarrow \mathrm{R}\) da seguinte forma: \[ J(\mathbf{\mathbf{\beta}}) = \langle \mathbf{X}\mathbf{\beta} - \mathbf{y},\mathbf{X}\mathbf{\beta} - \mathbf{y} \rangle \ ; \ \mathbf{\beta} \in \mathrm{R}^{p}. \] Então, o problema de minimização: encontrar \(\hat \beta \in \mathrm{R}^{p}\) tal que \[ J(\mathbf{\hat\beta}) = min\left\{J(\mathbf{\beta}) \ ; \ \mathbf{\beta} \in \mathrm{R}^{p} \right\} \] é equivalente ao Sistema Normal \[ \mathbf{X}^t\mathbf{X}\mathbf{\beta} = \mathbf{X}^t \mathbf{y}, \] que resulta em \[ \mathbf{\hat \beta} = \left(\mathbf{X}^t\mathbf{X}\right)^{-1}\mathbf{X}^t \mathbf{y}. \]

Demonstração Para encontrar o mínimo de \(J(\beta)\) devemos calcular suas derivadas parciais: \[\begin{eqnarray} \nabla J(\mathbf{\beta}) &=& \nabla \left(\mathbf{X}\mathbf{\beta} - \mathbf{y})^t (\mathbf{X}\mathbf{\beta} - \mathbf{y}\right)\nonumber\\ &=& \nabla \left(\mathbf{\beta}^t\mathbf{X}^t\mathbf{X}\mathbf{\beta} - \mathbf{\beta}^t\mathbf{X}^t\mathbf{X}\mathbf{y} - \mathbf{y}^t\mathbf{X}\mathbf{\beta} + \mathbf{y}^t\mathbf{y} \right)\nonumber\\ &=& \nabla tr \left(\mathbf{\beta}^t\mathbf{X}^t\mathbf{X}\mathbf{\beta} - \mathbf{\beta}^t\mathbf{X}^t\mathbf{y} - \mathbf{y}^t\mathbf{X}\mathbf{\beta} + \mathbf{y}^t\mathbf{y} \right)\nonumber\\ &=& \nabla \left(tr \mathbf{\beta}^t\mathbf{X}^t\mathbf{X}\mathbf{\beta} - 2 tr \mathbf{y}^t\mathbf{X}\mathbf{\beta} \right)\nonumber\\ &=& \nabla \left(\mathbf{X}^t \mathbf{X}\mathbf{\beta} + \mathbf{X}^t \mathbf{X}\mathbf{\beta} - 2\mathbf{X}^t \mathbf{y} \right)\nonumber\\ &=& 2 \mathbf{X}^t \mathbf{X}\mathbf{\beta} - 2 \mathbf{X}^t\mathbf{y}\nonumber \end{eqnarray}\] Dizemos que \(\mathbf{\hat \beta}\) é um ponto crítico de \(J\) se, e somente se, \[\begin{eqnarray} \nabla J(\mathbf{\beta})(v) = 2 \left \langle \mathbf{X}^t\mathbf{X} \mathbf{\hat \beta} - \mathbf{X}^t \mathbf{y} , \mathbf{v} \right \rangle = 0 \ , \textrm{para todo } \mathbf{v} \in \mathrm{R}^{p} \nonumber \end{eqnarray}\] \(\nabla J(\mathbf{\beta})(v)\) é derivada direcional de \(J\) no ponto \(\mathbf{\hat \beta}\) na direção de \(\mathbf{v} \in \mathrm{R}^{p}\). Portanto, o único ponto crítico do funcional \(J\) é caracterizado como: \[ \mathbf{\hat \beta} = \left(\mathbf{X}^t\mathbf{X}\right)^{-1}\mathbf{X}^t \mathbf{y}. \]
\(\blacksquare\)

Gradiente Descendente (GD)


No caso descrito, obtivemos \(\mathbf{\hat \beta}\) analiticamente. Entretanto, com frequência nos deparamos com situações do tipo: dificuldade em se obter a solução do sistema na forma fechada (ou ela nem existir), quando \(n\) é grande (p. ex., \(>100000\)) o cálculo da inversa ser muito caro computacionalmente etc. O Gradiente descendente (GD) minimiza esses (e outros problemas). O método obtém a solução do sistema de forma iterativa, simples e barata. Além de servir de base para vários algoritmos de segunda ordem que aceleram essa convergência (Método de Newton, Método do Gradiente Conjugado etc.). Em termos do problema de minimização da função \(J(\mathbf{\beta})\), cada passo é obtido da forma \[ \mathbf{\beta^{(k+1)}} = \underset{\mathbf{\beta}}{argmin}\ J(\mathbf{\beta}^{(k)}) + \nabla J(\mathbf{\beta}^{(k)})^t \left(\mathbf{\beta} - \mathbf{\beta}^{(k)}\right) + \frac{1}{2\alpha}\left\lVert{\mathbf{\beta} - \mathbf{\beta^{(k)}}}\right\rVert^2 . \]

Trata-se de uma aproximação quadrática ao trocar \(\frac{1}{2\alpha}I\) por \(\nabla^2 J(\mathbf{\beta}^{(k)})\). Note que não se fez qualquer menção à função \(J\) (assumimos apenas que a primeira derivada exista). Ou seja, tem-se uma garantia de convergência para um mínimo local, independente da complexidade do modelo. Derivando com relação a \(\mathbf{\beta}\) e igualando a zero, temos: \[ 0 = \nabla J(\mathbf{\beta}^{(k)}) +\frac{1}{\alpha}(\mathbf{\beta} - \mathbf{\beta}^{(k)}), \] e chega-se em: \[ \mathbf{\beta}^{(k+1)} = \mathbf{\beta}^{(k)} - \alpha \nabla J(\mathbf{\beta}^{(k)}). \]

Algoritmo: Escolha um chute inicial, \(\mathbf{\beta}^{(0)} \in \mathrm{R}^{p}\), repita: \[ \mathbf{\beta}^{(k+1)} = \mathbf{\beta}^{(k)} - \alpha_k \nabla J(\mathbf{\beta}^{(k)}), \ k= 0,1, \dots \] pare quando atingir convergência.

Em resumo, os \(\beta's\) são atualizados de acordo com o negativo do gradiente. O tamanho do passo em cada iteração é determinados de acordo com a taxa de aprendizagem, \(\alpha>0\). Note que

  1. Se \(\nabla J(\mathbf{\beta}^{(k)}) < 0\), temos \(\mathbf{\beta}^{(k+1)} > \mathbf{\beta}^{(k)}\), pois o mínimo está à direita;
  2. Se \(\nabla J(\mathbf{\beta}^{(k)}) > 0\), temos \(\mathbf{\beta}^{(k+1)} < \mathbf{\beta}^{(k)}\), pois o mínimo está à esquerda.
Abaixo, temos a reta de uma regressão linear simples obtida iterativamente. As linhas azuis representam cada passo do algoritmo, e a vermelha a reta resultante.



 Figura 1: Erro de predição segundo a flexibilidade do modelo.

Exemplo simulado


Considere os dados simulados a partir do modelo Michaelis Menten: \[ f(x_{i}, \theta) = \frac{\theta_a x_i }{\theta_v + x_i}, \quad \theta = (\theta_a, \theta_v)'. \]

A linha azul representa a função \(f\) quando (\(\theta_a, \theta_v\)) = (10, 2).

# Carregando pacotes.
rm(list = objects())
library(lattice)
library(latticeExtra)

#-----------------------------------------------------------------------
# Simulando dados do modelo Michaelis Menten.

# Retorna os valores determinados pelo modelo para a média de Y|x.
eta <- function(x, theta) {
    theta[1] * x/(theta[2] + x)
}

# Valores para os parâmetros.
pars <- c(theta_a = 10,
          theta_v = 2)

# Retorna um rúido normal padrão.
epsilon <- function(x, ...) {
    rnorm(n = length(x), mean = 0, ...)
}

# Simulando do modelo.
set.seed(321)
da <- data.frame(x = sort(runif(n = 50, min = 3, max = 20)))
da$y <- eta(x = da$x, theta = pars) + epsilon(x = da$x, sd = 0.1)

plot(y ~ x, data = da)
curve(eta(x, theta = pars),
      add = TRUE,lwd=2,
      col = "dodgerblue3")

O nosso interesse será encontrar tais valores através do algoritmo GD, minimizando a função objetivo: \[ L(\theta) = \sum_{i = 1}^{n} [y_{i} - g(x_{i}, \theta)]^2. \]

Por simplicidade, vamos considerar \(g(x_{i}, \theta) =\frac{\theta_a x_i }{\theta_v + x_i}\), pois queremos apenas elucidar a utilização do método. No gráfico da esquerda temos as curvas de nível do \(log\) da função objetivo (aplicou-se essa transformação para facilitar a visualização), com destaque para o ponto que desejamos encontrar numericamente. Na direita, o comportamento do \(log(rss)\) de acordo com os parâmetros.

# Função perda (custo ou objetivo).
loss <- function(theta, y, x) {
    sum((y - eta(x, theta = theta))^2)
}

# Grid de valores dos parâmetros.
grid <- expand.grid(theta_a = seq(0, 20, length.out = 41),
                    theta_v = seq(-3, 8, length.out = 41))

# Avalia a função objetivo em cada ponto do grid.

grid$rss <- mapply(theta_a = grid$theta_a,
                   theta_v = grid$theta_v,
                   FUN = function(theta_a, theta_v) {
                       loss(theta = c(theta_a, theta_v),
                            y = da$y,
                            x = da$x)
                   })

# Define escala de cores para visualização.
colr <- colorRampPalette(rev(brewer.pal(11, "Spectral")),
                         space = "rgb")

# Exibe a superfície da função objetivo (com log para mais detalhes).

plot1 <- levelplot(log(rss) ~ theta_a + theta_v,
          data = grid,
          col.regions = colr,
          contour = TRUE) +
    layer(panel.abline(v = pars["theta_a"],
                       h = pars["theta_v"],
                       lty = 2))

plot2 <- wireframe(log(rss) ~ theta_a + theta_v,
          data = grid,
          drape = TRUE,
          col.regions = colr(101),
          #panel.3d.wireframe = wzRfun::panel.3d.contour,
          type = c("bottom"))

grid.arrange(plot1,plot2, ncol=2)

Aplicamos o método do Gradiente Descendente com critérios de parada: 10000 iterações ou \(|\mathbf{\theta}^{(k+1)} - \mathbf{\theta}^{(k)}|< 0.0001\) e taxa de aprendizado \(\alpha = 0,001\). A figura abaixo ilustra os passos do algoritmo em cada iteração até a convergência para os pontos \((\theta_a,\theta_v) = (10,03; 2,05)\). Satisfazendo o critério de parada na iteração 2412.

#-----------------------------------------------------------------------
# Gradiente.

# Derivadas analíticas com o R.
# loss_expr <- expression((y - theta_a * x/(theta_v + x))^2)
# D(expr = loss_expr, name = "theta_a")
# D(expr = loss_expr, name = "theta_v")

# Função gradiente.
grad <- function(theta, y, x) {
    theta_a <- theta[1]
    theta_v <- theta[2]
    part <-
        cbind(partial_theta_a =
                  -(2 * (x/(theta_v + x) * (y - theta_a * x/(theta_v + x)))),
              partial_theta_v =
                  2 * (theta_a * x/(theta_v + x)^2 * (y - theta_a * x/(theta_v + x))))
    # Derivada da soma = soma das derivadas.
    colSums(part)
}

# Avaliando a função objetivo.
# loss(theta = c(10, 2), y = da$y, x = da$x)

# Avaliando a função gradiente.
# grad(theta = c(10, 2),
#     y = da$y,
#     x = da$x)

#-----------------------------------------------------------------------
# Algortimo de Gradiente Descendente.

# Tolerância: parar quando a mudança for irrelevante ou exceder número
# máximo de iterações.
tol <- 1E-5
niter <- 10000

# Define taxa de aprendizado (parâmetro de tunning).
alpha <- 0.001

# Cria matriz vazia para ir preenchendo.
th <- matrix(0, ncol = 2, nrow = niter)

# Inclui valor inicial e dispara o loop.
k <- 1
th[k, ] <- c(theta_a = 18, theta_v = -1)
repeat {
    delta <- alpha * grad(theta = th[k, ],
                          y = da$y,
                          x = da$x)
    th[k + 1, ] <- th[k, ] - delta
    # print(th[k, ])
    if (k > niter | all(abs(delta) < tol)) break
    k <- k + 1
}

th <- th[1:k, ]

#-----------------------------------------------------------------------
# Traço do algoritmo.

lth1 <- apply(th,
             MARGIN = 1,
             FUN = loss,
             y = da$y,
             x = da$x)

levelplot(rss ~ theta_a + theta_v,
          data = grid,
          col.regions = colr,
          contour = TRUE) +
    layer(panel.abline(v = pars["theta_a"],
                       h = pars["theta_v"],
                       lty = 2)) +
    layer(panel.points(x = th[, 1],
                       y = th[, 2],
                       type = "o",
                       cex = 0.1))


Gradiente Descendente Estocástico


O método do Gradiente Descendente apresentado acima (também chamado Gradiente Descendente Batch) utiliza – em cada iteração – toda a amostra de treino. Trata-se de um algoritmo útil para matrizes bem condicionadas e problemas fortemente convexos. Entretanto, frequentemente traz problemas, pois problemas interessantes não são fortemente convexos ou bem condicionados. Além disso, tal abordagem apresenta outros aspectos que merecem ser ressaltados:

  1. Pode ser muito custoso computacionalmente, pois para cada atualização do conjunto de treinamento, precisamos ensinar o algoritmo desde o início. Inviabilizando, inclusive, qualquer método online (expediente em que os exemplos chegam em um fluxo contínuo).

  2. Ademais, considerando que o erro padrão da média estimada é dado por \(\sigma/\sqrt{n}\), quando aumentamos o número de exemplos, temos um retorno (ganho) menor que o linear. Por exemplo, compare as estimativas advindas de uma amostra de tamanho 100 e outra de tamanho 10000. O segundo caso requer 100 vezes mais computação, mas reduz o erro padrão da média em apenas 10 vezes.

  3. Outro aspecto é com relação à observações redundantes. Muitas vezes o conjunto completo de treinamento apresenta exemplos que pouco contribuem para atualização dos parâmetros.

O Gradiente Descendente Estocástico (Stochastic Gradient Descent) minimiza esses problemas, pois utiliza em cada iteração apenas uma observação. Considere o par \((\mathbf{x}_i,y_i)\) amostrado do treinamento, a atualização dos parâmetros é dado como segue.

Algoritmo: Escolha um chute inicial, \(\mathbf{\beta}^{(0)} \in \mathrm{R}^{p}\), repita: \[ \mathbf{\beta}^{(k+1)} = \mathbf{\beta}^{(k)} - \alpha_k \nabla J(\mathbf{\beta}^{(k)}; \mathbf{x}_i,y_i), \ k= 0,1, \dots \] pare quando atingir convergência.

Um meio termo entre o Batch e Estocástico é o Gradiente Descendente Mini-batch, no qual em cada passo obtêm-se uma amostra do conjunto treino. Neste caso, o tamanho da amostra (chamado mini-batch size) será mais um parâmetro de estudo (essa é uma desvantagem do método), tipicamente, utiliza-se valores entre 2 e 100. A vantagem de se usar tal abordagem sobre o Batch é a velocidade em se atualizar os parâmetros. Já sobre o Estocástico, está na implementação vetorizada (paralelizando os cálculos).


Exemplo simulado (continuação)


Aplicamos o método do Gradiente Descendente Mini-Batch ao exemplo apresentado anteriormente, com mini-batch size = 3. Pela figura abaixo, pode-se observar que o algoritmo convergiu para o mesmo ponto.

#-----------------------------------------------------------------------
# Algortimo de Gradiente Descendente Estocástico.

# Para quando exceder número máximo de iterações.
niter <- 2500

# Define taxa de aprendizado (parâmetro de tunning).
alpha <- 0.01

# Número de elementos amostrados do conjunto de treinamento.
index <- 1:nrow(da)
nele <- 3

# Cria matriz vazia para ir preenchendo.
th <- matrix(0, ncol = 2, nrow = niter)

# Inclui valor inicial e dispara o loop.
k <- 1
th[k, ] <- c(theta_a = 18, theta_v = -1)
while (k < niter) {
    db <- da[sample(index, size = nele), ]
    delta <- alpha * grad(theta = th[k, ],
                          y = db$y,
                          x = db$x)
    th[k + 1, ] <- th[k, ] - delta
    # print(th[k, ])
    k <- k + 1
}

#-----------------------------------------------------------------------
# O traço do algorítmo.

th <- th[1:k, ]
levelplot(rss ~ theta_a + theta_v,
          data = grid,
          col.regions = colr,
          contour = TRUE) +
    layer(panel.abline(v = pars["theta_a"],
                       h = pars["theta_v"],
                       lty = 2)) +
    layer(panel.points(x = th[, 1],
                       y = th[, 2],
                       type = "l"))

Mas, como era de se esperar, o método Mini-Batch utiliza mais iterações para atingir tal convergência. É digno de nota que, devido à trajetória estocástica do algoritmo, não utilizamos critérios de parada baseados em diferenças de valores consecutivos dos parâmetros, pois isso acarretaria paradas antecipadas (por exemplo, veja o comportamento das 500 últimas iterações). Dessa forma, é comum executar o algoritmo até exceder o número máximo de iterações, neste caso 2500.

# A descida do Estocástico
lth2 <- apply(th,
             MARGIN = 1,
             FUN = loss,
             y = da$y,
             x = da$x)

layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE))

plot(log(lth1), type = "o", cex = 0.3, main = "Batch", 
     xlab = "Iteração", ylab = "log(rss)", col = "dodgerblue3")
plot(log(lth2), type = "o", cex = 0.3, main = "Mini-Batch", 
     xlab = "Iteração", ylab = "log(rss)",col = "dodgerblue3")

# Vendo apenas as últimas iterações.
plot((tail(lth2, n = 500)), type = "o", cex = 0.3, main = "Últimas 500 iterações GD Mini-Batch", 
     xlab = "Iteração", ylab = "log(rss)",col = "dodgerblue3")


Escolha da taxa de aprendizagem


A taxa de aprendizagem, \(\alpha\), é um parâmetro crucial, especialmente no Gradiente Descendente Mini-Batch (e Estocástico). Pode-se considerá-la fixa, mas idealmente a diminuímos ao longo do tempo. Isso ocorre, porque estamos introduzindo uma fonte de variabilidade (devida à amostragem) que não desaparece mesmo quando nos aproximamos do mínimo. Se considerarmos \(\alpha\) pequeno, o método fica lento e muitas vezes não consegue escapar de um mínimo local. Se \(\alpha\) for muito grande o método pode divergir. Atualmente, existem algoritmos que atualizam \(\alpha\) em cada iteração (p. ex., line search). Assim, pode-se começar com valores maiores, e a medida que se aproxima do mínimo, diminui-se o tamanho do passo.

 Figura 1: Erro de predição segundo a flexibilidade do modelo.

Gradiente Descendente Boosting


Vimos no capítulo sobre Métodos baseados em árvores que o princípio básico do Boosting é propor um modelo básico (weak learner) e o aprimorar (ou “ensinar”) em cada iteração. O processo consiste em filtrar os resultados corretos, e concentrar-se naqueles que o modelo não soube lidar.

Drawing


Considere que estamos interessados em explicar \(Y\) através da função \(h(x)\), ou seja: \[\begin{eqnarray} \label{eq1} Y = h(x) + \textrm{resíduo} \end{eqnarray}\] Se o \(\textrm{resíduo}\) não for um ruído branco (mas algo relacionado com \(Y\)), podemos obter novas informações, além daquelas obtidas com o primeiro modelo, da seguinte forma: \[\begin{eqnarray} \label{eq2} \textrm{resíduo} = g(x) + \textrm{resíduo2} \end{eqnarray}\] Note que \(g(x)\) representa uma função não contemplada em \(h(x)\), que ajuda a explicar \(Y\). Podendo ser, inclusive, uma função que corrige possíveis erros de \(h(x)\). Combinando (\ref{eq1}) e (\ref{eq2}) \[\begin{eqnarray} Y = h(x) + g(x) + \textrm{resíduo2} \nonumber \end{eqnarray}\]

Pode-se dizer que \(h(x)\) foi atualizada com uma parte do \(\textrm{resíduo}\), ou seja \[ h(x)^{(2)} = h(x)^{(1)} + g(x) \]

Na ocasião, aprendemos sobre o AdaBoost, em que o resíduo decorrente do modelo anterior era utilizado na construção das novas árvores, como ilustrado na figura abaixo.

#-----------------------------------------------------------------------
# Gerando dados.

# Retorna os valores determinados pelo modelo.
eta <- function(x, beta = c(10, 0.05, 2)) {
    beta[1] + beta[2] * x + beta[3] * sin(x/10)
}

# Simulando do modelo.
set.seed(567)
da <- data.frame(x = sort(runif(200, 0, 100)))
da$y <- eta(da$x) + epsilon(da$x, sd = 1)

# A função h(x) será uma função degrau.

# Modelo inicial com predito sendo a média das observações.
m0 <- lm(y ~ 1, data = da)

# Função degrau. Retorna ajustado e resíduo.
h <- function(beta, y, x) {
    dico <- ifelse(x < beta, 0, 1)
    fit <- ave(y, dico, FUN = mean)
    list(dico = dico, fit = fit, res = y - fit)
}

# Função custo.
loss <- function(beta, y, x) {
    fit <- h(beta, y, x)$fit
    sum((y - fit)^2)
}

# Verifica o comportamento da função.
# bseq <- seq(0, 100, 1)
# lseq <- sapply(bseq, FUN = loss, y = da$y, x = da$x)
# plot(lseq ~ bseq, type = "o")

# Usando a optim().
# optim(par = c(50),
#       fn = loss,
#       y = da$y,
#       x = da$x)$par

#-----------------------------------------------------------------------
# O algorítmo.

# Cria uma lista para armazenar os artefatos produzidos.
n <- nrow(da)
k <- 12
resul <- replicate(k,
                   data.frame(fit = numeric(n),
                              res = numeric(n)),
                   simplify = FALSE)

# Inicia a lista com resultados do primeiro ajuste.
resul[[1]]$fit <- fitted(m0)
resul[[1]]$res <- residuals(m0)
bval <- c(NA, numeric(k - 1))

# Executa o gradiente boosting.
for (i in 2:k) {
    # Encontra o beta minimizando a função perda.
    beta <- optimize(f = loss,
                     y = resul[[i - 1]]$res,
                     x = da$x,
                     interval = c(0, 100))$minimum
    # Exibe o progresso.
    cat(i, "\t", beta, "\n")
    bval[i] <- beta
    # Extrai valores ajustados e resíduos.
    hfit <- h(beta,
              y = resul[[i - 1]]$res,
              x = da$x)
    # Atualiza os valores ajustados e resíduos.
    resul[[i]]$fit <- resul[[i - 1]]$fit + hfit$fit
    resul[[i]]$res <- da$y - resul[[i]]$fit
}

#-----------------------------------------------------------------------
# Exibe os resultados.

# Lista replicando a mesma fórmula.
form <- replicate(k, y ~ x)
names(form) <- 1:k

xyplot.list(form,
            data = da,
            type = "n",
            as.table = TRUE) +
    layer({
        i <- which.packet()
        fit <- resul[[i]]$fit
        panel.points(x = x,
                     y = y,
                     col = "gray")
        panel.lines(da$x,
                    fit,
                    col = "dodgerblue3",
                    lwd = 3)
        panel.abline(v = bval[i], lty = 2)
    })

Mas agora, como este método está relacionado com Gradiente Descendente? Note que durante o treinamento do modelo, estamos interessados em minimizar uma função custo, por exemplo, \[ J(y_i,h(\mathbf{x})) = \frac{1}{2n}\sum_{i=1}^{n}[y_i - h(\mathbf{x}_i)]^2 \]

Mas, derivando com relação a \(h(\mathbf{x}_i)\) temos \[ \frac{\partial J(\mathbf{y},h(\mathbf{x}))}{\partial h(\mathbf{x}_i)} = h(\mathbf{x}_i) - y_i. \]

Dessa forma, concluímos que algoritmo utiliza o negativo do gradiente como informação para os resíduos: \[\begin{eqnarray} residuos &=& y_i - h(\mathbf{x}_i) \nonumber\\ &=& -\frac{\partial J(y_i,h(\mathbf{x}))}{\partial h(\mathbf{x}_i)} \nonumber \end{eqnarray}\]
\[\begin{eqnarray} \textrm{resíduo} &\Leftrightarrow& \textrm{negativo do gradiente} \nonumber\\ \textrm{Atualizar} \ \ h(\mathbf{x}_i) \ \ \textrm{com o resíduo} &\Leftrightarrow& \textrm{Atualizar} \ \ h(\mathbf{x}_i) \ \ \textrm{com o negativo do gradiente} \nonumber \end{eqnarray}\]

Algoritmo: Escolha um chute inicial, \(h(\mathbf{x}_i)^{(0)}\), faça:

  1. Calcule \(\displaystyle -\frac{\partial J(y_i, h(\mathbf{x})^{(k)})}{\partial h(\mathbf{x}_i)^{(k)}}\)
  2. Ajuste um modelo de regressão \(g(\mathbf{x}_i)^{(k)}\) baseado no negativo do gradiente; \[ h(\mathbf{x}_i)^{(k+1)} = h(\mathbf{x}_i)^{(k)} + \rho g(\mathbf{x}_i)^{(k)}, \ k= 0,1, \dots \] pare quando atingir convergência.

Temos a implementação do algoritmo Gradiente Descendente Boosting no software R, disponível na biblioteca [mboost], cujo autor também generaliza o método para ajuste de modelos lineares generalizados, regressão quantílica, modelos de sobrevivência, dentre outros.


Extreme Gradient Boosting


Extreme Gradient Boosting (XGBoost) utiliza a essência do gradiente boosting, mas vai muito além em termos de qualidade. Ele é o algoritmo mais popular de ML atualmente. Desde a sua criação, em 2014, tornou-se o “true love” dos competidores do Kaggle.

                                  Modelos de Regressão                           Extreme Gradient Boosting

Drawing


Algumas características que o torna tão superior são:

  • Computação em paralelo: por default, ele utiliza todos os cores da máquina;
  • Regularização: controla o trade-off entre vício e variância, permitindo selecionar variáveis etc.
  • Validação cruzada: já possui internamente esse recurso. Não necessitando funções externas.
  • Missing Values: se existir alguma tendência nos dados faltantes, o algoritmo captará;
  • Disponibilidade: atualmente encontra-se disponível em R, Python, Java, Julia e Scala;
  • Save and Reload: permite salvar a matriz de dados e o modelo, podendo recarregá-lo em outro momento.

Certamente, trata-se de um algoritmo extremamente potente. Justamente por isso, otimizá-lo (“tunar” os parâmetros) é uma tarefa complexa. Neste material, vamos apresentar apenas os parâmetros de tuning mais frequentes. A lista completa encontra-se na documentação oficial. Dividiremos em três categorias:

  1. Parâmetros gerais

    • Booster[default=gbtree]
      • Determina o tipo de booster (gbtree, gblinear ou dart). Para classificação podemos utilizar os dois primeiros. Para regressão, qualquer dos três;
  2. Parâmetros de booster

    • nrounds[default=100]
      • Controla o número de iterações. Para classificação, é similar ao número de árvores.
    • eta[default=0.3][range: (0,1)]
      • Controla a taxa de aprendizado. Tipicamente, utilizamos valores entre 0,01 e 0,3.
    • max_depth[default=6][range: (0,Inf)]
      • Controla o tamanho de cada árvore. Geralmente, utilizamos árvores menores, para evitar o superajuste. P. ex., com duas partições.
    • gamma[default=0][range: (0,Inf)]
      • Controla a regularização. Na prática, comece com 0 e verifique se a taxa de erro de validação é elevada. O parâmetro traz bons resultados, especialmente, em árvores pequenas.
    • subsample[default=1][range: (0,1)]
      • Controla o tamanho da amostra em cada árvore. Tipicamente, são valores entre 0,5 e 0,8.
    • colsample_bytree[default=1][range: (0,1)]
      • Controla o número de variáveis apresentadas à árvore. Tipicamente, são valores entre 0,5 e 0,9.
    • lambda[default=0]
      • Controla o penalty \(\ell_2\) nos pesos (equivalente ao Ridge). Utilizado para evitar superajuste.
    • alpha[default=1]
      • Controla o penalty \(\ell_1\) nos pesos (equivalente ao Lasso). Além de encurtar o parâmetro, funciona como selecionador de variáveis. Assim, é mais útil para dados em alta dimensão.
  3. Learning Task Parameters

    • eval_metric [no default, depends on objective selected]
      • Métrica utilizada para avaliar o modelo:
        1. mae - Mean Absolute Error (usada em regressão);
        2. Logloss - Negative loglikelihood (usada em classificação);
        3. AUC - Area under curve (usada em classificação);
        4. RMSE - Root mean square error (usada em regressão);
        5. error - Binary classification error rate (usada em classificação);
        6. mlogloss - multiclass logloss (usada em classificação).