"""This tutorial introduces Contractive auto-encoders (cA) using Theano.

They are based on auto-encoders as the ones used in Bengio et

al. 2007. An autoencoder takes an input x and first maps it to a

hidden representation y = f_{\theta}(x) = s(Wx+b), parameterized by

\theta={W,b}. The resulting latent representation y is then mapped

back to a "reconstructed" vector z \in [0,1]^d in input space z =

g_{\theta'}(y) = s(W'y + b'). The weight matrix W' can optionally be

constrained such that W' = W^T, in which case the autoencoder is said

to have tied weights. The network is trained such that to minimize

the reconstruction error (the error between x and z). Adding the

squared Frobenius norm of the Jacobian of the hidden mapping h with

respect to the visible units yields the contractive auto-encoder:

- \sum_{k=1}^d[ x_k \log z_k + (1-x_k) \log( 1-z_k)]

+ \| \frac{\partial h(x)}{\partial x} \|^2

References :

- S. Rifai, P. Vincent, X. Muller, X. Glorot, Y. Bengio: Contractive

Auto-Encoders: Explicit Invariance During Feature Extraction, ICML-11

- S. Rifai, X. Muller, X. Glorot, G. Mesnil, Y. Bengio, and Pascal

Vincent. Learning invariant features through local space

contraction. Technical Report 1360, Universite de Montreal

- Y. Bengio, P. Lamblin, D. Popovici, H. Larochelle: Greedy Layer-Wise

Training of Deep Networks, Advances in Neural Information Processing

Systems 19, 2007

"""

from __future__ import print_function

import os

import sys

import timeit

import numpy

import theano

import theano.tensor as T

from logistic_sgd import load_data

from utils import tile_raster_images

try:

import PIL.Image as Image

except ImportError:

import Image

class cA(object):

""" Contractive Auto-Encoder class (cA)

The contractive autoencoder tries to reconstruct the input with an

additional constraint on the latent space. With the objective of

obtaining a robust representation of the input space, we

regularize the L2 norm(Froebenius) of the jacobian of the hidden

representation with respect to the input. Please refer to Rifai et

al.,2011 for more details.

If x is the input then equation (1) computes the projection of the

input into the latent space h. Equation (2) computes the jacobian

of h with respect to x. Equation (3) computes the reconstruction

of the input, while equation (4) computes the reconstruction

error and the added regularization term from Eq.(2).

.. math::

h_i = s(W_i x + b_i) (1)

J_i = h_i (1 - h_i) * W_i (2)

x' = s(W' h + b') (3)

L = -sum_{k=1}^d [x_k \log x'_k + (1-x_k) \log( 1-x'_k)]

+ lambda * sum_{i=1}^d sum_{j=1}^n J_{ij}^2 (4)

"""

def __init__(self, numpy_rng, input=None, n_visible=784, n_hidden=100,

n_batchsize=1, W=None, bhid=None, bvis=None):

"""Initialize the cA class by specifying the number of visible units

(the dimension d of the input), the number of hidden units (the

dimension d' of the latent or hidden space) and the contraction level.

The constructor also receives symbolic variables for the input, weights

and bias.

:type numpy_rng: numpy.random.RandomState

:param numpy_rng: number random generator used to generate weights

:type theano_rng: theano.tensor.shared_randomstreams.RandomStreams

:param theano_rng: Theano random generator; if None is given

one is generated based on a seed drawn from `rng`

:type input: theano.tensor.TensorType

:param input: a symbolic description of the input or None for

standalone cA

:type n_visible: int

:param n_visible: number of visible units

:type n_hidden: int

:param n_hidden: number of hidden units

:type n_batchsize int

:param n_batchsize: number of examples per batch

:type W: theano.tensor.TensorType

:param W: Theano variable pointing to a set of weights that should be

shared belong the dA and another architecture; if dA should

be standalone set this to None

:type bhid: theano.tensor.TensorType

:param bhid: Theano variable pointing to a set of biases values (for

hidden units) that should be shared belong dA and another

architecture; if dA should be standalone set this to None

:type bvis: theano.tensor.TensorType

:param bvis: Theano variable pointing to a set of biases values (for

visible units) that should be shared belong dA and another

architecture; if dA should be standalone set this to None

"""

self.n_visible = n_visible

self.n_hidden = n_hidden

self.n_batchsize = n_batchsize

# note : W' was written as `W_prime` and b' as `b_prime`

if not W:

# W is initialized with `initial_W` which is uniformely sampled

# from -4*sqrt(6./(n_visible+n_hidden)) and

# 4*sqrt(6./(n_hidden+n_visible))the output of uniform if

# converted using asarray to dtype

# theano.config.floatX so that the code is runable on GPU

initial_W = numpy.asarray(

numpy_rng.uniform(

low=-4 * numpy.sqrt(6. / (n_hidden + n_visible)),

high=4 * numpy.sqrt(6. / (n_hidden + n_visible)),

size=(n_visible, n_hidden)

),

dtype=theano.config.floatX

)

W = theano.shared(value=initial_W, name='W', borrow=True)

if not bvis:

bvis = theano.shared(value=numpy.zeros(n_visible,

dtype=theano.config.floatX),

borrow=True)

if not bhid:

bhid = theano.shared(value=numpy.zeros(n_hidden,

dtype=theano.config.floatX),

name='b',

borrow=True)

self.W = W

# b corresponds to the bias of the hidden

self.b = bhid

# b_prime corresponds to the bias of the visible

self.b_prime = bvis

# tied weights, therefore W_prime is W transpose

self.W_prime = self.W.T

# if no input is given, generate a variable representing the input

if input is None:

# we use a matrix because we expect a minibatch of several

# examples, each example being a row

self.x = T.dmatrix(name='input')

else:

self.x = input

self.params = [self.W, self.b, self.b_prime]

def get_hidden_values(self, input):

""" Computes the values of the hidden layer """

return T.nnet.sigmoid(T.dot(input, self.W) + self.b)

def get_jacobian(self, hidden, W):

"""Computes the jacobian of the hidden layer with respect to

the input, reshapes are necessary for broadcasting the

element-wise product on the right axis

"""

return T.reshape(hidden * (1 - hidden),

(self.n_batchsize, 1, self.n_hidden)) * T.reshape(

W, (1, self.n_visible, self.n_hidden))

def get_reconstructed_input(self, hidden):

"""Computes the reconstructed input given the values of the

hidden layer

"""

return T.nnet.sigmoid(T.dot(hidden, self.W_prime) + self.b_prime)

def get_cost_updates(self, contraction_level, learning_rate):

""" This function computes the cost and the updates for one trainng

step of the cA """

y = self.get_hidden_values(self.x)

z = self.get_reconstructed_input(y)

J = self.get_jacobian(y, self.W)

# note : we sum over the size of a datapoint; if we are using

# minibatches, L will be a vector, with one entry per

# example in minibatch

self.L_rec = - T.sum(self.x * T.log(z) +

(1 - self.x) * T.log(1 - z),

axis=1)

# Compute the jacobian and average over the number of samples/minibatch

self.L_jacob = T.sum(J ** 2) // self.n_batchsize

# note : L is now a vector, where each element is the

# cross-entropy cost of the reconstruction of the

# corresponding example of the minibatch. We need to

# compute the average of all these to get the cost of

# the minibatch

cost = T.mean(self.L_rec) + contraction_level * T.mean(self.L_jacob)

# compute the gradients of the cost of the `cA` with respect

# to its parameters

gparams = T.grad(cost, self.params)

# generate the list of updates

updates = []

for param, gparam in zip(self.params, gparams):

updates.append((param, param - learning_rate * gparam))

return (cost, updates)

def test_cA(learning_rate=0.01, training_epochs=20,

dataset='mnist.pkl.gz',

batch_size=10, output_folder='cA_plots', contraction_level=.1):

"""

This demo is tested on MNIST

:type learning_rate: float

:param learning_rate: learning rate used for training the contracting

AutoEncoder

:type training_epochs: int

:param training_epochs: number of epochs used for training

:type dataset: string

:param dataset: path to the picked dataset

"""

datasets = load_data(dataset)

train_set_x, train_set_y = datasets[0]

# compute number of minibatches for training, validation and testing

n_train_batches = train_set_x.get_value(borrow=True).shape[0] // batch_size

# allocate symbolic variables for the data

index = T.lscalar() # index to a [mini]batch

x = T.matrix('x') # the data is presented as rasterized images

if not os.path.isdir(output_folder):

os.makedirs(output_folder)

os.chdir(output_folder)

####################################

# BUILDING THE MODEL #

####################################

rng = numpy.random.RandomState(123)

ca = cA(numpy_rng=rng, input=x,

n_visible=28 * 28, n_hidden=500, n_batchsize=batch_size)

cost, updates = ca.get_cost_updates(contraction_level=contraction_level,

learning_rate=learning_rate)

train_ca = theano.function(

[index],

[T.mean(ca.L_rec), ca.L_jacob],

updates=updates,

givens={

x: train_set_x[index * batch_size: (index + 1) * batch_size]

}

)

start_time = timeit.default_timer()

############

# TRAINING #

############

# go through training epochs

for epoch in range(training_epochs):

# go through trainng set

c = []

for batch_index in range(n_train_batches):

c.append(train_ca(batch_index))

c_array = numpy.vstack(c)

print('Training epoch %d, reconstruction cost ' % epoch, numpy.mean(

c_array[0]), ' jacobian norm ', numpy.mean(numpy.sqrt(c_array[1])))

end_time = timeit.default_timer()

training_time = (end_time - start_time)

print(('The code for file ' + os.path.split(__file__)[1] +

' ran for %.2fm' % ((training_time) / 60.)), file=sys.stderr)

image = Image.fromarray(tile_raster_images(

X=ca.W.get_value(borrow=True).T,

img_shape=(28, 28), tile_shape=(10, 10),

tile_spacing=(1, 1)))

image.save('cae_filters.png')

os.chdir('../')

if __name__ == '__main__':

test_cA()