Construyendo un autoencoder en Tensorflow para superar PCA

Hinton y Salakhutdinov en la reducción de la dimensionalidad de los datos con redes neuronales, Science 2006 propusieron un PCA no lineal mediante el uso de un autoencoder profundo. He intentado construir y entrenar un codificador automático de PCA con Tensorflow varias veces, pero nunca he podido obtener un mejor resultado que el PCA lineal.

¿Cómo puedo entrenar eficientemente un autoencoder?

(Edición posterior de @amoeba: la versión original de esta pregunta contenía código de Python Tensorflow que no funcionaba correctamente. Se puede encontrar en el historial de edición).

Aquí está la figura clave del artículo de Ciencia de 2006 de Hinton y Salakhutdinov:

Muestra la reducción de dimensionalidad del conjunto de datos MNIST ( $28\times 28$ imágenes en blanco y negro de un solo dígito) desde las dimensiones originales del 784 a dos.

El código está copiado de un cuaderno Jupyter. En Python 3.6 necesita instalar matplotlib (para pylab), NumPy, seaborn, TensorFlow y Keras. Cuando se ejecuta en Python Shell, es posible que deba agregar plt.show()para mostrar los gráficos.

Inicialización

%matplotlib notebook

import pylab as plt
import numpy as np
import seaborn as sns; sns.set()

import keras
from keras.datasets import mnist
from keras.models import Sequential, Model
from keras.layers import Dense
from keras.optimizers import Adam

(x_train, y_train), (x_test, y_test) = mnist.load_data()
x_train = x_train.reshape(60000, 784) / 255
x_test = x_test.reshape(10000, 784) / 255

PCA

mu = x_train.mean(axis=0)
U,s,V = np.linalg.svd(x_train - mu, full_matrices=False)
Zpca = np.dot(x_train - mu, V.transpose())

Rpca = np.dot(Zpca[:,:2], V[:2,:]) + mu    # reconstruction
err = np.sum((x_train-Rpca)**2)/Rpca.shape[0]/Rpca.shape[1]
print('PCA reconstruction error with 2 PCs: ' + str(round(err,3)));

Esto produce:

PCA reconstruction error with 2 PCs: 0.056

Entrenando al autoencoder

m = Sequential()
m.add(Dense(512,  activation='elu', input_shape=(784,)))
m.add(Dense(128,  activation='elu'))
m.add(Dense(2,    activation='linear', name="bottleneck"))
m.add(Dense(128,  activation='elu'))
m.add(Dense(512,  activation='elu'))
m.add(Dense(784,  activation='sigmoid'))
m.compile(loss='mean_squared_error', optimizer = Adam())
history = m.fit(x_train, x_train, batch_size=128, epochs=5, verbose=1, 
                validation_data=(x_test, x_test))

encoder = Model(m.input, m.get_layer('bottleneck').output)
Zenc = encoder.predict(x_train)  # bottleneck representation
Renc = m.predict(x_train)        # reconstruction

Esto lleva ~ 35 segundos en mi escritorio de trabajo y resultados:

Train on 60000 samples, validate on 10000 samples
Epoch 1/5
60000/60000 [==============================] - 7s - loss: 0.0577 - val_loss: 0.0482
Epoch 2/5
60000/60000 [==============================] - 7s - loss: 0.0464 - val_loss: 0.0448
Epoch 3/5
60000/60000 [==============================] - 7s - loss: 0.0438 - val_loss: 0.0430
Epoch 4/5
60000/60000 [==============================] - 7s - loss: 0.0423 - val_loss: 0.0416
Epoch 5/5
60000/60000 [==============================] - 7s - loss: 0.0412 - val_loss: 0.0407

así que ya puedes ver que superamos la pérdida de PCA después de solo dos épocas de entrenamiento.

(Por cierto, es instructivo cambiar todas las funciones de activación activation='linear'y observar cómo la pérdida converge precisamente a la pérdida de PCA. Esto se debe a que el autoencoder lineal es equivalente a PCA).

Trazar la proyección PCA lado a lado con la representación del cuello de botella

plt.figure(figsize=(8,4))
plt.subplot(121)
plt.title('PCA')
plt.scatter(Zpca[:5000,0], Zpca[:5000,1], c=y_train[:5000], s=8, cmap='tab10')
plt.gca().get_xaxis().set_ticklabels([])
plt.gca().get_yaxis().set_ticklabels([])

plt.subplot(122)
plt.title('Autoencoder')
plt.scatter(Zenc[:5000,0], Zenc[:5000,1], c=y_train[:5000], s=8, cmap='tab10')
plt.gca().get_xaxis().set_ticklabels([])
plt.gca().get_yaxis().set_ticklabels([])

plt.tight_layout()

Reconstrucciones

Y ahora veamos las reconstrucciones (primera fila - imágenes originales, segunda fila - PCA, tercera fila - autoencoder):

plt.figure(figsize=(9,3))
toPlot = (x_train, Rpca, Renc)
for i in range(10):
    for j in range(3):
        ax = plt.subplot(3, 10, 10*j+i+1)
        plt.imshow(toPlot[j][i,:].reshape(28,28), interpolation="nearest", 
                   vmin=0, vmax=1)
        plt.gray()
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)

plt.tight_layout()

Se pueden obtener resultados mucho mejores con una red más profunda, cierta regularización y capacitación más larga. Experimentar. ¡El aprendizaje profundo es fácil!

Enormes apoyos a @amoeba por hacer este gran ejemplo. Solo quiero mostrar que el procedimiento de reconstrucción y entrenamiento del codificador automático descrito en esa publicación también se puede hacer en R con similar facilidad. El codificador automático a continuación está configurado para que emule el ejemplo de ameba lo más cerca posible: el mismo optimizador y la arquitectura general. Los costos exactos no son reproducibles debido a que el backend TensorFlow no se sembró de manera similar.

Inicialización

library(keras)
library(rARPACK) # to use SVDS
rm(list=ls())
mnist   = dataset_mnist()
x_train = mnist$train$x
y_train = mnist$train$y
x_test  = mnist$test$x
y_test  = mnist$test$y

# reshape & rescale
dim(x_train) = c(nrow(x_train), 784)
dim(x_test)  = c(nrow(x_test), 784)
x_train = x_train / 255
x_test = x_test / 255

PCA

mus = colMeans(x_train)
x_train_c =  sweep(x_train, 2, mus)
x_test_c =  sweep(x_test, 2, mus)
digitSVDS = svds(x_train_c, k = 2)

ZpcaTEST = x_test_c %*% digitSVDS$v # PCA projection of test data

Autoencoder

model = keras_model_sequential() 
model %>%
  layer_dense(units = 512, activation = 'elu', input_shape = c(784)) %>%  
  layer_dense(units = 128, activation = 'elu') %>%
  layer_dense(units = 2,   activation = 'linear', name = "bottleneck") %>%
  layer_dense(units = 128, activation = 'elu') %>% 
  layer_dense(units = 512, activation = 'elu') %>% 
  layer_dense(units = 784, activation='sigmoid')

model %>% compile(
  loss = loss_mean_squared_error, optimizer = optimizer_adam())

history = model %>% fit(verbose = 2, validation_data = list(x_test, x_test),
                         x_train, x_train, epochs = 5, batch_size = 128)

# Unsurprisingly a 3-year old laptop is slower than a desktop
# Train on 60000 samples, validate on 10000 samples
# Epoch 1/5
#  - 14s - loss: 0.0570 - val_loss: 0.0488
# Epoch 2/5
#  - 15s - loss: 0.0470 - val_loss: 0.0449
# Epoch 3/5
#  - 15s - loss: 0.0439 - val_loss: 0.0426
# Epoch 4/5
#  - 15s - loss: 0.0421 - val_loss: 0.0413
# Epoch 5/5
#  - 14s - loss: 0.0408 - val_loss: 0.0403

# Set the auto-encoder
autoencoder = keras_model(model$input, model$get_layer('bottleneck')$output)
ZencTEST = autoencoder$predict(x_test)  # bottleneck representation  of test data

Trazar la proyección PCA lado a lado con la representación del cuello de botella

par(mfrow=c(1,2))
myCols = colorRampPalette(c('green',     'red',  'blue',  'orange', 'steelblue2',
                            'darkgreen', 'cyan', 'black', 'grey',   'magenta') )
plot(ZpcaTEST[1:5000,], col= myCols(10)[(y_test+1)], 
     pch=16, xlab = 'Score 1', ylab = 'Score 2', main = 'PCA' ) 
legend( 'bottomright', col= myCols(10), legend = seq(0,9, by=1), pch = 16 )

plot(ZencTEST[1:5000,], col= myCols(10)[(y_test+1)], 
     pch=16, xlab = 'Score 1', ylab = 'Score 2', main = 'Autoencoder' ) 
legend( 'bottomleft', col= myCols(10), legend = seq(0,9, by=1), pch = 16 )

Reconstrucciones

Podemos hacer la reconstrucción de los dígitos de la manera habitual. (La fila superior son los dígitos originales, la fila central las reconstrucciones de PCA y la fila inferior las reconstrucciones del codificador automático).

Renc = predict(model, x_test)        # autoencoder reconstruction
Rpca = sweep( ZpcaTEST %*% t(digitSVDS$v), 2, -mus) # PCA reconstruction

dev.off()
par(mfcol=c(3,9), mar = c(1, 1, 0, 0))
myGrays = gray(1:256 / 256)
for(u in seq_len(9) ){
  image( matrix( x_test[u,], 28,28, byrow = TRUE)[,28:1], col = myGrays, 
         xaxt='n', yaxt='n')
  image( matrix( Rpca[u,], 28,28, byrow = TRUE)[,28:1], col = myGrays , 
         xaxt='n', yaxt='n')
  image( matrix( Renc[u,], 28,28, byrow = TRUE)[,28:1], col = myGrays, 
         xaxt='n', yaxt='n')
}

$k$$0.0356$$0.0359$

Aquí está mi cuaderno jupyter donde trato de replicar su resultado, con las siguientes diferencias:

• en lugar de usar tensorflow directamente, lo uso para ver keras
• relu permeable en lugar de relu para evitar la saturación (es decir, la salida codificada es 0)
  • Esto podría ser una razón para el bajo rendimiento de AE
• la entrada del autoencoder es datos escalados a [0,1]
  • Creo que leí en alguna parte que los codificadores automáticos con relu funcionan mejor con datos [0-1]
  • ejecutar mi computadora portátil con la entrada de autoencoders siendo la media = 0, std = 1 dio MSE para AE> 0.7 para todas las reducciones de dimensionalidad, así que quizás este sea uno de sus problemas
• La entrada de PCA se mantiene como datos con media = 0 y estándar = 1
  • Esto también puede significar que el resultado MSE de PCA no es comparable al resultado MSE de PCA
  • Tal vez volveré a ejecutar esto más tarde con [0-1] datos para PCA y AE
• La entrada PCA también se escala a [0-1]. PCA también funciona con datos (media = 0, estándar = 1), pero el MSE sería incomparable a AE

Mi MSE resulta para PCA de reducción de dimensionalidad de 1 a 6 (donde la entrada tiene 6 columnas) y para AE de tenue. rojo. de 1 a 6 están abajo:

Con la entrada PCA (media = 0, estándar = 1) mientras que la entrada AE es [0-1] rango - 4e-15: PCA6 - .015: PCA5 - .0502: AE5 - .0508: AE6 - .051: AE4 - .053: AE3 - .157: PCA4 - .258: AE2 - .259: PCA3 - .377: AE1 - .483: PCA2 - .682: PCA1

• 9e-15: PCA6
• .0094: PCA5
• .0502: AE5
• .0507: AE6
• .0514: AE4
• .0532: AE3
• .0772: PCA4
• .1231: PCA3
• .2588: AE2
• .2831: PCA2
• .3773: AE1
• .3885: PCA1

El PCA lineal sin reducción de dimensionalidad puede alcanzar 9e-15 porque solo puede empujar lo que no pudo encajar en el último componente.