rbm.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:f62d84a1bb24e9288b9c0f7f4a40758e9a204668ca350cce30d7700af002107d"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Training layer one\n",
      "\n",
      "This notebook shows how to train a single-layer RBM on MNIST using Numpy.\n",
      "\n",
      "Based off the demo code at\n",
      "    http://www.cs.toronto.edu/~hinton/MatlabForSciencePaper.html"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import plotting\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Load the data\n",
      "\n",
      "First, we download and open the MNIST data set."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import gzip\n",
      "import cPickle as pickle\n",
      "import urllib\n",
      "\n",
      "filename = 'mnist.pkl.gz'\n",
      "\n",
      "if not os.path.exists(filename):\n",
      "    url = 'http://deeplearning.net/data/mnist/mnist.pkl.gz'\n",
      "    urllib.urlretrieve(url, filename=filename)\n",
      "\n",
      "with gzip.open(filename, 'rb') as f:\n",
      "    train, valid, test = pickle.load(f)\n",
      "    \n",
      "# use the training set images\n",
      "images = train[0]\n",
      "test_images = test[0]\n",
      "image_shape = (28, 28)\n",
      "assert np.prod(image_shape) == images.shape[1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### RBM parameters\n",
      "\n",
      "Next, we define the core parameters for our RBM. These include the number of visual and hidden nodes, and the initial weights `W` and biases `c` and `b`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rng = np.random\n",
      "\n",
      "n_vis = images.shape[1]\n",
      "n_hid = 500\n",
      "\n",
      "Wmag = 4 * np.sqrt(6. / (n_vis + n_hid))\n",
      "W = rng.uniform(low=-Wmag, high=Wmag, size=(n_vis, n_hid))\n",
      "c = np.zeros(n_hid)\n",
      "b = np.zeros(n_vis)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This RBM will have sparse weights, i.e. each neuron has a limited receptive field. To accomplish this, we create a mask that denotes the non-zero weights."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rf_shape = (9, 9)\n",
      "M, N = image_shape\n",
      "m, n = rf_shape\n",
      "\n",
      "# find random positions for top-left corner of each RF\n",
      "i = rng.randint(low=0, high=M-m+1, size=n_hid)\n",
      "j = rng.randint(low=0, high=N-n+1, size=n_hid)\n",
      "\n",
      "mask = np.zeros((M, N, n_hid), dtype='bool')\n",
      "for k in xrange(n_hid):\n",
      "    mask[i[k]:i[k]+m, j[k]:j[k]+n, k] = True\n",
      "\n",
      "mask = mask.reshape(n_vis, n_hid)\n",
      "W = W * mask  # make initial W sparse"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### RBM propagation functions\n",
      "\n",
      "We need to define the propagation functions to allow us to get the probabilities of the hidden layer given the visual (`probHgivenV`), the visual layer given the hidden (`probVgivenH`), and to sample the hidden layer (`sampHgivenV`). If we wanted to do more than one iteration of Gibbs sampling, we would need to also define `sampVgivenH` to sample the visual layer given the hidden."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def sigmoid(x):\n",
      "    return 1. / (1 + np.exp(-x))\n",
      "\n",
      "def probHgivenV(vis):\n",
      "    x = np.dot(vis, W) + c\n",
      "    return sigmoid(x)\n",
      "\n",
      "def probVgivenH(hid):\n",
      "    x = np.dot(hid, W.T) + b\n",
      "    return sigmoid(x)\n",
      "\n",
      "def sampHgivenV(vis):\n",
      "    hidprob = probHgivenV(vis)\n",
      "    hidsamp = hidprob > rng.uniform(size=hidprob.shape)\n",
      "    return hidprob, hidsamp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Train the RBM using CD\n",
      "\n",
      "We can now train the RBM using contrastive divergence (CD). We define the number of epochs, where each epoch goes through the training set once. We also divide the training set into mini-batches, where the gradient is updated after each mini-batch.\n",
      "\n",
      "We define the learning rate, a cost on large-magnitude weights, and the momentum. The momentum allows smoothing over the gradient estimates of many mini-batches, and may also help avoid local minima. To implement momentum, we need to track the last increments of the weights and biases (`Winc`, `cinc`, `binc`)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n_epochs = 5\n",
      "\n",
      "batch_size = 100\n",
      "batches = images.reshape(-1, batch_size, images.shape[1])\n",
      "\n",
      "rate = 0.1\n",
      "weightcost = 2e-4\n",
      "momentum = 0.5\n",
      "\n",
      "Winc = np.zeros_like(W)\n",
      "cinc = np.zeros_like(c)\n",
      "binc = np.zeros_like(b)\n",
      "\n",
      "for epoch in range(n_epochs):\n",
      "\n",
      "    # train on each mini-batch\n",
      "    costs = []\n",
      "    for batch in batches:\n",
      "        data = batch\n",
      "\n",
      "        # compute positive phase\n",
      "        poshidprob, poshidsamp = sampHgivenV(data)\n",
      "\n",
      "        posprods = np.dot(data.T, poshidprob) / batch_size\n",
      "        posvisact = np.mean(data, axis=0)\n",
      "        poshidact = np.mean(poshidprob, axis=0)\n",
      "\n",
      "        # compute negative phase\n",
      "        negdata = probVgivenH(poshidsamp)\n",
      "        neghidprob = probHgivenV(negdata)\n",
      "        negprods = np.dot(negdata.T, neghidprob) / batch_size\n",
      "        negvisact = np.mean(negdata, axis=0)\n",
      "        neghidact = np.mean(neghidprob, axis=0)\n",
      "\n",
      "        # compute error\n",
      "        rmse = np.sqrt(np.mean((data - negdata)**2, axis=1))\n",
      "        err = np.mean(rmse)\n",
      "\n",
      "        # compute updates\n",
      "        Winc = momentum * Winc + rate * (\n",
      "            (posprods - negprods) - weightcost * W)\n",
      "        cinc = momentum * cinc + rate * (poshidact - neghidact)\n",
      "        binc = momentum * binc + rate * (posvisact - negvisact)\n",
      "\n",
      "        W += Winc * mask\n",
      "        c += cinc\n",
      "        b += binc\n",
      "\n",
      "        costs.append(err)\n",
      "\n",
      "    print \"Epoch %d: %0.3f\" % (epoch, np.mean(costs))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To evaluate our RBM, we can plot its reconstructions of the testing data, as well as the filters (\"features\") that it has learned."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# --- plot reconstructions on test set\n",
      "plt.figure(2)\n",
      "plt.clf()\n",
      "codes = probHgivenV(test_images)\n",
      "recons = probVgivenH(codes)\n",
      "plotting.compare(\n",
      "    [test_images.reshape(-1, *image_shape),\n",
      "     recons.reshape(-1, *image_shape)],\n",
      "    rows=3, cols=10)\n",
      "\n",
      "# --- plot filters\n",
      "plt.figure(3)\n",
      "plt.clf()\n",
      "filters = W.T[mask.T].reshape(n_hid, *rf_shape)\n",
      "plotting.filters(filters, rows=5, cols=10)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}