2 Gradient Descent.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thinking in tensors, writing in PyTorch\n",
    "\n",
    "A hands-on course by [Piotr Migdał](https://p.migdal.pl) (2019-2022).\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/stared/thinking-in-tensors-writing-in-pytorch/blob/master/2%20Gradient%20Descent.ipynb\" target=\"_parent\">\n",
    "    <img src=\"https://colab.research.google.com/assets/colab-badge.svg\"/>\n",
    "</a>\n",
    "\n",
    "## Notebook 2: Gradient descent\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/7/79/Topographic_map_example.png/310px-Topographic_map_example.png)\n",
    "\n",
    "[Topographic map of Stowe, Vermont from Wikipedia](https://en.wikipedia.org/wiki/Topographic_map)\n",
    "\n",
    "> X: I want to learn Deep Learning.  \n",
    "> Me: Do you know what is gradient?  \n",
    "> X: Yes.  \n",
    "> Me: Then, it an easy way downhill!\n",
    "\n",
    "In this notebook, we explore:\n",
    "\n",
    "* What is gradient?\n",
    "* What is gradient descent?\n",
    "* Why does it (sort of) work?\n",
    "\n",
    "## See also:\n",
    "\n",
    "* [Gradient Descent Algorithm and Its Variants](https://towardsdatascience.com/gradient-descent-algorithm-and-its-variants-10f652806a3) by Imad Dabbura\n",
    "* [An overview of gradient descent optimization algorithms](https://ruder.io/optimizing-gradient-descent/) by Sebastian Ruder\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parabola\n",
    "\n",
    "Let's start with something easy, a [parabola](https://en.wikipedia.org/wiki/Parabola):\n",
    "\n",
    "$$y = x^2$$\n",
    "\n",
    "If we want to see how the function changes, we can look at its [derivative](https://en.wikipedia.org/wiki/Derivative), i.e. its slope at a given point. It is defined by:\n",
    "\n",
    "$$ \\frac{\\partial y}{\\partial x} = \\lim_{h \\to 0} \\frac{y(x + h) - y(x)} {h}$$\n",
    "\n",
    "For $y=x^2$ we can calculate it:\n",
    "\n",
    "$$\\frac{\\partial y}{\\partial x}\n",
    "= \\lim_{h \\to 0} \\frac{x^2 + 2 x h + h^2 - x^2}{h}\n",
    "= \\lim_{h \\to 0} \\left( 2 x + h \\right)  = 2x$$\n",
    "\n",
    "Now, let's see how does it look."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we approximate with 100 points, from -4 to 4\n",
    "X = np.linspace(-4, 4, num=100)\n",
    "Y = X**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax0, ax1) = plt.subplots(nrows=2, ncols=1,\n",
    "                               sharex=True, figsize=(7, 4))\n",
    "\n",
    "ax0.plot(X, Y)\n",
    "ax0.set(title='', xlabel='', ylabel='y')\n",
    "\n",
    "ax1.plot(X, 2 * X)\n",
    "ax1.set(title='', xlabel='x', ylabel='dy/dx')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical derivative in NumPy\n",
    "\n",
    "If you don't know how to differentiate a function, we can do it numerically by picking a step and dividing the changes of y by the changes of x."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dX = X[1:] - X[:-1]  # np.diff(X)\n",
    "dY = Y[1:] - Y[:-1]  # np.diff(Y)\n",
    "\n",
    "plt.plot(X[:-1], dY / dX)\n",
    "plt.title(\"dy/dx for h = {:.3f}\".format(dX[0]));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Symbolic derivative in PyTorch\n",
    "\n",
    "One of the main features of every neural network package is automatic differentiation. \n",
    "In PyTorch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = 10.  #@param {type:\"number\"}\n",
    "\n",
    "x = torch.tensor(x0, requires_grad=True)\n",
    "print(\"    x = {}\".format(x))\n",
    "\n",
    "y = x.pow(2)\n",
    "print(\"    y = {}\".format(y))\n",
    "\n",
    "# we calculate derivative of y with respect to everything\n",
    "y.backward()\n",
    "print(\"dy/dx = {}\".format(x.grad))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = 0.2  #@param {type:\"number\"}\n",
    "x0 = 4.   #@param {type:\"number\"}\n",
    "\n",
    "xs = [x0]\n",
    "x = torch.tensor(x0, requires_grad=True)\n",
    "\n",
    "for i in range(10):\n",
    "    y = x.pow(2)\n",
    "    y.backward()\n",
    "    x.data.add_(- lr * x.grad.data)\n",
    "    x.grad.data.zero_()\n",
    "    xs.append(x.item())\n",
    "\n",
    "# and plotting that\n",
    "    \n",
    "points_X = np.array(xs)\n",
    "points_Y = points_X**2\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 4))\n",
    "ax.plot(X, Y)\n",
    "ax.plot(points_X, points_Y, '-')\n",
    "ax.plot(points_X, points_Y, 'r.')\n",
    "ax.set(title='Gradient descent', xlabel='x', ylabel='y');\n",
    "    \n",
    "xs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Try other learning rates, e.g.:\n",
    "\n",
    "* 0.1\n",
    "* 0.5\n",
    "* 0.75\n",
    "* 1.\n",
    "* 1.5\n",
    "* -0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The exact solution\n",
    "\n",
    "In the case of $y = x^2$ it is possible to give an exact solution.\n",
    "\n",
    "We start with a point $x_0$ and with each step, we modify position \n",
    "\n",
    "$$ x_{t+1} = x_{t} - \\varepsilon \\left.\\frac{\\partial y}{\\partial x}\\right|_{x=x_t}\n",
    "= x_{t} - 2 \\varepsilon x_{t} = \\left( 1 - 2 \\varepsilon \\right) x_{t} = \\left( 1 - 2 \\varepsilon \\right)^t x_{0} $$\n",
    "\n",
    "That it is, is a [geometric sequence](https://en.wikipedia.org/wiki/Geometric_progression). They are found in the radioactive decay or contagious infections. Depending on the coefficient, they have a property to decay to zero, or explode. \n",
    "\n",
    "If, and only if, $-1 < (1 - 2 \\varepsilon) < 1$ then $x_{t}$ converges to zero.\n",
    "\n",
    "### Questions\n",
    "\n",
    "* For which learning rate ($\\varepsilon$) the position converges to zero?\n",
    "* If we had another function, e.g. $y = 10 x^2$, would this criterion be the same? (Why? Or why not?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slightly more complicated functions\n",
    "\n",
    "For $y = x^2$ we know that the lowest value is in 0, and we can easily solve it analytically. But... does this technique work for more complicated functions?\n",
    "\n",
    "Let's make things complicated, but only a tiny bit.\n",
    "Calculating derivatives is simple. Though, if you want to do it automatically, for example [with SymPy package](https://docs.sympy.org/latest/tutorial/calculus.html) or online ((https://www.wolframalpha.com/examples/mathematics/calculus-and-analysis/))[with Wolfram Alpha]. Polynomials are the simplest - you can give it a try write $(x+h)^n - x^n$ and see what happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return x - 4 * x**2 + 0.25 * x**4\n",
    "\n",
    "def df(x):\n",
    "    return 1 - 8 * x + x**3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.linspace(-4, 4, num=100)\n",
    "\n",
    "fig, (ax0, ax1) = plt.subplots(nrows=2, ncols=1,\n",
    "                               sharex=True, figsize=(7, 4))\n",
    "\n",
    "ax0.plot(X, f(X))\n",
    "ax0.set(title='', xlabel='', ylabel='y')\n",
    "\n",
    "ax1.plot(X, df(X))\n",
    "ax1.hlines(y=0, xmin=-4, xmax=4, linestyles='dashed')\n",
    "ax1.set(title='', xlabel='x', ylabel='dy/dx')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = 0.1  #@param {type:\"number\"}\n",
    "x0 = 4.   #@param {type:\"number\"}\n",
    "\n",
    "xs = [x0]\n",
    "x = torch.tensor(x0, requires_grad=True)\n",
    "\n",
    "for i in range(10):\n",
    "    y = f(x)\n",
    "    y.backward()\n",
    "    x.data.add_(- lr * x.grad.data)\n",
    "    x.grad.data.zero_()\n",
    "    xs.append(x.item())\n",
    "\n",
    "xs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points_X = np.array(xs)\n",
    "points_Y = f(points_X)\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 4))\n",
    "ax.plot(X, f(X))\n",
    "ax.plot(points_X, points_Y, '-')\n",
    "ax.plot(points_X, points_Y, 'r.')\n",
    "ax.set(title='Gradient descent', xlabel='x', ylabel='y');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient in 2d\n",
    "\n",
    "Gradients make sense for more dimensions. For mountains, gradient would be a vector directed along to the biggest slope.\n",
    "\n",
    "\n",
    "\n",
    "For example, let's have a function:\n",
    "\n",
    "$$y = g(x_1, x_2) = x_1^2 + \\sin(x_2)$$\n",
    "\n",
    "In this case, gradient is a vector. To calculate gradient we use [partial derivative](https://en.wikipedia.org/wiki/Partial_derivative).\n",
    "\n",
    "$$\\nabla g = \\left[ \\frac{\\partial g}{\\partial x_1}, \\frac{\\partial g}{\\partial x_2} \\right] = \\left[ 2 x_1, \\cos(x_2) \\right] $$\n",
    "\n",
    "\n",
    "Gradient symbol:\n",
    "\n",
    "$$\\nabla = \\left[ \\frac{\\partial }{\\partial x_1}, \\frac{\\partial }{\\partial x_2} \\right]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X0_ = np.linspace(-4, 4, num=100)\n",
    "X1_ = np.linspace(-4, 4, num=100)\n",
    "X0, X1 = np.meshgrid(X0_, X1_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# purely technically, so that we have the same code\n",
    "# for NumPy and PyTorch\n",
    "\n",
    "def sin(x):\n",
    "    if type(x) == torch.Tensor:\n",
    "        return x.sin()\n",
    "    else:\n",
    "        return np.sin(x)\n",
    "    \n",
    "def cos(x):\n",
    "    if type(x) == torch.Tensor:\n",
    "        return x.cos()\n",
    "    else:\n",
    "        return np.cos(x)\n",
    "\n",
    "# now the functions and their gradients\n",
    "# (calculated by hand)\n",
    "    \n",
    "def g(x0, x1):\n",
    "    return 0.25 * x0**2 + sin(x1)\n",
    "\n",
    "def dg_dx0(x0, x1):\n",
    "    return 0.5 * x0\n",
    "\n",
    "def dg_dx1(x0, x1):\n",
    "    return cos(x1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, let's draw a [contour plot](https://en.wikipedia.org/wiki/Contour_line), well known from topographic maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax0, ax1) = plt.subplots(nrows=1, ncols=2, figsize=(10, 5),\n",
    "                               sharex=True)\n",
    "\n",
    "cs = ax0.contour(X0, X1, g(X0, X1), cmap='coolwarm')\n",
    "ax0.clabel(cs, inline=1, fontsize=10)\n",
    "ax0.set_title(\"g\")\n",
    "ax0.set_xlabel(\"x0\")\n",
    "ax0.set_ylabel(\"x1\")\n",
    "\n",
    "X0_less = X0[::5, ::5]\n",
    "X1_less = X1[::5, ::5]\n",
    "ax1.set_title(r'$\\nabla g$')\n",
    "ax1.quiver(X0_less, X1_less,\n",
    "           dg_dx0(X0_less, X1_less), dg_dx1(X0_less, X1_less),\n",
    "           dg_dx0(X0_less, X1_less)**2 + dg_dx1(X0_less, X1_less)**2,\n",
    "           units='width', cmap='coolwarm')\n",
    "ax1.set_xlabel(\"x0\")\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Gradient descent in 2d\n",
    "\n",
    "Usually neural networks use more variables than one. For example - millions.\n",
    "But let's keep it simple and reduce to 2. Gradient works... and much alike water flowing the steepest descent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = 0.25  #@param {type:\"number\"}\n",
    "v = [-3.5, 1.0]\n",
    "\n",
    "xs = [v[0]]\n",
    "ys = [v[1]]\n",
    "v = torch.tensor(v, requires_grad=True)\n",
    "\n",
    "for i in range(20):\n",
    "    y = g(v[0], v[1])\n",
    "    y.backward()\n",
    "    v.data.add_(- lr * v.grad.data)\n",
    "    v.grad.data.zero_()\n",
    "    \n",
    "    xs.append(v[0].item())\n",
    "    ys.append(v[1].item())\n",
    "\n",
    "\n",
    "cs = plt.contour(X0, X1, g(X0, X1), cmap='coolwarm')\n",
    "plt.clabel(cs, inline=1, fontsize=10)\n",
    "plt.plot(xs, ys, '-')\n",
    "plt.plot(xs, ys, 'o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More funny stuff\n",
    "\n",
    "> TIL that changing random stuff until your program works is \"hacky\" and \"bad coding practice\" but if you do it fast enough it's \"#MachineLearning\" and pays 4x your current salary - [@smaine](https://twitter.com/smaine/status/994723834434502658)\n",
    "\n",
    "From [this tweet by @jebbery](https://twitter.com/jebbery/status/995491957559439360):\n",
    "\n",
    "![](https://pbs.twimg.com/media/DdCyqMXW4AAZJol?format=jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pytorch-macos-m1-gpu",
   "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.13"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "032715a8d60757555a2623145b534a689e8b5634c5398e8bbec1014c4a8ede12"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}