08_Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Research involves lots of parmeter tuing. When there are just a few paramters, we often tune them by hand using our intuition or run a grid search. But when the number of parameters is large or it is difficult to get any intuition, we need a systematic method for optimization.\n",
    "\n",
    "Optimization is in general considered as minimization or maximization of a certain objective function $f(x)$ where $x$ is a parameter vector. There are different cases:\n",
    "\n",
    "* If the mathematical form of the objective function $f(x)$ is known:\n",
    "\n",
    "    - put derivatieves $\\frac{\\partial f(x)}{\\partial x}=0$ and solve for $x$.\n",
    "    - check the signs of second order derivatives $\\frac{\\partial^2 f(x)}{\\partial x^2}$\n",
    "    \n",
    "        - if all positive, that is a minimum\n",
    "        - if all negative, that is a maximum\n",
    "        - if mixed, that is a saddle point\n",
    "\n",
    "\n",
    "* If analytic solution of $\\frac{\\partial f(x)}{\\partial x}=0$ is hard to derive:\n",
    "\n",
    "    - gradient descent/ascent\n",
    "    - Newton-Raphson method\n",
    "    - conjugate gradient method\n",
    "\n",
    "\n",
    "* If the derivatives of $f(x)$ is hard to derive:\n",
    "\n",
    "    - genetic/evolutionary algorithms\n",
    "    - sampling methods (next week)\n",
    "\n",
    "\n",
    "* If $f(x)$ needs to be optimized under constraints, such as $g(x)\\le 0$ or $h(x)=0$:\n",
    "\n",
    "    - penalty function\n",
    "    - Lagrange multiplyer method\n",
    "    - linear programming if $f(x)$ is linear\n",
    "    - quadratic programming if $f(x)$ is quadratic\n",
    "\n",
    "References:\n",
    "* Jan A. Snyman: Practical Mathematial Optimization. Springer, 2005.\n",
    "* SciPy Lecture Notes: 1.5.5 Optimization and fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Example\n",
    "For the sake of visualization, consider a function in 2D space $x=(x_1,x_2)$\n",
    "\n",
    "$$ f(x) = x_1^4 - \\frac{8}{3}x_1^3 - 6x_1^2 + x_2^4$$\n",
    "\n",
    "The gradient is\n",
    "\n",
    "$$ \\nabla f(x) = \\frac{\\partial f(x)}{\\partial x} = \n",
    "    \\left(\\begin{array}{c} 4x_1^3 - 8x_1^2 - 12x_1\\\\\n",
    "    4x_2^3\\end{array}\\right).$$\n",
    "\n",
    "By putting $\\nabla f(x)=0$, we have\n",
    "\n",
    "$$ x_1(x_1^2 - 2x_1 - 3) = x_1(x_1 + 1)(x_1 - 3) = 0 $$\n",
    "\n",
    "$$ x_2^3 = 0, $$\n",
    "\n",
    "so there are three points with zero gradient: $(-1,0)$, $(0,0)$, $(3,0)$.\n",
    "\n",
    "You can check the second-order derivative, or *Hessian*, to see if they are a minimum, a saddle point, or a maximum.\n",
    "\n",
    "$$ \\nabla^2 f(x) = \\frac{\\partial^2 f(x)}{\\partial x^2} = \n",
    "    \\left(\\begin{array}{cc} 12x_1^2 - 16x_1 - 12 & 0\\\\\n",
    "    0 & 12x_2^2\\end{array}\\right).$$\n",
    "\n",
    "As $\\frac{\\partial^2 f(x)}{\\partial x_1^2}$ is positive for $x_1=-1$ and $x_1=3$ and negative for $x_1=0$, $(-1,0)$ and $(3,0)$ are minima and $(0,0)$ is a saddle point.\n",
    "\n",
    "Let us visualize this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dips(x):\n",
    "    \"\"\"a function to minimize\"\"\"\n",
    "    f = x[0]**4 - 8/3*x[0]**3 - 6*x[0]**2 + x[1]**4\n",
    "    return(f)\n",
    "\n",
    "def dips_grad(x):\n",
    "    \"\"\"gradient of dips(x)\"\"\"\n",
    "    df1 = 4*x[0]**3 - 8*x[0]**2 - 12*x[0]\n",
    "    df2 = 4*x[1]**3\n",
    "    return(np.array([df1, df2]))\n",
    "\n",
    "def dips_hess(x):\n",
    "    \"\"\"hessian of dips(x)\"\"\"\n",
    "    df11 = 12*x[0]**2 - 16*x[0] - 12\n",
    "    df12 = 0\n",
    "    df22 = 12*x[1]**2\n",
    "    return(np.array([[df11, df12], [df12, df22]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x2 = np.meshgrid(np.linspace(-5, 5, 20), np.linspace(-5, 5, 20))\n",
    "fx = dips([x1, x2])\n",
    "# 3D plot\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.plot_surface(x1, x2, fx, cmap='viridis')\n",
    "ax.set_xlabel('x1')\n",
    "ax.set_ylabel('x2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.contour(x1, x2, fx)\n",
    "dfx = dips_grad([x1, x2])\n",
    "plt.quiver(x1, x2, dfx[0], dfx[1])\n",
    "plt.axis('square')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Descent/Ascent\n",
    "*Gradient descent/ascent* is the most basic method of min/maximization of a function using its gradient.\n",
    "\n",
    "From an initial state $x_0$ and a coefficient $\\eta>0$, repeat\n",
    "\n",
    "$$ x_{i+1} = x_i - \\eta\\nabla f(x_i) $$\n",
    "\n",
    "for minimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_descent(f, df, x0, eta=0.01, eps=1e-6, imax=1000):\n",
    "    \"\"\"Gradient descent\"\"\"\n",
    "    xh = np.zeros((imax+1, len(np.ravel([x0]))))  # history\n",
    "    xh[0] = x0\n",
    "    f0 = f(x0)  # initialtization\n",
    "    for i in range(imax):\n",
    "        x1 = x0 - eta*df(x0)\n",
    "        f1 = f(x1)\n",
    "        # print(x1, f1)\n",
    "        xh[i+1] = x1\n",
    "        if(f1 <= f0 and f1 > f0 - eps):  # small decrease\n",
    "            return(x1, f1, xh[:i+2])\n",
    "        x0 = x1\n",
    "        f0 = f1\n",
    "    print(\"Failed to converge in \", imax, \" iterations.\")\n",
    "    return(x1, f1, xh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, fmin, xhist = grad_descent(dips, dips_grad, [1,2], 0.02)\n",
    "print(xmin, fmin)\n",
    "plt.contour(x1, x2, fx)\n",
    "plt.plot(xhist[:,0], xhist[:,1], '.-')\n",
    "#plt.axis([1, 4, -1, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton-Raphson Method\n",
    "A problem with the gradient descence/ascent is the choice of the coefficient $\\eta$. If the second-order derivative, called the *Hessian*, \n",
    "\n",
    "$$\\nabla^2f(x)=\\frac{\\partial^2 f}{\\partial x^2}$$ \n",
    "\n",
    "is available, we can use the Newton method to find the solution for \n",
    "\n",
    "$$\\nabla f(x)=\\frac{\\partial f}{\\partial x}=0$$\n",
    "\n",
    "by repeating\n",
    "\n",
    "$$ x_{i+1} = x_i - \\nabla^2f(x_i)^{-1} \\nabla f(x_i). $$\n",
    "\n",
    "This is called Newton-Raphson method. It works efficiently when the Hessian is positive definite ($f(x)$ is like a parabolla), but can be unstable when the Hessian has a negative eigenvalue (near the saddle point)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton_raphson(f, df, d2f, x0, eps=1e-6, imax=1000):\n",
    "    \"\"\"Newton-Raphson method\"\"\"\n",
    "    xh = np.zeros((imax+1, len(np.ravel([x0]))))  # history\n",
    "    xh[0] = x0\n",
    "    f0 = f(x0)  # initialtization\n",
    "    for i in range(imax):\n",
    "        x1 = x0 - np.linalg.inv(d2f(x0)) @ df(x0)\n",
    "        f1 = f(x1)\n",
    "        #print(x1, f1)\n",
    "        xh[i+1] = x1\n",
    "        if( f1 <= f0 and f1 > f0 - eps):  # decreasing little\n",
    "            return(x1, f1, xh[:i+2])\n",
    "        x0 = x1\n",
    "        f0 = f1\n",
    "    print(\"Failed to converge in \", imax, \" iterations.\")\n",
    "    return(x1, f1, xh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, fmin, xhist = newton_raphson(dips, dips_grad, dips_hess, [4,2])\n",
    "print(xmin, fmin)\n",
    "plt.contour(x1, x2, fx)\n",
    "plt.plot(xhist[:,0], xhist[:,1], '.-')\n",
    "#plt.axis([1, 4, -1, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## scipy.optimize\n",
    "To address those issues, advanced optimization algorithms have been developed and implemented in `scipy.optimize` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The default method for unconstrained minimization is 'BFGS' (Broyden-Fletcher-Goldfarb-Shanno) method, a variant of gradient descent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = minimize(dips, [-1,2], jac=dips_grad, options={'disp': True})\n",
    "print( result.x, result.fun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the gradient function is not specified, it is estimated by finite difference method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = minimize(dips, [2,2], options={'disp': True})\n",
    "print( result.x, result.fun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 'Newton-CG' (Newton-Conjugate-Gradient) is a variant of Newton-Raphson method using linear search in a *conjugate* direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = minimize(dips, [2,2], method='Newton-CG', \n",
    "          jac=dips_grad, hess=dips_hess, options={'disp': True})\n",
    "print( result.x, result.fun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 'Nelder-Mead' is a *simplex* method that uses a set of $n+1$ points to estimate the gradient and select a new point by flipping the simplex.\n",
    "    - note that it is totally different from the *simplex* method for linear programming."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = minimize(dips, [2,2], method='Nelder-Mead', options={'disp': True})\n",
    "print( result.x, result.fun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constrained Optimization\n",
    "Often we want to minimize/maximize $f(x)$ under constraints on $x$, e.g.\n",
    "* inequality constraints\n",
    "\n",
    "$$g_j(x)\\le 0 \\quad (j=1,...,m)$$\n",
    "\n",
    "* equality constraints\n",
    "\n",
    "$$h_j(x)=0 \\quad (j=1,...,r)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Penalty function\n",
    "Define a function with penalty terms:\n",
    "\n",
    "$$ P(x,\\rho) = f(x) + \\sum_j \\beta_j(x) g_j(x)^2 + \\sum_j \\rho h_j(x)^2 $$\n",
    "\n",
    "$$ \\beta_j(x) = \\left\\{ \\begin{array}{ccc}\n",
    "    0 & \\mbox{if} & g_j(x)\\le 0 \\\\\n",
    "    \\rho & \\mbox{if} & g_j(x)>0 \\end{array}\\right.$$\n",
    "\n",
    "and increase $\\rho$ to a large value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lagrange multiplyer method\n",
    "For minimization of $f(x)$ with equality constraints\n",
    "$h_j(x)=0$, $(j=1,...,r)$, define a *Lagrangian function*\n",
    "\n",
    "$$ L(x,\\lambda) = f(x) + \\sum_j \\lambda_j h_j(x). $$\n",
    "\n",
    "The necessary condition for a minimum is:\n",
    "\n",
    "$$ \\frac{\\partial L(x,\\lambda)}{\\partial x_i} = 0 \\qquad (i=1,...,n) $$\n",
    "\n",
    "$$ \\frac{\\partial L(x,\\lambda)}{\\partial \\lambda_j} = 0 \\qquad (j=1,...,r) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy implements `SLSQP` (Sequential Least SQuares Programming) method. Constraints are defined in a sequence of dictionaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# h(x) = - x[0] + x[1] - 0.6 = 0\n",
    "def h(x):\n",
    "    return -x[0] + x[1] - 0.6\n",
    "def h_grad(x):\n",
    "    return np.array([-1, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With equality constraint $h(x)=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cons = ({'type':'eq', 'fun':h, 'jac':h_grad})\n",
    "result = minimize(dips, [1,-3], jac=dips_grad,\n",
    "            method='SLSQP', constraints=cons, options={'disp': True})\n",
    "print( result.x, result.fun)\n",
    "plt.contour(x1, x2, fx)\n",
    "plt.plot([-4,4], [-3.4,4.4])  # h(x) = 0\n",
    "plt.plot(result.x[0], result.x[1], 'o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With inequality constraint $h(x)>0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cons = ({'type': 'ineq', 'fun': h, 'jac':h_grad})\n",
    "result = minimize(dips, [1,-3], jac=dips_grad,\n",
    "            method='SLSQP', constraints=cons, options={'disp': True})\n",
    "print( result.x, result.fun)\n",
    "plt.contour(x1, x2, fx)\n",
    "plt.plot([-4,4], [-3.4,4.4])  # h(x) = 0\n",
    "plt.plot(result.x[0], result.x[1], 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Genetic/Evolutionaly Algorithms\n",
    "For objective functions with many local minima/maxima, stochastic search methods are preferred. They are called *genetic algorithm (GA)* or *evolutionay algorithm (EA)*, from an analogy with mutation and selection in genetic evolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evol_min(f, x0, sigma=0.1, imax=100):\n",
    "    \"\"\"simple evolutionary algorithm\n",
    "    f: function to be minimized\n",
    "    x0: initial population (p*n)\n",
    "    sigma: mutation size\"\"\"\n",
    "    p, n = x0.shape  # population, dimension\n",
    "    x1 = np.zeros((p, n))\n",
    "    xh = np.zeros((imax, n))  # history\n",
    "    for i in range(imax):\n",
    "        f0 = f(x0.T)  # evaluate the current population\n",
    "        fmin = min(f0)\n",
    "        xmin = x0[np.argmin(f0)]\n",
    "        #print(xmin, fmin)\n",
    "        xh[i] = xmin  # record the best one\n",
    "        # roulette selection\n",
    "        fitness = max(f0) - f0  # how much better than the worst\n",
    "        prob = fitness/sum(fitness)  # selection probability\n",
    "        #print(prob)\n",
    "        for j in range(p):  # pick a parent for j-th baby\n",
    "            parent = np.searchsorted(np.cumsum(prob), np.random.random())\n",
    "            x1[j] = x0[parent] + sigma*np.random.randn(n)\n",
    "        x0 = x1\n",
    "    return(xmin, fmin, xh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.random.rand(20, 2)*10 - 5\n",
    "xmin, fmin, xhist = evol_min(dips, x0, 0.1)\n",
    "print(xmin, fmin)\n",
    "plt.contour(x1, x2, fx)\n",
    "plt.plot(xhist[:,0], xhist[:,1], '.-')\n",
    "plt.plot(x0[:,0], x0[:,1], '*')\n",
    "plt.plot(xhist[-1,0], xhist[-1,1], 'o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "For more advanced genetic/evolutionary algorithms, you can use `deap` package: https://github.com/DEAP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}