unit07_optim/demo1_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": 3,
   "metadata": {},
   "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 below computes `f(w)` its gradient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f     = 260.000000\n",
      "fgrad = [132 192]\n"
     ]
    }
   ],
   "source": [
    "def feval(w):\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])\n",
    "    \n",
    "    return f, fgrad\n",
    "\n",
    "# Point to evaluate \n",
    "w = np.array([2,4])\n",
    "f, fgrad = feval(w)\n",
    "\n",
    "print('f     = %f' % f)\n",
    "print('fgrad = ' + str(fgrad))\n"
   ]
  },
  {
   "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": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Jgrad = [15.63967135 -2.77296801]\n"
     ]
    }
   ],
   "source": [
    "def Jeval(w):\n",
    "    \n",
    "    # Unpack vector\n",
    "    a = w[0]\n",
    "    b = w[1]\n",
    "    \n",
    "    # Compute the loss function\n",
    "    yerr = y-a*np.exp(-b*x)\n",
    "    J = 0.5*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])\n",
    "    return J, Jgrad\n",
    "\n",
    "# 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",
    "w = np.array([a,b])\n",
    "\n",
    "J, Jgrad = Jeval(w)\n",
    "print('Jgrad = ' + str(Jgrad))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3:  Loss Function for a Log-Linear Model\n",
    "\n",
    "Now suppose we have a logarithmic linear model:\n",
    "\n",
    "    yhat = log(z),  z = w_0 + \\sum_j w_jx_j\n",
    "    \n",
    "If we have data `X, y`, the prediction on sample `i` is:\n",
    "\n",
    "    yhat_i = log(z_i),  z_i = w_0 + \\sum_j w_j*X_{ij}\n",
    " \n",
    "    \n",
    "Suppose we use MSE loss:\n",
    "\n",
    "    J(w) = \\sum_i (y_i - yhat_i )**2\n",
    "    \n",
    "To compute the components `dJ/dw_j`, first write `z = Aw` where `A=[1 X]`, the matrix with ones on the first column.\n",
    "Then, `dz_i/dw_j = A_{ij}`.  Also, `dyhat_i / dz_i = 1/z_i`, so with the multi-variable chain rule: \n",
    "\n",
    "    dJ/dw_j = \\sum_i dJ/dyhat_i * dyhat_i / dw_j \n",
    "            = \\sum_i (dJ/dyhat_i) * (dyhat_i / dz_i) * (dz_i / dw_j) = \n",
    "            = 2*\\sum_i (yhat_i - y_i) * 1/z_i* A_{ij}\n",
    "    \n",
    "We can implement the loss and gradient computation as follows:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Jeval(w,X,y):\n",
    "\n",
    "    # Create matrix A=[1 X]\n",
    "    n = X.shape[0]\n",
    "    A = np.column_stack((np.ones(n), X))\n",
    "\n",
    "    # Compute function\n",
    "    z = A.dot(w)\n",
    "    yhat = np.log(z)\n",
    "    J = np.sum((y-yhat)**2)\n",
    "    \n",
    "    # Compute gradient\n",
    "    dJ_dz = 2*(yhat-y)/z\n",
    "    Jgrad = A.T.dot(dJ_dz)\n",
    "    \n",
    "    return J, Jgrad  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whenever you implement a gradient, you should always check the gradients.  Errors in the gradient are the number \n",
    "\n",
    "*  Take two points `w0` and `w1` close to one another\n",
    "*  Verify that `J(w1)-J(w0)` is close to `Jgrad(w0).dot(w1-w0)`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Actual difference:      -1.1895e-03\n",
      "Estimated difference:   -1.1896e-03\n"
     ]
    }
   ],
   "source": [
    "# Generate random positive data\n",
    "n = 100\n",
    "d = 5\n",
    "X = np.random.uniform(0,1,(n,d))\n",
    "w0 = np.random.uniform(0,1,(d+1,))\n",
    "y = np.random.uniform(0,2,(n,))\n",
    "\n",
    "# Compute function and gradient at point w0\n",
    "J0, Jgrad0 = Jeval(w0,X,y)\n",
    "\n",
    "# Take a small perturbation\n",
    "step = 1e-4\n",
    "w1 = w0 + step*np.random.normal(0,1,(d+1,))\n",
    "\n",
    "# Evaluate the function at perturbed point\n",
    "J1, Jgrad1 = Jeval(w1,X,y)\n",
    "\n",
    "dJ = J1-J0\n",
    "dJ_est = Jgrad0.dot(w1-w0)\n",
    "print('Actual difference:     %12.4e' % dJ)\n",
    "print('Estimated difference:  %12.4e' % dJ_est)"
   ]
  },
  {
   "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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feval(W,a,b):\n",
    "    # Function\n",
    "    f = a.dot(W.dot(b))\n",
    "\n",
    "    # Gradient -- Use python broadcasting\n",
    "    fgrad = a[:,None]*b[None,:]\n",
    "    \n",
    "    return f, fgrad\n",
    "    \n",
    "# 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",
    "f, fgrad = feval(W,a,b)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## In-Class Exercise:  An Exponential Model\n",
    "\n",
    "Consider a model,\n",
    "\n",
    "    yhat = w[0]*exp(-w[1]*(x-w[2])**2/2)\n",
    "   \n",
    "where the parameter `w[2] > 0` is positive.\n",
    "\n",
    "Now, suppose that, given data `x` and `y`, we want to minimize the MSE loss function,\n",
    "\n",
    "    J = mean( (y[i] - yhat[i])**2 ) \n",
    "   \n",
    "Complete the following function to compute `J` and its gradient for parameters `w` and data (`x,y`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Jeval(w,X,y):\n",
    "    # TODO    \n",
    "    return J, Jgrad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test the gradient:\n",
    "\n",
    "* Generate random `x` and `y` vectors with length `n=100`.\n",
    "* Compute `J0, Jgrad0`, the loss function and its gradient, at the `w0` value below\n",
    "* Compute some `w1`close to `w0`. \n",
    "* Compute `J1, Jgrad1`, the loss function and its gradient, at `w1`.\n",
    "* Check that `J1-J0` is approximately `Jgrad0.dot(w1-w0)`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "w0 = np.array([0.4, 2, 2])"
   ]
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}