Gradiente descendente

Nesse post farei uma breve introdução sobre o método do gradiente descendente (gradient descent ou steepest descent) e algumas variações importantes para o contexto de aprendizado de máquina.

Esse método é utilizado para otimizar (maximizar/minimizar) algumas funções de modo que, a partir de um ponto inicial, caminha-se em direção ao ponto ótimo controlando o tamanho do passo. Ao decorrer do post ficará mais claro como isso é feito.

É importante ressaltar que, no contexto de aprendizado de máquina, muitas vezes buscamos parâmetros de forma a otimizar (minimizar/maximizar) determinadas funções de interesse. Daí a importância desse método nessa área.

As implementações apresentadas aqui, longe de conter uma programação otimizada e recomendável em uma aplicação, servem para facilitar a compreensão do método. O código apresentado ao longo do texto é implementado em R. No entanto, apresentarei uma versão em Python ao final do post.

Antes de começar, vamos definir o que é o gradiente. O gradiente de uma função $f$ em um ponto $x$ nos fornece a direção em que se obtém o maior incremento na função. Fazendo uma analogia com a subida de um morro, seria para onde daríamos o próximo passo que nos levaria ao ponto mais alto mais rapidamente. Em notação matemática, é dado pelo vetor de derivadas parciais. Considerando os casos uni-variado e bi-variado, teremos os gradientes da função $f$ dados por:

$$\mathrm{\nabla }f(x)=[\frac{\partial }{\partial x}f(x)]\text{e}\mathrm{\nabla }f(x,y)=\left[\begin{array}{c}\frac{\partial }{\partial x}f(x)\\ \frac{\partial }{\partial y}f(y)\end{array}\right].$$

Para uma representação visual do conceito, assista esse vídeo da Khan Academy.

Tenha em mente que a aplicação desse método aqui terá como objetivo otimizar funções (no nosso caso, minimizar funções de perda). Por exemplo, minimizar o erro quadrático médio em função de alguns parâmetros.

Regressão passando pela origem

Começaremos com um exemplo simples, considerando um único parâmetro a ser aprendido/estimado, com base em três observações:

x	y
1	0.519
2	5.577
3	5.043

Considerando a notação usual, y será a variável resposta/target/valor a ser predito e x a variável explicativa/preditora. Trabalharemos com o seguinte modelo:

$$y=\theta x.$$

Note que, esse é um modelo linear passando pela origem, ou seja, pelo ponto $(0,0)$. A princípio, podemos considerar infinitos valores de $\theta$ para esse modelo. Veja alguns exemplos na figura a seguir:

Entre essas infinitas possibilidades, como escolher um valor único para $\theta$? Como comparar os resultados obtidos com diferentes $\theta$´s? Uma alternativa seria considerar uma função de perda (loss function) $L(\theta )$ que penalize valores preditos muito distantes dos valores observados para cada observação e determinar $\theta$ de forma a minimizar essa função de perda. Podemos considerar diferentes funções para essa comparação. Por exemplo,

$$\text{perda quadrática: }L(\theta )=\frac{1}{n}\sum _{i=1}^n(y_i-\theta x_i{)}^2$$

$$\text{perda absoluta: }L(\theta )=\frac{1}{n}\sum _{i=1}^n|y_i-\theta x_i|.$$

Nesse exemplo utilizaremos a perda quadrática. Assim, para esse caso particular com 3 observações e desenvolvendo a expressão, teremos:

$$L(\theta )={\theta }^2(\frac{1}{3}\sum _{i=1}^3{x}_i^2)-\theta (\frac{2}{3}\sum _{i=1}^3{y}_i{x}_i)+\frac{1}{3}\sum _{i=1}^3{y}_i^2.$$ Note que essa função de perda é uma função quadrática em $\theta$. Assim, substituindo os valores de $\mathbf{x}$ e $\mathbf{y}$ pelos valores observados na equação anterior, podemos fazer o gráfico da função de perda em função de $\theta$.

Lembre-se que estamos buscando o valor de $\theta$ que minimiza a função de perda. Dessa forma, queremos encontrar o ponto de mínimo da função apresentada no gráfico acima. Vamos considerar dois chutes iniciais para o valor de $\theta$ que minimiza essa função:

• ponto A: o gradiente é negativo, pois a função está aumentando ao se deslocar para esquerda no eixo horizontal. Portanto, se dermos um passo no sentido do gradiente (negativo), andaremos para esquerda no eixo horizontal e aumentaremos o valor da função de perda. Dessa forma, daremos um passo no sentido contrário ao gradiente da curva.

• ponto B: o gradiente da curva é positivo, pois a função está aumentando ao se avançar no eixo horizontal. Portanto, se dermos um passo no sentido do gradiente (positivo), andaremos para direita no eixo horizontal e aumentaremos o valor da função de perda. Dessa forma, daremos um passo no sentido contrário ao gradiente da curva.

Note que, nas duas situações, o sentido do passo foi oposto ao gradiente. Dessa forma, para nos aproximar do ponto de mínimo, seguiremos uma sequência como:

$${\theta }^{\text{(próximo passo)}}={\theta }^{\text{(passo anterior)}}-\text{fator}\times \text{gradiente}.$$

No gráfico a seguir é apresentada uma descrição do processo apresentado anteriormente partindo do ponto A.

Uma definição mais genérica do gradiente descendente é dada a seguir.

Considere um cenário com um vetor de parâmetros $\theta =({\theta }_1,\dots ,{\theta }_p)$ e uma função $L:{\mathbb{R}}^p\to \mathbb{R}$ que se deseja minimizar.

• Inicialize com um valor ${\theta }^{(0)}\in {\mathbb{R}}^p$

• Repita até convergência ou algum critério de parada com $k=0,1,2\dots$

  $${\theta }^{(k+1)}={\theta }^{(k)}-\eta {\mathrm{\nabla }}_{\theta }L({\theta }^{(k)})$$

• Retorne os parâmetros

Aqui $\eta$ é chamado de taxa de aprendizado e controla o tamanho do passo em cada iteração. O valor dessa taxa é crítico pois valores muito pequenos ou muito altos podem fazer com que o processo demore muito tempo ou não convirja para um valor razoável. É possível considerar valores de $\eta$ variando ao longo das iterações e um grid de valores. Ainda, ${\mathrm{\nabla }}_{\theta }L({\theta }^{(k)})$ indica o gradiente da função $L$ no ponto ${\theta }^{(k)}$.

Na figura a seguir apresentamos um caso com um valor alto para $\eta$.

Em aplicaões, pode se considerar um número grande de iterações com interrupção do processo quando a diferença na função de perda entre duas iterações for muito pequena ou quando a norma do gradiente estiver abaixo de um valor de tolerância fixado.

Implementando o método

Agora que já entendemos o que o algoritmo deve fazer, vamos implementar essa função considerando o modelo $y=\theta x$ com a função de perda quadrática e os dados apresentados anteriormente no primeiro exemplo.

# inicializa os dados 

dados <- tibble(x = c(1, 2, 3),
                y = c(0.519, 5.577, 5.043))

Assim, para o caso unidimensional, o gradiente é dado por:

$$\frac{\partial L(\theta )}{\partial \theta }=\frac{\partial }{\partial \theta }[\frac{1}{3}\sum _{i=1}^3(y_i-\theta x_i{)}^2]=-\frac{2}{3}\sum _{i=1}^3(y_i-\theta x_i)x_i$$

Para nos auxiliar nesse processo, criaremos algumas funções . É importante ressaltar que estamos fazendo uma implementação didática. Dessa forma, o esforço aqui é para tornar os passos o mais explícito possível e não fornecer uma versão otimizada desse método.

A primeira função será uma função que calcula o gradiente.

gradiente <- function(theta) -(2/3) * sum((dados$y - theta*dados$x)*dados$x)

A seguir criaremos a função gd (gradiente descendente) que recebe como argumento theta_0 que será o chute inicial para $\theta$, eta que será a taxa de aprendizagem e iteracoes que definirá o número de iterações realizada.

gd <- function(theta_0, eta, iteracoes) {
  
  theta_i <- theta_0
  
  for (i in 1:iteracoes) theta_i <- theta_i - eta * gradiente(theta_i)
  
  return(theta_i)
  
}

Aplicando a função com theta_0 = -20, eta = 0.1 e iteracoes = 10, chegamos ao seguinte resultado:

gd(-20, .1, 10)

[1] 1.914429

Note que esse valor já é muito próximo ao valor calculado analiticamente ($\approx 1.91$).

Exemplo - Modelo de Regressão Linear

Agora abordaremos uma situação mais próxima do processo de modelagem. Considere dados gerados de acordo com o seguinte modelo:

$$y={\theta }_1{x}_1+{\theta }_2{x}_2+{ϵ}_i\text{ em que }{ϵ}_i\sim \text{N}(0,{\sigma }^2)$$

Nesse exemplo utilizaremos ${10}^5$ observações com ${\theta }_1=5$, ${\theta }_2=2$ e $\sigma =10$.

set.seed(321)

n <- 10^5

dados <- tibble(x_1 = runif(n, -10, 10), 
                x_2 = runif(n, -10, 10), 
                y = 5*x_1 + 2*x_2 + rnorm(n, 0, 10)) 

Utilizaremos, novamente, a perda quadrática $L(\theta )$ dada por:

$$L(\theta )=\frac{1}{n}\sum _{i=1}^n(y_i-({\theta }_1{x}_1+{\theta }_2{x}_2){)}^2.$$

Portanto, o gradiente é dado por

$${\mathrm{\nabla }}_{\theta }L(\theta )=\left(\begin{array}{c}-\frac{2}{n}\sum _{i=1}^n{x}_{1i}(y_i-({\theta }_1{x}_{1i}+{\theta }_2{x}_{2i}))\\ -\frac{2}{n}\sum _{i=1}^n{x}_{2i}(y_i-({\theta }_1{x}_{1i}+{\theta }_2{x}_{2i}))\end{array}\right).$$

Assim, o método pode ser implementado considerando

$$\left(\begin{array}{c}{\theta }_1^{(k+1)}\\ {\theta }_2^{(k+1)}\end{array}\right)=\left(\begin{array}{c}{\theta }_1^{(k)}\\ {\theta }_2^{(k)}\end{array}\right)-\eta \left(\begin{array}{c}-\frac{2}{n}\sum _{i=1}^n{x}_{1i}(y_i-({\theta }_1^{(k)}{x}_{1i}+{\theta }_2^{(k)}{x}_{2i}))\\ -\frac{2}{n}\sum _{i=1}^n{x}_{2i}(y_i-({\theta }_1^{(k)}{x}_{1i}+{\theta }_2^{(k)}{x}_{2i}))\end{array}\right)$$ para $k=0,1,2,\dots$. Note que, nesse método, utilizamos todas as observações em cada iteração. Essa forma de cálculo é denominada gradiente descendente batch. Nesse caso em particular, com um número reduzido de amostras e parâmetros, podemos fazer essas contas facilmente. No entanto, esse procedimento fica impráticavel no contexto de modelos de deep learning, em que podemos considerar milhões de observações e parâmetros.

A seguir abordaremos a primeira variação desse método.

Gradiente Descendente Estocástico (Stochastic Gradient Descent - SGD)

Nessa variação executamos o procedimento e atualização em apenas uma observação por vez. De forma geral, podemos descrever da seguinte forma:

$$\theta -\eta {\mathrm{\nabla }}_{\theta }L(\theta ;{\mathbf{x}}_i,y_i)\text{para }i=1,2,\dots$$

em que ${\mathbf{x}}_i$ e $y_i$ correspondem a $i$-ésima observação sorteada para essa iteração.

Suponha que sorteamos, como primeira unidade para execução do método, a sétima observação. Assim,

$$\left(\begin{array}{c}{\theta }_1^{(1)}\\ {\theta }_2^{(1)}\end{array}\right)=\left(\begin{array}{c}{\theta }_1^{(0)}\\ {\theta }_2^{(0)}\end{array}\right)-\eta \left(\begin{array}{c}-\frac{2}{1}{x}_{1,7}(y_7-({\theta }_1^{(0)}{x}_{1,7}+{\theta }_2^{(0)}{x}_{2,7}))\\ -\frac{2}{1}{x}_{2,7}(y_7-({\theta }_1^{(0)}{x}_{1,7}+{\theta }_2^{(0)}{x}_{2,7}))\end{array}\right).$$

Se no segundo passo a terceira observação for sorteada, na sequência, teremos

$$\left(\begin{array}{c}{\theta }_1^{(2)}\\ {\theta }_2^{(2)}\end{array}\right)=\left(\begin{array}{c}{\theta }_1^{(1)}\\ {\theta }_2^{(1)}\end{array}\right)-\eta \left(\begin{array}{c}-\frac{2}{1}{x}_{1,3}(y_3-({\theta }_1^{(1)}{x}_{1,3}+{\theta }_2^{(1)}{x}_{2,3}))\\ -\frac{2}{1}{x}_{2,3}(y_3-({\theta }_1^{(1)}{x}_{1,3}+{\theta }_2^{(1)}{x}_{2,3}))\end{array}\right).$$

Após executar esse procedimento com todas as observações, dizemos que foi executada uma época/epoch. Podemos repetir esse procedimento para um número arbitrário de épocas.

Na próxima seção apresentaremos uma variação muito utilizada em modelos de deep learning. O gradiente descendente mini batch (gradient descent mini batch).

Gradiente Descendente mini-batch

Nesse método, o conjunto de dados é dividido em pequenos lotes (mini-batches) e o processo é executado para cada lote separadamente. Por exemplo, considere que os dados observados são separados em $B$ lotes. Assim, o processo é executado da seguinte forma

$$\theta -\eta {\mathrm{\nabla }}_{\theta }L(\theta ;{\mathbf{x}}_b,y_b)\text{ para }b=1,\dots ,B$$ em que ${\mathbf{x}}_b$ e $y_b$ indicam que as observações pertencentes ao $b$-ésimo mini-batch. Após percorrer todos os lotes/batches dizemos que o método executou uma época/epoch.

Começamos criando uma função auxiliar para, dado um número de amostras (nrow(dados)), retornar um vetor com as quantidades de elementos de cada lote.

vetor_batches <- function(batch_size) {

 # define o n de batches e o número de observações em cada cada batch
 # A funcao %% dá o resto da divisão: 7 %% 2 retorna 1
  
  if (nrow(dados) %% batch_size == 0) {
    
    n_batch <- nrow(dados)/batch_size
    
    batches <- rep(1:n_batch, batch_size)
  
  } else { 
    
    n_batch <- round((nrow(dados)/ batch_size)) + 1
    
    batches <- c(rep(1:(n_batch - 1), batch_size), 
                 rep(n_batch, nrow(dados) %% batch_size))
    
  }
  
  return(batches)
}

Note que poderíamos considerar uma atribuição aleatória dos batches na função anterior. Por exemplo, fazer a função retornar sample(batches) no lugar de batches. Como estamos trabalhando com dados simulados de forma completamente independente, não adicionei esse passo na função.

A seguir definiremos uma função para calcular, dado um valor de $\theta$ e um conjunto de observações contido em um data frame df, a função de perda:

funcao_perda <- function(theta, df) 0.5*mean((df$y - (theta[[1]]*df$x_1 + theta[[2]]*df$x_2))^2)

O gradiente será calculado pela função de mesmo nome:

gradiente <- function(theta, df) {
  
  return(c(-2*mean(df$x_1*(df$y - (theta[[1]]*df$x_1 + theta[[2]]*df$x_2))), 
           -2*mean(df$x_2*(df$y - (theta[[1]]*df$x_1 + theta[[2]]*df$x_2))))) 
 
}

A partir das funções auxiliares definidas anteriormente, podemos escrever uma função para execução do gradiente descendente mini-batch. Essa função receberá como argumentos:

• theta_incial: um vetor com dois elementos contendo o chute inicial para $\theta$

• eta: taxa de aprendizado $\eta$

• epochs: o número de vezes que o método deverá percorrer todos os lotes

• batch_size: número de elementos em cada lote

gd_mb <- function(theta_inicial, eta, epochs, batch_size) {
  
  batches <- vetor_batches(batch_size)
  
  n_batch <- max(batches)
  
  theta <- tibble(theta_1 = c(theta_inicial[[1]], rep(NA, epochs*n_batch)), 
                  theta_2 = c(theta_inicial[[2]], rep(NA, epochs*n_batch)))
    
  custo <- vector("numeric", length = nrow(theta))
  
  custo[1] <- funcao_perda(theta[1,], dados)
  
  idx <- 2
  
  for (e in 1:epochs) {
    
    for (b in 1:n_batch) {
      
      theta[idx,] <- theta[idx - 1,] - eta * gradiente(theta[idx - 1,], dados[batches == b,])

      custo[idx] <- funcao_perda(theta[idx,], dados[batches = b,])
      
      idx <- idx + 1
    
    }
    
  }
  
  return(bind_cols(theta, custo = custo))
    
}

Vamos inicializar com theta_inicial = c(3, 1), eta = 0.001, epochs = 5 e batch_size = 32.

fit <- lm(y ~ ., dados)

Considerando esses parâmetros, foram obtidas as seguintes estimativas: $\widehat{{\theta }_1}$ = 5.0418137 e $\widehat{{\theta }_2}$ = 1.9349076. Note que esses valores são relativamente próximos dos valores ótimos dados por 5.0052234 e 1.9969684.

A seguir apresentamos um subconjunto da função $L(\theta )$ em função do número da iteração.

No gráfico a seguir é apresentado o caminho percorrido pela execução do método. Note que há uma pequena variação na região central da figura.

Para entender o impacto do tamanho dos lotes, reexecutamos o processo considerando agora batch_size = 320 e apresentamos o resultado no gráfico abaixo.

Note a diferença de variabilidade entre os dois cenários. Sob uma ótica estatística, podemos entender isso pela variabilidade do estimador da função $L(\theta )$ que, nesse caso, é dado por uma média. Podemos escrever $r_i=(y_i-({\theta }_1{x}_{1,i}+{\theta }_2{x}_{2,i}{)}^2$ e interpretar

$$L(\theta )=\frac{1}{n}\sum _{i=1}^n(y_i-({\theta }_1{x}_{1,i}+{\theta }_2{x}_{2,i}{)}^2=\frac{1}{n}\sum _{i=1}^n{r}_i$$

como uma média amostral. Assim, para os dois cenários anteriores com batches de 32 e 320 e considerando $\text{var}(r_i)={\sigma }^2$, temos:

$$\text{Var}(L(\theta ))=\frac{{\sigma }^2}{32}\text{ e }\text{Var}(L(\theta ))=\frac{{\sigma }^2}{320}.$$

Note que a variância do segundo cenário corresponde a 10% da variância do primeiro cenário. Embora fazer o procedimento com tamanhos de amostras maiores seja estatisticamente melhor em termos de variabilidade, não é a melhor escolha computacional considerando tempo de processamento e outros aspectos. A seguir é apresentado um post com um toque de humor do pesquisador Yann LeCun:

Para acessar o artigo citado no post, utilize o link https://arxiv.org/abs/1804.07612.

Conclusão

Como citado no início do post, o objetivo é apresentar de forma simples e didática esse método de otimização de funções. É importante ressaltar que a aplicação mais famosa se dá em modelos de deep learning, mas é um método que pode ser utilizado em diferentes contextos.

A tabela a seguir, retirada do curso do Andrew Ng, ilustra um resumo das variações de métodos apresentadas:

• gradiente descendente batch: utiliza todas as observações em cada iteração

• gradiente descendente estocástico: utilizar uma observação em cada iteração

• gradiente descendente mini-batch: utiliza o número de observações de cada lote em cada iteração

Em um próximo post farei a aplicação desse método no contexto de redes neurais.

Caso tenha alguma crítica, sugestão ou comentário, me envie uma mensagem!

Implementação em Python

Primeiro, será considerada a implementação em python do [Gradiente Descendente mini-batch]. Começaremos carregando as bibliotecas necessárias e gerando um conjunto de dados.

import numpy as np
import pandas as pd

np.random.seed(15)

n = 10**5

# cria dados simulados
dados = pd.DataFrame({'x_1': np.random.uniform(-10, 10, n), 
                      'x_2': np.random.uniform(-10, 10, n)})

dados['y'] = 5*dados['x_1'] + 2*dados['x_2'] + np.random.normal(0, 10, n)

A seguir definiremos as funções auxiliares vetor_batches, funcao_perda e gradiente.

# função para definir o vetor com a identificação dos batches

def vetor_batches(batch_size):
  
    if dados.shape[0] % batch_size == 0:
    
        n_batch = dados.shape[0] / batch_size
    
        batches = np.repeat(np.arange(1, n_batch + 1), batch_size)
  
    else:
        
        n_batch = round(dados.shape[0] / batch_size) + 1
    
        batches = np.concatenate((np.repeat(np.arange(1, n_batch), batch_size), 
                                  np.repeat(n_batch, dados.shape[0] % batch_size)))
    
    return np.array(batches, dtype = 'int')


# cálculo da função de perda de acordo com valores de theta e banco de dados

def funcao_perda(theta, df):
    
    return 0.5 * np.mean((df['y'] - (theta[0]*df['x_1'] + theta[1]*df['x_2'])**2))


# cálculo do gradiente de acordo com theta e banco de dados

def gradiente(theta, df):
    
    return np.array([-2*np.mean(df['x_1']*(df['y'] - (theta[0]*df['x_1'] + theta[1]*df['x_2']))), 
                     -2*np.mean(df['x_2']*(df['y'] - (theta[0]*df['x_1'] + theta[1]*df['x_2'])))])

Por fim, definiremos a função gd_mb e faremos a chamada do método.

# função para o gradiente descendente mini-batch

def gd_mb(theta_inicial, eta, epochs, batch_size):
    
    batches = vetor_batches(batch_size)
  
    n_batch = np.max(batches)
  
    theta = pd.DataFrame({'theta_1': np.concatenate((np.array([theta_inicial[0]]), 
                                                     np.repeat(None, epochs*n_batch))),
                          'theta_2': np.concatenate((np.array([theta_inicial[0]]), 
                                                     np.repeat(None, epochs*n_batch)))})
  
    custo = np.repeat(None, theta.shape[0])
  
    custo[0] = funcao_perda(theta.iloc[0, 0:2], dados) 
  
    idx = 1
    
    for e in range(1, epochs + 1):
        
        for b in range(1, n_batch + 1):
            
            theta.iloc[idx,] = theta.iloc[idx - 1,] - eta * gradiente(theta.iloc[idx - 1,], dados.iloc[batches == b,])
      
            custo[idx] = funcao_perda(theta.iloc[idx, 0:2], dados.iloc[batches == b, ])
      
            idx = idx + 1
      
    theta['custo'] = custo
      
    return theta


# executa a tarefa

gd_mb([3, 1], 0.001, 5, 32)

Referências

• Géron, Aurélien. (2019) Hands-on Machine Learning with Scikit-Learn, Keras and TensorFlow: Concepts, Tools, and Techniques to Build Intelligent Systems. 2nd ed. CA 95472: O’Reilly.

• Goodfellow, I., Bengio, Y. and Courville, A. (2016) Deep Learning. MIT Press.

• Izbicki, R. e Santos, T. M. dos. (2020) Aprendizado de máquina: uma abordagem estatística. 1ᵃ edição. 272 páginas. ISBN: 978-65-00-02410-4.