modelos_gerativos

aula de hoje:

comparação de um autoencoder com PCA

baseada em: https://www.r-bloggers.com/pca-vs-autoencoders-for-dimensionality-reduction/

variational autoencoders (VAE)

referência: Deep Learning in R, Chollet & Allaire, 8.4-generating-images-with-vaes.Rmd

1. comparação de um autoencoder com PCA

Autoencoders (AE) são um tipo de rede neural que podem ser usadas tanto para redução de dimensionalidade quanto como componentes de redes gerativas.

Um AE objetiva reproduzir seu input mas, para conseguir reduzir a dimensionalidade do input, ele tem uma “gargalo” no meio.

Essa gargalo tem dimensão menor que o input e assim, comprime a informação do input em um número menor de componentes: as “variáveis latentes”.

Um AE tem duas partes:

-o encoder - que transforma o input em variáveis latentes

-o decoder - que usa as variáveis latentes para reconstruir o input

Um AE poder ser treinado com gradiente descendente, minimizando o erro quadrático médio.

Vamos comparar a capacidade de redução de dimensionalidade do PCA com um AE:

# inicialização
t00 = Sys.time()     # para o tempo de processamento desse script
set.seed(1234)     # para reprodutibilidade

# carregando as bibliotecas:
library(keras)
library(ggplot2)
library(scatterplot3d)

os dados:

vamos analisar 1000 espectros ópticos sintéticos de galáxias, amostrados entre 3500 e 8500 A de 50 em 50 A

espectros = read.table(file='spec1kd50.txt',header=FALSE)
lambda = read.table(file='lambda1kd50.txt',header=FALSE)

# vamos ilustrar com 4 espectros
par(mfcol=c(2,2)) 
x = lambda[,1]
y = espectros[1,]
plot(x,y,type='l',xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo normalizado',col='red')
y = espectros[2,]
plot(x,y,type='l',xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo normalizado',col='red')
y = espectros[3,]
plot(x,y,type='l',xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo normalizado',col='red')
y = espectros[4,]
plot(x,y,type='l',xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo normalizado',col='red')

Principal Components Analysis

PCA descreve o espaço dos dados com um sistema de coordenadas ortogonais em direções de maior variância, as PCs, e reduz a dimensionalidade representando os dados num sub-espaço definido pelas PCs que cobrem a maior parte da variância.

pré-processamento dos dados:

# pré-processamento: vamos colocar os dados (fluxos em cada comprimento de onda) entre 0 e 1
minmax <- function(x) (x - min(x))/(max(x) - min(x))
x_train <- apply(espectros, 2, minmax)

PCA:

t0 = Sys.time()  
pca <- prcomp(x_train)

# fração da variância cumulativa em função do número de componentes
par(mfcol=c(1,1))
x = 1:dim(x_train)[2]
y = cumsum(pca$sdev)/sum(pca$sdev)
qplot(x[1:10], y[1:10], geom = "line")

# contribuição à variância das 10 primeiras componentes
y[1:10]

##  [1] 0.6602135 0.7627732 0.8445692 0.8975840 0.9207894 0.9399514 0.9540042
##  [8] 0.9614353 0.9683823 0.9735721

#A figura mostra que 4 componentes explicam ~90% da variância.

# visualização de PC1 x PC2:
ggplot(as.data.frame(pca$x), aes(x = PC1, y = PC2)) + geom_point()

# visualização de PC1 x PC2 X PC3:
scatterplot3d(pca$x[,1],pca$x[,2],pca$x[,3],main="PCA",asp=1)

Sys.time()-t0

## Time difference of 1.226029 secs

Vamos agora examinar como se compara a compressão do PCA com a de um AE usando o keras:

t0 = Sys.time()  
# preparando os dados:
x_train <- as.matrix(x_train)

# modelo do AE: camadas densas
# ativação tanh: as variáveis estão entre -1 e 1

# espaço latente com 2 variáveis
model <- keras_model_sequential()
model %>%
  layer_dense(units = 6, activation = "tanh", input_shape = ncol(x_train)) %>%
  # vamos dar um nome à camada gargalo:
  layer_dense(units = 2, activation = "tanh", name = "bottleneck") %>%
  layer_dense(units = 6, activation = "tanh") %>%
  layer_dense(units = ncol(x_train))

summary(model)

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #     
## ===========================================================================
## dense (Dense)                    (None, 6)                     546         
## ___________________________________________________________________________
## bottleneck (Dense)               (None, 2)                     14          
## ___________________________________________________________________________
## dense_1 (Dense)                  (None, 6)                     18          
## ___________________________________________________________________________
## dense_2 (Dense)                  (None, 90)                    630         
## ===========================================================================
## Total params: 1,208
## Trainable params: 1,208
## Non-trainable params: 0
## ___________________________________________________________________________

# compilação:
model %>% compile(
  loss = "mean_squared_error", 
  optimizer = "adam"               #um método de otimização estocástica
)

# treinamento (fit)
model %>% fit(
  x = x_train, 
  y = x_train, 
  epochs = 2000,
  verbose = 0
)

# resultado
mse.ae2 <- evaluate(model, x_train, x_train)
mse.ae2

##         loss 
## 0.0003825233

# vamos extrair as variáveis latentes da gargalo:
intermediate_layer_model <- keras_model(inputs = model$input, outputs = get_layer(model, "bottleneck")$output)
intermediate_output <- predict(intermediate_layer_model, x_train)

ggplot(data.frame(AE1 = intermediate_output[,1], AE2 = intermediate_output[,2]), aes(x = AE1, y = AE2)) + geom_point()

#Vamos agora considerar o espaço latente com 3 variáveis

model3 <- keras_model_sequential()
model3 %>%
  layer_dense(units = 6, activation = "tanh", input_shape = ncol(x_train)) %>%
  layer_dense(units = 3, activation = "tanh", name = "bottleneck") %>%
  layer_dense(units = 6, activation = "tanh") %>%
  layer_dense(units = ncol(x_train))


summary(model3)

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #     
## ===========================================================================
## dense_3 (Dense)                  (None, 6)                     546         
## ___________________________________________________________________________
## bottleneck (Dense)               (None, 3)                     21          
## ___________________________________________________________________________
## dense_4 (Dense)                  (None, 6)                     24          
## ___________________________________________________________________________
## dense_5 (Dense)                  (None, 90)                    630         
## ===========================================================================
## Total params: 1,221
## Trainable params: 1,221
## Non-trainable params: 0
## ___________________________________________________________________________

# compilação:
model3 %>% compile(
  loss = "mean_squared_error", 
  optimizer = "adam"
)

# treinamento
model3 %>% fit(
  x = x_train, 
  y = x_train, 
  epochs = 2000,
  verbose = 0
)

# resultado
evaluate(model3, x_train, x_train)

##         loss 
## 0.0001703455

# gargalo
intermediate_layer_model <- keras_model(inputs = model3$input, outputs = get_layer(model3, "bottleneck")$output)
intermediate_output <- predict(intermediate_layer_model, x_train)


# plot das variáveis latentes
ae3 <- data.frame(node1 = intermediate_output[,1], node2 = intermediate_output[,2], node3 = intermediate_output[,3])

scatterplot3d(intermediate_output[,1], intermediate_output[,2], intermediate_output[,3], main="AE",asp=1)

Sys.time()-t0

## Time difference of 2.642732 mins

Vamos agora comparar o erro na reconstrução = média da diferença quadrática entre o input e o output

reconstrução via PCA:

t0 = Sys.time()  
pca.recon <- function(pca, x, k){
  mu <- matrix(rep(pca$center, nrow(pca$x)), nrow = nrow(pca$x), byrow = T)
  recon <- pca$x[,1:k] %*% t(pca$rotation[,1:k]) + mu
  mse <- mean((recon - x)^2)
  return(list(x = recon, mse = mse))
}

# erro
xhat <- rep(NA, 10)
for(k in 1:10){
  xhat[k] <- pca.recon(pca, x_train, k)$mse
}
# erros até a componente 5:
xhat[1:5]

## [1] 1.342324e-03 6.794064e-04 2.577391e-04 8.060598e-05 4.666837e-05

Sys.time()-t0

## Time difference of 0.1454558 secs

reconstrução via AE (até 5 variáveis latentes):

t0 = Sys.time()  
ae.mse <- rep(NA, 5)
for(k in 1:5){
  modelk <- keras_model_sequential()
  modelk %>%
    layer_dense(units = 6, activation = "tanh", input_shape = ncol(x_train)) %>%
    layer_dense(units = k, activation = "tanh", name = "bottleneck") %>%
    layer_dense(units = 6, activation = "tanh") %>%
    layer_dense(units = ncol(x_train))

  modelk %>% compile(
    loss = "mean_squared_error", 
    optimizer = "adam"
  )

  modelk %>% fit(
    x = x_train, 
    y = x_train, 
    epochs = 5000,
    verbose = 0
  )

# erro:
  ae.mse[k] <- unname(evaluate(modelk, x_train, x_train))
}
# erros até a componente 5:
ae.mse[1:5]

## [1] 1.134849e-03 3.636937e-04 1.630001e-04 8.191094e-05 4.481221e-05

t0 = Sys.time()

comparação:

df <- data.frame(k = c(1:10, 1:5), mse = c(xhat, ae.mse), method = c(rep("pca", 10), rep("autoencoder", 5)))
ggplot(df, aes(x = k, y = mse, col = method)) + geom_line()

Vemos que o AE é melhor que o PCA quando o número de componentes é pequeno, e converge para o mesmo resultado que o PCA conforme o número de componentes aumenta.

2. VAE: variational autoencoders

VAE são uma variante mais poderosa de AE:

em vez de compactar o input em um “código” fixo no espaço latente, ela transforma o input nos parâmetros de uma distribuição estatística: uma média e uma variância.

Autoencoder

VAE

isso significa que estamos assumindo que a imagem de entrada foi gerada por um processo estatístico e que a aleatoriedade desse processo deve ser levada em consideração durante a codificação e decodificação.

O VAE usa os parâmetros de média e variância para amostrar aleatoriamente um elemento da distribuição e decodifica esse elemento tentando obter o input original.

O caráter estocástico deste processo melhora a robustez e força o espaço latente a codificar representações significativas.

Eis um resumo de como um VAE funciona:

o encoder transforma o input em dois parâmetros no espaço latente, que designaremos como “z_mean” e “z_log_variance”.

amostramos um ponto “z” de uma distribuição normal latente via z = z_mean + exp(z_log_variance) * epsilon onde epsilon é um tensor aleatório com valores pequenos. Como epsilon é pequeno, todo ponto próximo no espaço latente vai ser decodificado próximo ao input

o decoder produz um output que busca replicar o input.

Uma VAE usa duas funções de custo:

custo de reconstrução- força o output a reproduzir o input

custo de regularização: ajuda a estruturar o espaço latente

Esquematicamente, uma VAE no Keras tem essa cara:

codifica o input com parâmetros de média e variância:

c(z_mean, z_log_variance) %<% encoder(input_img)

amostra um ponto latente com epsilon pequeno:

z <- z_mean + exp(z_log_variance) * epsilon

decodifica z de volta para uma imagem:

reconstructed_img <- decoder(z)

cria um modelo:

model <- keras_model(input_img, reconstructed_img)

treina o modelo com as duas funções de custo, de reconstrução e de regularização

encoder: uma convnet que mapeia o input em 2 vetores: z_mean e z_log_variance

# inicialização
t00 = Sys.time()     # para o tempo de processamento desse script
set.seed(1234)     # para reprodutibilidade

library(keras)

img_shape <- c(28, 28, 1)
batch_size <- 16
latent_dim <- 2L  # Dimensionality of the latent space: a plane

input_img <- layer_input(shape = img_shape)

x <- input_img %>% 
  layer_conv_2d(filters = 32, kernel_size = 3, padding = "same", 
                activation = "relu") %>% 
  layer_conv_2d(filters = 64, kernel_size = 3, padding = "same", 
                activation = "relu", strides = c(2, 2)) %>%
  layer_conv_2d(filters = 64, kernel_size = 3, padding = "same", 
                activation = "relu") %>%
  layer_conv_2d(filters = 64, kernel_size = 3, padding = "same", 
                activation = "relu") 

shape_before_flattening <- k_int_shape(x)

x <- x %>% 
  layer_flatten() %>% 
  layer_dense(units = 32, activation = "relu")

z_mean <- x %>% 
  layer_dense(units = latent_dim)

z_log_var <- x %>% 
  layer_dense(units = latent_dim)

amostra um ponto do espaço latente, z:

sampling <- function(args) {
  c(z_mean, z_log_var) %<-% args
  epsilon <- k_random_normal(shape = list(k_shape(z_mean)[1], latent_dim),
                             mean = 0, stddev = 1)
  z_mean + k_exp(z_log_var) * epsilon
}

z <- list(z_mean, z_log_var) %>% 
  layer_lambda(sampling)

decoder: reformata-se o vetor z como uma imagem que entra numa convnet com saída com as mesmas dimensões da imagem de input

# This is the input where we will feed `z`.
decoder_input <- layer_input(k_int_shape(z)[-1])

x <- decoder_input %>% 
  # Upsample to the correct number of units
  layer_dense(units = prod(as.integer(shape_before_flattening[-1])),
              activation = "relu") %>% 
  # Reshapes into an image of the same shape as before the last flatten layer
  layer_reshape(target_shape = shape_before_flattening[-1]) %>% 
  # Applies and then reverses the operation to the initial stack of 
  # convolution layers
  layer_conv_2d_transpose(filters = 32, kernel_size = 3, padding = "same",
                          activation = "relu", strides = c(2, 2)) %>%  
  layer_conv_2d(filters = 1, kernel_size = 3, padding = "same",
                activation = "sigmoid")  
  # We end up with a feature map of the same size as the original input.

# This is our decoder model.
decoder <- keras_model(decoder_input, x)

# We then apply it to `z` to recover the decoded `z`.
z_decoded <- decoder(z)

Aqui define-se uma função para implementar as duas funções de custo:

library(R6)

CustomVariationalLayer <- R6Class("CustomVariationalLayer",
                                  
  inherit = KerasLayer,
  
  public = list(
    
    vae_loss = function(x, z_decoded) {
      x <- k_flatten(x)
      z_decoded <- k_flatten(z_decoded)
      xent_loss <- metric_binary_crossentropy(x, z_decoded)
      kl_loss <- -5e-4 * k_mean(
        1 + z_log_var - k_square(z_mean) - k_exp(z_log_var), 
        axis = -1L
      )
      k_mean(xent_loss + kl_loss)
    },
    
    call = function(inputs, mask = NULL) {
      x <- inputs[[1]]
      z_decoded <- inputs[[2]]
      loss <- self$vae_loss(x, z_decoded)
      self$add_loss(loss, inputs = inputs)
      x
    }
  )
)

layer_variational <- function(object) { 
  create_layer(CustomVariationalLayer, object, list())
} 

# Call the custom layer on the input and the decoded output to obtain
# the final model output
y <- list(input_img, z_decoded) %>% 
  layer_variational()

inicializa e treina

não se passa ‘target’, só x_train (que é o target)

vae <- keras_model(input_img, y)

vae %>% compile(
  optimizer = "rmsprop",
  loss = NULL
)

# Trains the VAE on MNIST digits
mnist <- dataset_mnist() 
c(c(x_train, y_train), c(x_test, y_test)) %<-% mnist

x_train <- x_train / 255
x_train <- array_reshape(x_train, dim =c(dim(x_train), 1))

x_test <- x_test / 255
x_test <- array_reshape(x_test, dim =c(dim(x_test), 1))

vae %>% fit(
  x = x_train, y = NULL,
  epochs = 10,
  batch_size = batch_size,
  validation_data = list(x_test, NULL)
)

agora que o modelo está treinado, pode-se usar o decoder para transformar vetores arbitrários do espaço latente em imagens:

n <- 15            # Number of rows / columns of digits
digit_size <- 28   # Height / width of digits in pixels

# Transforms linearly spaced coordinates on the unit square through the inverse
# CDF (ppf) of the Gaussian to produce values of the latent variables z,
# because the prior of the latent space is Gaussian
grid_x <- qnorm(seq(0.05, 0.95, length.out = n))
grid_y <- qnorm(seq(0.05, 0.95, length.out = n))

op <- par(mfrow = c(n, n), mar = c(0,0,0,0), bg = "black")
for (i in 1:length(grid_x)) {
  yi <- grid_x[[i]]
  for (j in 1:length(grid_y)) {
    xi <- grid_y[[j]]
    z_sample <- matrix(c(xi, yi), nrow = 1, ncol = 2)
    z_sample <- t(replicate(batch_size, z_sample, simplify = "matrix"))
    x_decoded <- decoder %>% predict(z_sample, batch_size = batch_size)
    digit <- array_reshape(x_decoded[1,,,], dim = c(digit_size, digit_size))
    plot(as.raster(digit))
  }
}

par(op)

Exercício:

Verifique o desempenho de um autoencoder com uma camada densa a mais no encoder e no decoder. Considere que esta camada tem 6 unidades e função de ativação relu.

Sys.time()-t00

## Time difference of 30.29975 mins

modelos_gerativos

lsj

November 9, 2019

aula de hoje:

1. comparação de um autoencoder com PCA

2. VAE: variational autoencoders

Exercício: