unit05_optim/demo05_1_computing_gradients.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demo:  Computing Gradients\n",
    "\n",
    "Most numerical optimization methods require that we compute gradients of the loss function that we are attempting to minimize.  In this demo, we illustrate how to compute gradients efficiently in python for a few simple examples.  As much as possible, we avoid for loops for fast implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1:  A Simple Vector-Input Function\n",
    "\n",
    "Suppose `f(w) = w_0^2 + 2w_0w_1^3`.  Then the function and gradient at `w=[2,4]` can be computed as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Point to evaluate \n",
    "w = np.array([2,4])\n",
    "\n",
    "# Function\n",
    "f = w[0]**2 + 2*w[0]*(w[1]**3)\n",
    "\n",
    "# Gradient\n",
    "df0 = 2*w[0]+2*(w[1]**3)\n",
    "df1 = 6*w[0]*(w[1]**2)\n",
    "fgrad = np.array([df0, df1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2:  Non-Linear Least Squares for an Exponential Model\n",
    "\n",
    "Consider an exponential model \n",
    "\n",
    "    yhat = a*exp(-b*x)\n",
    "    \n",
    "for parameters `w=[a,b]`.  Given training data `(x[i],y[i])` a natural loss function is given by\n",
    "\n",
    "    J(w) := \\sum_i (y[i] - yhat[i])**2,   yhat[i] = a*exp(-b*x[i])\n",
    "    \n",
    "The following code computes the the loss function `J(w)` and its gradient `dJ/dw`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate some random data\n",
    "ny = 100\n",
    "y = np.random.randn(ny)\n",
    "x = np.random.rand(ny)\n",
    "\n",
    "# Some arbitrary parameters \n",
    "# to compute the gradient at\n",
    "a = 1\n",
    "b = 2\n",
    "\n",
    "# Compute the loss function\n",
    "yerr = y-a*np.exp(-b*x)\n",
    "J = np.sum(yerr**2)\n",
    "\n",
    "# Compute the gradient\n",
    "dJ_da = -np.sum( yerr*np.exp(-b*x))\n",
    "dJ_db = np.sum( yerr*a*x*np.exp(-b*x))\n",
    "Jgrad = np.array([dJ_da, dJ_db])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3:  A Function of a Matrix.\n",
    "\n",
    "Suppose `f(W) = a'*W*b`.  Then, `fgrad(W) = a*b.T`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Some random data\n",
    "m = 4\n",
    "n = 3\n",
    "W = np.random.randn(m,n)\n",
    "a = np.random.randn(m)\n",
    "b = np.random.randn(n)\n",
    "\n",
    "# Function\n",
    "f = a.dot(W.dot(b))\n",
    "\n",
    "# Gradient -- Use python broadcasting\n",
    "fgrad = a[:,None]*b[None,:]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}