Tutorial 4.A: Redes Neurais Artificiais com o pacote Torch (R)

Este tutorial é adaptado de um script produzido pelos cientistas de dados Daniel Falbel e Sigrid Keydana, responsáveis peos pacotes que usaremos a seguir, e é baseado nos exemplos do livro An Introduction to Statistical Learning, de Gareth James, Daniela Witten, Trevor Hastie e Rob Tibshirani (https://www.statlearning.com/).

Rede de Camada Única

Começamos carregando os dados Hitters, que contém informações sobre jogadores de beisebol da Major League Baseball (MLB) das temporadas de 1986 e 1987, diponível no pacote ISLR2. Excluimos NAs (com a função na.omit()) e separamos os dados remanescentes em um conjunto de treino e um conjunto de teste:

library(ISLR2)
Gitters <- na.omit(Hitters)
n <- nrow(Gitters)
set.seed(13)
ntest <- trunc(n / 3)
testid <- sample(1:n, ntest)

Primeiro, ajustaremos um modelo linear. Observe o uso do comando with() abaixo: o primeiro argumento é um data frame, e o segundo é uma expressão que pode se referir aos elementos/variáveis do data frame pelo nome. Neste caso, o data frame corresponde aos dados de teste e a expressão calcula o erro absoluto médio de predição nesses dados (rode esses elementos idividualmente no R se quiser entender o que eles representam):

lfit <- lm(Salary ~ ., data = Gitters[-testid, ])
lpred <- predict(lfit, Gitters[testid, ])
with(Gitters[testid, ], mean(abs(lpred - Salary)))
[1] 254.6687

Em seguida, normalizamos as variáveis independentes criando um novo objeto chamado x e mantemos os valores originais de nossa variável dependente em um novo objeto chamado y:

x <- scale(model.matrix(Salary ~ . - 1, data = Gitters))
y <- Gitters$Salary

A primeira linha do script acima utiliza a função model.matrix(), que produz a mesma matriz utilizada por lm() (o termo -1 omite o intercepto). Note que essa função converte automaticamente fatores em variáveis dummies (i.e., com \(0\)s e \(1\)s). A função scale() padroniza a matriz para que cada coluna tenha média \(0\) e variância \(1\).

A seguir, carregamos os pacotes necessários para ajuste de redes neurais artificiais que utilizaremos nesse tutorial:

library(torch)
library(luz) # interface para o pacote `torch`
library(torchvision) # para conjuntos de dados e transformações de imagens
library(torchdatasets) # para os conjuntos de dados
library(zeallot)
torch_manual_seed(13)

Para ajustar uma rede neural, primeiro definimos sua estrutura (note que os scripts a seguir são bastante complexos e provavelmente tomará bastante tempo e prática para entendê-los totalmente):

modnn <- nn_module(
  initialize = function(input_size) {
    self$hidden <- nn_linear(input_size, 50)
    self$activation <- nn_relu()
    self$dropout <- nn_dropout(0.4)
    self$output <- nn_linear(50, 1)
  },
  forward = function(x) {
    x %>% 
      self$hidden() %>% 
      self$activation() %>% 
      self$dropout() %>% 
      self$output()
  }
)

Criamos um simples objeto em R (modelo) chamado modnn definindo as funções initialize e forward, dentro da função nn_module(). A função initialize é responsável por inicializar os submódulos que serão usados pela rede neural. Em forward, implementamos o que acontece quando o modelo é chamado com dados de entrada. Neste caso, usamos as camadas definidas em initialize(), seguindo a ordem específicada.

No R, self é um objeto especial, similar a uma lista, usado para compartilhar informações entre os métodos dentro do nn_module(). Quando você atribui um objeto a self em initialize(), ele poderá ser acessado posteriormente pela função forward().

O operador pipe (i.e., %>%) passa o termo anterior como o primeiro argumento da próxima função, e retorna o resultado. Para entender melhor, vamos ilustrar o uso do operador pipe em um exemplo simples. Anteriormente, criamos x usando o comando:

x <- scale(model.matrix(Salary ~ . - 1, data = Gitters))

Primeiro criamos uma matriz e, em seguida, normalizamos cada uma das variáveis. No entanto, podemos obter o mesmo resultado usando o operador pipe (%>%):

x <- model.matrix(Salary ~ . - 1, data = Gitters) %>% scale()

Ou seja, usar o operador pipe facilita o acompanhamento da sequência de operações.

Agora voltamos para nossa rede neural e vamos interpretá-la:

  • O objeto modnn possui \(1\) única camada oculta (hidden layer), definida como self$hidden, com \(50\) unidades (nodes) ou neurônios (neurons) e uma função de ativação (self$activation) do tipo ReLU (nn_relu()).

  • Em seguida, temos uma “camada” de self$dropout (dropout layer), onde aleatórios \(40\%\) dos \(50\) neurônios da camada anterior são definidos como zero durante cada iteração do algoritmo de descida de gradiente estocástico (stochastic gradient descent, ou SGD). Note que a dropout layer não é propriamente uma camada típica de cálculo: em geral, ela é uma técnica usada para evitar o overfitting (sobreajuste) em redes neurais. Durante o treinamento, nem todos os neurônios são usados. Como cada iteração do treinamento pode ter um conjunto diferente de neurônios “ativos”, o modelo é forçado a aprender dependências/relações mais distribuídas entre eles, e não confiar em um único neurônio ou conjunto de neurônios. Isso ajuda o modelo a generalizar melhor quando aplicado a dados novos.

  • Finalmente, a camada de saída (self$output) possui apenas \(1\) neurônio (ou seja, convergindo de \(50\) para \(1\): nn_linear(50, 1)), sem função de ativação, indicando que o modelo fornece uma saída quantitativa única.

Mas ainda não estamos pronto para ajustar nossa rede neural. Precisamos adicionar alguns “detalhes” ao modnn que controlam o algoritmo de ajuste. Vamos minimizar o erro quadrático médio (MSE, em inglês), que representa a função de perda da nossa rede (loss = nn_mse_loss()). Essa função de perda é uma métrica comum para quantificar a diferença entre os valores previstos pelo modelo e os valores reais em tarefas de regressão, onde a saída esperada é um valor contínuo. O algoritmo acompanhará o erro absoluto médio (MAE, em inglês), nos dados de treinamento e nos dados de validação (metrics = list(luz_metric_mae())):

modnn <- modnn %>% 
  setup(
    loss = nn_mse_loss(),
    optimizer = optim_rmsprop,
    metrics = list(luz_metric_mae())
  ) %>% 
  set_hparams(input_size = ncol(x))

No script acima, o operador pipe (%>%) passa modnn como o primeiro argumento para setup(). A função setup() embute essas especificações em um novo objeto de modelo. Também usamos set_hparam() para especificar os argumentos que devem ser passados para o método initialize() de modnn—um pouco confuso, não é mesmo?

Agora ajustaremos o modelo com o script abaixo. Fornecemos os dados de treinamento (através da linha list(x[-testid, ], matrix(y[-testid], ncol = 1))) e o número de epochs. Por padrão, a cada passo do SGD, o algoritmo seleciona aleatoriamente \(32\) observações de treinamento para o cálculo do gradiente. Lembre-se que uma época (epoch) corresponde ao número de passos/iterações necessários para processar/passar \(n\) observações pela rede. Como nosso conjunto de treinamento tem \(n=176\), uma epoch corresponde a \(176/32=5.5\) passos de SGD. (Por sinal, antes de prosseguir, tente entender o que list(x[-testid, ], matrix(y[-testid], ncol = 1)) faz com os dados de treinamento e de onde vem \(n=176\).)

A função fit() tem um argumento valid_data. Os dados alocados ao valid_data não são usados no ajuste e podem ser usados para acompanhar o progresso do modelo (neste caso, reportando o erro absoluto médio—ou MAE, em inglês). Aqui, fornecemos os dados de teste para que possamos ver o MAE tanto dos dados de treinamento quanto dos dados de teste conforme as epoch avançam (para ver mais opções de ajuste, use ?fit.luz_module_generator):

fitted <- modnn %>% 
  fit(
    data = list(x[-testid, ], matrix(y[-testid], ncol = 1)),
    valid_data = list(x[testid, ], matrix(y[testid], ncol = 1)),
    epochs = 5
  )

Note que aqui utilizamos um número de epochs bastante baixo para diminuir o tempo de computação, porém podemos ver que o modelo claramente melhoraria com mais iterações… Na verdade, é comum ver redes rodando com milhares de epochs, tudo depende de como a acurácia/erro continuam a mudar durante cada iteração do modelo (tanto no conjunto de treinamento quanto no conjunto de validação).

Podemos plotar o modelo fitted para exibir o erro absoluto médio para os dados de treinamento e teste:

plot(fitted)

Por fim, fazemos previsões a partir do modelo final e avaliamos seu desempenho nos dados de teste (devido ao uso do SGD, os resultados podem variar ligeiramente a cada ajuste):

npred <- predict(fitted, x[testid, ])
mean(abs(y[testid] - as.matrix(npred)))
[1] 349.7425

Vamos plotar um gráfico de dispersão, colocando no eixo \(x\) os valores verdadeiros y[testid] (as respostas reais do conjunto de teste) e no eixo \(y\) as predições npred feitas pelo modelo nos mesmos dados de teste. Se o modelo for perfeito, os pontos ficarão alinhados na diagonal (linha reta de 45º). Quanto mais espalhados os pontos em relação a essa linha ideal, pior a predição:

plot(y[testid], as.matrix(npred), abline(a = 0, b = 1, col = "red", lwd = 2))

Note o detalhe que tivemos que converter o objeto npred para uma matriz, pois o método predict atual retorna um objeto da classe torch_tensor.

class(npred)
[1] "torch_tensor" "R7"          

Rede Multicamadas

O pacote torchvision vem com diversos conjuntos de dados de exemplo, incluindo o conjunto de dados MNIST(Modified National Institute of Standards and Technology database). Esse conjunto de dados é muito famoso em machine learning e deep learning, e é especialmente usado para treinar e testar algoritmos de classificação de imagens.

Nosso primeiro passo para ajustar uma rede multicamadas é carregar o conjunto de dados MNIST. A função mnist_dataset() é fornecida para esse propósito. Essa função retorna uma estrutura de dados implementada no torch, incluindo o processo de aquisição e armazenamento dos dados no R:

train_ds <- mnist_dataset(root = ".", train = TRUE, download = TRUE)
test_ds <- mnist_dataset(root = ".", train = FALSE, download = TRUE)

str(train_ds[1])
List of 2
 $ x: int [1:28, 1:28] 0 0 0 0 0 0 0 0 0 0 ...
 $ y: int 6
str(test_ds[2])
List of 2
 $ x: int [1:28, 1:28] 0 0 0 0 0 0 0 0 0 0 ...
 $ y: int 3
length(train_ds)
[1] 60000
length(test_ds)
[1] 10000

\(60.000\) imagens nos dados de treinamento e \(10.000\) nos dados de teste. As imagens são de tamanho \(28\times 28\), e armazenadas como matrizes de pixels.

Primeiramente, precisamos transformar cada imagem em um vetor. Redes neurais são um tanto sensíveis à escala das entradas. Aqui, as entradas são valores armazenados em \(8\) bits, ou seja, variando entre \(0\) e \(255\). Então, fazemos o reescalonamento para o intervalo unitário, de \(0\) a \(1\). (Curiosidade: \(8\) bits significa \(2^8\), que é igual a \(256\). Como a convenção é começar do \(0\), os valores possíveis variam de \(0\) a \(255\).)

Para aplicar essas transformações, vamos definir train_ds e test_ds passando o argumento transform, que aplicará uma transformação a cada uma das imagens de entrada:

transform <- function(x) {
  x %>% 
    torch_tensor() %>% 
    torch_flatten() %>% 
    torch_div(255)
}
train_ds <- mnist_dataset(
  root = ".", 
  train = TRUE, 
  download = TRUE, 
  transform = transform
)
test_ds <- mnist_dataset(
  root = ".", 
  train = FALSE, 
  download = TRUE,
  transform = transform
)

Agora estamos prontos para ajustar nossa rede neural multicamada. Já que essa parte do script pode ser bastante confusa, ela vem com comentários sobre cada linha:

# Criando um modelo neural com três camadas lineares e duas camadas de `dropout`:
modelnn <- nn_module(
  initialize = function() {
    # Definindo a primeira camada linear (`input layer`) com 784 neurônios de entrada (das imagens 28x28) e 256 neurônio de saídas:
    self$linear1 <- nn_linear(in_features = 784, out_features = 256)
    
    # Definindo a segunda camada linear com 256 neurônios de entrada (da camada anterior) e 128 neurônios de saída:
    self$linear2 <- nn_linear(in_features = 256, out_features = 128)
    
    # Definindo a terceira camada linear com 128 neurônios de entrada (da camada anterior) e 10 neurônios de saída (equivalente a 10 classes consideradas para a classifição dos pixels nas images 28x28 originais):
    self$linear3 <- nn_linear(in_features = 128, out_features = 10)
    
    # Definindo uma camada de `dropout` com 40% dos neurônios sendo descartados:
    self$drop1 <- nn_dropout(p = 0.4)
    
    # Definindo uma segunda camada de dropout com 30% dos neurônios sendo descartados:
    self$drop2 <- nn_dropout(p = 0.3)
    
    # Definindo a função de ativação ReLU, a ser usada entre as camadas lineares:
    self$activation <- nn_relu()
  },
  
  # Função que descreve o que acontece quando o modelo é chamado com uma entrada (`forward pass`):
  forward = function(x) {
    x %>% 
      
      # Passando a entrada pela primeira camada linear (`sef$linear1`) seguida da função de ativação ReLU e `dropout`:
      self$linear1() %>% 
      self$activation() %>% 
      self$drop1() %>% 
      
      # Passando a saída da primeira camada pela segunda camada linear (`sef$linear2`), ativação ReLU e `dropout`:
      self$linear2() %>% 
      self$activation() %>% 
      self$drop2() %>% 
      
      # Passando a saída da segunda camada pela terceira camada linear (`sef$linear3`):
      self$linear3()
  }
)

No script acima, definimos as funções initialize() e forward() do nn_module():

  • Em initialize, especificamos todas as camadas que são usadas no modelo. Por exemplo, self$linear1 <- nn_linear(784, 256) define uma camada densa que vai de \(28\times28=784\) neurônios de entrada para uma camada oculta de 256 neurônios. O modelo terá \(3\) dessas camadas, cada uma diminuindo o número de neurônios de saída. A última terá 10 neurônios de saída, porque cada unidade será associada a uma classe diferente, portanto temos um problema de classificação com 10 classes.

  • Também definimos camadas de dropout usando nn_dropout(). Essas serão usadas para realizar a regularização por dropout.

  • Finalmente, definimos a camada de ativação usando nn_relu(). Note que a saída da última camada (self$linear3()) não é seguida por nenhuma função de ativação explícita na função forward. Ou seja, a saída da última camada não passa por nenhuma função de ativação. No entanto, em redes de classificação multiclasse, os algoritmos tendem a aplicar “internamente” a função Softmax à saída da última camada. (Para problemas de classificação binária, a função de ativação mais comum na camada de saída é a Sigmoid.)

  • Em forward(), definimos a ordem em que essas camadas são chamadas. Chamamos elas em blocos (como linear1() + activation() + drop1()), exceto pela última camada, que não utiliza uma função de ativação ou dropout.

Por fim, usamos print para resumir o modelo e garantir que fizemos tudo certo (tente entender de onde vem os números de parâmetros listados abaixo):

print(modelnn())
An `nn_module` containing 235,146 parameters.

── Modules ─────────────────────────────────────────────────────────────────────
• linear1: <nn_linear> #200,960 parameters
• linear2: <nn_linear> #32,896 parameters
• linear3: <nn_linear> #1,290 parameters
• drop1: <nn_dropout> #0 parameters
• drop2: <nn_dropout> #0 parameters
• activation: <nn_relu> #0 parameters

Os parâmetros para cada camada incluem um termo de “viés” (bias), o que resulta em um total de 235.146 parâmetros. Por exemplo, a primeira camada oculta envolve \((784+1)\times 256=200.960\) parâmetros.

Em seguida, adicionamos detalhes ao modelo para especificar seu algoritmo de ajuste: usaremos nn_cross_entropy_loss() dessa vez. A função de perda nn_cross_entropy_loss() é especificamente utilizada para problemas de classificação multiclasse. Curiosidade: a função de entropia cruzada (nn_cross_entropy_loss()) não exige que o alvo (variável dependente) seja codificado em one-hot encoding (se você não está familiarizado com esse termo, faça uma breve pesquisa):

modelnn <- modelnn %>% 
  setup(
    loss = nn_cross_entropy_loss(),
    optimizer = optim_rmsprop, 
    metrics = list(luz_metric_accuracy())
  )

Agora estamos prontos para começar. O passo final é fornecer os dados de treinamento e ajustar o modelo (aumente o número de epochs para melhorar o ajuste):

system.time(
   fitted <- modelnn %>%
      fit(
        data = train_ds, 
        epochs = 5, 
        valid_data = 0.2,
        dataloader_options = list(batch_size = 256),
        verbose = TRUE
      )
 )
Epoch 1/5
Train metrics: Loss: 1.9187 - Acc: 0.7463
Valid metrics: Loss: 0.2833 - Acc: 0.9227
Epoch 2/5
Train metrics: Loss: 0.4122 - Acc: 0.8826
Valid metrics: Loss: 0.2287 - Acc: 0.9344
Epoch 3/5
Train metrics: Loss: 0.345 - Acc: 0.8992
Valid metrics: Loss: 0.2165 - Acc: 0.9333
Epoch 4/5
Train metrics: Loss: 0.3105 - Acc: 0.9122
Valid metrics: Loss: 0.1706 - Acc: 0.9504
Epoch 5/5
Train metrics: Loss: 0.2842 - Acc: 0.9207
Valid metrics: Loss: 0.1847 - Acc: 0.9476
   user  system elapsed 
 234.61  373.79  114.93 
plot(fitted)

A saída é um relatório de progresso sobre o ajuste do modelo, agrupado por época. Isso é muito útil, já que em conjuntos de dados grandes o processo de ajuste pode levar bastante tempo.

Aqui, especificamos uma divisão de dados de validação de \(20\%\). Ou seja, o treinamento da rede foi realizado em \(80\%\) das \(60.000\) imagens disponíveis nos dados de treinamento. Esta é uma alternativa ao fornecimento “real” de dados de validação (veja ?fit.luz_module_generator para todos os argumentos opcionais de ajuste). O SGD usou lotes (batches) de \(256\) observações para calcular o gradiente. Ou seja, fazendo as contas, vemos que uma epoch corresponde a \(188\) iterações.

Para obter o erro da rede baseado nos dados de teste, primeiro escrevemos uma função simples accuracy() que compara as classes previstas e reais, e depois a usamos para avaliar nossas previsões:

# Define a função `accuracy` que recebe duas entradas:
# `pred`: as previsões feitas pelo modelo (valores previstos)
# `truth`: os valores reais ou verdadeiros (valores de classe verdadeira)
accuracy <- function(pred, truth) {
   # Compara cada valor de 'pred' com o valor correspondente em 'truth' e retorna a média de acertos, ou seja, a proporção de previsões corretas:
   mean(pred == truth)
}

# obtém as classes verdadeiras de todas as observações em `test_ds`: 
truth <- sapply(seq_along(test_ds), function(x) test_ds[x][[2]])

fitted %>% 
  predict(test_ds) %>% 
  torch_argmax(dim = 2) %>%  # a classe escolhida é a que tem a maior probabilidade segundo o output da rede
  as_array() %>% # convertemos para um objeto R
  accuracy(truth)
[1] 0.9495

Algo interessante sobre redes neurais é que podemos utilizá-las para ajustar modelos lineares e logísticos também, basicamente omitindo as camadas ocultas (hidden layers).

Por exemplo, embora pacotes como o glmnet consigam lidar com regressão logística multiclasse, eles podem ser bastante lentos com um grande conjunto de dados. É muito mais rápido ajustar tal regressão usando redes neurais. Nesse caso, nós apenas precisamos de apenas \(1\) camada de entrada e \(1\) camada de saída, sem camadas ocultas:

modellr <- nn_module(
  initialize = function() {
    self$linear <- nn_linear(784, 10)
  },
  forward = function(x) {
    self$linear(x)
  }
)
print(modellr())
An `nn_module` containing 7,850 parameters.

── Modules ─────────────────────────────────────────────────────────────────────
• linear: <nn_linear> #7,850 parameters

Ajustamos o modelo da mesma forma que antes (aumente o número de epochs para melhorar o ajuste):

fit_modellr <- modellr %>%
  setup(
    loss = nn_cross_entropy_loss(), # usada para classificações
    optimizer = optim_rmsprop,
    metrics = list(luz_metric_accuracy())
  ) %>%
  fit(
    data = train_ds,
    epochs = 5,
    valid_data = 0.2,
    dataloader_options = list(batch_size = 128)
  )

fit_modellr %>%
  predict(test_ds) %>%
  torch_argmax(dim = 2) %>%  
  as_array() %>% 
  accuracy(truth)
[1] 0.9237
# alternativamente, pode-se usar a função `evaluate` para obter os resultados no test_ds (muito mais fácil!)
evaluate(fit_modellr, test_ds)
A `luz_module_evaluation`
── Results ─────────────────────────────────────────────────────────────────────
loss: 0.2941
acc: 0.9237

Parabéns, você sobreviveu a este tutorial complicadíssimo sobre redes neurais em R com o pacote torch! Agora você pode começar o Tutorial 4.B: Redes Neurais Artificiais com o pacote Keras (Python)—não se preocupe, ele é bem menos complexo que este!