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:

\[ \nabla f(x) = \left[\frac{\partial}{\partial x}f(x) \right] \quad \text{e} \quad \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:

\[{\Large 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\Bigg(\frac{1}{3}\sum_{i = 1}^{3}x_{i}^{2}\Bigg) - \theta\Bigg(\frac{2}{3}\sum_{i = 1}^{3}y_{i}x_i\Bigg) + \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 \(\mathbf{\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 \(\mathbf{\theta}^{(0)} \in \mathbb{R}^{p}\)

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

  \[\mathbf{\theta}^{(k + 1)} = \mathbf{\theta}^{(k)} - \eta \nabla_\mathbf{\theta}L(\mathbf{\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, \(\nabla_\mathbf{\theta}L(\mathbf{\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} \left[ \frac{1}{3} \sum_{i=1}^{3}(y_i - \theta x_i)^2 \right] = -\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 \(\mathbf{\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_1x_1 + \theta_2x_2 + \epsilon_i \qquad \text{ em que } \epsilon_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(\mathbf{\theta})\) dada por:

\[L(\theta) = \frac{1}{n}\sum_{i=1}^{n}(y_i - (\theta_1x_1 + \theta_2x_2))^2.\]

Portanto, o gradiente é dado por

\[\nabla_\mathbf{\theta}L(\mathbf{\theta}) = \left( \begin{array}{c} -\frac{2}{n}\sum_{i=1}^{n}x_{1i}(y_i - (\theta_1x_{1i} + \theta_2x_{2i}))\\ -\frac{2}{n}\sum_{i=1}^{n}x_{2i}(y_i - (\theta_1x_{1i} + \theta_2x_{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:

\[\mathbf{\theta} -\eta \nabla_\mathbf{\theta}L(\mathbf{\theta}; \mathbf{x}_i, y_i) \qquad \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

\[\mathbf{\theta} -\eta \nabla_\mathbf{\theta}L(\mathbf{\theta}; \mathbf{x}_{b}, y_b) \qquad \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: \(\hat{\theta_1}\) = 5.0418137 e \(\hat{\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(\mathbf{\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(\mathbf{\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(\mathbf{\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(\mathbf{\theta})) = \frac{\sigma^2}{32} \qquad \text{ e } \qquad \text{Var}(L(\mathbf{\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.