<commit_before>.. _SdA:

Stacked Denoising Autoencoders (SdA)
====================================

.. note::
  This section assumes the reader has already read through :doc:`logreg`
  and :doc:`mlp`. Additionally it uses the following Theano functions
  and concepts : TODO


The Stacked Denoising Autoencoder (SdA) is an extension of the stacked 
autoencoder [Bengio07]_ and it was introduced in [Vincent08]_. We will start the 
tutorial with a short discussion on :ref:`autoencoders`
and then move on to how classical
autoencoders are extended to denoising autoencoders (:ref:`DA`).

.. _autoencoders:

Autoencoders
+++++++++++++

See section 4.6 of [Bengio09]_ for an overview of auto-encoders.
An autoencoder takes an input :math:`\mathbf{x} \in [0,1]^d` and first 
maps it (with an *encoder*) to a hidden representation :math:`\mathbf{y} \in [0,1]^{d'}` 
through a deterministic mapping, e.g.:

.. math::
  
  \mathbf{y} = s(\mathbf{W}\mathbf{x} + \mathbf{b})

Where :math:`s` is a non-linearity such as the sigmoid.
The latent representation :math:`\mathbf{y}`, or **code** is then mapped back (with a *decoder*) into a
**reconstruction** :math:`\mathbf{z}` of same shape as
:math:`\mathbf{x}` through a similar transformation, e.g.:

.. math::

  \mathbf{z} = s(\mathbf{W'}\mathbf{y} + \mathbf{b'})

where ' does not indicate transpose, and
:math:`\mathbf{z}` should be seen as a prediction of :math:`\mathbf{x}`, given the code :math:`\mathbf{y}`.
The weight matrix :math:`\mathbf{W'}` of the reverse mapping may be
optionally constrained by :math:`\mathbf{W'} = \mathbf{W}^T`, which is
an instance of *tied weights*. The parameters of this model (namely 
:math:`\mathbf{W}`, :math:`\mathbf{b}`, 
:math:`\mathbf{b'}` and, if one doesn't use tied weights, also 
:math:`\mathbf{W'}`) are optimized such that the average reconstruction 
error is minimized. The reconstruction error can be measured in many ways, depending
on the appropriate distributional assumptions on the input given the code, e.g., using the 
traditional *squared error* :math:`L(\mathbf{x}, \mathbf{z}) = || \mathbf{x} - \mathbf{z} ||^2`, 
or if the input is interpreted as either bit vectors or vectors of 
bit probabilities by the reconstruction *cross-entropy* defined as : 

.. math::

  L_{H} (\mathbf{x}, \mathbf{z}) = - \sum^d_{k=1}[\mathbf{x}_k \log
          \mathbf{z}_k + (1 - \mathbf{x}_k)\log(1 - \mathbf{z}_k)] 

The hope is that the code :math:`\mathbf{y}` is a distributed representation
that captures the coordinates along the main factors of variation in the data
(similarly to how the projection on principal components captures the main factors
of variation in the data).
Because :math:`\mathbf{y}` is viewed as a lossy compression of :math:`\mathbf{x}`, it cannot
be a good compression (with small loss) for all :math:`\mathbf{x}`, so learning
drives it to be one that is a good compression in particular for training
examples, and hopefully for others as well, but not for arbitrary inputs.
That is the sense in which an auto-encoder generalizes: it gives low reconstruction
error to test examples from the same distribution as the training examples,
but generally high reconstruction error to uniformly chosen configurations of the
input vector.

If there is one linear hidden layer (the code) and
the mean squared error criterion is used to train the network, then the :math:`k`
hidden units learn to project the input in the span of the first :math:`k`
principal components of the data. If the hidden
layer is non-linear, the auto-encoder behaves differently from PCA,
with the ability to capture multi-modal aspects of the input
distribution. The departure from PCA becomes even more important when
we consider *stacking multiple encoders* (and their corresponding decoders)
when building a deep auto-encoder [Hinton06]_.

We want to implement an auto-encoder using Theano, in the form of a class,
that could be afterwards used in constructing a stacked autoencoder. The
first step is to create shared variables for the parameters of the 
autoencoder ( :math:`\mathbf{W}`, :math:`\mathbf{b}` and 
:math:`\mathbf{b'}`, since we are using tied weights in this tutorial ): 



.. code-block:: python
 
    class AutoEncoder(object):

        def __init__(self, n_visible= 784, n_hidden= 500, input = None):

          # initial values for weights and biases
          # note : W' was written as `W_prime` and b' as `b_prime`

          # W is initialized with `initial_W` which is uniformly sampled
          # from -6./sqrt(n_visible+n_hidden) and 6./sqrt(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.random.uniform( \
              low = -numpy.sqrt(6./(n_visible+n_hidden)), \
              high = numpy.sqrt(6./(n_visible+n_hidden)), \
              size = (n_visible, n_hidden)), dtype = theano.config.floatX)
          initial_b       = numpy.zeros(n_hidden)
          initial_b_prime= numpy.zeros(n_visible)
     
    
          # theano shared variables for weights and biases
          self.W       = theano.shared(value = initial_W,       name = "W")
          self.b       = theano.shared(value = initial_b,       name = "b")
          # tied weights, therefore W_prime is W transpose
          self.W_prime = W.T 
          self.b_prime = theano.shared(value = initial_b_prime, name = "b'")

Note that we pass the symbolic ``input``  to the autoencoder as a
parameter. This is such that later we can concatenate layers of 
autoencoders to form a deep network: the symbolic output (the :math:`\mathbf{y}` above, self.y 
in the code below) of
the k-th layer will be the symbolic input of the (k+1)-th.

Now we can express the computation of the latent representation and of the reconstructed
signal: 

.. code-block:: python

    self.y    = T.nnet.sigmoid(T.dot(self.x, self.W      ) + self.b)
    self.z    = T.nnet.sigmoid(T.dot(self.y, self.W_prime) + self.b_prime)
    # note : we sum over the size of a data point; if we are using minibatches,
    #        L will  be a vector, with one entry per example in minibatch 
    self.L    = - T.sum( self.x*T.log(self.z) + (1-self.x)*T.log(1-self.z), axis=1 ) 
    # 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 average cost over the minibatch 
    self.cost = T.mean(self.L)

Training the autoencoder can be done by iteratively updating the parameters ``W``, 
``b`` and ``b_prime``, e.g., by stochastic gradient descent, such that the 
reconstruction cost is approximately minimized.

.. code-block:: python

  train = theano.function( [x], cost, updates = { \
            self.W       : self.W       - T.grad(self.cost, self.W      )*learning_rate, 
            self.b       : self.b       - T.grad(self.cost, self.b      )*learning_rate,
            self.b_prime : self.b_prime - T.grad(self.cost, self.b_prime)*learning_rate})

Note that the training script for the stacked denoising autoencoder will be
slightly different. The above algorithm was used as a building block
(for unsupervised pre-training of each layer) in [Bengio07]
to build deep networks. The main idea of unsupervised pre-training is
the following: each layer of a deep network is first trained using a
local unsupervised criterion (such as minimizing reconstruction error
at each level);
afterwards, all the encoders parameters are used to initialize a deep 
supervised neural network, and all the parameters are jointly
fine-tuned to minimize the supervised cost function.

One serious potential issue with auto-encoders is that if there is no other
constraint besides minimizing the reconstruction error, 
then an auto-encoder with :math:`n` inputs and an
encoding of dimension at least :math:`n` could potentially just learn
the identity function, for which many encodings would be useless (e.g.,
just copying the input), i.e., the autoencoder would not differentiate
test examples (from the training distribution) from other input configurations.
Surprisingly, experiments reported in [Bengio07]_ nonetheless 
suggest that in practice, when trained with
stochastic gradient descent, non-linear auto-encoders with more hidden units
than inputs (called overcomplete) yield useful representations
(in the sense of classification error measured on a network taking this
representation in input). A simple explanation is based on the 
observation that stochastic gradient
descent with early stopping is similar to an L2 regularization of the
parameters. To achieve perfect reconstruction of continuous
inputs, a one-hidden layer auto-encoder with non-linear hidden units
(exactly like in the above code)
needs very small weights in the first (encoding) layer (to bring the non-linearity of
the hidden units in their linear regime) and very large weights in the
second (decoding) layer.
With binary inputs, very large weights are
also needed to completely minimize the reconstruction error. Since the
implicit or explicit regularization makes it difficult to reach
large-weight solutions, the optimization algorithm finds encodings which
only work well for examples similar to those in the training set, which is
what we want. It means that the representation is exploiting statistical
regularities present in the training set, rather than learning to
replicate the identity function.

There are different ways that an auto-encoder with more hidden units
than inputs could be prevented from learning the identity, and still
capture something useful about the input in its hidden representation.
One is the addition of sparsity (forcing many of the hidden units to
be zero or near-zero), and it has been exploited very successfully
by many [Ranzato07]_ [Lee08]_. Another is to add randomness in the transformation from
input to reconstruction. This is exploited in Restricted Boltzmann
Machines (discussed later in this tutorial), as well as in
Denoising Auto-Encoders, discussed below. 

.. _DA:

Denoising Autoencoders
++++++++++++++++++++++

The idea behind denoising autoencoders is simple. In order to force
the hidden layer to discover more robust features and prevent it
from simply learning the identity, we train the
autoencoder to *reconstruct the input from a corrupted version of it*.

The denoising auto-encoder is a stochastic version of the auto-encoder.
Intuitively, a denoising auto-encoder does two things: try to encode the
input (preserve the information about the input), and try to undo the
effect of a corruption process stochastically applied to the input of the
auto-encoder. The latter can only be done by capturing the statistical
dependencies between the inputs. The denoising
auto-encoder can be understood from different perspectives 
( the manifold learning perspective, 
stochastic operator perspective, 
bottom-up -- information theoretic perspective, 
top-down -- generative model perspective ), all of which are explained in 
[Vincent08]. 
See also section 7.2 of [Bengio09]_ for an overview of auto-encoders.

In [Vincent08], the stochastic corruption process
consists in randomly setting some of the inputs (as many as half of them)
to zero. Hence the denoising auto-encoder is trying to *predict the corrupted (i.e. missing)
values from the uncorrupted (i.e., non-missing) values*, for randomly selected subsets of
missing patterns. Note how being able to predict any subset of variables
from the rest is a sufficient condition for completely capturing the
joint distribution between a set of variables (this is how Gibbs
sampling works).

To convert the autoencoder class into a denoising autoencoder class, all we 
need to do is to add a stochastic corruption step operating on the input. The input can be
corrupted in many ways, but in this tutorial we will stick to the original 
corruption mechanism of randomly masking entries of the input by making
them zero. The code below 
does just that : 

.. code-block:: python

	from theano.tensor.shared_randomstreals import RandomStreams

	theano_rng = RandomStreams()
	corrupted_x = x * theano.rng.binomial(x.shape, 1, 0.9)

The final denoising autoencoder class becomes : 

.. code-block:: python

   class dA(object):

       def __init__(self, n_visible= 784, n_hidden= 500, input= None):

          self.n_visible = n_visible
          self.n_hidden  = n_hidden
    
          # create a Theano random generator that gives symbolic random values
          theano_rng = RandomStreams()
          # create a numpy random generator
          numpy_rng = numpy.random.RandomState()
    
     
          # initial values for weights and biases
          # note : W' was written as `W_prime` and b' as `b_prime`

          initial_W = numpy.asarray( numpy.random.uniform( \
              low = -numpy.sqrt(6./(n_visible+n_hidden)), \
              high = numpy.sqrt(6./(n_visible+n_hidden)), \
              size = (n_visible, n_hidden)), dtype = theano.config.floatX)
          initial_b       = numpy.zeros(n_hidden)
          initial_b_prime= numpy.zeros(n_visible)
     
    
          # theano shared variables for weights and biases
          self.W       = theano.shared(value = initial_W,       name = "W")
          self.b       = theano.shared(value = initial_b,       name = "b")
          # tied weights, therefore W_prime is W transpose
          self.W_prime = self.W.T 
          self.b_prime = theano.shared(value = initial_b_prime, name = "b'")

          # if no input is given, generate a variable representing the input
          if input == 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.tilde_x  = theano_rng.binomial( self.x.shape,  1,  0.9) * self.x
          self.y   = T.nnet.sigmoid(T.dot(self.tilde_x, self.W      ) + self.b)
          self.z        = T.nnet.sigmoid(T.dot(self.y, self.W_prime) + self.b_prime)
          self.L = - T.sum( self.x*T.log(self.z) + (1-self.x)*T.log(1-self.z), axis=1 ) 
          # 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
          self.cost = T.mean(self.L)
          # note : y is computed from the corrupted `tilde_x`. Later on, 
          #        we will need the hidden layer obtained from the uncorrupted 
          #        input when for example we will pass this as input to the layer 
          #        above
          self.hidden_values = T.nnet.sigmoid( T.dot(self.x, self.W) + self.b)



.. _stacked_autoencoders:

Stacked Autoencoders
++++++++++++++++++++

The denoising autoencoders can be stacked to form a deep network by
feeding the latent representation (output code)
of the denoising auto-encoder found on the layer 
below as input to the current layer. The **unsupervised pre-training** of such an 
architecture is done one layer at a time. Each layer is trained as 
a denoising auto-encoder by minimizing the reconstruction of its input
(which is the output code of the previous layer).
Once the first :math:`k` layers 
are trained, we can train the :math:`k+1`-th layer because we can now 
compute the code or latent representation from the layer below. 
Once all layers are pre-trained, the network goes through a second stage
of training called **fine-tuning**. Here we consider **supervised fine-tuning**
where we want to minimize prediction error on a supervised task.
For this we first add a logistic regression 
layer on top of the network (more precisely on the output code of the
output layer). We then
train the entire network as we would train a multilayer 
perceptron. At this point, we only consider the encoding parts of
each auto-encoder.
This stage is supervised, since now we use the target class during
training (see the :ref:`mlp` for details on the multilayer perceptron).

This can be easily implemented in Theano, using the class defined
before for a denoising autoencoder :

.. code-block:: python

  class StackedAutoencoder():
   
    def __init__(self, input, n_ins, hidden_layers_sizes, n_outs):
        """ This class is made to support a variable number of layers. 

        :param input: symbolic variable describing the input of the SdA

        :param n_ins: dimension of the input to the sdA

        :param n_layers_sizes: intermediate layers size, must contain 
        at least one value

        :param n_outs: dimension of the output of the network
        """

Next step, we create an denoising autoencoder for each layer and link them
together:

.. code-block:: python

        self.layers =[]

        if len(hidden_layers_sizes) < 1 :
            raiseException (' You must have at least one hidden layer ')

        # add first layer:
        layer = dA(n_ins, hidden_layers_sizes[0], input = input)
        self.layers += [layer]
        # add all intermediate layers
        for i in xrange( 1, len(hidden_layers_sizes) ):
            # input size is that of the previous layer
            # input is the output of the last layer inserted in our list 
            # of layers `self.layers`
            layer = dA( hidden_layers_sizes[i-1],             \
                        hidden_layers_sizes[i],               \
                        input = self.layers[-1].hidden_values )
            self.layers += [layer]
        

        self.n_layers = len(self.layers)


Note that during the second stage of training (fine-tuning) we need to 
use the weights of the autoencoders to define a multilayer perceptron.
This is already given by the above lines of code, in the sense that 
the ``hidden_values`` of the last denoising autoencoder already computes what
should be the input of the logistic regression layer that sits at the 
top of the MLP. All we need now is to add the logistic layer. We will 
use the ``LogisticRegression`` class introduced in :ref:`logreg`. 

.. code-block:: python

        # add a logistic layer on top
        self.logLayer = LogisticRegression(\
                         input = self.layers[-1].hidden_values,\
                         n_in = hidden_layers_sizes[-1], n_out = n_outs)


The negative log likelihood of this MLP (formed from reusing the weights 
of the denoising autoencoders) is given by the negative log likelihood
function of the logistic layer : 

.. code-block:: python

    def negative_log_likelihood(self, y):
        """Return the mean of the negative log-likelihood of the prediction
        of this model under a given target distribution. In our case this 
        is given by the logistic layer.

        :param y: corresponds to a vector that gives for each example the
        :correct label
        """
        return self.logLayer.negative_log_likelihood(y)

    def errors(self, y):
        """Return a float representing the number of errors in the minibatch 
        over the total number of examples of the minibatch 
        """
        
        return self.logLayer.errors(y)



Putting it all together
+++++++++++++++++++++++

The few lines of code below constructs the stacked denoising
autoencoder : 

.. code-block:: python 


    # construct the logistic regression class
    classifier = SdA( input=x, n_ins=28*28, \
                      hidden_layers_sizes = [500, 500, 500], n_outs=10)
    

There are two stages in training this network, a layer wise pre-training and 
fine-tuning afterwards. 

For the pre-training stage, we will loop over all the layers of the
network. For each layer we will compile a theano function that
implements a SGD step towards optimizing the weights for reducing 
the reconstruction cost of that layer. This function will be applied 
to the training set for a fixed number of epochs given by
``pretraining_epochs``


.. code-block:: python


    ## Pre-train layer-wise 
    for i in xrange(classifier.n_layers):
        # compute gradients of layer parameters
        gW       = T.grad(classifier.layers[i].cost, classifier.layers[i].W)
        gb       = T.grad(classifier.layers[i].cost, classifier.layers[i].b)
        gb_prime = T.grad(classifier.layers[i].cost, \
                                               classifier.layers[i].b_prime)
        # updated value of parameters after each step
        new_W       = classifier.layers[i].W      - gW      * pretraining_lr
        new_b       = classifier.layers[i].b      - gb      * pretraining_lr
        new_b_prime = classifier.layers[i].b_prime- gb_prime* pretraining_lr
        cost = classifier.layers[i].cost
        layer_update = theano.function([index], [cost], \
          updates = { 
              classifier.layers[i].W       : new_W \
            , classifier.layers[i].b       : new_b \
            , classifier.layers[i].b_prime : new_b_prime },
          givens = {
              x :train_set_x[index*batch_size:(index+1)*batch_size]})
        # go through pretraining epochs 
        for epoch in xrange(pretraining_epochs):
            # go through the training set
            for batch_index in xrange(n_train_batches):
                c = layer_update(batch_index)
            print 'Pre-training layer %i, epoch %d'%(i,epoch)
 

The fine-tuning loop is very similar with the one in the :ref:`mlp`, we just
have a slightly more complex training function. The reason is that now we 
need to update all parameters of the network in one call of the training 
function ( this includes the weight and biases of the denoising
autoencoders plus those of the logistic regression layer). To create 
this function, we will loop over the layers and create an update list 
containing pairs of the form (parameter before the SGD step, parameter after 
the SGD step). The new value of a parameter can be easily computed by calling 
``T.grad`` to compute the corresponding gradient, multiply it with the 
learning rate and subtract the result from the old value of the
parameter: 

.. code-block:: python 

    # Fine-tune the entire model
    # the cost we minimize during training is the negative log likelihood of 
    # the model
    cost = classifier.negative_log_likelihood(y) 

    # compute the gradient of cost with respect to theta and add them to the 
    # updates list
    updates = []
    for i in xrange(classifier.n_layers):        
        g_W   = T.grad(cost, classifier.layers[i].W)
        g_b   = T.grad(cost, classifier.layers[i].b)
        new_W = classifier.layers[i].W - learning_rate * g_W
        new_b = classifier.layers[i].b - learning_rate * g_b
        updates += [ (classifier.layers[i].W, new_W) \
                   , (classifier.layers[i].b, new_b) ]
    # add the gradients of the logistic layer
    g_log_W   = T.grad(cost, classifier.logLayer.W)
    g_log_b   = T.grad(cost, classifier.logLayer.b)
    new_log_W = classifier.logLayer.W - learning_rate * g_log_W
    new_log_b = classifier.logLayer.b - learning_rate * g_log_b
    updates += [ (classifier.logLayer.W, new_log_W) \
               , (classifier.logLayer.b, new_log_b) ]

    # compiling a theano function `train_model` that returns the cost, but  
    # in the same time updates the parameter of the model based on the rules 
    # defined in `updates`
    train_model = theano.function([index], cost, updates=updates,
          givens = {
            x: train_set_x[index*batch_size:(index+1)*batch_size],
            y: train_set_y[index*batch_size:(index+1)*batch_size]})

Now we pass this ``training_model`` (together with a ``validate_model`` and 
a ``test_model`` generated as in the other tutorials) to the 
early stopping loop and we are done. 




Running the Code
++++++++++++++++

pretraining duration : 79.4971 minutes


282.195 min , 144epochs 1.36 validation , 1.6 test
TODO

References
++++++++++

.. [Hinton06] `Reducing the Dimensionality of Data with Neural Networks <http://www.cs.toronto.edu/~rsalakhu/papers/science.pdf>`_, <NAME> and <NAME>,  Science, 28 July 2006, Vol. 313. no. 5786, pp. 504 - 507.

.. [Bengio07] `Greedy Layer-Wise Training of Deep Networks <http://www.iro.umontreal.ca/~lisa/publications2/index.php/publications/show/190>`_, <NAME>, <NAME>, <NAME> and <NAME>, in Advances in Neural Information Processing Systems 19 (NIPS'06), pages  153-160, MIT Press 2007.

.. [Ranzato07] `Efficient Learning of Sparse Representations with an Energy-Based Model <http://yann.lecun.com/exdb/publis/pdf/ranzato-06.pdf>`_, <NAME>, <NAME>, <NAME> and <NAME>, in J. Platt et al. (Eds), Advances in Neural Information Processing Systems (NIPS 2006), MIT Press, 2007.

.. [Lee08] `Sparse deep belief net model for visual area V2 <http://www.stanford.edu/~hllee/nips07-sparseDBN.pdf>`_, <NAME>, <NAME>, and <NAME>., in Advances in Neural Information Processing Systems (NIPS) 20, 2008.
<commit_msg>added some results ; they should be improved though (1.6% testing error)<commit_after>.. _SdA:

Stacked Denoising Autoencoders (SdA)
====================================

.. note::
  This section assumes the reader has already read through :doc:`logreg`
  and :doc:`mlp`. Additionally it uses the following Theano functions
  and concepts : TODO


The Stacked Denoising Autoencoder (SdA) is an extension of the stacked 
autoencoder [Bengio07]_ and it was introduced in [Vincent08]_. We will start the 
tutorial with a short discussion on :ref:`autoencoders`
and then move on to how classical
autoencoders are extended to denoising autoencoders (:ref:`DA`).

.. _autoencoders:

Autoencoders
+++++++++++++

See section 4.6 of [Bengio09]_ for an overview of auto-encoders.
An autoencoder takes an input :math:`\mathbf{x} \in [0,1]^d` and first 
maps it (with an *encoder*) to a hidden representation :math:`\mathbf{y} \in [0,1]^{d'}` 
through a deterministic mapping, e.g.:

.. math::
  
  \mathbf{y} = s(\mathbf{W}\mathbf{x} + \mathbf{b})

Where :math:`s` is a non-linearity such as the sigmoid.
The latent representation :math:`\mathbf{y}`, or **code** is then mapped back (with a *decoder*) into a
**reconstruction** :math:`\mathbf{z}` of same shape as
:math:`\mathbf{x}` through a similar transformation, e.g.:

.. math::

  \mathbf{z} = s(\mathbf{W'}\mathbf{y} + \mathbf{b'})

where ' does not indicate transpose, and
:math:`\mathbf{z}` should be seen as a prediction of :math:`\mathbf{x}`, given the code :math:`\mathbf{y}`.
The weight matrix :math:`\mathbf{W'}` of the reverse mapping may be
optionally constrained by :math:`\mathbf{W'} = \mathbf{W}^T`, which is
an instance of *tied weights*. The parameters of this model (namely 
:math:`\mathbf{W}`, :math:`\mathbf{b}`, 
:math:`\mathbf{b'}` and, if one doesn't use tied weights, also 
:math:`\mathbf{W'}`) are optimized such that the average reconstruction 
error is minimized. The reconstruction error can be measured in many ways, depending
on the appropriate distributional assumptions on the input given the code, e.g., using the 
traditional *squared error* :math:`L(\mathbf{x}, \mathbf{z}) = || \mathbf{x} - \mathbf{z} ||^2`, 
or if the input is interpreted as either bit vectors or vectors of 
bit probabilities by the reconstruction *cross-entropy* defined as : 

.. math::

  L_{H} (\mathbf{x}, \mathbf{z}) = - \sum^d_{k=1}[\mathbf{x}_k \log
          \mathbf{z}_k + (1 - \mathbf{x}_k)\log(1 - \mathbf{z}_k)] 

The hope is that the code :math:`\mathbf{y}` is a distributed representation
that captures the coordinates along the main factors of variation in the data
(similarly to how the projection on principal components captures the main factors
of variation in the data).
Because :math:`\mathbf{y}` is viewed as a lossy compression of :math:`\mathbf{x}`, it cannot
be a good compression (with small loss) for all :math:`\mathbf{x}`, so learning
drives it to be one that is a good compression in particular for training
examples, and hopefully for others as well, but not for arbitrary inputs.
That is the sense in which an auto-encoder generalizes: it gives low reconstruction
error to test examples from the same distribution as the training examples,
but generally high reconstruction error to uniformly chosen configurations of the
input vector.

If there is one linear hidden layer (the code) and
the mean squared error criterion is used to train the network, then the :math:`k`
hidden units learn to project the input in the span of the first :math:`k`
principal components of the data. If the hidden
layer is non-linear, the auto-encoder behaves differently from PCA,
with the ability to capture multi-modal aspects of the input
distribution. The departure from PCA becomes even more important when
we consider *stacking multiple encoders* (and their corresponding decoders)
when building a deep auto-encoder [Hinton06]_.

We want to implement an auto-encoder using Theano, in the form of a class,
that could be afterwards used in constructing a stacked autoencoder. The
first step is to create shared variables for the parameters of the 
autoencoder ( :math:`\mathbf{W}`, :math:`\mathbf{b}` and 
:math:`\mathbf{b'}`, since we are using tied weights in this tutorial ): 



.. code-block:: python
 
    class AutoEncoder(object):

        def __init__(self, n_visible= 784, n_hidden= 500, input = None):

          # initial values for weights and biases
          # note : W' was written as `W_prime` and b' as `b_prime`

          # W is initialized with `initial_W` which is uniformly sampled
          # from -6./sqrt(n_visible+n_hidden) and 6./sqrt(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.random.uniform( \
              low = -numpy.sqrt(6./(n_visible+n_hidden)), \
              high = numpy.sqrt(6./(n_visible+n_hidden)), \
              size = (n_visible, n_hidden)), dtype = theano.config.floatX)
          initial_b       = numpy.zeros(n_hidden)
          initial_b_prime= numpy.zeros(n_visible)
     
    
          # theano shared variables for weights and biases
          self.W       = theano.shared(value = initial_W,       name = "W")
          self.b       = theano.shared(value = initial_b,       name = "b")
          # tied weights, therefore W_prime is W transpose
          self.W_prime = W.T 
          self.b_prime = theano.shared(value = initial_b_prime, name = "b'")

Note that we pass the symbolic ``input``  to the autoencoder as a
parameter. This is such that later we can concatenate layers of 
autoencoders to form a deep network: the symbolic output (the :math:`\mathbf{y}` above, self.y 
in the code below) of
the k-th layer will be the symbolic input of the (k+1)-th.

Now we can express the computation of the latent representation and of the reconstructed
signal: 

.. code-block:: python

    self.y    = T.nnet.sigmoid(T.dot(self.x, self.W      ) + self.b)
    self.z    = T.nnet.sigmoid(T.dot(self.y, self.W_prime) + self.b_prime)
    # note : we sum over the size of a data point; if we are using minibatches,
    #        L will  be a vector, with one entry per example in minibatch 
    self.L    = - T.sum( self.x*T.log(self.z) + (1-self.x)*T.log(1-self.z), axis=1 ) 
    # 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 average cost over the minibatch 
    self.cost = T.mean(self.L)

Training the autoencoder can be done by iteratively updating the parameters ``W``, 
``b`` and ``b_prime``, e.g., by stochastic gradient descent, such that the 
reconstruction cost is approximately minimized.

.. code-block:: python

  train = theano.function( [x], cost, updates = { \
            self.W       : self.W       - T.grad(self.cost, self.W      )*learning_rate, 
            self.b       : self.b       - T.grad(self.cost, self.b      )*learning_rate,
            self.b_prime : self.b_prime - T.grad(self.cost, self.b_prime)*learning_rate})

Note that the training script for the stacked denoising autoencoder will be
slightly different. The above algorithm was used as a building block
(for unsupervised pre-training of each layer) in [Bengio07]
to build deep networks. The main idea of unsupervised pre-training is
the following: each layer of a deep network is first trained using a
local unsupervised criterion (such as minimizing reconstruction error
at each level);
afterwards, all the encoders parameters are used to initialize a deep 
supervised neural network, and all the parameters are jointly
fine-tuned to minimize the supervised cost function.

One serious potential issue with auto-encoders is that if there is no other
constraint besides minimizing the reconstruction error, 
then an auto-encoder with :math:`n` inputs and an
encoding of dimension at least :math:`n` could potentially just learn
the identity function, for which many encodings would be useless (e.g.,
just copying the input), i.e., the autoencoder would not differentiate
test examples (from the training distribution) from other input configurations.
Surprisingly, experiments reported in [Bengio07]_ nonetheless 
suggest that in practice, when trained with
stochastic gradient descent, non-linear auto-encoders with more hidden units
than inputs (called overcomplete) yield useful representations
(in the sense of classification error measured on a network taking this
representation in input). A simple explanation is based on the 
observation that stochastic gradient
descent with early stopping is similar to an L2 regularization of the
parameters. To achieve perfect reconstruction of continuous
inputs, a one-hidden layer auto-encoder with non-linear hidden units
(exactly like in the above code)
needs very small weights in the first (encoding) layer (to bring the non-linearity of
the hidden units in their linear regime) and very large weights in the
second (decoding) layer.
With binary inputs, very large weights are
also needed to completely minimize the reconstruction error. Since the
implicit or explicit regularization makes it difficult to reach
large-weight solutions, the optimization algorithm finds encodings which
only work well for examples similar to those in the training set, which is
what we want. It means that the representation is exploiting statistical
regularities present in the training set, rather than learning to
replicate the identity function.

There are different ways that an auto-encoder with more hidden units
than inputs could be prevented from learning the identity, and still
capture something useful about the input in its hidden representation.
One is the addition of sparsity (forcing many of the hidden units to
be zero or near-zero), and it has been exploited very successfully
by many [Ranzato07]_ [Lee08]_. Another is to add randomness in the transformation from
input to reconstruction. This is exploited in Restricted Boltzmann
Machines (discussed later in this tutorial), as well as in
Denoising Auto-Encoders, discussed below. 

.. _DA:

Denoising Autoencoders
++++++++++++++++++++++

The idea behind denoising autoencoders is simple. In order to force
the hidden layer to discover more robust features and prevent it
from simply learning the identity, we train the
autoencoder to *reconstruct the input from a corrupted version of it*.

The denoising auto-encoder is a stochastic version of the auto-encoder.
Intuitively, a denoising auto-encoder does two things: try to encode the
input (preserve the information about the input), and try to undo the
effect of a corruption process stochastically applied to the input of the
auto-encoder. The latter can only be done by capturing the statistical
dependencies between the inputs. The denoising
auto-encoder can be understood from different perspectives 
( the manifold learning perspective, 
stochastic operator perspective, 
bottom-up -- information theoretic perspective, 
top-down -- generative model perspective ), all of which are explained in 
[Vincent08]. 
See also section 7.2 of [Bengio09]_ for an overview of auto-encoders.

In [Vincent08], the stochastic corruption process
consists in randomly setting some of the inputs (as many as half of them)
to zero. Hence the denoising auto-encoder is trying to *predict the corrupted (i.e. missing)
values from the uncorrupted (i.e., non-missing) values*, for randomly selected subsets of
missing patterns. Note how being able to predict any subset of variables
from the rest is a sufficient condition for completely capturing the
joint distribution between a set of variables (this is how Gibbs
sampling works).

To convert the autoencoder class into a denoising autoencoder class, all we 
need to do is to add a stochastic corruption step operating on the input. The input can be
corrupted in many ways, but in this tutorial we will stick to the original 
corruption mechanism of randomly masking entries of the input by making
them zero. The code below 
does just that : 

.. code-block:: python

	from theano.tensor.shared_randomstreals import RandomStreams

	theano_rng = RandomStreams()
	corrupted_x = x * theano.rng.binomial(x.shape, 1, 0.9)

The final denoising autoencoder class becomes : 

.. code-block:: python

   class dA(object):

       def __init__(self, n_visible= 784, n_hidden= 500, input= None, corruption_level = 0.1):

          self.n_visible = n_visible
          self.n_hidden  = n_hidden
    
          # create a Theano random generator that gives symbolic random values
          theano_rng = RandomStreams()
          # create a numpy random generator
          numpy_rng = numpy.random.RandomState()
    
     
          # initial values for weights and biases
          # note : W' was written as `W_prime` and b' as `b_prime`

          initial_W = numpy.asarray( numpy.random.uniform( \
              low = -numpy.sqrt(6./(n_visible+n_hidden)), \
              high = numpy.sqrt(6./(n_visible+n_hidden)), \
              size = (n_visible, n_hidden)), dtype = theano.config.floatX)
          initial_b       = numpy.zeros(n_hidden)
          initial_b_prime= numpy.zeros(n_visible)
     
    
          # theano shared variables for weights and biases
          self.W       = theano.shared(value = initial_W,       name = "W")
          self.b       = theano.shared(value = initial_b,       name = "b")
          # tied weights, therefore W_prime is W transpose
          self.W_prime = self.W.T 
          self.b_prime = theano.shared(value = initial_b_prime, name = "b'")

          # if no input is given, generate a variable representing the input
          if input == 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
          # keep 90% of the inputs the same and zero-out randomly selected subset of 10% of the inputs
          # note : first argument of theano.rng.binomial is the shape(size) of 
          #        random numbers that it should produce
          #        second argument is the number of trials 
          #        third argument is the probability of success of any trial
          #
          #        this will produce an array of 0s and 1s where 1 has a 
          #        probability of 1 - ``corruption_level`` and 0 with
          #        ``corruption_level``
          self.tilde_x  = theano_rng.binomial( self.x.shape,  1,  1 - corruption_level) * self.x
          self.y   = T.nnet.sigmoid(T.dot(self.tilde_x, self.W      ) + self.b)
          self.z        = T.nnet.sigmoid(T.dot(self.y, self.W_prime) + self.b_prime)
          self.L = - T.sum( self.x*T.log(self.z) + (1-self.x)*T.log(1-self.z), axis=1 ) 
          # 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
          self.cost = T.mean(self.L)
          # note : y is computed from the corrupted `tilde_x`. Later on, 
          #        we will need the hidden layer obtained from the uncorrupted 
          #        input when for example we will pass this as input to the layer 
          #        above
          self.hidden_values = T.nnet.sigmoid( T.dot(self.x, self.W) + self.b)



.. _stacked_autoencoders:

Stacked Autoencoders
++++++++++++++++++++

The denoising autoencoders can be stacked to form a deep network by
feeding the latent representation (output code)
of the denoising auto-encoder found on the layer 
below as input to the current layer. The **unsupervised pre-training** of such an 
architecture is done one layer at a time. Each layer is trained as 
a denoising auto-encoder by minimizing the reconstruction of its input
(which is the output code of the previous layer).
Once the first :math:`k` layers 
are trained, we can train the :math:`k+1`-th layer because we can now 
compute the code or latent representation from the layer below. 
Once all layers are pre-trained, the network goes through a second stage
of training called **fine-tuning**. Here we consider **supervised fine-tuning**
where we want to minimize prediction error on a supervised task.
For this we first add a logistic regression 
layer on top of the network (more precisely on the output code of the
output layer). We then
train the entire network as we would train a multilayer 
perceptron. At this point, we only consider the encoding parts of
each auto-encoder.
This stage is supervised, since now we use the target class during
training (see the :ref:`mlp` for details on the multilayer perceptron).

This can be easily implemented in Theano, using the class defined
before for a denoising autoencoder :

.. code-block:: python

  class StackedAutoencoder():
   
    def __init__(self, input, n_ins, hidden_layers_sizes, n_outs, corruption_levels):
        """ This class is made to support a variable number of layers. 

        :param input: symbolic variable describing the input of the SdA

        :param n_ins: dimension of the input to the sdA

        :param n_layers_sizes: intermediate layers size, must contain 
        at least one value

        :param n_outs: dimension of the output of the network
        
        :param corruption_levels: amount of corruption to use for each 
        layer
 
        """

Next step, we create an denoising autoencoder for each layer and link them
together:

.. code-block:: python

        self.layers =[]

        if len(hidden_layers_sizes) < 1 :
            raiseException (' You must have at least one hidden layer ')

        # add first layer:
        layer = dA(n_ins, hidden_layers_sizes[0], input = input, \
                   corruption_level = corruption_levels[0])
        self.layers += [layer]
        # add all intermediate layers
        for i in xrange( 1, len(hidden_layers_sizes) ):
            # input size is that of the previous layer
            # input is the output of the last layer inserted in our list 
            # of layers `self.layers`
            layer = dA( hidden_layers_sizes[i-1],             \
                        hidden_layers_sizes[i],               \
                        input = self.layers[-1].hidden_values,\
                        corruption_level = corruption_levels[i])
            self.layers += [layer]
        

        self.n_layers = len(self.layers)


Note that during the second stage of training (fine-tuning) we need to 
use the weights of the autoencoders to define a multilayer perceptron.
This is already given by the above lines of code, in the sense that 
the ``hidden_values`` of the last denoising autoencoder already computes what
should be the input of the logistic regression layer that sits at the 
top of the MLP. All we need now is to add the logistic layer. We will 
use the ``LogisticRegression`` class introduced in :ref:`logreg`. 

.. code-block:: python

        # add a logistic layer on top
        self.logLayer = LogisticRegression(\
                         input = self.layers[-1].hidden_values,\
                         n_in = hidden_layers_sizes[-1], n_out = n_outs)


The negative log likelihood of this MLP (formed from reusing the weights 
of the denoising autoencoders) is given by the negative log likelihood
function of the logistic layer : 

.. code-block:: python

    def negative_log_likelihood(self, y):
        """Return the mean of the negative log-likelihood of the prediction
        of this model under a given target distribution. In our case this 
        is given by the logistic layer.

        :param y: corresponds to a vector that gives for each example the
        :correct label
        """
        return self.logLayer.negative_log_likelihood(y)

    def errors(self, y):
        """Return a float representing the number of errors in the minibatch 
        over the total number of examples of the minibatch 
        """
        
        return self.logLayer.errors(y)



Putting it all together
+++++++++++++++++++++++

The few lines of code below constructs the stacked denoising
autoencoder : 

.. code-block:: python 


    # construct the logistic regression class
    classifier = SdA( input=x, n_ins=28*28, \
                      hidden_layers_sizes = [1000, 1000, 1o00], n_outs=10)
    

There are two stages in training this network, a layer wise pre-training and 
fine-tuning afterwards. 

For the pre-training stage, we will loop over all the layers of the
network. For each layer we will compile a theano function that
implements a SGD step towards optimizing the weights for reducing 
the reconstruction cost of that layer. This function will be applied 
to the training set for a fixed number of epochs given by
``pretraining_epochs``


.. code-block:: python


    ## Pre-train layer-wise 
    for i in xrange(classifier.n_layers):
        # compute gradients of layer parameters
        gW       = T.grad(classifier.layers[i].cost, classifier.layers[i].W)
        gb       = T.grad(classifier.layers[i].cost, classifier.layers[i].b)
        gb_prime = T.grad(classifier.layers[i].cost, \
                                               classifier.layers[i].b_prime)
        # updated value of parameters after each step
        new_W       = classifier.layers[i].W      - gW      * pretraining_lr
        new_b       = classifier.layers[i].b      - gb      * pretraining_lr
        new_b_prime = classifier.layers[i].b_prime- gb_prime* pretraining_lr
        cost = classifier.layers[i].cost
        layer_update = theano.function([index], [cost], \
          updates = { 
              classifier.layers[i].W       : new_W \
            , classifier.layers[i].b       : new_b \
            , classifier.layers[i].b_prime : new_b_prime },
          givens = {
              x :train_set_x[index*batch_size:(index+1)*batch_size]})
        # go through pretraining epochs 
        for epoch in xrange(pretraining_epochs):
            # go through the training set
            for batch_index in xrange(n_train_batches):
                c = layer_update(batch_index)
            print 'Pre-training layer %i, epoch %d'%(i,epoch)
 

The fine-tuning loop is very similar with the one in the :ref:`mlp`, we just
have a slightly more complex training function. The reason is that now we 
need to update all parameters of the network in one call of the training 
function ( this includes the weight and biases of the denoising
autoencoders plus those of the logistic regression layer). To create 
this function, we will loop over the layers and create an update list 
containing pairs of the form (parameter before the SGD step, parameter after 
the SGD step). The new value of a parameter can be easily computed by calling 
``T.grad`` to compute the corresponding gradient, multiply it with the 
learning rate and subtract the result from the old value of the
parameter: 

.. code-block:: python 

    # Fine-tune the entire model
    # the cost we minimize during training is the negative log likelihood of 
    # the model
    cost = classifier.negative_log_likelihood(y) 

    # compute the gradient of cost with respect to theta and add them to the 
    # updates list
    updates = []
    for i in xrange(classifier.n_layers):        
        g_W   = T.grad(cost, classifier.layers[i].W)
        g_b   = T.grad(cost, classifier.layers[i].b)
        new_W = classifier.layers[i].W - learning_rate * g_W
        new_b = classifier.layers[i].b - learning_rate * g_b
        updates += [ (classifier.layers[i].W, new_W) \
                   , (classifier.layers[i].b, new_b) ]
    # add the gradients of the logistic layer
    g_log_W   = T.grad(cost, classifier.logLayer.W)
    g_log_b   = T.grad(cost, classifier.logLayer.b)
    new_log_W = classifier.logLayer.W - learning_rate * g_log_W
    new_log_b = classifier.logLayer.b - learning_rate * g_log_b
    updates += [ (classifier.logLayer.W, new_log_W) \
               , (classifier.logLayer.b, new_log_b) ]

    # compiling a theano function `train_model` that returns the cost, but  
    # in the same time updates the parameter of the model based on the rules 
    # defined in `updates`
    train_model = theano.function([index], cost, updates=updates,
          givens = {
            x: train_set_x[index*batch_size:(index+1)*batch_size],
            y: train_set_y[index*batch_size:(index+1)*batch_size]})

Now we pass this ``training_model`` (together with a ``validate_model`` and 
a ``test_model`` generated as in the other tutorials) to the 
early stopping loop and we are done. 




Running the Code
++++++++++++++++

The user can run the code by calling:

.. code-block:: bash
  
  python code/SdA.py

By default the code runs 15 pre-training epochs for each layer, with 
a corruption level of 0.1 and a learning rate of 0.1. Pre-training takes
79.4971 minutes. Fine-tuning is completed after 144 epochs in 282.195 
minutes and results  in a validation score of 1.36 %, with a test 
performace of 1.6 %.


References
++++++++++

.. [Hinton06] `Reducing the Dimensionality of Data with Neural Networks <http://www.cs.toronto.edu/~rsalakhu/papers/science.pdf>`_, <NAME> and <NAME>,  Science, 28 July 2006, Vol. 313. no. 5786, pp. 504 - 507.

.. [Bengio07] `Greedy Layer-Wise Training of Deep Networks <http://www.iro.umontreal.ca/~lisa/publications2/index.php/publications/show/190>`_, <NAME>, <NAME>, <NAME> and <NAME>, in Advances in Neural Information Processing Systems 19 (NIPS'06), pages  153-160, MIT Press 2007.

.. [Ranzato07] `Efficient Learning of Sparse Representations with an Energy-Based Model <http://yann.lecun.com/exdb/publis/pdf/ranzato-06.pdf>`_, <NAME>, <NAME>, <NAME> and <NAME>, in J. Platt et al. (Eds), Advances in Neural Information Processing Systems (NIPS 2006), MIT Press, 2007.

.. [Lee08] `Sparse deep belief net model for visual area V2 <http://www.stanford.edu/~hllee/nips07-sparseDBN.pdf>`_, <NAME>, <NAME>, and <NAME>., in Advances in Neural Information Processing Systems (NIPS) 20, 2008.

{'fraction_non_alphanumeric': 0.08285564386063773, 'fraction_numerical': 0.007247925694188057, 'mean_word_length': 3.079587792310741, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1715383', 'n_tokens_mistral': 14963, 'n_tokens_neox': 14090, 'n_words': 6154}