5_strategies/CMAES.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/logo.png\" width=\"20%\" align=\"right\" style=\"margin:0px 20px\">\n",
    "\n",
    "\n",
    "# Evolutionary Algorithms\n",
    "\n",
    "## Covariance Matrix Adaptation Evolutionary Strategy\n",
    "\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-sa/4.0/\"><img alt=\"Creative Commons License\" align=\"left\" src=\"https://i.creativecommons.org/l/by-sa/4.0/80x15.png\" /></a>&nbsp;| Dennis G. Wilson | <a href=\"https://d9w.github.io/evolution/\">https://d9w.github.io/evolution/</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Covariance Matrix Adaptation Evolutionary Strategy, or CMA-ES [1, 2], is one of the most well-known evolutionary algorithms in general and is a state-of-the-art algorithm for continuous optimization. The strength of this method is that it adapts the distribution it uses to generate the next population based on the current distribution of individuals. Until now, we were limited to a Normal distribution with a static $\\sigma$. The adaptive distribution of CMA-ES means it will cross search spaces faster and narrow in more exactly on optimal points.\n",
    "\n",
    "[1] Hansen, Nikolaus, and Andreas Ostermeier. \"Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation.\" Proceedings of IEEE international conference on evolutionary computation. IEEE, 1996.\n",
    "\n",
    "[2] Hansen, Nikolaus, and Andreas Ostermeier. \"Completely derandomized self-adaptation in evolution strategies.\" Evolutionary computation 9.2 (2001): 159-195."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Specifically, the things that CMA-ES improves over the previous Evolutionary Strategies we've seen is that it:\n",
    "+ combines information from multiple individuals (as in $(\\mu, \\lambda)$ and Canonical ES)\n",
    "+ combines information from multiple generations\n",
    "+ transforms the distribution of the new population to match the search space\n",
    "+ adapts the step size to prevent premature convergence\n",
    "\n",
    "One improvement in CMA-ES which we've already seen is the use of **fitness rank** instead of direct fitness for the update. However, instead of directly updating a centroid, CMA-ES updates a covariance matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "CMA-ES uses two principles to achieve this: maximum likelihood estimation and step-size control. We'll start with maximum likelihood estimation, which consists of increasing the probability of successful points. CMA-ES will update the population center by taking the weighted average of high-fitness individuals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 1</h3>\n",
    "    \n",
    "CMA-ES has been widely used in many applications. Discuss one of the following application papers in your group. What is the application? How is the problem encoded in a way that CMA-ES can solve it? Is CMA-ES compared to other methods and if so, how does it do?\n",
    "    \n",
    "+ [Gagné, C., Beaulieu, J., Parizeau, M., & Thibault, S. (2008). Human-competitive lens system design with evolution strategies. Applied Soft Computing, 8(4), 1439-1452.](http://www.genetic-programming.org/hc2007/06-Gagne/Gagne-RT-LVSN-2007-01.pdf)\n",
    "+ [Bayer, P., & Finkel, M. (2007). Optimization of concentration control by evolution strategies: Formulation, application, and assessment of remedial solutions. Water resources research, 43(2).](https://www.academia.edu/download/42649750/Optimization_of_concentration_control_by20160213-31306-8zy30l.pdf)\n",
    "+ [Kämpf, J. H., & Robinson, D. (2009). A hybrid CMA-ES and HDE optimisation algorithm with application to solar energy potential. Applied Soft Computing, 9(2), 738-745.](https://www.sciencedirect.com/science/article/pii/S1568494608001336?casa_token=w27xSjQ0zY8AAAAA:6XzGZ73KLEsXGbj5m7YE4LgZ1zCm5sQhQP8I2Zx6XH-HPP5PmfF5VO3d1pu7fXSn2jh97cN9AdU)\n",
    "+ [Maki, A., Sakamoto, N., Akimoto, Y., Nishikawa, H., & Umeda, N. (2020). Application of optimal control theory based on the evolution strategy (CMA-ES) to automatic berthing. Journal of Marine Science and Technology, 25(1), 221-233.](https://link.springer.com/article/10.1007/s00773-019-00642-3)\n",
    "+ [Loshchilov, I., & Hutter, F. (2016). CMA-ES for hyperparameter optimization of deep neural networks. arXiv preprint arXiv:1604.07269.](https://openreview.net/pdf/xnrA4qzmPu1m7RyVi38Z.pdf)\n",
    "+ [Fukagata, K., Kern, S., Chatelain, P., Koumoutsakos, P., & Kasagi, N. (2008). Evolutionary optimization of an anisotropic compliant surface for turbulent friction drag reduction. Journal of Turbulence, (9), N35.](http://www.cse-lab.ethz.ch/wp-content/papercite-data/pdf/fukagata2008a.pdf)\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import ticker, cm\n",
    "import matplotlib.colors as colors\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.animation as animation\n",
    "from IPython import display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def himmelblau(x, y):\n",
    "    return (x**2 + y - 11)**2 + (x + y**2 - 7)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def rosenbrock(x, y, a=1, B=100):\n",
    "      return (a-x)**2 + B*((y-x**2))**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def rastrigin(X, Y):\n",
    "    Z = (X**2 - 10 * np.cos(2 * np.pi * X)) + \\\n",
    "      (Y**2 - 10 * np.cos(2 * np.pi * Y)) + 20\n",
    "    return Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def rastrigin_shifted(X, Y, a=0.3, b=1.2, c=0.6, d=2.9):\n",
    "    X = a*X + b\n",
    "    Y = c*Y + d\n",
    "    Z = (X**2 - 10 * np.cos(2 * np.pi * X)) + \\\n",
    "      (Y**2 - 10 * np.cos(2 * np.pi * Y)) + 20\n",
    "    return Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "optf = himmelblau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "X = np.arange(-5, 5, 0.1)\n",
    "Y = np.arange(-5, 5, 0.1)\n",
    "X, Y = np.meshgrid(X, Y)\n",
    "Z = optf(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "norm = colors.Normalize(vmin=Z.min(), vmax=Z.max())\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral', alpha=0.4, norm=norm)\n",
    "fig.colorbar(cs)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Covariance Matrix Update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "First, we select the top $\\mu$ individuals based on their fitness values. We'll then weight their importance by their fitness rank, using a static weight vector $w$. This weight vector $w$ is a hyperparameter in CMA-ES, but there is a standard normalized logarithmic weight scale used by most. Note that the weight is not their fitness values but rather just the rank. This has been demonstrated to aid in search in ill-formed search spaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "lam = 70\n",
    "x = np.random.rand(2) * 5 - 2.5\n",
    "N = np.random.normal(size=(lam,2)) * 0.5\n",
    "x_t = x + N\n",
    "fits = optf(x_t[:, 0], x_t[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mu = 10\n",
    "sorted_ids = np.argsort(fits)\n",
    "w = np.log(mu + 1/2) - np.log(np.arange(1, mu+1))\n",
    "w /= np.sum(w)\n",
    "weighted = N[sorted_ids[:mu]] * w.reshape(mu, 1)\n",
    "s = np.mean(weighted, axis=0)\n",
    "s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This $s$ direction vector is similar to the gradient estimation from the last exercise. However, in CMA-ES, only the top $\\mu$ individuals are used to calculate $s$, and instead of directly using fitness values, $s$ is calculated using weight based on fitness rank. We can compare the two vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = (fits - np.mean(fits)) / np.std(fits)\n",
    "G = -np.dot(A, N) / lam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "norm = colors.Normalize(vmin=Z.min(), vmax=Z.max())\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral', alpha=0.4, norm=norm)\n",
    "plt.scatter(x_t[:, 0], x_t[:, 1], c=fits, cmap='Spectral', norm=norm)\n",
    "plt.quiver(x[0], x[1], G[0], G[1], color='b', scale=4, label='fitness')\n",
    "plt.quiver(x[0], x[1], s[0], s[1], color='g', scale=1, label='rank')\n",
    "plt.scatter(x[0], x[1], c='k')\n",
    "fig.colorbar(cs)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "While we're only using the fitness information for a ranked weighting, this direction vector $s$ can be seen as a gradient approximation. As before, we can use this gradient information to update the center of our distribution. This is the first step of maximum likelihood estimation, aligning the distribution center with the weighted average of the best solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = x + s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "A fundamental idea behind CMA-ES is to estimate the covariance between the selected samples, and to use this covariance matrix to shape the next distribution. A covariance matrix measures the joint variability between the different dimensions of a vector. In CMAES, we measure the covariance of the gradient estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We start CMA-ES with the identity matrix as the initial covariance matrix, assuming each variable is independent and has a standard deviation of 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "C = np.eye(2)\n",
    "C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "At each step, we update $C = C + ss^T$. This will use the weighted sample average as an estimate of the covariance. In CMA-ES, it is referred to as the rank one covariance update. For problems of higher dimension, there is an additional matrix update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "C = C + np.outer(s, (s).T)\n",
    "C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In practice, we will use a time constant to update the covariance matrix at different speeds.\n",
    "\n",
    "$C = (1 - c_i)C + c_i s s^T$\n",
    "\n",
    "where $c_i$ is the rank one covariance matrix update learning rate. $c_i = 2/n^2$ where $n$ is the number of dimensions is a default parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The eigen decomposition of the covariance matrix is then used to transform the Normal distribution to match the search space. This is the second part of maximum likelihood estimation: by pulling samples from a distribution fit to performant parts of the search space, there is a higher chance of sampling good individuals. The mean of the distribution $x$ and the shape given by $C$ together maximize this likelihood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "N = np.random.normal(size=(lam, 2))\n",
    "x_N = x + N\n",
    "eigenvalues, eigenvectors = np.linalg.eig(C)\n",
    "x_C = x + np.array([np.dot(eigenvectors, np.sqrt(eigenvalues) * N[i, :]) for i in range(lam)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "norm = colors.Normalize(vmin=Z.min(), vmax=Z.max())\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral', alpha=0.4, norm=norm)\n",
    "plt.scatter(x_N[:, 0], x_N[:, 1], c='k', alpha=0.3)\n",
    "plt.scatter(x_C[:, 0], x_C[:, 1], c='b')\n",
    "plt.scatter(x[0], x[1], c='r')\n",
    "fig.colorbar(cs)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 2</h3>\n",
    "    Simulate multiple covariance matrix updates from the same starting point, but using x_C (the transformed distribution) as the new population. Are the gradient directions sampled similar? How does the distribution change as the covariance matrix is updated multiple times?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Step size control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In CMA-ES, the standard deviation $\\sigma$ of the population distribution also updates automatically. A simple rule used before CMA-ES which demonstrates this concept is the one-fifth rule: \n",
    "\n",
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/cmaes_onefifth.png\" width=\"40%\" height=\"auto\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The rule used in CMA-ES is a bit more complex; it estimates the recent evolutionary path and uses this estimate to \"slow down\" or \"speed up\" search.\n",
    "\n",
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/cmaes_pathlength.png\" width=\"40%\" height=\"auto\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/cmaes_pathlength2.png\" width=\"40%\" height=\"auto\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "An effect of these two different updates, which are often referred to as two \"evolutionary paths\", is that CMA-ES will continue to search even after finding a local optima. $\\sigma$ updates faster than $C$, narrowing in on selected areas or expanding on others, while the shape of the sampled distribution is refined over time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## CMAES in practice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/cmaes.png\" width=\"40%\" height=\"auto\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To see CMA-ES in practice, we'll use the pycma package (https://github.com/CMA-ES/pycma). If you haven't installed it, do so here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "`pycma` is actively maintained by the CMA-ES author and has many more algorithm tricks and features than we've discussed here. [4] provides a good review of different CMA-ES modifications.\n",
    "\n",
    "[4] Hansen, Nikolaus. \"The CMA evolution strategy: a comparing review.\" Towards a new evolutionary computation. Springer, Berlin, Heidelberg, 2006. 75-102."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import cma\n",
    "es = cma.CMAEvolutionStrategy(2 * [0], 0.1, {'popsize': 20, 'verb_disp': 1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "solutions = np.array(es.ask())\n",
    "es.tell(solutions, [optf(x[0], x[1]) for x in solutions])\n",
    "es.disp_annotation()\n",
    "es.disp();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6))\n",
    "ax = plt.subplot()\n",
    "norm = colors.Normalize(vmin=Z.min(), vmax=Z.max())\n",
    "cs = plt.contour(X, Y, Z, levels=100, cmap='Spectral', alpha=0.4, norm=norm)\n",
    "plt.scatter(solutions[:, 0], solutions[:, 1])\n",
    "fig.colorbar(cs)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def generation(i, es, ax):\n",
    "    solutions = np.array(es.ask())\n",
    "    es.tell(solutions, [optf(x[0], x[1]) for x in solutions])\n",
    "    \n",
    "    ax.clear()\n",
    "    ax.contour(X, Y, Z, levels=100, cmap='Spectral', alpha=0.4, norm=norm)\n",
    "    ax.scatter(solutions[:, 0], solutions[:, 1])\n",
    "    ax.set_ylim(-5, 5)\n",
    "    ax.set_xlim(-5, 5)\n",
    "    ax.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "es = cma.CMAEvolutionStrategy(2 * [0], 0.1, {'popsize': 20})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "ax = plt.subplot()\n",
    "fig.tight_layout()\n",
    "frames = 20\n",
    "animator = animation.FuncAnimation(fig, generation, fargs=(es, ax), frames=frames, interval=250, blit=False)\n",
    "display.display(display.HTML(animator.to_html5_video()))\n",
    "plt.close();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 3</h3>\n",
    "Change to the shifted rosenbrock function. Does CMA-ES get stuck in a local minimum? If so, are there parameters which allow it to find the global optimum?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "celltoolbar": "Slideshow",
  "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}