WakeSleep.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we start, let's get some Python plumbing out of the way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from itertools import islice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I usually program in languages where this is built in :)\n",
    "# https://stackoverflow.com/questions/6822725/rolling-or-sliding-window-iterator-in-python\n",
    "\n",
    "def window(seq, n=2):\n",
    "    \"Returns a sliding window (of width n) over data from the iterable\"\n",
    "    \"   s -> (s0,s1,...s[n-1]), (s1,s2,...,sn), ...                   \"\n",
    "    it = iter(seq)\n",
    "    result = tuple(islice(it, n))\n",
    "    if len(result) == n:\n",
    "        yield result\n",
    "    for elem in it:\n",
    "        result = result[1:] + (elem,)\n",
    "        yield result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... back to regular scheduled programming\n",
    "\n",
    "The report starts here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The wake-sleep algorithm for unsupervised neural networks\n",
    "\n",
    "Paper by Geoffrey E Hinton, Peter Dayan, Brendan J Frey and Radford M Neal\n",
    "\n",
    "Read, interpreted and Pythonized by Vadim Liventsev for _Graphical Models of Statistical Inference_ course at [Skoltech](skoltech.ru/en/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is this all about?\n",
    "\n",
    "[Artificial neural networks](https://en.wikipedia.org/wiki/Artificial_neural_network) are typically used for supervised learning: given some training data with inputs and outputs create a model that maps inputs to outputs as closely as possible. This paper introduces _wake-sleep algorithm_ for unsupervised learning (training data contains only inputs) with neural networks and explains its theoretical underpinnings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How does it work?\n",
    "\n",
    "Wake-sleep algorithm mostly comes down to 2 clever tricks:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clever trick 1: Probability - Entropy approach\n",
    "\n",
    "A neural network can be seen as a coding scheme: the first layer takes an input vector and the last one outputs its representation (a code). If a reference set of \"correct\" representations is known, one can minimize mean squared deviation of actual representations from reference representations. But it's not necessary. Instead, one can take another approach inspired by coding theory: minimize the complexity of representations. Set the weights of the connections in a neural network so that inputs are encoded in the most compact way possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sounds simple, but developing this idea into an actual algorithm requires some preparation:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Meet Binary Stochastic Neuron\n",
    "\n",
    "Binary Stochastic Neuron is a neuron that takes a vector of binary values and outputs 0 or 1 with probability defined by the logistic function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Stochastic Neuron](ws-clipart/1-stochastic-neuron.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where _s<sub>v</sub>_ is the output of neuron _v_, _b<sub>v</sub>_ is the bias of neuron _v_, _s<sub>u</sub>_ is the output of neuron _u_ and _w<sub>uv</sub>_ is the weight of the connection from neuron _u_ to neuron _v_\n",
    "\n",
    "or, in Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class StochasticNeuron():\n",
    "    def __init__(self, input_weights = np.array([]), bias = 0):\n",
    "        self.input_weights = np.array(input_weights)\n",
    "        self.bias = bias\n",
    "        \n",
    "    def activation_probability(self, inputs):\n",
    "        inputs = np.array(inputs)\n",
    "        return 1 / (1 + np.exp(- self.bias - np.sum(self.input_weights * inputs)))\n",
    "    \n",
    "    def activation_function(self, inputs):\n",
    "        return np.random.binomial(1, self.activation_probability(inputs))\n",
    "    \n",
    "    def __call__(self, inputs):\n",
    "        return self.activation_function(inputs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = StochasticNeuron([0, -0.5, 1], 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3775406687981454"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n.activation_probability([False, True, False])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n([False, True, False])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6224593312018546"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n.activation_probability([False, True, True])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n([False, True, True])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can see that the neuron is stochastic by calculating it's output several times in a row "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[n([False, True, False]) for i in range(15)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic neural networks\n",
    "\n",
    "A layer of a neural network is just a tuple of neurons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Unbiased, unadjusted layer\n",
    "def create_layer(input_size, size):\n",
    "    return [StochasticNeuron(np.array([1.0 for i in range(input_size)]), 0) for i in range(size)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 1, 0, 0, 0, 0, 0, 0, 0, 0]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "layer = create_layer(3, 10)\n",
    "inp = [1, 2, -4]\n",
    "[neuron(inp) for neuron in layer]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again, let's infer the ouput of the layer several times and note the stochasticity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 0, 1, 0, 0, 0, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 1, 0, 0, 0, 0],\n",
       " [0, 0, 0, 0, 0, 1, 1, 0, 0, 0],\n",
       " [0, 0, 0, 0, 1, 1, 0, 0, 0, 0],\n",
       " [0, 1, 0, 1, 0, 0, 0, 0, 0, 0]]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Once again, it is stochastic\n",
    "[[neuron(inp) for neuron in layer] for i in range(5)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A multi-layer neural network is built by stacking several layers together with outputs of layer _n_ as inputs of layer _n+1_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_network(layer_sizes):\n",
    "    return [create_layer(prev_size, size) for prev_size, size in window(layer_sizes, 2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's build the network displayed below:\n",
    "\n",
    "![JustNN](ws-clipart/just-nn.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "net = create_network([3, 4, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The state of the entire network can be inferred by calculating neuron outputs layer-by-layer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def infer(layers, inp):\n",
    "    state = [inp]\n",
    "    for layer in layers:\n",
    "        state.append([neuron(state[-1]) for neuron in layer])\n",
    "    return state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try to set the input vector to 1,0,0 and infer the state of the network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[True, False, False], [0, 1, 1, 1], [1, 1]]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state = infer(net, [True, False, False])\n",
    "state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### OK, but what does it have to do with graphical models?\n",
    "\n",
    "Great question! The image below is (probably) the most cited graphical model:\n",
    "\n",
    "![Sprinkler](ws-clipart/bayes.png)\n",
    "\n",
    "It is a [Bayesian network](https://en.wikipedia.org/wiki/Bayesian_network) describing the reasoning of a certain British detective who woke up late on a foggy British day to notice that his British lawn is wet. The arcs (edges, arrows) of the graph define conditional probabilities of lawn being wet because of rain or because of a sprinkler.\n",
    "\n",
    "Notice that there is little meaningful difference between nodes in a Bayesian network and our stochastic binary neurons. It is easy to show:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "sprinkler_neuron = StochasticNeuron([-10], 0)\n",
    "grass_neuron = StochasticNeuron([5,5], 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rain: True\n",
      "Sprinkler active: False\n",
      "Grass wet: True\n"
     ]
    }
   ],
   "source": [
    "rain = True\n",
    "print(\"Rain: \" + str(bool(rain)))\n",
    "sprinkler = sprinkler_neuron([rain])\n",
    "print(\"Sprinkler active: \" + str(bool(sprinkler)))\n",
    "grass_wet = grass_neuron([rain, sprinkler])\n",
    "print(\"Grass wet: \" + str(bool(grass_wet)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, essentially, stochastic neural networks and graphical models are 2 languages used to discuss the same mathematical structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### But didn't we set out to solve learning, not inference?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes! Statistical inference and learning are in a sense complementary problems: inference discusses \"how would this neural network/graphical model behave?\" and learning - \"how to build a neural network/graphical model that behaves the way we want?\". Hinton et al defined \"the way we want\" as minimizing _the description length_ of the net's output given the input.\n",
    "\n",
    "Consider the following problem: you have _n_ objects, some of which are more likely to occur than others, defined by some probability distribution. You have to come up with a code for every object so that the expected length of the code will be as short as possible. This problem is less theoretical than it might seem: languages (natural and programming languages alike) attempt to solve whis exact problem by describing concepts that occur often with short words (_dog_) and rare ones with long words (_concatenation_). Languages are far from perfect at it, but [if a perfect language were to exist](https://en.wikipedia.org/wiki/Lojban), every word in it would be of length _-logP_ where _P_ is the probability of this word\n",
    "\n",
    "_-logP_ is called _description length_ (Shannon coding theorem) and it is applicable to any probability distribution as every probability distribution, in theory, defines a coding scheme. Our stochastic neural network is no exception: when an input vector and weights of all connections are fixed, every possible output of the net has a description length."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is derived as follows:\n",
    "\n",
    "The cost _C_ of describing a single neuron is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Neuron DL](ws-clipart/2-neuron-dl.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where _α_ is network state (aggregate of states of all enurons), _s<sub>j</sub><sup>α</sup>_ is the state of neuron _j_ given network state _α_ and _p<sub>j</sub><sup>α</sup>_ is the activation probability of neuron _j_ (P(_s<sub>j</sub><sup>α</sup>_)=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the cost of describing state _α_ in its entirety is simply the cost of describing all the hidden states in all the hidden layers plus the cost of describing the input vector given the hidden states:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Full cost](ws-clipart/3-full-cost.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where _d_ is the input vector. Note that decomposition of _C(α,d)_ into 2 addants is very meaningful. _C(d|α)_ indicates how well the state of the network represents the input vector while _C(α)_ corresponds to the complexity of said state\n",
    "\n",
    "And if we were to minimize this cost (we are) using gradient descent, we would (we will) use the following delta rule derived from the above formula:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Delta rule](ws-clipart/4-gen-delta.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Long time, no Python! Let's implement it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def learn(layers, state, rate):\n",
    "    for layer, (inputs, outputs) in zip(layers, window(state, 2)):\n",
    "        for neuron, output in zip(layer, outputs):\n",
    "            # Yes, I could greatly optimize it by caching probabilities\n",
    "            # As soon as I have some spare time\n",
    "            p = neuron.activation_probability(inputs)\n",
    "            neuron.input_weights -= np.array(\n",
    "                [rate * inpt * (output - p) for inpt, weight in zip(inputs, neuron.input_weights)]\n",
    "                ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clever trick 2: using 2 reverse neural nets to train each other"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several issues that haven't been addressed yet:\n",
    "\n",
    "- Our delta rules requires state _α_ to already exist. And using some dummy state is a bad idea since that will jeopardize the quality of learning\n",
    "- We have a network that outputs represenations of the input vector, but we don't have a good way to reconstruct the input vector from its representations\n",
    "\n",
    "Hinton et al proposed the following solution:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![2-way network](ws-clipart/net.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This architecture can be thought of as either a network where each pair of adjacent layers is connected by 2 sets of connections (one in each direction) or as 2 networks with shared state: one is responsible for recognition: building a representation of input data, another for generation: reconstructing the input vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Learning happens iteratively, each iteration consists of:\n",
    "1. *Recognition inference*. Calculating state _α_.\n",
    "2. *Generative learning*. Changing generative weights so that _d_ is better described by _α_ using the delta rule mentioned above.\n",
    "3. *Generative inference*. Calculating a new (fantasy) _d_\n",
    "4. *Recognition learning*. Changing recognition weights so that _α_ is better described by _d_\n",
    "\n",
    "English -> Python translation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wakesleep(layer_sizes, inputs, learning_rate, iterations):\n",
    "    recognition_connections = create_network(layer_sizes)\n",
    "    generative_connections = create_network(reversed(layer_sizes))\n",
    "    state = [inputs]\n",
    "    \n",
    "    for i in range(iterations):\n",
    "        state = infer(recognition_connections, state[-1]) # That reverses state\n",
    "        learn(generative_connections, reversed(state), learning_rate)\n",
    "        state = infer(generative_connections, state[-1]) # That reverses it once again\n",
    "        learn(recognition_connections, reversed(state), learning_rate)\n",
    "\n",
    "    return reversed(state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1, 1],\n",
       " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n",
       " [0,\n",
       "  0,\n",
       "  1,\n",
       "  1,\n",
       "  1,\n",
       "  1,\n",
       "  1,\n",
       "  1,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  1,\n",
       "  1,\n",
       "  1,\n",
       "  1,\n",
       "  1,\n",
       "  0,\n",
       "  1,\n",
       "  1,\n",
       "  1,\n",
       "  0,\n",
       "  0,\n",
       "  1,\n",
       "  0,\n",
       "  0,\n",
       "  0,\n",
       "  1,\n",
       "  1,\n",
       "  1])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tuple(wakesleep([2, 10, 30], [-1, 1], 1, 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voila!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}