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/).
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"
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
Há \(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!