Entendendo o Gradiente Descent (GD)

O Gradient Descent (GD), também conhecido como steepest descent ou até como Stochastic Gradient Descent (SGD)¹, é um algoritmo fundamental em otimização que desempenha um papel crucial no ajuste de parâmetros de modelos de machine learning. Basicamente, o GD utiliza a informação do gradiente de uma função, que se quer otimizar, para encontrar o “caminho” em que ela decaí mais rapidamente.

No caso mais típico, da minimização de uma função de perda, o GD ou SGD usa o gradiente da função para encontrar os valores dos parâmetros que minimizam esta função. O algoritmo opera iterativamente e, idealmente, aproxima-se gradualmente da solução correta.

Assim, a utilidade geral do GD reside na sua capacidade de iterativamente encontrar o mínimo de funções usando relativamente pouca informação.

Regressão Simples

Um pouco de matemática

Para tornar as coisas mais concretas vamos considerar o caso de uma regressão linear simples. Temos uma variável de resposta \(y\) que será “explicada” por uma variável independente (feature) \(x\) . A relação é modelada de maneira linear como:

\[ y = \beta_{0} + \beta_{1}x + \varepsilon \]

A equação acima descreve o nosso modelo, que será estimado posteriormente. Com este modelo podemos calcular os valores preditos \(\hat{y_{i}}\) que serão comparados com os valores observados (reais) \(y_{i}\). Nosso objetivo é ajustar os parâmetros \(\beta_{0}\) e \(\beta_{1}\) para minimizar a função de custo, que pode ser representada como o Erro Quadrático Médio (MSE, do inglês Mean Squared Error):

\[ \text{MSE} = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y_{i}})^2= \frac{1}{n} \sum_{i=1}^{n} (y_i - (\beta_0 + \beta_1 X_i))^2 \]

Para atualizar os parâmetros \(\beta_0\) e \(\beta_1\) usando GD, precisamos calcular o gradiente da função de custo em relação a esses parâmetros. A primeira derivada da função de custo em relação a \(\beta_0\) e \(\beta_1\) nos dá as direções em que devemos ajustar esses parâmetros para minimizar a função de custo.

A primeira derivada da função de custo em relação a \(\beta_0\) é dada por:

\[ \frac{\partial \text{MSE}}{\partial \beta_0} = \frac{-2}{n} \sum_{i=1}^{n} (y_i - (\beta_0 + \beta_1 x_i)) \]

E a primeira derivada da função de custo em relação a \(\beta_1\) é dada por:

\[ \frac{\partial \text{MSE}}{\partial \beta_1} = \frac{-2}{n} \sum_{i=1}^{n} x_i (y_i - (\beta_0 + \beta_1 x_i)) \]

Essas derivadas nos fornecem a direção em que devemos ajustar os parâmetros para minimizar a função de custo.

Implementando o algoritmo

O algortimo de GD funciona iterativamente O algoritmo de SGD atualiza o parâmetro \(\beta^{t}\) a cada iteração t, onde \(\beta^{0}\) é dado, usando o gradiente \(\nabla_{\beta}^t\) da seguinte maneira:

\[ \beta^{t+1} = \beta^{t} - \gamma\nabla_{\beta}^{t} \]

onde \(\gamma\) é um número real não-negativo, tipicamente próximo de 0.001, chamado “learning rate”. Quanto maior for o valor de \(\gamma\) maiores serão os “passos” no processo de atualização; inversamente, quanto menor for o valor de \(\gamma\) menores serão os “passos”no processo iterativo.

Para implementar o passo-a-passo do algoritmo vamos usar a base mtcars. O primeiro passo é ajustar os dados usando a função scale. Em seguida, separamos alguns objetos úteis para facilitar a exposição.

#> Regularizar os dados
mtcars_scaled <- as.data.frame(scale(mtcars))

y <- mtcars_scaled$mpg
x <- mtcars_scaled$wt
N <- nrow(mtcars_scaled)

Vamos fazer a regressão de mpg (milhas por galão), uma medida da eficiência do veículo, contra wt (peso), o peso do veículo. Visualmente, parece haver uma relação linear decrescente entre as variáveis.

plot(y ~ x)

Antes de fazer o loop, vamos decompor o algoritmo em etapas. Primeiro, precisa-se de valores iniciais para os parâmetros \(\beta_{0}\) e \(\beta_{1}\). Por simplicidade, vamos sortear números aleatórios entre 0 e 1 a partir de uma distribuição uniforme. Com estes valores será possível calcular o valor de \(\hat{y}_{i}^0\), onde 0 indica que estamos na iteração de valor 0.

b0 <- runif(1)
b1 <- runif(1)

(yhat <- b0 + b1 * x)

 [1] -0.17140720  0.05168499 -0.43386860  0.34914125  0.54598730  0.56348473
 [7]  0.65972058  0.32726947  0.29227461  0.54598730  0.54598730  1.09715625
[13]  0.79969999  0.84344356  2.12950443  2.28173204  2.21261721 -0.53885316
[19] -1.05065289 -0.85818120 -0.30701225  0.61597701  0.54161295  0.89593584
[25]  0.90031019 -0.77069407 -0.59134544 -1.13988977  0.30977204 -0.04017650
[31]  0.65972058 -0.03142778

Não é necessário, mas é instrutivo calcular a função de perda.

(mse <- mean((y - yhat)^2))

[1] 3.240961

Agora, calculamos o valor do gradiente neste ponto.

(gb0 <- sum(y - yhat) * (-2 / N))

[1] 0.7022194

(gb1 <- sum((y - yhat) * x) * (-2 / N))

[1] 3.339637

Por fim, o valor dos parâmetros é atualizado segundo a fórmula matemática do algoritmo. Utiliza-se \(g = 0.01\) como valor para a learning-rate.

g = 0.01

(b0_new <- b0 - g * gb0)

[1] 0.3440875

(b1_new <- b1 - g * gb1)

[1] 0.822628

Este processo será repetido \(T\) vezes até que se atinja algum critério de convergência. Em geral, estabelece-se

Um número máximo de iterações. (10.000 iterações, por exemplo).
Um valor mínimo de mudança na estimativa dos parâmetros. Isto é, quando o valor das estimativas para de mudar significativamente, entende-se que ele convergiu para um valor satisfatório.

Para deixar o loop abaixo mais simples, vou simplesmente estabelecer um número máximo de 5000 iterações. O código segue abaixo. Note que algumas partes do código acima foram repetidas por conveniência da leitura.

num_iterations <- 1000
#> Learning-rate
g <- 0.01

y <- mtcars_scaled$mpg
x <- mtcars_scaled$wt
N <- nrow(mtcars_scaled)

#> Valores iniciais para estimativas dos parâmetros
b0 <- runif(1)
b1 <- runif(1)

for (i in seq_len(num_iterations)) {
  if (i %% 100 == 0) {
    cat("Iteração: ", i, "\n")
  }

  #> Calcula o valor previsto
  yhat <- b0 + b1 * x

  #> Calcula a "função de perda"
  error <- y - yhat
  mse <- mean(error^2)

  if (i %% 100 == 0) {
    cat("Valor da perda: ", as.numeric(mse), "\n")
  }

  #> Calcula o gradiente nos pontos atuais
  gb0 <- sum(y - yhat) * (-2 / N)
  gb1 <- sum((y - yhat) * x) * (-2 / N)
  #> Atualiza o valor dos parâmetros usando o gradiente
  b0_new <- b0 - g * gb0
  b1_new <- b1 - g * gb1

  b0 <- b0_new
  b1 <- b1_new

  if (i %% 100 == 0) {
    cat("Betas: ", c(b0, b1), "\n\n")
  }
}

Iteração:  100 
Valor da perda:  0.3013997 
Betas:  0.1179619 -0.6505911 

Iteração:  200 
Valor da perda:  0.2406465 
Betas:  0.01564405 -0.8369768 

Iteração:  300 
Valor da perda:  0.2394667 
Betas:  0.002074708 -0.8633224 

Iteração:  400 
Valor da perda:  0.2394437 
Betas:  0.0002751468 -0.8670463 

Iteração:  500 
Valor da perda:  0.2394432 
Betas:  3.648985e-05 -0.8675727 

Iteração:  600 
Valor da perda:  0.2394432 
Betas:  4.839267e-06 -0.8676471 

Iteração:  700 
Valor da perda:  0.2394432 
Betas:  6.417815e-07 -0.8676576 

Iteração:  800 
Valor da perda:  0.2394432 
Betas:  8.511277e-08 -0.8676591 

Iteração:  900 
Valor da perda:  0.2394432 
Betas:  1.128762e-08 -0.8676593 

Iteração:  1000 
Valor da perda:  0.2394432 
Betas:  1.49696e-09 -0.8676594

Para recuperar o valor final dos betas.

c(b0, b1)

[1]  1.496960e-09 -8.676594e-01

Para efeito didático, vamos comparar estas estimativas finais contra os valores estimados pela função lm. Note que os valores estão muito similares

summary(model_lm <- lm(mpg ~ wt, data = mtcars_scaled))


Call:
lm(formula = mpg ~ wt, data = mtcars_scaled)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.75381 -0.39236 -0.02077  0.23388  1.14033 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.040e-15  8.934e-02   0.000        1    
wt          -8.677e-01  9.077e-02  -9.559 1.29e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5054 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Visualizando o processo

Sempre que possível, é útil visualizar o funcionamento do algoritmo em gráficos. Como estamos trabalhando com um exemplo simples, pode-se plotar os resultados gradativamente num gráfico de dispersão. O código abaixo mostra como a linha de ajuste (linha de regressão) vai se alterando à medida que se aumenta o número de amostras.

Code

num_iterations <- 1000
#> Learning-rate
alpha <- 0.01

y <- dat$mpg
x <- dat$wt
N <- nrow(dat)

#> Valores iniciais para estimativas dos parâmetros
b0 <- runif(1)
b1 <- runif(1)

betas <- matrix(ncol = 2, nrow = num_iterations)
for (i in seq_len(num_iterations)) {
  #> Calcula o valor previsto
  yhat <- b0 + b1 * x

  #> Calcula a "função de perda"
  error <- y - yhat
  mse <- mean(error^2)

  #> Calcula o gradiente nos pontos atuais
  gb0 <- sum(y - yhat) * (-2 / N)
  gb1 <- sum((y - yhat) * dat$wt) * (-2 / N)
  #> Atualiza o valor dos parâmetros usando o gradiente
  b0_new <- b0 - alpha * gb0
  b1_new <- b1 - alpha * gb1

  b0 <- b0_new
  b1 <- b1_new

  betas[i, 1] <- b0_new
  betas[i, 2] <- b1_new
}

sel <- c(1:10, seq(20, 100, 10), seq(100, 1000, 50))

library(ggplot2)
library(gganimate)
library(tibble)

tbl_betas <- tibble(
  iter = sel,
  beta0 = betas[sel, 1],
  beta1 = betas[sel, 2]
)

ggplot() +
  geom_point(
    data = dat,
    aes(x = wt, y = mpg),
    shape = 21
  ) +
  geom_abline(
    data = tbl_betas,
    aes(intercept = beta0, slope = beta1),
    color = "#CB181D",
    lwd = 0.8
  ) +
  geom_text(
    data = tbl_betas,
    aes(x = 1.8, y = 1.8, label = paste("Iteration:", iter)),
    size = 5
  ) +
  transition_states(iter) +
  enter_fade() +
  exit_shrink() +
  theme_bw()

Regressão Múltipla

Agora vamos considerar o caso de regressão múltipla, onde temos várias variáveis independentes. A equação do modelo é generalizada para:

\[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]

Isto é, agora temos \(k\) variáveis independentes, \(x_1, x_2, ..., x_k\), e temos \(k\) coeficientes, \(\beta_0, \beta_1, \beta_2, ..., \beta_k\), a ser estimados. A função de custo torna-se:

\[ \text{MSE} = \frac{1}{n} \sum_{i=1}^{n} (y_i - (\beta_0 + \beta_{1} x_{i1} + \beta_2 x_{i2} + ... + \beta_k X_{ik}))^2 \]

Agora, as derivadas parciais da função de custo em relação a cada parâmetro \(\beta\) são calculadas e utilizadas para atualizar os coeficientes durante o GD.

\[ \nabla_{\beta} = (\frac{\partial\text{MSE}}{\partial\beta_{0}}, \frac{\partial\text{MSE}}{\partial\beta_{1}}, ..., \frac{\partial\text{MSE}}{\partial\beta_{k}}) \]

Regressão Múltipla com Matrizes

Na regressão múltipla, podemos representar os dados de entrada \(X\) e os parâmetros do modelo \(\beta\) como matrizes. Esta forma de representação é mais prática quando temos muitas variáveis e permite dispensar o uso de somatórios.

Suponha que tenhamos \(n\) observações e \(k\) variáveis independentes.

As observações de entrada podem ser organizadas em uma matriz \(X\) de dimensão \(n \times (k+1)\), onde a primeira coluna é composta por \(1\)s para representar o intercepto do modelo. Assim, a matriz \(X\) é dada por:

\[ X = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{nk} \end{bmatrix} \]

Os parâmetros do modelo \(\beta\) podem ser representados como um vetor de coeficientes de dimensão \((k+1) \times 1\). Assim, o vetor \(\beta\) é dado por:

\[ \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_k \end{bmatrix} \]

A resposta \(y\) pode ser representada como um vetor de dimensão \(n \times 1\).

\[ y = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \end{bmatrix} \]

Note que o problema de minimização acima é equivalente a minimizar a soma do erro ao quadrado do modelo. Isto acontece pois o erro, no caso mais simples, é simplesmente

\[ \varepsilon = y_i - \hat{y_i} = y_i - \beta_{0} - \beta_{1}x_i \]

e o problema de minimizar:

\[ \text{min } \varepsilon^2 = (y_i - \hat{y_i})^2 \]

Em termos matriciais temos:

\[ e = y - \hat{y} = y - X\beta \] e agora o problema de minimizar

\[ \text{min } e^te = (y - X\beta)^t(y-XB) \]

Para encontrar o gradiente da função de custo \(e^te\), em relação aos parâmetros \(\beta\), podemos usar cálculo matricial.

O gradiente \(\nabla_{\beta} (e'e)\) é dado por:

\[ \nabla_{\beta} (e'e) = -2X^T(y - X\beta) \]

Onde \(X^T\) representa a transposta da matriz X. Este gradiente nos fornece a direção em que devemos ajustar os (múltiplos) parâmetros \(\beta\) para minimizar a função de custo \(e'e\).

Implementando o algoritmo

Desta vez, vamos implementar tanto o gradiente como a função de perda como functions no R.

grad <- function(beta) {

  (2/N) * t(X) %*% (X %*% beta - y)

}

loss <- function(beta) {

  e = y - X %*% beta

  t(e) %*% e

}

O modelo

Nosso modelo de regressão agora terá a forma:

\[ \text{mpg} = \beta_{0} + \beta_{1}\text{wt} + \beta_{2}\text{qsec}+ \beta_{3}\text{am} \]

onde mpg e wt tem as mesmas definições dadas acima; já qsec é uma medida de velocidade do veículo e am é uma variável binária que indica se o câmbio do veículo é manual ou automático.

Os dados

Desta vez, o preparo dos dados será feito usando o pacote dplyr.

library(dplyr)

dat <- mtcars |>
  select(c("mpg", "wt", "qsec", "am")) |>
  mutate(across(everything(), ~as.numeric(scale(.x))))

y <- dat$mpg
X <- as.matrix(dat[, c("wt", "qsec", "am")])
X <- cbind(1, X)
colnames(X)[1] <- c("coef")
N <- nrow(X)

Primeira iteração

Novamente, para ganhar um pouco de intuição vamos rodar o

#> Valor inicial para os betas
beta <- runif(ncol(X))

# Opcional
#> Computa a "predição" do modelo
yhat <- X %*% beta
#> Calcula o valor da função de perda
l <- loss(beta)

#> Atualiza o valor dos beta
beta_new <- beta - alpha * grad(beta)

O loop completo

O código abaixo mostra o loop completo. Fora algumas pequenas modificações, ele é exatamente igual ao loop anterior. Neste segundo exemplo, eu reduzo o valor da learning-rate e aumento o número de iterações.

Code

beta <- runif(ncol(X))
num_iterations <- 10000
alpha <- 0.001

for (i in seq_len(num_iterations)) {

  if (i %% 1000 == 0) cat("Iteração: ", i, "\n")

  #> Calcula o valor previsto
  yhat <- X %*% beta

  #> Calcula a "função de perda"
  vl_loss <- loss(beta)

  if (i %% 1000 == 0) {
    cat("Valor da perda: ", as.numeric(vl_loss), "\n")
  }

  #> Calcula o gradiente nos pontos atuais
  grad_current <- grad(beta)
  #> Atualiza o valor dos parâmetros usando o gradiente
  beta_current <- beta - alpha * grad_current

  beta <- beta_current

  if (i %% 1000 == 0) {
    cat("Betas: ", beta, "\n\n")
  }

}

Iteração:  1000 
Valor da perda:  7.843567 
Betas:  0.1035636 -0.1690523 0.5729491 0.6570115 

Iteração:  2000 
Valor da perda:  5.904957 
Betas:  0.01398776 -0.3430823 0.5113749 0.53897 

Iteração:  3000 
Valor da perda:  5.220956 
Betas:  0.001889251 -0.4398704 0.4639382 0.4426585 

Iteração:  4000 
Valor da perda:  4.914118 
Betas:  0.0002551707 -0.5040486 0.4311833 0.377357 

Iteração:  5000 
Valor da perda:  4.775269 
Betas:  3.446451e-05 -0.5471717 0.4090413 0.3334177 

Iteração:  6000 
Valor da perda:  4.712415 
Betas:  4.654933e-06 -0.5761808 0.3941314 0.3038562 

Iteração:  7000 
Valor da perda:  4.683963 
Betas:  6.287163e-07 -0.5956981 0.3840982 0.2839671 

Iteração:  8000 
Valor da perda:  4.671083 
Betas:  8.491727e-08 -0.6088296 0.3773476 0.2705854 

Iteração:  9000 
Valor da perda:  4.665252 
Betas:  1.146931e-08 -0.6176647 0.3728056 0.2615821 

Iteração:  10000 
Valor da perda:  4.662613 
Betas:  1.549098e-09 -0.6236091 0.3697497 0.2555244

Por fim, temos o valor final dos betas estimados.

beta_current

              [,1]
coef  1.549098e-09
wt   -6.236091e-01
qsec  3.697497e-01
am    2.555244e-01

Novamente, podemos comparar estas estimativas com aquelas calculadas pela função lm. Note que, neste caso, mesmo após 10.000 iterações ainda há algumas pequenas divergências entre os valores.

summary(model_lm <- lm(mpg ~ wt + qsec + am, data = dat))


Call:
lm(formula = mpg ~ wt + qsec + am, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.5776 -0.2581 -0.1204  0.2341  0.7734 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.485e-15  7.212e-02   0.000 1.000000    
wt          -6.358e-01  1.155e-01  -5.507 6.95e-06 ***
qsec         3.635e-01  8.559e-02   4.247 0.000216 ***
am           2.431e-01  1.168e-01   2.081 0.046716 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.408 on 28 degrees of freedom
Multiple R-squared:  0.8497,    Adjusted R-squared:  0.8336 
F-statistic: 52.75 on 3 and 28 DF,  p-value: 1.21e-11

Conclusão

O Gradiente Descendente é uma ferramenta poderosa para otimização em aprendizado de máquina e outros campos. Neste post, explicamos a matemática por trás do GD, mostramos como derivar o gradiente tanto para regressão linear simples quanto múltipla, e como entender o funcionamento deste algoritmo fundamental.

Posts relacionados

Footnotes

Apesar dos métodos serem formalmente diferentes, a sua essência é idêntica a ponto de ser comum confundi-los. O SGD é uma “aproximação” do GD, onde apenas uma parte (amostra) dos dados é utilizada para calcular o gradiente. Neste sentido, o SGD é muito mais eficiente do ponto de vista computacional.↩︎