Chapter02_Optimization_Basics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ba7c1da4-4d8c-4acb-83ae-63a06f599959",
   "metadata": {},
   "source": [
    "<h1>Gradient and Intro to Gradient Based Optimization</h1>\n",
    "\n",
    "<img src=\"https://img.icons8.com/external-itim2101-lineal-color-itim2101/344/external-professor-life-style-avatar-itim2101-lineal-color-itim2101.png\" alt=\"Instructor\" width=60>  In this section we briefly reviw the framework and mathematics of optimization theory (30 mins). <br>We review the general form of the problem talk about categorization of solutions here. \n",
    "<center>\n",
    "<img alt=\"File:Max paraboloid.svg\" src=\"//upload.wikimedia.org/wikipedia/commons/thumb/7/72/Max_paraboloid.svg/700px-Max_paraboloid.svg.png\" decoding=\"async\" width=\"500\" height=\"560\" srcset=\"//upload.wikimedia.org/wikipedia/commons/thumb/7/72/Max_paraboloid.svg/1050px-Max_paraboloid.svg.png 1.5x, //upload.wikimedia.org/wikipedia/commons/thumb/7/72/Max_paraboloid.svg/1400px-Max_paraboloid.svg.png 2x\" data-file-width=\"700\" data-file-height=\"560\"></center>\n",
    "<ul><li>In optimization we have a function and we want to say at what input the value of the function is maximized (or minimized).\n",
    "\n",
    " </li>\n",
    " <li>It is convenient to establish a general form for all optimization problems.</li>\n",
    " <li>We focus on minimization problem. Every maximization problem can be converted to a minimization (by multiplying it with -1).</li>\n",
    "    <li>In optimization theory: we have two types of minimums: local minimum and global minimum.\n",
    "    <li> We have to notation which is used in optimization: <i>min</i> and <i>argmin</i>.\n",
    "</ul>\n",
    "<table>\n",
    "    <tr>  \n",
    "      <td><img src=\"images/argmin.png\" width=300></td><td><img src=\"https://i.stack.imgur.com/5x9UM.png\"></td>\n",
    "    </tr>\n",
    "</table>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64e91a0c-15ac-4c4b-b5ca-ed780be855e2",
   "metadata": {},
   "source": [
    "<img src=\"https://img.icons8.com/external-soft-fill-juicy-fish/344/external-maths-school-soft-fill-soft-fill-juicy-fish.png\" alt=\"Math Tip\" width=50>\n",
    "<h1> How can we find the minimum?</h1>\n",
    "    <img src=\"https://img.icons8.com/external-itim2101-lineal-color-itim2101/344/external-professor-life-style-avatar-itim2101-lineal-color-itim2101.png\" alt=\"Instructor\" width=50> \n",
    "    <ul>\n",
    "        <li> Search Approach: evaluating function at different values and find the minimum value.</li>\n",
    "        <li> Using Fermat's theory in calculus (<a href=\"https://en.wikipedia.org/wiki/Fermat%27s_theorem_(stationary_points)\"> Visit this</a>)\n",
    "    </ul>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a910787-b7c3-4b97-a63c-6e8a61b6a248",
   "metadata": {},
   "source": [
    "<h2> Fermat's Theorem</h2>\n",
    "<table>\n",
    "    <tr>\n",
    "        <td><img src=\"https://upload.wikimedia.org/wikipedia/commons/f/f3/Pierre_de_Fermat.jpg\" alt=\"Fermat\" width=150></td>\n",
    "        <td><img src=\"images/Fermat.png\"></td>\n",
    "    </tr>\n",
    "</table>\n",
    "<ul>\n",
    "<li>Fermat's theorem states that the (local) exterma points of a function happens in its <b>critical points</b>.</li> \n",
    "<li>A Critical point is a point that either derivative does not exist or it is zero</li>\n",
    "<li> Not every critical point is a local exterma of the function, but it is a candidate. So, you need to check them.</li>\n",
    "</ul>\n",
    "Let's check the following example:\n",
    "<center>\n",
    "<img src=\"https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcSfzLelxmWiK2tGyQmn-j4LXOyEgsldFq2x6Q&usqp=CAU\" width=500>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f102c17-94f3-4cf7-b302-04452b46b489",
   "metadata": {},
   "source": [
    "For multi-variate functions, the derivative is replaced by <b>Gradient</b>. Gradient is vector and each element corresponds to derivative for an input:\n",
    "<table>\n",
    "    <tr>\n",
    "        <td> <img src=\"https://wikimedia.org/api/rest_v1/media/math/render/svg/d632a346cd0677aef80d9fa32f476a5b5bf4dc58\" width=300></td>\n",
    "        <td> <img src=\"images/gradient_vis.png\" width= 600></td>\n",
    "    </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "578d190b-9760-4559-9686-dd4b2c4dc6a0",
   "metadata": {},
   "source": [
    "<h3>Using Gradient in Optimization</h3>\n",
    "For using Fermat's theorem, we need to calculate the gradient and solve the equation f'(x) = 0, to find the critical points. But, solving this equation would be very hard. Instead of solving this equation we take a numerical approach to the problem. We know the gradient is a vector pointing to the direction that the function increases. So, if we move to the opposite direction of gradient, the value of the function decreases. If we continue taking steps in the opposite direction we can reach to a local minima of the function. \n",
    "<table>\n",
    "    <tr>\n",
    "        <td>\n",
    "            <img src=\"https://builtin.com/sites/www.builtin.com/files/styles/ckeditor_optimize/public/inline-images/national/gradient-descent-range.png\">\n",
    "        </td>\n",
    "        <td><img src=\"https://i2.wp.com/ucanalytics.com/blogs/wp-content/uploads/2016/03/bowl-and-ball.jpg\" width =500>\n",
    "        </td>\n",
    "    </tr>\n",
    "</table>\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "197526a8-538b-43b5-847c-b4a6edeb18dd",
   "metadata": {},
   "source": [
    "<h1> Optimization Methods</h1>\n",
    "As pointed before there are two groups of methods for optimizations:\n",
    "<h4>Gradient-Based Methods</h4>\n",
    "In this category of methods, we take small steps in the opposite direction of gradient to get to a local minima. It starts with a random starting point. It calculates the gradient of the function many times.\n",
    "<h4> Gradient-Free Methods</h4>\n",
    "In this category, there is no need to calculate the gradient and the value of the function is evaluated many times. The goal of these approach is to find a better guess for the next point to evaluate the function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6802339e-a8dd-4f4e-bb7b-476457eb823f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import axes3d\n",
    "import numpy as np\n",
    "\n",
    "def f(x,y):\n",
    "    return  (x**2 + y**2)\n",
    "\n",
    "def grad_f(x,y):\n",
    "    return 2 * x, 2 * y\n",
    "\n",
    "\n",
    "x, y = np.meshgrid(np.linspace(-5, 5, 1000), np.linspace(-5, 5, 1000))\n",
    "z = f(x,y)\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(24,12))\n",
    "ax = fig.add_subplot(1,2,1, projection='3d')\n",
    "ax.plot_wireframe(x,y,z, rstride=10, cstride=10)\n",
    "ax.set_title(\"f(x,y)\")\n",
    "ax = fig.add_subplot(1,2,2)\n",
    "ax.contourf(x, y, z,levels=np.linspace(z.min(), z.max(), 20))\n",
    "x, y = np.meshgrid(np.linspace(-5, 5, 10), np.linspace(-5, 5, 10))\n",
    "gx,gy = grad_f(x,y)\n",
    "ax.quiver(x,y, gx, gy)\n",
    "ax.set_title(\"Gradient\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e4428b4-5bb3-4c16-8050-90d7e4930701",
   "metadata": {},
   "source": [
    "<h2>Basic Gradient Decent Algorithm</h2>\n",
    "In Gradient Decent(GD), we start at a random point and then we take small steps in opposite direction of gradient until we reach to optimal point. GD is an iterative algorithm. At each step the solution(θ) is updated using the following rule: <br>\n",
    "<center><img src=\"images/vanilla_gd.png\" alt=\"x_new = x_old - s. Grad(f(x_old))\" width = 250></center>\n",
    "Here, f is the objective function (the function we are minimizing) and x = [x1, ..., xn] are inputs.\n",
    "<br> At the optimal point the graident become zero and no updates are possible. But it is unlikely to hit that point in practice. So, some sort of stop condition is needed.  We will study GD in further details later. Here, we develop a basic implementation of method here to provide context for coming topics.\n",
    "\n",
    "<br><br>\n",
    "The function f(x) = x^2 has one global minimum at zero. Let's find it using GD. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1acd32ab-1093-4b56-a588-57058ec1adc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "grad = lambda x: 8 *x\n",
    "f = lambda x: 4.*x**2\n",
    "x_0 = 11\n",
    "alpha = .1\n",
    "x_opt = x_0\n",
    "num_iter = 50\n",
    "for i in range(num_iter):\n",
    "    x_opt = x_opt - alpha * grad(x_opt)\n",
    "print( f\"Optimal point: {x_opt}, Objective function value:{f(x_opt)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6a6b644-bc10-477c-ab1d-04565b423be5",
   "metadata": {},
   "source": [
    "Let's add some visualization to understand what is going on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0709308c-802a-4d8d-9b76-616d7389b80c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "# Function Definitions\n",
    "grad = lambda x: 8 *x\n",
    "f = lambda x:  4*x**2\n",
    "x_0 = 15\n",
    "alpha = .03\n",
    "x_opt = x_0\n",
    "num_iter = 50\n",
    "\n",
    "x_opt_hist = [x_opt]\n",
    "f_opt_hist = [f(x_opt)]\n",
    "# Optimization loop\n",
    "for i in range(num_iter):\n",
    "    x_opt = x_opt - alpha * grad(x_opt)\n",
    "    x_opt_hist.append(x_opt)\n",
    "    f_opt_hist.append(f(x_opt))\n",
    "# Plotting part\n",
    "animate = lambda i: l.set_data(x_opt_hist[:i], f_opt_hist[:i])\n",
    "x = np.linspace(-20,20, 4000)\n",
    "f_x = f(x)\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, f_x)\n",
    "l, = ax.plot(x_opt, f(x_opt), 'o')\n",
    "for i in range(num_iter):\n",
    "    animate(i)\n",
    "    clear_output(wait=True)\n",
    "    display(fig)\n",
    "    \n",
    "clear_output(wait=True)\n",
    "print( f\"Optimal point: {x_opt}, Objective function value:{f(x_opt)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72f078f2-7bc2-4284-9de1-f7c69eb8bbf4",
   "metadata": {},
   "source": [
    "<img src= \"https://img.icons8.com/external-flaticons-flat-flat-icons/344/external-question-100-most-used-icons-flaticons-flat-flat-icons.png\" width=70> __Question__:<br>\n",
    "Apply the following changes to the above program and check the outcome for each step:<br>\n",
    "<ul>\n",
    "    <li> num_iter = 100</li>\n",
    "    <li> alpha = .001</li>\n",
    "    <li> alpha = 0.005 </li>\n",
    "    <li> alpha = 0.22 </li>\n",
    "    <li> alpha = 0.25</li>\n",
    "    <li> alpha = 0.26</li>\n",
    "</ul>\n",
    "Can you explain the outcome? (Using the animation cell can be helpful!)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b0df2df-ec7a-4620-a6ba-f48da1b2eb5b",
   "metadata": {},
   "source": [
    "## Calculations of Gradient\n",
    "In the above example we define the gradient as function ourself. But, is there a way to ask python to do it?<br>\n",
    "The answer is yes. Basically, there are three categorise of solutions.\n",
    "<ul>\n",
    "    <li>The symbolic derivative of a function.</li>\n",
    "    <li>Compute numerical derivatives of a function defined only by a sequence of data points.</li>\n",
    "    <li>Compute numerical derivatives of a analytically supplied function.</li>\n",
    "</ul>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaa5641e-bb2b-4fb8-8976-72231b9e130c",
   "metadata": {},
   "source": [
    "<img src=\"https://www.sympy.org/static/images/logo.png\" width=120> <font size=20><b>Sympy</b></font> <br> Sympy is a <a href=\"https://en.wikipedia.org/wiki/Computer_algebra\">symbolic computation</a> package for Python. Sympy is writen in Python and defines the symbolic language as Python itself. It is very light and powerful. <br>\n",
    "Use this command to install it from Anaconda ```conda install -c anaconda sympy```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c95483d-b3ef-4812-80a2-459ee4fbf12a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75db0723-fc34-4016-8f2a-656ee298b71a",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_s = sympy.sqrt(2)\n",
    "x_m = math.sqrt(2)\n",
    "print(\"Symbolic: \", x_s)\n",
    "print(\"Math: \", x_m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fac9d349-3803-494a-98b9-53219832560a",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Symbolic: \", x_s**2)\n",
    "print(\"Math: \", x_m**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61abd5ad-1222-412c-8152-de81dad606d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sympy.symbols(\"x\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4278b7b2-5fb1-4cf9-b133-cd9419e76257",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = x**2\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5eae8626-de7f-45c8-a38b-e262a1c2e751",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = f + sympy.sin(x)\n",
    "g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1db26d0-f303-4b93-a8db-a42ae8168199",
   "metadata": {},
   "outputs": [],
   "source": [
    "val = g.subs([(x,sympy.pi)])\n",
    "val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aac53c08-4000-4984-8c8f-7fab8338a090",
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.N(val) # pi * pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4de7a296-bd2d-4d33-b11a-c2108255c692",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = sympy.symbols(\"x y\")\n",
    "exper = (x + 2)*(y+4)\n",
    "exper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f5e75ee-9ed2-4935-9b15-d48ebe23d41c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.expand(exper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99d70c62-2993-4e7c-9514-c5f4cb910397",
   "metadata": {},
   "outputs": [],
   "source": [
    "exper = x*x + 2*x*y + y*y\n",
    "exper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "767ac42b-4b4c-459b-87ca-cd999e265303",
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.factor(exper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd55245d-7892-45cc-8507-8774988d774b",
   "metadata": {},
   "outputs": [],
   "source": [
    "exper = x*x + 6*x*y + y*y + 4*x*y + x*y*y*y - 3*y*y\n",
    "sympy.simplify(exper)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff6a1654-8bcc-40d6-bdbc-5053853d67df",
   "metadata": {},
   "source": [
    "#### Calculus Using Symbolic Operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23202c38-d098-4929-908d-c5a587e6ad36",
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.diff(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae47793f-706a-42d1-8c47-1b6a8929f69c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.integrate(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2b4c31f-3f8d-4691-b564-35b12000e071",
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.diff(g, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ea389c7-4345-4224-8a12-aa3a0ebc738a",
   "metadata": {},
   "outputs": [],
   "source": [
    "x,y = sympy.symbols(\"x y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8928b823-0d74-452e-a4a6-2903528a7a4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = x**2 + sympy.cos(y)\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd944c64-4d54-49ff-94ee-66b77e08dbe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.diff(f, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9c6031b-b760-473e-b1fa-dc5ea9679174",
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.diff(f, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e037de67-3620-49a6-a6ec-7f440e072518",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculating gradient\n",
    "from sympy.tensor.array import derive_by_array\n",
    "f = x**2 + y**2\n",
    "derive_by_array(f, (x,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29d3ff33-5579-4e40-b492-5dadb5f779c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A complex one\n",
    "f = x**3*sympy.cos(y*sympy.sin(x**2+y**3 + 2*x*y**2))\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccf26130-85ca-4480-b992-50f4a699237e",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_prime = derive_by_array(f, (x,y))\n",
    "f_prime"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44609751-cb40-487f-aa8d-e40456db627b",
   "metadata": {},
   "source": [
    "<img src=\"images/autograd_logo.png\" width=200> <font size=20><b>Autograd</b></font><br>\n",
    "Autograd is an <a href= \"https://en.wikipedia.org/wiki/Automatic_differentiation\"> automatic differentiation</a> package for python. Each computer program is built of simple operations like add, multiply, loops and so on. In AD, a graph of necessary operations for function is built and we know the deravative of building blocks. So, it is possible to calculate the derivative using this graph. Autograd is used in some of DeepLearning packages like Torch. \n",
    "<br>\n",
    "For installing autograd use: ```conda install -c conda-forge autograd``` <br>\n",
    "<br>\n",
    "If you are curios how Autograd/Jax works check <a href=\"https://github.com/HIPS/autograd/blob/master/docs/tutorial.md\">this</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f8bbe27-3d04-4076-a008-f98cd5718edf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from autograd import numpy as np\n",
    "from autograd import grad, value_and_grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8209383-0746-4318-b276-f92699887a73",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The numpy in JAX is the same as numpy but JAX is aware of details.\n",
    "np.cos(np.pi/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a1fb1f7-9da0-4f5e-8a37-1ae50be02b1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_prime = grad(np.cos)  # -sin(x)\n",
    "f_prime(np.pi/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4acde913-dce7-4540-896d-5f0c56c1f329",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x: x[0]**3 + x[1]**4\n",
    "f(np.array([1,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63b361d4-4e37-49eb-a622-b2ba37b1fe99",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_prime = grad(f)\n",
    "f_prime(np.array([1.0,1.0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "503117a3-3ff8-4e05-a282-07aa2bc4b934",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Beware of types (type casting is not differentiable).\n",
    "# The following line raise an exception\n",
    "f_prime(np.array([1,1])) # Here int object is not differentiable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53b912cb-d07a-43b7-9be0-494e98e1eab8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the dtype\n",
    "f_prime(np.array([1,1], dtype=np.float32))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9deab203-2f9c-4062-9edb-b0e1a007d4ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's work on something more interesting\n",
    "def rosenbrock(x):\n",
    "    return 100*(x[1] - x[0]**2)**2 + (1 - x[0])**2\n",
    "\n",
    "# Build a function that also returns gradients using autograd.\n",
    "rosenbrock_with_grad = value_and_grad(rosenbrock)\n",
    "rosenbrock_with_grad(np.array([0,0], dtype=np.float32))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77543076-45a8-486e-bbcb-53b477336c91",
   "metadata": {},
   "source": [
    "<img src= \"https://img.icons8.com/external-flaticons-flat-flat-icons/344/external-question-100-most-used-icons-flaticons-flat-flat-icons.png\" width=70> __Question__:<br>\n",
    "Write a program to find the minmum of the function f(x,y,z) = x^2 + y^2 + z^2 using gradient decent. Calculate the gradient using autograd library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc45bda5-83c0-4df1-8e51-813bc026e402",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gradient of a more complex function (Autograd example)\n",
    "def taylor_sine(x):  # Taylor approximation to sine function\n",
    "    ans = currterm = x\n",
    "    i = 0\n",
    "    while np.abs(currterm) > 0.001: # There is a while loop here !!!\n",
    "        currterm = -currterm * x**2 / ((2 * i + 3) * (2 * i + 2))\n",
    "        ans = ans + currterm\n",
    "        i += 1\n",
    "    return ans\n",
    "\n",
    "grad_sine = grad(taylor_sine)\n",
    "print (\"Gradient of sin(pi) is\", grad_sine(np.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7afa599f-4f67-4936-83a8-42212fcc9d7b",
   "metadata": {},
   "source": [
    "<img src=\"https://img.icons8.com/external-kosonicon-lineal-color-kosonicon/344/external-lab-tool-back-to-school-kosonicon-lineal-color-kosonicon.png\" alt=\"LAB\" width=80 > <b>Training Linear Regression</b>\n",
    "In linear regression a continuous function dependent variable is written as a linear combination of independent variables. <br>\n",
    "<center><img src=\"images/lr.png\" width=500></center>\n",
    "One can find the weights (β_i) using Mean Square Error objective through an optimizaiton by GD. <br>\n",
    "\n",
    "Let's assume we are getting data from a black box which generates data like (x,y) where y = 2x + 1 + n where n is a random noise. Write a program for following steps:\n",
    "<ul>\n",
    "    <li> Write a function to generate samples from model according to above equation</li>\n",
    "    <li>Draw 100 samples </li>\n",
    "    <li> Write a function to calculate the MSE error for given parameters β0, β1</li>\n",
    "    <li> Use auto grad to calculate the gradient of MSE.</li>\n",
    "    <li> and use GD to find β0, β1 </li>\n",
    "    <li> Plot data points and fitted line</li>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7a557ff-90cb-45e9-99c9-c9b4ee68b0b6",
   "metadata": {},
   "source": [
    "<img src=\"https://img.icons8.com/external-flaticons-lineal-color-flat-icons/344/external-coffee-cup-bakery-flaticons-lineal-color-flat-icons.png\" alt=\"icon\" width=80> <b>Takehome Question:</b><br>\n",
    "\n",
    "Repeat the lab experiment for a logistic regression problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20814668-b8ba-45ad-92b0-189a9f8bc4aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can use auto grad for calculation of higher order deravatives\n",
    "from autograd import elementwise_grad as egrad\n",
    "\n",
    "def tanh(x):\n",
    "    return (1.0 - np.exp(-x))  / (1.0 + np.exp(-x))\n",
    "\n",
    "fig = plt.figure(figsize=(6,6))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "\n",
    "x = np.linspace(-7, 7, 200)\n",
    "_=ax.plot(x, tanh(x),\n",
    "         x, egrad(tanh)(x),                                     # first derivative\n",
    "         x, egrad(egrad(tanh))(x),                              # second derivative\n",
    "         x, egrad(egrad(egrad(tanh)))(x),                       # third derivative\n",
    "         x, egrad(egrad(egrad(egrad(tanh))))(x),                # fourth derivative\n",
    "         x, egrad(egrad(egrad(egrad(egrad(tanh)))))(x),         # fifth derivative\n",
    "         x, egrad(egrad(egrad(egrad(egrad(egrad(tanh))))))(x))  # sixth derivative\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2147785c-f7ef-4708-8758-bc6c94bb054b",
   "metadata": {},
   "source": [
    "<img src=\"https://img.icons8.com/color/344/light.png\" width=70> \n",
    "<b> <a href=\"https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant\">Jacobian </a></b> is extension of definition of gradient to a vector function. A vector function is vector of scalar functions. For finding the jacobian you can iterate over functions in vector or you can you use ready jacobian function in autograd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bb39993-f2c7-4c93-8acf-9f2e2957b467",
   "metadata": {},
   "outputs": [],
   "source": [
    "from autograd import jacobian\n",
    "scalar_funs = [\n",
    "        lambda x: np.sum(x**3),\n",
    "        lambda x: np.prod(np.sin(x) + np.sin(x)),\n",
    "        lambda x: grad(lambda y: np.exp(y) * np.tanh(x[0]))(x[1])]\n",
    "\n",
    "vector_fun = lambda x: np.array([f(x) for f in scalar_funs])\n",
    "\n",
    "x = np.random.randn(5)\n",
    "jac = jacobian(vector_fun)(x)\n",
    "grads = [grad(f)(x) for f in scalar_funs]\n",
    "\n",
    "np.allclose(jac, np.vstack(grads))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a87f474a-77e6-48f3-8002-b5ef3a14f656",
   "metadata": {},
   "source": [
    "<img src=\"https://img.icons8.com/color/344/light.png\" width=70>\n",
    "Also, the second derivative of a scalar function is called <a href=\"https://en.wikipedia.org/wiki/Hessian_matrix\"><b>Hessian</a></b>. \n",
    "You can caluclate the hessian matrix by iteration or use the corresponding autograd function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81958189-f58e-4297-bab2-8a302efaf9a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from autograd import hessian\n",
    "fun = lambda x: x[0]**2 + x[1]**3\n",
    "H = hessian(fun)\n",
    "H_indirect = jacobian(grad(fun))\n",
    "x = np.array([1,1],dtype=np.float32)\n",
    "\n",
    "np.allclose(H(x), H_indirect(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40b6249d-ecff-43a2-93b1-6bb892eda5d0",
   "metadata": {},
   "source": [
    "<img src=\"https://img.icons8.com/external-flaticons-lineal-color-flat-icons/344/external-coffee-cup-bakery-flaticons-lineal-color-flat-icons.png\" alt=\"icon\" width=80> <b>Takehome Question:</b><br>\n",
    "    Check the following example. If you want autograd not looking at a function, you should use @primitive decorative. Also, you need to define the gradient function separately. Check the code below and try to understand what is going on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abd16b8f-49b6-4616-b2dc-d17e24b07f0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Working interactively with autograd\n",
    "from autograd.extend import primitive, defvjp\n",
    "from autograd.test_util import check_grads\n",
    "# The @primitive tells autograd not look at the function body. So, if the gradient is needed you have to provide it separately.\n",
    "@primitive\n",
    "def logsumexp(x):\n",
    "    \"\"\"Numerically stable log(sum(exp(x))), also defined in scipy.special\"\"\"\n",
    "    max_x = np.max(x)\n",
    "    return max_x + np.log(np.sum(np.exp(x - max_x)))\n",
    "\n",
    "def logsumexp_vjp(ans, x):\n",
    "    # If you want to be able to take higher-order derivatives, then all the\n",
    "    # code inside this function must be itself differentiable by Autograd.\n",
    "    \n",
    "    # This closure multiplies g with the Jacobian of logsumexp (d_ans/d_x).\n",
    "    # Because Autograd uses reverse-mode differentiation, g contains\n",
    "    # the gradient of the objective w.r.t. ans, the output of logsumexp.\n",
    "    # This returned VJP function doesn't close over `x`, so Python can\n",
    "    # garbage-collect `x` if there are no references to it elsewhere.\n",
    "    x_shape = x.shape\n",
    "    return lambda g: np.full(x_shape, g) * np.exp(x - np.full(x_shape, ans))\n",
    "\n",
    "# Now we tell Autograd that logsumexmp has a gradient-making function.\n",
    "defvjp(logsumexp, logsumexp_vjp)\n",
    "\n",
    "def example_func(y):\n",
    "    z = y**2\n",
    "    lse = logsumexp(z)\n",
    "    return np.sum(lse)\n",
    "\n",
    "grad_of_example = grad(example_func)\n",
    "print(\"Gradient: \\n\", grad_of_example(np.random.randn(10)))\n",
    "\n",
    "# Check the gradients numerically, just to be safe.\n",
    "check_grads(example_func, modes=['rev'])(np.random.randn(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd304dc0-8caa-400c-8c75-7faa779e0f74",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/google/jax/main/images/jax_logo_250px.png\" width=100><br>\n",
    "JAX is NumPy on the CPU, GPU, and TPU, with great automatic differentiation for high-performance machine learning research. JAX is the combination of Autograd with XLA. XLA (Accelerated Linear Algebra) itself is a domain-specific compiler for linear algebra. JAX and XLA are supported by Google for Tensorflow. You can work with Autograd directly, or you can use its faster implementation JAX. <br>\n",
    "For getting the latest jaxlib use: ```conda install -c conda-forge jaxlib```<br>\n",
    "For installing JAX use: ```conda install -c conda-forge jax```<br> \n",
    "For using JAX, everything is pretty much the same as autograd. It just differs in terms of computations. Specifically, if you are using GPUs and TPUs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bbcab70-95fa-46b2-b0d2-c22164040c84",
   "metadata": {},
   "source": [
    "### numdifftools\n",
    "numdifftools is a package for <a href=\"https://en.wikipedia.org/wiki/Numerical_differentiation\"> numerical gradient</a> calculation (estimation). Based on circumstances,  numdifftools can be faster or slower than automatic differentiation. Also, since it estimates the gradient, it might have some error (usually very small).\n",
    "For installation use ```conda install -c conda-forge numdifftools``` <br>\n",
    "\n",
    "If you are interested to know how it works check this <a href=\"https://numdifftools.readthedocs.io/en/latest/src/numerical/derivest.html\"> page</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60367bbf-351f-4616-b1c0-5dbd0a2a932b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numdifftools as nd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d98f9c5-d61c-4a08-b566-1f1359092038",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np # numdifftools does not have an implementation of numpy\n",
    "f = np.sin\n",
    "f_prime = nd.Derivative(f)\n",
    "f_prime(np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb70fd07-174b-41cb-b2ac-d307d96dd4e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# higher order derivatives\n",
    "x = np.linspace(-2, 2, 100)\n",
    "for i in range(10):\n",
    "    df = nd.Derivative(np.tanh, n=i)\n",
    "    y = df(x)\n",
    "    h = plt.plot(x, y/np.abs(y).max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c74a4f7f-b91f-4cb4-935d-fe3ccd643552",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hessian\n",
    "rosen = lambda x: (1-x[0])**2 + 105.*(x[1]-x[0]**2)**2\n",
    "h = nd.Hessian(rosen) \n",
    "x_0 = np.array([1,1])\n",
    "h(x_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dc9f317-7be1-4c9a-837b-e1ffb150df43",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.rand(5,3)\n",
    "b = np.random.rand(5)\n",
    "fun = lambda x: np.dot(x, A.T) - b\n",
    "x = np.random.rand(3)\n",
    "jac = nd.Jacobian(fun)(x)\n",
    "jac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "970aab45-e93b-41c5-be53-04062d3073cd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fun = lambda xy: np.r_[xy[0]**2, np.cos(xy[0] - xy[1])] # np.r_ makes a row out of slicings. Here it is used to make an array object. you can replace it with np.array.\n",
    "jac = nd.Jacobian(fun)([-2, -3])\n",
    "jac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b62d998-36f3-4f07-9b8a-84dc6fb72008",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For simple functions autograd is faster. If you think autograd is not fast for your app, try numdifftools to see what it can give you.\n",
    "import numpy as np\n",
    "x = np.linspace(0, np.pi/2, 100)\n",
    "f_prime_1 = nd.Derivative(np.sin)\n",
    "from autograd import numpy as np\n",
    "f_prime_2 = grad(np.sin)\n",
    "%timeit [f_prime_1(t) for t in x]\n",
    "%timeit [f_prime_2(t) for t in x]\n",
    "del x"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}