DAE_theano_with_comments.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This tutorial is meant to be done after MLP_theano_with_comments.ipynb\n",
    "\n",
    "# We are working with MNIST again, this time no labels are necessary - \n",
    "# the denoising autoencoder (DAE) is an unsupervised model that tries to reconstruct the original input.\n",
    "\n",
    "# All imports up here this time\n",
    "import pickle\n",
    "import numpy\n",
    "import numpy.random as rng\n",
    "import theano\n",
    "import theano.tensor as T\n",
    "import theano.sandbox.rng_mrg as RNG_MRG\n",
    "from utils import tile_raster_images\n",
    "from PIL import Image as pil_img\n",
    "from IPython.display import Image\n",
    "\n",
    "# Load our data \n",
    "# Download and unzip pickled version from here: http://www.iro.umontreal.ca/~lisa/deep/data/mnist/mnist.pkl.gz\n",
    "(train_x, _), (valid_x, _), (test_x, _) = pickle.load(open('data/mnist.pkl', 'r'))\n",
    "print \"Shapes:\"\n",
    "print train_x.shape\n",
    "print valid_x.shape\n",
    "print test_x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The DAE data flow looks like this: input -> input_(add noise) -> hiddens -> input\n",
    "# This can be repeated by sampling from the reconstructed input. Repeating like that is\n",
    "# a pseudo Gibbs sampling process. We can define how many times we want to repeat (known as walkbacks).\n",
    "\n",
    "# We can specify any hyperparameters to play with up here:\n",
    "input_size = 784  # 28x28 images\n",
    "hidden_size = 1000\n",
    "w_mean = 0.0\n",
    "w_std = 0.05\n",
    "noise = 0.3\n",
    "walkbacks = 3\n",
    "learning_rate = 0.1\n",
    "batch_size = 100\n",
    "epochs = 100\n",
    "check_frequency = 10\n",
    "\n",
    "# To make the organization better, lets define all the variables and parameters here.\n",
    "# Just like with the MLP, we need a symbolic matrix to input the images\n",
    "x = T.matrix('x')\n",
    "# Next, we need the weights matrix W_x. This will be used to go both from input -> hidden and\n",
    "# hidden -> input (by using its transpose). This is called tied weights.\n",
    "# Again, initialization has a lot of literature, but we are just goint to stick with gaussian at the moment.\n",
    "W_x = numpy.asarray(rng.normal(loc=w_mean, scale=w_std, size=(input_size, hidden_size)), dtype=theano.config.floatX)\n",
    "# (Don't forget to make parameters into shared variables so they can be updated!)\n",
    "W_x = theano.shared(W_x, \"W_x\")\n",
    "# Because we are outputting back into the input space, we also need a bias vector for both the input\n",
    "# and hidden layers.\n",
    "b_x = numpy.zeros((input_size,), dtype=theano.config.floatX)\n",
    "b_h = numpy.zeros((hidden_size,), dtype=theano.config.floatX)\n",
    "b_x = theano.shared(b_x, \"b_x\")\n",
    "b_h = theano.shared(b_h, \"b_h\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Now for the most important part of a denoising autoencoder - making the input noisy!\n",
    "# Noise acts as regularization so the autoencoder doesn't just memorize the training set.\n",
    "# This makes it more effective for test data by reducing overfitting.\n",
    "\n",
    "# We deal with adding noise during training but not testing in Theano with a switch variable!\n",
    "# Switches can be turned on or off to direct data flow in the computation graph - \n",
    "# so we turn on for train and off for test!\n",
    "# You guessed it - we need a shared variable to represent the switch condition so we can change it at runtime.\n",
    "noise_switch = theano.shared(1, \"noise_switch\")\n",
    "\n",
    "# One important thing to note - the type of noise has to correspond to the type of input.\n",
    "# i.e. we can't add real-value noise when the input is expected to be binary\n",
    "# So for these binary inputs, we will add salt-and-pepper masking noise!\n",
    "# This is a function so we can keep adding it during the computation chain when we alternate sampling\n",
    "# from input and reconstructing from hiddens.\n",
    "# Theano random number generator\n",
    "theano_rng = RNG_MRG.MRG_RandomStreams(1)\n",
    "def salt_and_pepper(variable):\n",
    "    mask = theano_rng.binomial(size=variable.shape, n=1, p=(1-noise), dtype=theano.config.floatX)\n",
    "    saltpepper = theano_rng.binomial(size=variable.shape, n=1, p=0.5, dtype=theano.config.floatX)\n",
    "    ones = T.eq(mask, 0) * saltpepper\n",
    "    # Randomly set some bits to 0 or 1 with equal probability.\n",
    "    noisy = variable*mask + ones\n",
    "    return T.switch(noise_switch,\n",
    "                    # true condition\n",
    "                    noisy,\n",
    "                    # false condition\n",
    "                    variable)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Now we are ready to create the computation graph!\n",
    "# Remember it is noisy_x -> hiddens -> x -> hiddens -> x .....\n",
    "\n",
    "inputs=[x]\n",
    "for walkback in range(walkbacks):\n",
    "    # First, we want to corrupt the input x\n",
    "    noisy_x = salt_and_pepper(inputs[-1])\n",
    "    # Now calculate the hiddens\n",
    "    h = T.tanh(\n",
    "        T.dot(noisy_x, W_x) + b_h\n",
    "    )\n",
    "    # From the hiddens, reconstruct x.\n",
    "    # We have to use an appropriate activation function for the type of inputs!\n",
    "    # In our case with MNIST images, it is binary so sigmoid works.\n",
    "    reconstruction = T.nnet.sigmoid(\n",
    "        T.dot(h, W_x.T) + b_x\n",
    "    )\n",
    "    # That is all for an autoencoder!\n",
    "    inputs.append(reconstruction)\n",
    "    \n",
    "# Remove the original input from our reconstructions list\n",
    "reconstructions = inputs[1:]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The output of our computation graph is the last reconstructed input in the Gibbs chain.\n",
    "output = reconstructions[-1]\n",
    "\n",
    "# Our cost function is now the reconstruction error between all reconstructions and the original input.\n",
    "# Again, because our input space is binary, using mean binary cross-entropy is a good analog for \n",
    "# reconstruction error.\n",
    "# For real-valued inputs, we could use mean square error.\n",
    "cost = numpy.sum([T.mean(T.nnet.binary_crossentropy(recon, x)) for recon in reconstructions])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Just like with the MLP, compute gradient updates for the parameters to use with training.\n",
    "parameters = [W_x, b_h, b_x]\n",
    "# Automagic differentiation! (Still love it)\n",
    "gradients = T.grad(cost, parameters)\n",
    "\n",
    "# Update the parameters for stochastic gradient descent!\n",
    "train_updates = [(param, param - learning_rate*gradient) for param, gradient in zip(parameters, gradients)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compile our training and testing function like before!\n",
    "# Train function updates the parameters and returns the total train cost to monitor.\n",
    "f_train = theano.function(\n",
    "    inputs=[x], \n",
    "    outputs=cost, \n",
    "    updates=train_updates, \n",
    "    allow_input_downcast=True\n",
    ")\n",
    "\n",
    "# Our test function will return the final reconstruction, and it needs to include updates from scan.\n",
    "f_test = theano.function(\n",
    "    inputs=[x], \n",
    "    outputs=output,\n",
    "    allow_input_downcast=True\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# That's it! Now perform SGD like before.\n",
    "# Main training loop\n",
    "\n",
    "train_batches = len(train_x) / batch_size\n",
    "\n",
    "try:\n",
    "    for epoch in range(epochs):\n",
    "        print epoch+1, \":\",\n",
    "\n",
    "        # Don't forget to turn on our noise switch for training! Just set the shared variable to 1 (True)\n",
    "        noise_switch.set_value(1.)\n",
    "\n",
    "        train_costs = []\n",
    "        for i in range(train_batches):\n",
    "            batch_x = train_x[i*batch_size:(i+1)*batch_size]\n",
    "\n",
    "            costs = f_train(batch_x)\n",
    "\n",
    "            train_costs.append(costs)\n",
    "        print \"cost:\", numpy.mean(train_costs)\n",
    "\n",
    "        if (epoch+1) % check_frequency == 0:\n",
    "            print \"Saving images...\"\n",
    "            # For validation, let's run a few images through and see the reconstruction \n",
    "            # (with the noise from training still added)\n",
    "            valid_recons = f_test(valid_x[:25])\n",
    "            # Use the tile_raster_image helper function to rearrange the matrix into a 5x10 image of digits\n",
    "            # (Two 5x5 images next to each other - the first the inputs, the second the reconstructions.)\n",
    "            valid_stacked = numpy.vstack(\n",
    "                [numpy.vstack([\n",
    "                        valid_x[i*5:(i+1)*5],\n",
    "                        valid_recons[i*5:(i+1)*5]\n",
    "                    ])\n",
    "                 for i in range(5)]\n",
    "            )\n",
    "            valid_image = pil_img.fromarray(\n",
    "                # helper from utils.py\n",
    "                tile_raster_images(valid_stacked, (28, 28), (5, 10), (1, 1))\n",
    "            )\n",
    "            valid_image.save(\"dae_valid_%d.png\"%(epoch+1))\n",
    "\n",
    "            # Now do the same for test, but don't add any noise\n",
    "            # This means set the noise switches to 0. (False)\n",
    "            noise_switch.set_value(0.)\n",
    "\n",
    "            test_recons = f_test(test_x[:25])\n",
    "            test_stacked = numpy.vstack(\n",
    "                [numpy.vstack([\n",
    "                        test_x[i*5:(i+1)*5],\n",
    "                        test_recons[i*5:(i+1)*5]\n",
    "                    ])\n",
    "                 for i in range(5)]\n",
    "            )\n",
    "            test_image = pil_img.fromarray(\n",
    "                tile_raster_images(test_stacked, (28, 28), (5, 10), (1, 1))\n",
    "            )\n",
    "            test_image.save(\"dae_test_%d.png\"%(epoch+1))\n",
    "except KeyboardInterrupt:\n",
    "    pass\n",
    "        \n",
    "# Let's finally save an image of the filters the DAE learned - this is simply the transpose of the weights!\n",
    "weight_filters = pil_img.fromarray(\n",
    "    tile_raster_images(\n",
    "        W_x.get_value(borrow=True).T,\n",
    "        img_shape=(28, 28),\n",
    "        tile_shape=(25, 40),\n",
    "        tile_spacing=(1, 1)\n",
    "    )\n",
    ")\n",
    "print \"Saving filters...\"\n",
    "weight_filters.save(\"dae_filters.png\")\n",
    "print \"Done!\"\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}