{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Regression with Gradient Descent Algorithm\n",
    "\n",
    "This notebook demonstrates the implementation of linear regression with gradient descent algorithm.  \n",
    "\n",
    "Consider the following implementation of the gradient descent loop with NumPy arrays based upon [1]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def gradient_descent_numpy(X, Y, theta, alpha, num_iters):\n",
    "    m = Y.shape[0]\n",
    "\n",
    "    theta_x = 0.0\n",
    "    theta_y = 0.0\n",
    "\n",
    "    for i in range(num_iters):\n",
    "        predict = theta_x + theta_y * X\n",
    "        err_x = (predict - Y)\n",
    "        err_y = (predict - Y) * X\n",
    "        theta_x = theta_x - alpha * (1.0 / m) * err_x.sum()\n",
    "        theta_y = theta_y - alpha * (1.0 / m) * err_y.sum()\n",
    "\n",
    "    theta[0] = theta_x\n",
    "    theta[1] = theta_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To speedup this implementation with Numba, we need to add the `@jit` decorator to annotate the function signature.  Then, we need to expand the NumPy array expressions into a loop.  The resulting code is shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numba import jit, f8, int32, void\n",
    "\n",
    "@jit(void(f8[:], f8[:], f8[:], f8, int32))\n",
    "def gradient_descent_numba(X, Y, theta, alpha, num_iters):\n",
    "    m = Y.shape[0]\n",
    "\n",
    "    theta_x = 0.0\n",
    "    theta_y = 0.0\n",
    "\n",
    "    for i in range(num_iters):\n",
    "        err_acc_x = 0.0\n",
    "        err_acc_y = 0.0\n",
    "        for j in range(X.shape[0]):\n",
    "            predict = theta_x + theta_y * X[j]\n",
    "            err_acc_x += predict - Y[j]\n",
    "            err_acc_y += (predict - Y[j]) * X[j]\n",
    "        theta_x = theta_x - alpha * (1.0 / m) * err_acc_x\n",
    "        theta_y = theta_y - alpha * (1.0 / m) * err_acc_y\n",
    "\n",
    "    theta[0] = theta_x\n",
    "    theta[1] = theta_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The rest of the code generates some artificial data to test our linear regression algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pylab\n",
    "from timeit import default_timer as timer\n",
    "\n",
    "def populate_data(N, slope, intercept, stdev=10.0):\n",
    "    noise = stdev*np.random.randn(N)\n",
    "    X = np.arange(N, dtype=np.float64)\n",
    "    Y = noise + (slope * X + intercept)\n",
    "    return X, Y\n",
    "\n",
    "def run(gradient_descent, X, Y, iterations=10000, alpha=1e-6):\n",
    "    theta = np.empty(2, dtype=X.dtype)\n",
    "\n",
    "    ts = timer()\n",
    "    gradient_descent(X, Y, theta, alpha, iterations)\n",
    "    te = timer()\n",
    "\n",
    "    timing = te - ts\n",
    "\n",
    "    print(\"x-offset = {}    slope = {}\".format(*theta))\n",
    "    print(\"time elapsed: {} s\".format(timing))\n",
    "\n",
    "    return theta, timing\n",
    "\n",
    "\n",
    "def plot(X, theta, c='r'):\n",
    "    result = theta[0] + theta[1] * X\n",
    "    pylab.plot(X, result, c=c, linewidth=2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will a benchmark with 50 elements to compare the pure python version against the numba version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = 10\n",
    "X, Y = populate_data(N, 3, 10)\n",
    "pylab.scatter(X, Y, marker='o', c='b')\n",
    "pylab.title('Linear Regression')\n",
    "\n",
    "print('NumPy'.center(30, '-'))\n",
    "theta_python, time_python = run(gradient_descent_numpy, X, Y)\n",
    "\n",
    "print('Numba'.center(30, '-'))\n",
    "theta_numba, time_numba  = run(gradient_descent_numba, X, Y)\n",
    "\n",
    "# make sure all method yields the same result\n",
    "assert np.allclose(theta_python, theta_numba)\n",
    "\n",
    "print('Summary'.center(30, '='))\n",
    "print('Numba speedup %.1fx' % (time_python / time_numba))\n",
    "\n",
    "plot(X, theta_numba, c='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## References\n",
    "\n",
    "[1] http://aimotion.blogspot.com/2011/10/machine-learning-with-python-linear.html"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}