handson-ml/gradient_descent_comparison.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division, unicode_literals\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib nbagg\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 100\n",
    "X = 2*np.random.rand(m, 1)\n",
    "X_b = np.c_[np.ones((m, 1)), X]\n",
    "y = 4 + 3*X + np.random.rand(m, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def batch_gradient_descent():\n",
    "    n_iterations = 1000\n",
    "    learning_rate = 0.05\n",
    "    thetas = np.random.randn(2, 1)\n",
    "    thetas_path = [thetas]\n",
    "    for i in range(n_iterations):\n",
    "        gradients = 2*X_b.T.dot(X_b.dot(thetas) - y)/m\n",
    "        thetas = thetas - learning_rate*gradients\n",
    "        thetas_path.append(thetas)\n",
    "\n",
    "    return thetas_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stochastic_gradient_descent():\n",
    "    n_epochs = 50\n",
    "    t0, t1 = 5, 50\n",
    "    thetas = np.random.randn(2, 1)\n",
    "    thetas_path = [thetas]\n",
    "    for epoch in range(n_epochs):\n",
    "        for i in range(m):\n",
    "            random_index = np.random.randint(m)\n",
    "            xi = X_b[random_index:random_index+1]\n",
    "            yi = y[random_index:random_index+1]\n",
    "            gradients = 2*xi.T.dot(xi.dot(thetas) - yi)\n",
    "            eta = learning_schedule(epoch*m + i, t0, t1)\n",
    "            thetas = thetas - eta*gradients\n",
    "            thetas_path.append(thetas)\n",
    "\n",
    "    return thetas_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mini_batch_gradient_descent():\n",
    "    n_iterations = 50\n",
    "    minibatch_size = 20\n",
    "    t0, t1 = 200, 1000\n",
    "    thetas = np.random.randn(2, 1)\n",
    "    thetas_path = [thetas]\n",
    "    t = 0\n",
    "    for epoch in range(n_iterations):\n",
    "        shuffled_indices = np.random.permutation(m)\n",
    "        X_b_shuffled = X_b[shuffled_indices]\n",
    "        y_shuffled = y[shuffled_indices]\n",
    "        for i in range(0, m, minibatch_size):\n",
    "            t += 1\n",
    "            xi = X_b_shuffled[i:i+minibatch_size]\n",
    "            yi = y_shuffled[i:i+minibatch_size]\n",
    "            gradients = 2*xi.T.dot(xi.dot(thetas) - yi)/minibatch_size\n",
    "            eta = learning_schedule(t, t0, t1)\n",
    "            thetas = thetas - eta*gradients\n",
    "            thetas_path.append(thetas)\n",
    "\n",
    "    return thetas_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_mse(theta):\n",
    "    return np.sum((np.dot(X_b, theta) - y)**2)/m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def learning_schedule(t, t0, t1):\n",
    "    return t0/(t+t1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta0, theta1 = np.meshgrid(np.arange(0, 5, 0.1), np.arange(0, 5, 0.1))\n",
    "r, c = theta0.shape\n",
    "cost_map = np.array([[0 for _ in range(c)] for _ in range(r)])\n",
    "for i in range(r):\n",
    "    for j in range(c):\n",
    "        theta = np.array([theta0[i,j], theta1[i,j]])\n",
    "        cost_map[i,j] = compute_mse(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "exact_solution = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)\n",
    "bgd_thetas = np.array(batch_gradient_descent())\n",
    "sgd_thetas = np.array(stochastic_gradient_descent())\n",
    "mbgd_thetas = np.array(mini_batch_gradient_descent())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "bgd_len = len(bgd_thetas)\n",
    "sgd_len = len(sgd_thetas)\n",
    "mbgd_len = len(mbgd_thetas)\n",
    "n_iter = min(bgd_len, sgd_len, mbgd_len)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 5))\n",
    "data_ax = fig.add_subplot(121)\n",
    "cost_ax = fig.add_subplot(122)\n",
    "\n",
    "cost_ax.plot(exact_solution[0,0], exact_solution[1,0], 'y*')\n",
    "cost_img = cost_ax.pcolor(theta0, theta1, cost_map)\n",
    "fig.colorbar(cost_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def animate(i):\n",
    "    data_ax.cla()\n",
    "    cost_ax.cla()\n",
    "\n",
    "    data_ax.plot(X, y, 'k.')\n",
    "\n",
    "    cost_ax.plot(exact_solution[0,0], exact_solution[1,0], 'y*')\n",
    "    cost_ax.pcolor(theta0, theta1, cost_map)\n",
    "\n",
    "    data_ax.plot(X, X_b.dot(bgd_thetas[i,:]), 'r-')\n",
    "    cost_ax.plot(bgd_thetas[:i,0], bgd_thetas[:i,1], 'r--')\n",
    "\n",
    "    data_ax.plot(X, X_b.dot(sgd_thetas[i,:]), 'g-')\n",
    "    cost_ax.plot(sgd_thetas[:i,0], sgd_thetas[:i,1], 'g--')\n",
    "\n",
    "    data_ax.plot(X, X_b.dot(mbgd_thetas[i,:]), 'b-')\n",
    "    cost_ax.plot(mbgd_thetas[:i,0], mbgd_thetas[:i,1], 'b--')\n",
    "\n",
    "    data_ax.set_xlim([0, 2])\n",
    "    data_ax.set_ylim([0, 15])\n",
    "    cost_ax.set_xlim([0, 5])\n",
    "    cost_ax.set_ylim([0, 5])\n",
    "\n",
    "    data_ax.set_xlabel(r'$x_1$')\n",
    "    data_ax.set_ylabel(r'$y$')\n",
    "    cost_ax.set_xlabel(r'$\\theta_0$')\n",
    "    cost_ax.set_ylabel(r'$\\theta_1$')\n",
    "\n",
    "    data_ax.legend(('Data', 'BGD', 'SGD', 'MBGD'))\n",
    "    cost_ax.legend(('Normal Equation', 'BGD', 'SGD', 'MBGD'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "animation = FuncAnimation(fig, animate, frames=n_iter)\n",
    "plt.show()"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}