chapter7-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "af61f11b",
   "metadata": {},
   "source": [
    "# 7.1 Golden-section search"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41ea9f97",
   "metadata": {},
   "source": [
    "Golden-section search (GSS) is a derivative-free univariate optimization method, which works in a similar fashion to the bisection method discussed in Chapter 2 for solving nonlinear equations. GSS starts with an initial interval where the minimizer (maximizer) is located and continuously reduces the range of the interval until it is small enough. GSS works for a particular type of cost function and we define the function below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a43490d",
   "metadata": {},
   "source": [
    "\\begin{definition}[Unimodal function]\\label{def:unimodal}\n",
    "A function $f(x)$ defined on $[a,b]$ is unimodal if there exists a point $m\\in [a,b]$ such that $f(x)$ is monotonically decreasing (resp., increasing) for $x\\le m$ and monotonically increasing (resp., decreasing) for $x\\ge m$.\n",
    "\\end{definition}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e39034a6",
   "metadata": {},
   "source": [
    "So a unimodal function has a unique minimum (or maximum) on the given interval. GSS applies to unimodal functions. To illustrate the algorithm, we use the following simple degree 4 polynomial:\n",
    "\\begin{equation}\\label{eq:golden_exa1}\n",
    "f(x) = 3x^4+2x^2-x-1\n",
    "\\end{equation}\n",
    "of which the plot is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40848379",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fe53657",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 3*x**4+2*x**2-x-1\n",
    "\n",
    "x = np.linspace(-1, 1, 1000)\n",
    "plt.plot(x, f(x))\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e814a092",
   "metadata": {},
   "source": [
    "It clear shows $f(x)$ is unimodal on the interval, e.g., $[-0.5, 0.5]$, and a unique minimum is located in this interval.\n",
    "\n",
    "We now consider how GSS works for minimization problems. Maximization problems can be converted to minimization problems by negating the cost functions. The following figure shows the plot of a cost function $f(x)$ that is unimodal on the interval $[a,b]$ with a minimizer at $x_0$.  Our goal now is to find another two points $x_1, x_2\\in [a,b]$, $x_1<x_2$, such that a new smaller interval can be constructed based on the four points. If $f(x_1)\\le f(x_2)$, then we know the minimizer must be between $a$ and $x_2$, and if $f(x_1) > f(x_2)$, then $x_0\\in [x_1, b]$. Both $[a, x_2]$ and $[x_1, b]$ are shorter compared to $[a,b]$. To guarantee an efficient reduction of the interval, we will ensure the following two conditions:\n",
    "1. The distance between $x_1$ and $a$ is the same as that between $x_2$ and $b$, so that the new intervals in both cases have the same length. Let $l_1=x_1-a$, $l_2=x_2-x_1$, and $l_3=b-x_2$. Then we have \n",
    "    \\begin{equation}\\label{eq:golden1}\n",
    "    l_1=l_3.\n",
    "    \\end{equation}\n",
    "2. If the interval $[a, x_2]$ becomes the new interval, the spacing of the three points $a, x_1, x_2$ is the same as that of $a, x_1, b$. Mathematically, it means\n",
    "   \\begin{equation}\\label{eq:golden2}\n",
    "   \\frac{l_1}{l_2+l_3} = \\frac{l_2}{l_1}.\n",
    "   \\end{equation}\n",
    "   Note that this also ensures if the interval $[x_1, b]$ becomes the new interval, the spacing of the three points $x_1, x_2, b$ is the same as that of $a, x_1, b$, in which case, we have\n",
    "   \\begin{equation}\\label{eq:golden3}\n",
    "   \\frac{l_1}{l_2+l_3} = \\frac{l_2}{l_3}.\n",
    "   \\end{equation}\n",
    "\n",
    "Rewriting $l_2$ and $l_1$ on the right hand side of Eqn. \\eqref{eq:golden2} by $(l_2+l_3)-l_3$ and $(l_2+l_3)-l_2$ respectively, and reorganizing, an equation of $\\frac{l_2+l_3}{l_1}$ can be obtained:\n",
    "\n",
    "\\begin{align}\\label{eq:golden_eqn}\n",
    "&\\frac{l_1}{l_2+l_3} = \\frac{(l_2+l_3)-l_3}{(l_2+l_3)-l_2} \\nonumber \\\\\n",
    "&\\Rightarrow (l_2+l_3)^2-l_3(l_2+l_3) = l_1(l_2+l_3)-l_1l_2 \\nonumber \\\\\n",
    "&\\Rightarrow (l_2+l_3)^2-2l_1(l_2+l_3) + l_1l_2 = 0 \\nonumber \\\\\n",
    "&\\Rightarrow (l_2+l_3)^2-2l_1(l_2+l_3) + l_1[(l_2+l_3)-l_3] = 0 \\nonumber \\\\\n",
    "&\\Rightarrow (l_2+l_3)^2-l_1(l_2+l_3) -l_1l_3 = 0 \\nonumber \\\\\n",
    "&\\Rightarrow\\left(\\frac{l_2+l_3}{l_1}\\right)^2 - \\frac{l_2+l_3}{l_1} - 1 = 0.\n",
    "\\end{align}\n",
    "So \n",
    "\\begin{equation}\\label{eq:golden_ratio}\n",
    "\\frac{l_2+l_3}{l_1} = r = \\frac{\\sqrt{5}+1}{2}\n",
    "\\end{equation}\n",
    "which is the **golden ratio**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a97e23ee",
   "metadata": {},
   "source": [
    "<img src=\"images/golden.png\" width=500 height=500 />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "358524e2",
   "metadata": {},
   "source": [
    "Therefore, once the interval $[a,b]$ is identified, we can find $x_1$ by \n",
    "\\begin{equation*}\n",
    "x_1=a+\\frac{l_1}{l_1+(l_2+l_3)}(b-a)=a+\\frac{2}{\\sqrt{5}+3}(b-a)=a+\\frac{3-\\sqrt{5}}{2}(b-a)\n",
    "\\end{equation*}\n",
    "and $x_2$ by \n",
    "\\begin{align*}\n",
    "x_2&=b-\\frac{l_3}{l_1+(l_2+l_3)}(b-a) = b-\\frac{l_1}{l_1+(l_2+l_3)}(b-a) \\\\\n",
    "    &=b-\\frac{2}{\\sqrt{5}+3}(b-a)=a+(b-a)+\\frac{\\sqrt{5}-3}{2}(b-a) = a+\\frac{\\sqrt{5}-1}{2}(b-a).\n",
    "\\end{align*}\n",
    "Then this process continues until a certain criterion is met, such as a specified number of steps is reached, or the length of the final interval is smaller than a prescribed threshold. The final estimated minimizer is either $a$ or $b$, depending on the values of $f(a)$ and $f(b)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2d8a07e",
   "metadata": {},
   "source": [
    "The following Python function **golden_section** implements golden-section search. The algorithm stops when the size of the interval is smaller than a user-provided tolerance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcc9c800",
   "metadata": {},
   "outputs": [],
   "source": [
    "def golden_section(f, a, b, tol):\n",
    "    \"\"\"\n",
    "    Perform golden-section search\n",
    "    input:\n",
    "    f: function, the function to be minimized\n",
    "    a: float, the initial left end\n",
    "    b: float, the initial right end\n",
    "    tol: float, the tolerance for the stopping criterion\n",
    "    \"\"\"\n",
    "    g = (np.sqrt(5)-1)/2\n",
    "    while b-a > tol:\n",
    "        if f(a+(1-g)*(b-a)) < f(a+g*(b-a)):\n",
    "            b = a + g*(b-a)\n",
    "        else:\n",
    "            a += (1-g)*(b-a)\n",
    "    if f(a) < f(b):\n",
    "        return a\n",
    "    else:\n",
    "        return b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1770e36",
   "metadata": {},
   "source": [
    "We now apply *golden_section* to locate the minimizer of the cost function in Eqn. \\eqref{eq:golden_exa1}."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2960f41f",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0.0\n",
    "b = 0.5\n",
    "tol = 1e-9\n",
    "p = golden_section(f, a, b, tol)\n",
    "print('The function f has a minimum of ', f(p), ' at ', p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4856772c",
   "metadata": {},
   "source": [
    "The minimum is achieved at around $0.22$, which is in accordance with what we see from the plot of $f$. We can animate the process of how the solution is found by golden-section search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a8ebf9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML\n",
    "\n",
    "# golden-section search for the animation\n",
    "def golden_section_nstep(f, a, b, nstep):\n",
    "    \"\"\"\n",
    "    Perform golden-section search for nstep\n",
    "    input:\n",
    "    f: function, the function to be minimized\n",
    "    a: float, the initial left end\n",
    "    b: float, the initial right end\n",
    "    nstep: int, the number of steps for which GSS is performed\n",
    "    \"\"\"\n",
    "    g = (np.sqrt(5)-1)/2\n",
    "    for i in range(nstep):\n",
    "        if f(a+(1-g)*(b-a)) < f(a+g*(b-a)):\n",
    "            b = a + g*(b-a)\n",
    "        else:\n",
    "            a += (1-g)*(b-a)\n",
    "    if f(a) < f(b):\n",
    "        return a, a, b\n",
    "    else:\n",
    "        return b, a, b\n",
    "\n",
    "# Animation for the golden-section search method\n",
    "a = -1\n",
    "b = 1\n",
    "true_sol = 0.21864333330456398\n",
    "x_plot = np.linspace(-1, 1, 100)\n",
    "y_plot = f(x_plot)\n",
    "fig, ax = plt.subplots(figsize=(12,5))\n",
    "def animate(i):\n",
    "    ax.cla()\n",
    "    ax.set_xlim(-1, 1)\n",
    "    ax.plot(x_plot, y_plot, label='f')\n",
    "    x, aa, bb = golden_section_nstep(f, a, b, i)\n",
    "    ax.plot(x, f(x), 'o', ms=6, label='Golden-section search')\n",
    "    ax.plot(aa, f(aa), 'x', ms = 6, label='bracket-left end')\n",
    "    ax.plot(bb, f(bb), 'x', ms = 6, label='bracket-right end')\n",
    "    ax.plot(true_sol, f(true_sol), 's', color='green', ms=6, label='True solution')\n",
    "    ax.plot()\n",
    "    ax.axhline()\n",
    "    ax.axvline()\n",
    "    ax.set_title('Golden-section search: iteration step '+str(i), size=16)\n",
    "    ax.set_xlabel('$x$', size=16)\n",
    "    ax.set_ylabel('$y$', size=16)\n",
    "    ax.legend(loc='upper right')\n",
    "plt.close()\n",
    "anim = FuncAnimation(fig, animate, frames=20, interval=2000)\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e195438",
   "metadata": {},
   "source": [
    "The animation shows the plot of the cost function, the true solution, the bracket of each step of golden-section search, and the estimated solution at each step. It can be seen that after around 10 steps, the estimated solution is already very close to the true minimizer."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8ed06c5",
   "metadata": {},
   "source": [
    "\\begin{example}\\label{example:golden_sin}\n",
    "Consider the negative sine function $f(x)=-\\sin{x}$ on the interval $[0, \\pi]$. Apply golden-section search to find the minimum.\n",
    "\\end{example}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e153e17e",
   "metadata": {},
   "source": [
    "We know the minimum is $-1$, achieved at $\\frac{\\pi}{2}$. Now we use golden-section search to obtain the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d294ac1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return -np.sin(x)\n",
    "\n",
    "a = 0.0\n",
    "b = np.pi\n",
    "tol = 1e-5\n",
    "p = golden_section(f, a, b, tol)\n",
    "print('The function f has a minimum of ', f(p), ' at ', p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42615066",
   "metadata": {},
   "source": [
    "The same results are obtained as expected."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2014459",
   "metadata": {},
   "source": [
    "\\begin{exercise}\\label{ex:7.1-1}\n",
    "Use golden-section search to find the maximum of the following function:\n",
    "\\begin{equation*}\n",
    "f(x) = -x^2+\\cos{(2x)}.\n",
    "\\end{equation*}\n",
    "\\end{exercise}"
   ]
  }
 ],
 "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.7"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}