{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "***\n", "# Lecture 3 - Gradient Descent\n", "***\n", "$\\newcommand{\\vct}[1]{\\mathbf{#1}}$\n", "$\\newcommand{\\mtx}[1]{\\mathbf{#1}}$\n", "$\\newcommand{\\e}{\\varepsilon}$\n", "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n", "$\\newcommand{\\minimize}{\\text{minimize}\\quad}$\n", "$\\newcommand{\\maximize}{\\text{maximize}\\quad}$\n", "$\\newcommand{\\subjto}{\\quad\\text{subject to}\\quad}$\n", "$\\newcommand{\\R}{\\mathbb{R}}$\n", "$\\newcommand{\\trans}{T}$\n", "$\\newcommand{\\ip}[2]{\\langle {#1}, {#2} \\rangle}$\n", "$\\newcommand{\\zerovct}{\\vct{0}}$\n", "$\\newcommand{\\diff}[1]{\\mathrm{d}{#1}}$" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Most modern optimization methods are iterative: they generate a sequence of points $\\vct{x}_0,\\vct{x}_1,\\dots$ in $\\R^n$\n", "in the hope that this sequences will converge to a local or global minimizer $\\vct{x}^*$ of a function $f(\\vct{x})$. A typical rule for generating such a sequence would be to start with a vector $\\vct{x}_0$, chosen by an educated guess, and then for $k\\geq 0$, move from step $k$ to $k+1$ by\n", "\n", "\\begin{equation*}\n", " \\vct{x}_{k+1} = \\vct{x}_k+\\alpha_k\\vct{p}_k,\n", "\\end{equation*}\n", "\n", "in a way that ensures that $f(\\vct{x}_{k+1})\\leq f(\\vct{x}_k)$.\n", "The parameter $\\alpha_k$ is called the **step length**, while $\\vct{p}_k$ is the **search direction**. In this lecture we discuss one such method, the method of gradient descent, or steepest descent, and discuss how to select the right step length." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Gradient descent\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the method of gradient descent, the search direction is chosen as\n", "\n", "\\begin{equation*}\n", " \\vct{p}_k = -\\nabla f(\\vct{x}_k).\n", "\\end{equation*}\n", "\n", "To see why this makes sense, let $\\vct{p}$ be a direction and consider the Taylor expansion\n", "\n", "\\begin{equation*}\n", " f(\\vct{x}_k+\\alpha \\vct{p}) = f(\\vct{x}_k)+\\alpha \\ip{\\vct{p}}{\\nabla f(\\vct{x}_k)}+O(\\alpha^2).\n", "\\end{equation*}\n", "\n", "Considering this as a function of $\\alpha$, the rate of change in direction $\\vct{p}$ at $\\vct{x}_k$ is the derivative of this function at $\\alpha=0$,\n", "\n", "\\begin{equation*}\n", " \\frac{\\diff{f}(\\vct{x}_k+\\alpha \\vct{p})}{\\diff{\\alpha}}|_{\\alpha=0} = \\ip{\\vct{p}}{\\nabla f(\\vct{x}_k)},\n", "\\end{equation*}\n", "\n", "also known as the **directional derivative** of $f$ at $\\vct{x}_k$ in the direction $\\vct{p}$.\n", "This formula indicates that the rate of change is *negative*, and we have a **descent direction**, if $\\ip{\\vct{p}}{\\nabla f(\\vct{x}_k)}<0$. \n", "\n", "The Cauchy-Schwarz inequality (see Preliminaries, Page 9) gives the bounds\n", "\n", "\\begin{equation*}\n", " -\\norm{\\vct{p}}_2 \\norm{\\nabla f(\\vct{x}_k)}_2 \\leq \\ip{\\vct{p}}{\\nabla f(\\vct{x}_k)} \\leq \\norm{\\vct{p}}_2 \\norm{\\nabla f(\\vct{x}_k)}_2.\n", "\\end{equation*}\n", "\n", "We see that the rate of change is the smallest when the first inequality is an equality, which happens if \n", "\n", "\\begin{equation*}\n", " \\vct{p} = -\\alpha \\nabla f(\\vct{x}_k)\n", "\\end{equation*}\n", "\n", "for some $\\alpha>0$. \n", "\n", "For a visual interpretation of what it means to be a descent direction, note that the **angle** $\\theta$ between a vector $\\vct{p}$ and the gradient $\\nabla f(\\vct{x})$ at a point $\\vct{x}$ is given by (see Preliminaries, Page 9)\n", "\n", "\\begin{equation*}\n", " \\ip{\\vct{x}}{\\nabla f(\\vct{x})} = \\norm{\\vct{p}}_2\\norm{\\nabla f(\\vct{x})}_2 \\cos(\\theta).\n", "\\end{equation*}\n", "\n", "This is negative if the vector $\\vct{p}$ forms and angle greater than $\\pi/2$ with the gradient. Recall that the gradient points in the direction of steepest ascent, and is orthogonal to the *level sets*. If you are standing on the slope of a mountain, walking along the level set lines will not change your elevation, the gradient points to the steepest upward direction, and the negative gradient to the steepest descent.\n", "\n", "![Descent](images/descent.png)\n", "\n", "Any multiple $\\alpha \\nabla f(\\vct{x}_k)$ points in the direction of steepest descent, but we have to choose a sensible parameter $\\alpha$ to ensure that we make sufficient progress, but at the same time don't overshoot. Ideally, we would choose the value $\\alpha_k$ that minimizes $f(\\vct{x}_k-\\alpha_k \\nabla f(\\vct{x}_k))$. While finding such a minimizer is in general not easy, for quadratic functions in can be given in closed form.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "### Linear least squares\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider a function of the form\n", "\n", "\\begin{equation*}\n", " f(\\vct{x}) = \\frac{1}{2}\\norm{\\mtx{A}\\vct{x}-\\vct{b}}_2^2.\n", "\\end{equation*}\n", "\n", "In Problem Sheet 1 you will show that that the Hessian is symmetric and positive semidefinite, with the gradient given by\n", "\n", "\\begin{equation*}\n", " \\nabla f(\\vct{x}) = \\mtx{A}^{\\trans}(\\mtx{A}\\vct{x}-\\vct{b}).\n", "\\end{equation*}\n", "\n", "The method of gradient descent proceeds as\n", "\n", "\\begin{equation*}\n", " \\vct{x}_{k+1} = \\vct{x}_k-\\alpha_k \\mtx{A}^{\\top}(\\mtx{A}\\vct{x}_k-\\vct{b}).\n", "\\end{equation*}\n", "\n", "To find the best $\\alpha_k$, we compute the minimum of the function\n", "\n", "[1]\n", "\\begin{equation*}\n", " \\alpha \\mapsto f(\\vct{x}_k-\\alpha \\mtx{A}^{\\trans}(\\mtx{A}\\vct{x}_k-\\vct{b})).\n", "\\end{equation*}\n", "\n", "If we set $\\vct{r}_k:=\\mtx{A}^{\\trans}(\\vct{b}-\\mtx{A}\\vct{x}_k) = -\\nabla f(\\vct{x}_k)$ and compute the minimum of [1] by differentiating, we get the step length\n", "\n", "\\begin{equation*}\n", " \\alpha_k = \\frac{\\vct{r}_k^{\\trans}\\vct{r}_k}{\\vct{r}_k^{\\trans}\\mtx{A}^\\trans\\mtx{A}\\vct{r}_k} = \\frac{\\norm{\\vct{r}_k}_2^2}{\\norm{\\mtx{A}\\vct{r}_k}_2^2}.\n", "\\end{equation*}\n", "\n", "(Verify this!) Note also that when we have $\\vct{r}_k$ and $\\alpha_k$, we can compute the next $\\vct{r}_k$ as\n", "\n", "\\begin{align*}\n", " r_{k+1} &= \\mtx{A}^{\\trans}(\\vct{b}-\\mtx{A}\\vct{x}_{k+1})\\\\\n", " &= \\mtx{A}^{\\trans}(\\vct{b}-\\mtx{A}(\\vct{x}_{k}+\\alpha_k \\vct{r}_k))\\\\\n", " &= \\mtx{A}^{\\trans}(\\vct{b}-\\mtx{A}\\vct{x}_k - \\alpha_k \\mtx{A}^{\\trans}\\mtx{A}\\vct{r}_k) = \\vct{r}_k - \\alpha_k \\mtx{A}^{\\trans}\\mtx{A}\\vct{r}_k.\n", "\\end{align*}\n", "\n", "The gradient descent algorithm for the linear least squares problem proceeds by first computing $\\vct{r}_0=\\mtx{A}^{\\trans}(\\vct{b}-\\mtx{A}\\vct{x}_0)$, and then at each step\n", "\n", "\\begin{align*}\n", " \\alpha_k &= \\frac{\\vct{r}_k^{\\trans}\\vct{r}_k}{\\vct{r}_k^{\\trans}\\mtx{A}^\\trans\\mtx{A}\\vct{r}_k}\\\\\n", " \\vct{x}_{k+1} &= \\vct{x}_k+\\alpha_k \\vct{r}_k\\\\\n", " \\vct{r}_{k+1} &= \\vct{r}_k - \\alpha_{k}\\mtx{A}^{\\trans}\\mtx{A}\\vct{r}_k.\n", "\\end{align*}\n", "\n", "Does this work? How do we know when to stop? It is worth noting that the residual satisfies $\\vct{r}=0$ if and only if $\\vct{x}$ is a stationary point, in our case, a minimizer. One criteria for stopping could then be to check whether $\\norm{\\vct{r}_k}_2\\leq \\e$ for some given tolerance $\\e>0$. One potential problem with this criterion is that the function can become *flat* long before reaching a minimum, so an alternative stopping method would be to stop when the difference between two successive points, $\\norm{\\vct{x}_{k+1}-\\vct{x}_k}_2$, becomes smaller than some $\\e>0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example.** We plot the trajectory of gradient descent with the data\n", "\\begin{equation*}\n", "\\mtx{A} = \\begin{pmatrix}\n", "2 & 0\\\\\n", "1 & 3\\\\\n", "0 & 1\n", "\\end{pmatrix}, \\quad\n", "\\vct{b} = \\begin{pmatrix}\n", "1\\\\ -1\\\\ 0\n", "\\end{pmatrix}.\n", "\\end{equation*}\n", "\n", "The Python code below implements the gradient descent algorithm above. The stopping criteria used is that the residual $\\vct{r}_k$ becomes smaller than the tolerance in the 2-norm." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXyR6ykYTsbCEQ9k2URVCDKCK2gliirdpq\n3X62tnVpq1Kp+OVbC/227l9bl/IFa1sKaJXiBoqpVQFF9h2SsGUjC9m3ycz5/XETDGGyzkzuLJ/n\n4zEPh+Ew98MQP7l533PPUVprhBBCeCc/swsQQgjhOtLkhRDCi0mTF0IILyZNXgghvJg0eSGE8GLS\n5IUQwos5pckrpeYopQ4ppY4opR5pZ8zzSqmjSqldSqkJzjiuEEKIjjnc5JVSfsCLwDXAaOC7SqkR\nbcZcC6RprYcB9wJ/cvS4QgghOueMM/nJwFGt9QmttQVYDcxrM2Ye8DqA1nobEKWUSnDCsYUQQnTA\nGU0+BTjV6tenm1/raEyenTFCCCGcTC68CiGEFwtwwnvkAQNb/bp/82ttxwzoZAwASilZTEcIIbpJ\na63sve6MM/mvgKFKqUFKqSDgZmB9mzHrge8DKKWmAuVa66IOiu3S44knnujyWLMet9xyi92/4y23\n3GJ6bd7yGXt6zZ5Wr9TsfvV2xOEmr7W2AvcDG4H9wGqt9UGl1L1KqXuax7wH5CqljgEvAz9y9Lie\nYunSpaSlpZ33WlpaGkuXLjWpIiGEL3FGXIPW+gNgeJvXXm7z6/udcSxPk5qayqZNm1i8eDGfffYZ\npaWlPP3006SmpppdmhDCB3j0hdeMjAyzS+iS1NRU3njjDVauXMlPfvITPv/8c7NL6jJP+Yxb87Sa\nPa1ekJp7g7PqVZ3lOb1NKaXdrSZn2rVrFwsWLCA7Oxul7F4nEUKIblFKoV144VV0w/jx4/H392fH\njh1mlyKE8AHS5HuZUoqFCxeyZs0as0sRQvgAiWtMIJGNEMKZJK5xMxLZCCF6izR5E0hkI4ToLRLX\nmEQiGyGEs0hc44YkshFC9AZp8iaRyEYI0RskrjGRRDZCCGeQuMZNSWQjhHA1afImkshGCOFqEteY\nTCIbIYSjJK5xYxLZCCFcSZq8ySSyEUK4ksQ1bkAiGyGEIySucXMS2QghXEWavBuQyEYI4SoS17gJ\niWyEED0lcY0HkMhGCOEK0uTdhEQ2QghXkLjGjUhkI4ToCYlrPIRENkIIZ5Mm70YkshFCOJvENW5G\nIhshRHdJXONBJLIRQjiTNHk3I5GNEMKZJK5xQxLZCCG6Q+IaDyORjRDCWaTJuyGJbIQQziJxjZuS\nyEYI0VUS13ggiWyEEM4gTd5NSWQjhHAGiWvcmEQ2QoiukLjGQ0lkI4RwlENNXikVrZTaqJQ6rJT6\nUCkV1c64PyulipRSexw5nq+RyEYI4ShHz+QfBT7SWg8HNgOPtTPu/4BrHDyWT8rMzGTt2rVIhCWE\n6AlHm/w8YFXz81XAfHuDtNafAWcdPJZPkshGCOEIR5t8vNa6CEBrXQjEO16SaE0iGyGEIwI6G6CU\n2gQktH4J0MDjdoY7JVNYsmTJuecZGRlkZGQ44209VmZmJgsWLGDZsmUyy0YIQVZWFllZWV0a69AU\nSqXUQSBDa12klEoEPtFaj2xn7CDgX1rrcZ28p0yhbENrTXp6OqtXr2bSpElmlyOEcDOunEK5Hri9\n+fkPgHc6qqP5IbpJIhshRE852uSXA1crpQ4Ds4BlAEqpJKXUhpZBSqm/AV8A6Uqpk0qpOxw8rs+R\nWTZCiJ6QO149hEQ2Qoj2yB2vXkAiGyFET8iZvAeRtWyEEPbImbyXkBujhBDdJU3eg0hkI4ToLolr\nPIxENkKItiSu8SIS2QghukOavIeRyEYI0R0S13ggiWyEEK1JXONlJLIRQnSVNHkPJJGNEKKrJK7x\nUBLZCCFaSFzjhSSyEUJ0hTR5DyWRjRCiKySu8WAS2QghQOIaryWRjRCiM9LkPZhENkKIzkhc4+Ek\nshFCSFzjxSSyEUJ0RJq8h5PIRgjREYlrvIBENkL4NolrvJxENkKI9kiT9wIS2Qgh2iNxjZeQyEYI\n3yVxjQ+QyEYIYY80eS8hkY0Qwh6Ja7yIRDZC+CaJa3yERDZCiLakyXsRiWyEEG1JXONlJLIRwvdI\nXONDJLIRQrQmTd7LSGQjhGhN4hovJJGNEL5F4hofI5GNEKKFNHkvJJGNEKKFxDVeSiIbIXyHy+Ia\npVS0UmqjUuqwUupDpVSUnTH9lVKblVL7lVJ7lVI/deSYomskshFCgONxzaPAR1rr4cBm4DE7Y5qA\nh7TWo4FpwI+VUiMcPK7ohEQ2QghwMK5RSh0CrtBaFymlEoEsrXWHDVwp9Tbwgtb643Z+X+IaJ5HI\nxnNoDTYbNNmM/9q08Wj9f4IfoBT4+YF/q4cQHcU1AQ6+d7zWughAa12olIrvpJDBwARgm4PHFV3Q\nOrKZNGmS2eX4lMYmqKqDmnqoaYDaBqhrNB71FmhohIYmY5ylCZqsRgNvadx+fqAwXmuh9TfN32oz\nHgCB/sYjKAACAyA40HiEBEJIEIQGQZ8g6BNsPMKCjXHCN3T6T62U2gQktH4J4wTjcTvD2z0FV0qF\nA+uAn2mtq7tZp+iB1pGNNHnna2yCsmoor4azNVBRC5W1UFkHViuEh0JYiNFUw4KN5/0ijeYbHAjB\nARAUCEH+EOBvNPbustqMbxCWJmi0QqPF+ObRYIH6RqizwNlqyG+A2sZvvun4+0F4iFFjRAhEhEJk\nH4gMhb5hRn3CO3Ta5LXWV7f3e0qpIqVUQqu45kw74wIwGvxftNbvdHbMJUuWnHuekZFBRkZGZ39E\ntCMzM5MFCxawbNkyiWwcUNsAZyqMR0kVlFYaTbNvGMSEQ98+MCQBopobZUjQ+WfhrtJy5t+dpqy1\n8U2gut54VNUZ35iKC4xvVBW1xnv2DTMe0eEQHQYxEcbfTb6MzJeVlUVWVlaXxjqayS8HyrTWy5VS\njwDRWutH7Yx7HSjRWj/UhfeUTN6JtNakp6ezevVqOZvvIq2bz37LIP8sFJUbTTEuCuKjIC7SOCOP\n7AN+XtjwtDa+gVXUGD+hnK02HmXVRtQUEw6xEdAvwvgc+kUYP5EI83SUyTva5GOANcAA4ASQqbUu\nV0olAa9qrb+llJoOfArsxYhzNLBIa/1BO+8pTd7JFi1ahNVqZfny5WaX4rZqG+BkCZxqfgT6Q3IM\nJEVDYrRxJitnsMY3u7JqKK2CkkoorjR+HR5ifPOLj4KEvsbzAH+zq/UdLmvyriBN3vlkls2FtDaa\nU04h5J4xsvT+sTAgDgbEGmfpomtsNuOMvyXOOlNhfLbRYUbDT4o2HhGhZlfqvaTJ+ziJbL5xthqO\nFEB2AVisRo4+JME4W5fpiM7TZDXO8gvLofAsFJw1Pt/kGEiOhpRY4/qFnHM4hzR54dORTYMFjuTD\noTzjQuOwJBiaBAlR0mR6i9bGBd38MuORV2ZMBU2Jaf4Jqp+c6TtCmrzwycimqBz2nYScIhjYD0b2\nh/79vPNiqafR2pjRk1cKp0rhdIkxI2lAP+PfKjnGmPcvukaavPCZyMamjZx9V64xQ2TMQBiZAqHB\nZlcmOqK1cSG35eL3mQojzx8cB4Pijamcon3S5AXg3ZGN1WrEMTtyjbs7J6RCaoKctXuqxiY4XQon\nzsDxYuOsfnA8pMYb10/k3/V80uQF4J2RjdVmNPftx4z525PSjB/1hffQ2riIm3sGcouM6a6D440L\n5gNiwV+makqTFwZvimy0NrL2rUeM5QKmDjPO8IT3q6w1/u1ziqCsCgbGGRfSB/bz3bn50uTFOd4Q\n2ZRUwn8OGrNmpo8wLtYJ31RTbzT77ELj62JwgjF7qn+sb02JlSYvzvHkyKbRAtuOwtECmDwMRg2Q\nbFZ8o6YejhUaXx8VNZCWCMNTILGv90+VlSYvzvHUyOZEMWTtM87QLh1hLJ8rRHsqa42b3o7kGddt\nLh5qTKH1VtLkxXk8KbJpbILPDhozLWaOkWhGdE/LRdsmq3dfkJcmL87jKZHNmQrYuMtY9+SyUXJz\njBDtcdlG3sIzufsm31rD3hOwYTtMSYdZ46TBC9FT0uR9kDtv8t1khc17Yf8puHGqMVNCCNFzEtf4\nKHeMbOoa4L0dxrz3WWNlH1IhukriGnEBd4tsymvgza3GErTXTJAGL4SzSJP3Ue4U2ZRUwtvbYGIq\nTE33/jnNQvQmiWt8mDtENmcqjAusl42S/F2InpK4RthldmRTUmk0+Iwx0uCFcBVp8j7MzMimvAb+\ntR0uH2WsJiiEcA1p8j4uMzOTtWvX0psRWW0D/OsrmDLMWD1QCOE60uR9XG9HNk1WeO9rSE8xFhgT\nQriWNHkf15uRjdbw7/0QHgqTh7r8cEIIpMkLei+y2X/KWCxq1liZJilEb5EmL3olsimphG1HYM5E\nudFJiN4kTV64PLJpssKm3TB9JPQNc8khhBDtkCYvANdGNl8dM5r78GSnv7UQohPS5AXgusimpBIO\nnoYrRksOL4QZpMkLwDWRjdaQtd9Yj6ZPsNPeVgjRDdLkxTnOjmyO5BuN3pv31hTC3UmTF+c4M7Jp\nssLWIzBjpMQ0QphJmrw4x5mRzf5T0C/S2J9VCGEeafLiPM6IbJqssCMHJg9zYmFCiB6RJi/O44zI\n5lAexEUaDyGEuaTJi/M4GtloDbuPG7s8CSHMJ01eXMCRyOZ0Kfj7QXKMCwoTQnSbQ01eKRWtlNqo\nlDqslPpQKRVlZ0ywUmqbUmqnUmqvUuoJR44pXM+RyGb/KRg9QGbUCOEuHD2TfxT4SGs9HNgMPNZ2\ngNa6AZiptZ4ITACuVUpNdvC4woV6GtnUW+BUCQyT5QuEcBuONvl5wKrm56uA+fYGaa1rm58GAwGA\n7NTt5noS2eQWQf9YCAl0YWFCiG5xtMnHa62LALTWhUC8vUFKKT+l1E6gENiktf7KweMKF+tJZHOs\nEIYmurAoIUS3ddrklVKblFJ7Wj32Nv/3ejvD7Z72aa1tzXFNf2CKUmqUg3ULF+tuZNPYBAVlMCjO\nxYUJIbql0+0btNZXt/d7SqkipVSC1rpIKZUInOnkvSqVUp8Ac4AD7Y1bsmTJuecZGRlkZGR0VqZw\ngczMTBYsWMCyZctQnVxJzSuF+L4QJFGNEC6XlZVFVlZWl8YqR+5sVEotB8q01suVUo8A0VrrR9uM\n6QdYtNYVSqlQ4ENgmdb6vXbeU7t6GzrRNVpr0tPTWb16NZMmTepw7H8OGCtNTkrrpeKEEOcopdBa\n2z0TczSTXw5crZQ6DMwCljUfMEkptaF5TBLwiVJqF7AN+LC9Bi/cS3cim/yzMjdeCHfk0Jm8K8iZ\nvHvZtWsXCxYsIDs7u93IxtIEKzbDXbPA37+XCxRCdHgmL1sqiw61nmXTXmRTUgXRYdLgPZnNpqlr\nsFFXb6OhUdPUpLHaNFobN7b5+ykCAhTBQYqQYD9Cg/3w95c73jyBNHnRodaRTbtNvtJYVli4r4ZG\nG6eLGskrbKSg2EJhiYXS8ibKKpqorLZSU2cjJMiPkGBFcJAfAQGKAH+FAmxaY7OBpUnTaDG+GdQ3\n2AgO8iMizI/IcH+iIwOIiQogLjqA+NhA4mMDSIoLIiLMr9OL9sK1JK4Rneossvn0AESFwnhZlKzX\n5ebmsnjxYvLy8khJSWHp0qUMGjSYkwWNHMyu48jxeo6dbOBMqYWkuED6JwaRFBdIQr9A4qIDiY7y\nJyrcn/Awf/z9ut6MW878q2psVFQ1cbbSSll5EyVnmzhTZqGoxEJ+sQU/Bf0Tg+ifGMSg5GAGJwcx\nuH8w4X3kxz5n6iiukSYvOtXZLJv1X8K4wTDY7q1wwlVyc3O5+uqryc7OPvdabNwgLrthJUnJgxk1\nNJThqSEMHRhC/8QgAgN694xaa01FtZXThY2cKmjkRH4jJ/IbOJ7XQFREAMMGBZM+OIThqaEMGRBM\ngMQ/PSaZvHBIZ5FNZR1E9TGhMB/30MM/P6/BA5QWn8C/7BVeevlvJlX1DaUUfSMC6BsRwJhh33yB\n2GyavDMWjh6v53BuHR9vqaSo1MKwQSGMGRbKmPQ+pA8OkabvJHImL7qkvchGa3h5I9w5CwLllMHl\nSkpKeP0vb/LnlX/nwN7/gLZdMGbmzJls3rzZhOp6rrrWysHsOvYfrWPv0ToKiy2MHhrKxFF9mDQ6\njPhYucuuI3ImLxzW3iybxibwU9LgXamkpIS3336b1/+ymq+++pKEgTO47tu3MTQ1jvXvrLtgfHKy\n5y0DGt7Hn0vGhnPJ2HAAKqut7D5Uy44DNax+r4yoCH8mjw1j8rhwhg0Klou53SBn8qLLFi1ahNVq\nZfny5edeq6iB9V/BbRnm1eWNWhr7mjVr2Lp1G2kjryAqeTb33jGP+bNTCA3xs5vJp6WlsWnTJlJT\nvecquNWmOXq8nq/21rBtTzX1DZppE8KZflE4w1NDpOEjF16Fk9iLbIorYPNeuGmGycV5gdaNfdu2\nbcyefQ2DR1zL6aqLmTsziRtnR9Mn9PxZKS2za06fPs2WLVv4+OOPmTHDu/8xThU08MXOav7zdRWN\nFs2MSRFkTI5gYFKw2aWZRpq8cAp7s2zyy2DrEVgw1eTiPFTbxn7NNdewcOFCJky6ilfWVREW6se9\nN8eTFBfU6XvdeeedjB49moceeqgXKjef1prjeY38Z3sVn26vIjrSnyunRnLZxRE+N0VTmrxwmraR\nzakS+DoH5steX13WXmOfO3cuYWFhfLKtkpX/LOHmuTHMuSyqy3HEBx98wJNPPsmWLVtc/DdwP1ab\nZs+hWj7eWsnOg7VcPCaMa2ZEMXKIb8Q50uSF07SNbE4Ww65cuF6afIc6a+xgNKq/vFPCtt01PHZv\nUrfjB4vFQmJiIjt37mTgwIGu+Gt4hMpqK1lfVvLhZxUEBijmXBZFxuRIQoIdXY/RfUmTF07TNrI5\nWQw7c2GeNPkLdKWxt7A0aZ5ZVUhllZVf3pVEZHjP4gZfi2w6orVmz+E63v+0nNi+Adyd6b1360mT\nF07VOrLJK4Uvj8ENU8yuyj10p7G3sDRpfr+iAG2DX9yZSGBgz884fTmy6YjNpvHrxrINnkaavHCq\n1pFNUbniPwdh4aVmV2WenjT2Flprnl1VRF2DjV/cmeTw0gMS2fgmV24aInxQ6xujggKM9eR9TUlJ\nCa+99hqzZ88mLS2NjRs3cvfdd5Ofn8+aNWtYuHBhpw0eYPV7ZRSWWHj4jkSnrC0TGBjI/PnzWbfu\nwpukhG+SJi+6rfVaNsGBUG8xu6Le4azG3mL7vho+2VrJY/ckERzkvP8VFy5cyNq1a532fsKzSVwj\neqQlsjlyNJtXNyn+3zXG5hLexpEopiPllU08tOwkD/8widFDQ51YsUQ2vkjiGuF0LZHN7l1GZFPX\naHZFzuPsM3Z7/vxmMRmTI53e4EEiG3E+afKiR1pHNuEhUF1vdkWO6Y3G3mLfkVoO59Zz01zX7Xwu\nkY1oIXGN6LGWyOaFtdkMT1YMTTK7ou5xVRTTEa01jz19mrmXR3H5Ja7bM1EiG98icY1wiZbIpvD4\nDipqza6ma3rzjN2e3YdqqamzMX1ShMuOARLZiG9Ikxc91hLZfPHRGsprzK6mfWY39tY2ZJUz78q+\n3dpPtackshEgTV44KDMzk4/eX0tplXtFbO7U2FuUlTdxOLeeyy527Vl8i1mzZnHkyBFOnjzZK8cT\n7kmavHDI+PHjCQz0Z9euHdgu3ImuV7ljY2/ti13VXDI2zKlz4jsikY0AafLCQUopMhcuZM/nazhr\nQmTj7o29tS/3VDN1fHivHlMiGyFNXjgsMzOTrz9dS+HZ3olsPKmxt2hotHH0RD1jhjl/XnxHJLIR\nsv2ycNj48eMJCvTnsy07GD1wUud/oAfsTXe8++67+ec//+l2Dd2e3NMNpMQHXbB9n6u1jmx6Y/lh\nq1VTU9tEda2VxkYbTVZtbPQe6EdIsB8R4QGEBPvWrk1mkyYvHKaU4oYFC9nwzhruvcl5Td7TG3tr\nOacaGDLAnD1IFy5cyJNPPunUJm+zaXJP1rLvcBVHc2s4mVdHwZkGzlZYCAn2I6xPAMFBfgT4K2xa\nY7Fo6uqtVNU04acUsTGBxMUEk5IYzIDkUIYM7MPQ1DBiozvf5lB0j9wMJZxix85dXDN3Abk52YSH\n9nx6oBk3KPWG19aeIT42kOuvjO71Y1ssFpKSkti5cycDBgzo8fs0NdnYvqeCT74o5cud5fQJ9Wfs\nyAhGpIUzqH8oSQkhxPYN7HA9fK01dfU2Ss82UlTSQH5hAyfz6sg5Wcux3BqCghQjh0YwZkQE40dF\nMiw1DH9/L1wUyclkPXnhclpr+g9K5/mXV3Pjtd07m/fWxt7aUy/nc+XUyF6/8NrizjvvZMyYMTz4\n4IPd/rOVVRb++UEhb39QRHJCMFdO78elF0eTlBDi1Bq11hSeaWD/kSr2Hqpi1/5KSs82cvG4vkyb\nFM20SX2JjAh06jG9hTR50Svuvn8RlTVW/vF/yzsd6wuNvbXHnj7FbfP6MSqtdy+8tujJjlGNFhtv\nvVfI39/J49JJ0dw8L5lB/fu4sMoLlZ5t5Mud5Xy+vYyd+yoZOTScjEtjuXxKjDT8VqTJi17x2dZd\n3HDDAgrzsu3e0elrjb21n/7mBA/fnsigFHNy+e5GNidO17Lk6SMk9Avmvh8MZlCKOd+cWqurt7Jt\nZzmffF7C9j0VXDK+L9fNimfSuCiv3tqvKzpq8nLhVThNcnwk1RXFXDTpEsaOHsHSpUuJiIjwmoun\njrDZtKnZcmBgIPPmzWPdunWdRjafbi3lD6/kcM8tA5l7ZTzKTTYKCA3xJ2NaLBnTYqmqaWLzZyW8\n/MYJauusXD87gblXxsvZvR1yJi+cIjc3l6uvvprs7Oxzr4WGhuLn58fcuXN95oy9PfcvPc4jdyUx\nIMmcM3noWmTz8WclvLTqOL99bATpQ8y5ftAdWmsOHavmn+8XsmXHWWZfHkfmt5NJiDPvczaDxDXC\n5W699Vb++te/XvD6TTfdxOrVq02oyL38fPlJ7r05nmGDnHuxsjs6i2y27y7nNy8c4w+/HsWQgb2b\nvTtDcWkD694t4P3NxUy/JJpbb+xPSqJ5n3dvctlSw0qpaKXURqXUYaXUh0qpqA7G+imldiil1jty\nTOGe8vLy7L5+5syZXq7EPYX18aOm1mpqDa0jm7bKKyw89cIxnnhwmEc2eIC42GDu+/5g3nhhAvH9\ngrnvsb0882oOpWe9aNuyHnB0WYNHgY+01sOBzcBjHYz9GXDAweMJN5WSkmL39eTk5F6uxD1FRwZQ\nVmFukwfO7ebV1vMrcrnqsn5MGN3ueZrHiIwI5I6bBvD6cxMIDvLjjgd389b7BWaXZRpHm/w8YFXz\n81XAfHuDlFL9gbnAaw4eT7ippUuXkpaWdt5rQ4aksXTpUpMqci9xMQGcKbWYXQazZs3i6NGjnDp1\n6txruadq2bm/kh/e3PMbpdxR38hAfvSDwbz8u7EMHeSb14LA8SYfr7UuAtBaFwLx7Yx7BvgFIGG7\nl0pNTWXTpk3ccsstzJw5k8tm38KKv28iNTXV7NLcQnJ8EPlnzI8N7EU2azcUcMOcRK9dUyYpPoRx\no1y31aK763QKpVJqE5DQ+iWMZv24neEXNHGl1HVAkdZ6l1Iqo/nPd2jJkiXnnmdkZJCRkdHZHxFu\nIDU1lTfeeAOAQ6fhWKHJBbmRwSnBvLXprNllAN+sZfPggw/SZNV8tq2M134/zmXHy83NZfHixeTl\n5ZGSksLSpUvlm7+DsrKyyMrK6tJYh2bXKKUOAhla6yKlVCLwidZ6ZJsxTwG3Ak1AKBABvKW1/n47\n7ymza7yAxQqrPoHM6RBp/n00pmuyam79RTYrfpPa6ytRttV6lk15TRTP//k4r/6Pa5q8vam1LctD\nJyQOpL7BhtWqCQn2o0+fAAJknZoeceXNUOuB24HlwA+Ad9oO0FovAhY1F3IF8HB7DV54j0B/GJ4M\n+0/CtOFmV2O+AH9F2sAQDufWM3GUuflw68hm0IibGTHUdfPhFy9efF6DB8jOzmb6zP/H2KlLCA3x\nw89PUd9go77BSmREAMkJoQwe0If0tDBGpUcydLAsUuYIR5v8cmCNUuqHwAkgE0AplQS8qrX+loPv\nLzzY2EHw1la4ZCgEeGfc2y2jh4ay90id6U0evolsvnvXfAYmu2YueV5eHl98sdXu76UPtrDxH5ee\n91qTVXO2vJHTBXUcP1nL4WPVvLkhn/JKC1MviiHj0n5cMjG617ZP9BZyM5RwqXe3w+AEGO1dEzd6\n5FBOHX/8+xme+9Ugs0s5F9n88CdvMfPyEVw7s705E92Tl5fHunXrWLt2Lbt278NGGHXV+ReMu+WW\nW85dv+lM4Zl6Pv+qjE8+Lyb3ZC0zL+3Ht69JYnia+9+R21tcdjOUEJ0Znwq7ckG+b8OwwSGUV1kp\nLDF/KmVLZLN35wcORyF5eXk899xzzJgxg7Fjx7Jz504umX4X19/6AR9/9NEFU2vT0ro3tTYxPoQb\nr0vmxafGs+KZicT3C2bRb/bzs8V72L7rLHJS2DFp8sKlUmIgKABy5cZX/P0UU8eH8cXOKrNLAYzI\n5uCe96mr6/5NWvYa+6JFiygsLOTWH/6OworRvPjbi5g2ZeS5qbWjRo0iKSmJTZt6PrU2IS6E72cO\nZPXLlzBnZgLPvJLNfY/sZs+Bih69ny+QuEa4XHYh7MiB70wDN1nQ0DT7jtby6ppinl000PTVHS0W\nCzGxCSx+6j1+ef/UTse3jmIOHDjA9ddfT2ZmJldddRVBQca2fU1NNr5733Z+/dAIxo48f256dXU1\nKSkp5ObmEhMT45S/g9Wq+fg/xbz8l+OMHxXJj+4YQr8Y39tCUOIaYaohCcaUypMlZldivlFpoTRY\nNMdONphdCoGBgcy4/Fo++vCCSXHndHTGvnLlSubOnXuuwQNkfVFCckLIBQ0eIDw8nKuuuop33mn/\neN3l76+YnRHPG/87ifi4YG7/2dds2FQoEU4rciYvesXRfNh9HG6Us3ne3FhGQbGF+29J6Hywiz3z\n3Cs8+si5og8wAAASe0lEQVQvuXTaxHM3KgUFBXV6xt6eP/zpGANTQln4bftrGa1evZrXX3+d9957\n79xr1TVNfP5lKXsPVlJQVEd9g43AQEVcTDD9k0MZNTyS0cMju3R/Qc6JGv7r6cMMSA7llz8eSkS4\nb6wvL0sNC9PZNPzjM2PO/GDnTOTwWBVVTfz4v07w0hODiQw3b25pbm4uV111FTk5OedeCwkJISgo\niBtuuKHLjb21n/xqDz/IHMDF4+1vWN42stn07zP84Y9HGD+6LxeP70tKUiihIf40WmwUlzZw8nQt\n+w5Vciy3hmkXxzBvThITx/btMOpqaLTxp1W5fLG9jOWPj2bwAM9cVbM7ZGcoYTo/BVOGwdYjMCjO\nt8/moyICmDYxnPf+Xc7N18WaVsfixYvPa/AA9fX1zJ8/n5UrV/boPf0UHTbg1pFNXP85vPbGcV5a\nPpEhzQuI2VsC4Ud3TKSquomN/y7iD388RnCwH7ffNIjLpsbaPVZwkB8/uzuN9LRwfvqrPTzx8xFM\nGte3R38fbyBNXvSa1ATjAuyRfBhu/6d5nzF/VjSLnjnNvFnRhIaYc2msvT0AioqKevyeEeEBnC3v\neCG2hQsXsuL/VqGjhvHnZy5iYIpxpm1vCYStW7eem41z43Up3HBtMp9/Wcqrbxxn/YcF/PL+dOL7\n2d8F6torE0iMC+bX/3OIRT9NZ9rFzrnY62nkwqvoNUoZcc22o9Bk/tLqpkpJCGLCiD689+9y82pw\nwR4AI4dFsP9wx1NEx4wZw+aPP2L/lp+y6JF7yM3NBWDRokV2l0C49ts/5qWVOZw4VYufn+Kyqf1Y\n8exFjB4ewZ0Pfs2/vyhu91gTx/Zl2a9GsXu/D0+x1Fq71cMoSXizDdu1/jrb7CrMd6qgQX//kWxd\nWd1kyvFzcnJ0Wlqaxlg9VgM6LS1N5+Tk9Pg9Dx2t1Jn3fKmtVluXj5mUlKS/853v6ICAgPNeb/0I\nCAjSsYkX6wce26hLyhrOvd/+wxX6htu36Dc3nO5xzd6guW/a7alyJi963fQRsDMHas2fRWiq/olB\nTBkXxpsflply/NZ7AFxxRQb90+bw6v+949AywOlp4YQE+7Fzn/0zZ3sLlhUUFJCTk8O3vtX+UldN\nTY2UFm7nuWXXsfD29WR9bpy9j0qP5MXfjufvb53mzXftx0++Tpq86HV9w2BEf9hy2OxKzHfzdbFs\n3lZp2oYiLXsAZGV9wrPPreC9LOXQHHOlFJnXp7DqHyftvk971wGioqJ4+umn8Q/oeF1qrS0c2/UI\nv//jUfYeNL6RJCeG8vxT41m5+gQHj1T2uHZvJU1emOKSNDhVAoXusY+GaWKiAph/VTQr3iw2/Qae\nedckUFjcwJavHftHuWZmAmXljXz+5YU/oXR0HSA1NZUBgzpfl/rUyRwef3AEjy87QE1tEwBJCSE8\ncM9Qlj59iCarTMFuTZq8MEVQIFw6ArL2g81mdjXm+vbMaIpKm9i6q9rUOgIC/HjgriE8+1ruuebZ\no/fxVzx831D+8KdjlFeevxibvb2AExIHnVuwbPKk8/YcatfUSTEM7t/nvAuqsy6LJyoykE+3yK3V\nrUmTF6YZlgShQbDnhNmVmCswQHHfzfG8tq6Y6lpzpx1NGhfFxeP78tIqx/5RJo7py6zL4vjt80ew\ntjqzbrsX8PwFNzF6yu/JzTOWDV627DcMTh3SpWMkJ4Zw/FTtea/d+K0U3t1U4FDt3kaavDCNUnDF\naPg6GyrrzK7GXKOGhjJlfDgr3mx/OmBv+fEPBrFzfwVZX5Q69D733jaY+norL63MPe/1lusAmzdv\n5p9vruZPz8zh+VeP8da7eaSmprL544+46abvEhwcYfd977//fjZsLGDL9jKunBF3/nsP7ENRsY9f\n0W9DmrwwVd8wY835f++TNedvm9ePA9n1bNtjbmwT1ieAXz+QzrOv5XAyr+fffQMD/fjvR0eybUcZ\nr6852e64Yanh/O/yCbz1bj5P/O4A0TH9Wb36b9TXV3LX3fedN3bA0AUcyM/kb2+d4vmnxpMYf+Gu\nVk1NPv6F1IasXSNMZ7XBui9g3GAY2d/sasx1MLuO371WwNOPDSQ60twb0jd8VMQ/1ufz0m/HEhHW\n81pKyhp5YPEeZk6P44ffbX+J5YYGK6+8cZyNWUXce1sqc2Ylnrexd0ODlZpaKyHBfoSG+tt9n5df\nz6Wu3soD9wztcb2eSBYoE26vuBL+9RVkXgrhHc+i83p/f7eUQzl1/PrHKfj7mbvIzwsrcsk5Wcvv\nfjWSwMCe/+B/tryRh5fsIz0tnIf+31CCOnivw8eqeOG1bArO1HPD3GSunZVIbHTni6S980E+K/52\ngpeWTyAlybe+iKTJC4/w1VEoLIdvXezbC5hZbZolL+QxemioqQuYgbEpx5I/HCEwUPGrnw5zaKvA\n2jor//3MYSqqLDz58xH0i7W/5kyLQ8eqeOvdfD7dUsLIYRFMGt+XoanhpCSFEN78k0VJWSPHT9by\n3keFFBbX8z+/Hkv/ZN9q8CBNXngIqw3e2gojU2CM+Xtdm6qsoolf/O4kP/peApNGh5laS0ODlUd/\ne4jEuGB+cV8afg78dGGzaV5fe4q33s3nwXvTmDk9rtM/U19vZduOMnbtryDnRA35hfXU1jahgdjo\nIPonhzJjSj9mXxHv0E8bnkyavPAYZ6uNRr9gKkSHm12NuQ5m17H8tQKeerA/yfHmbmlXW2flkd8c\nZGBKKA/dM8Thzb/3H67kv589QnpqGPffOYS4Ts7qRcekyQuPsu8k7D9p7Anrb96eGm5h42cVrP/k\nLMseHkB4H3M/jNo6K79afojoqEAeu3+ow2fN9Q1W/rL2FO98UEDm9SncNC+F4GAf/wfvIWnywqNo\nDe/vhMhQmNG1GyC92p/XFXOyoIHH70shMMDcixUNjTb+6+kjNDTaWPJw+rls3BH5hXW8tDKXA0eq\nuGlef749O7FLW/2Jb0iTFx6nvhH+8TlcPsrYbMSXWW2a5a8UEB7mx09uTehw56Xe0GTVvLgil90H\nK/ntYyNJjHNO1HL4WBV/WXeKXfsquO7qRObPSSIp4cJ58OJC0uSFRyo4C+/vMGKbSO/fprND9Q02\nfv38acam9+G2ef3MLgetNeveLeTvb+ex+IFhTBwT5bT3zi+sY92GfDZmnWHwwD7MnB7HZVNi290B\nqqv1mv3N0ZWkyQuPtTMXjhXAgimSz1dWW1n0zCmuvjSKebPsb5Td27bvKec3zx9j4XVJ3Dwv2aGZ\nN201Wmx8ueMs/95awhdflREbHcSkcX0ZOzKSoYPDSE4MISDA/nUBrTVFxQ3sO1TJpk+LSYoP4YF7\n0uyO9QbS5IXH0ho+2GksZJYxxuxqzFdy1sKiZ07zndkxzJ7hvLNnRxQVN/DkM0cI6+PPoz8e2qUb\nl7rLatUcOlbFzr0V7D9cSc7JWkpKG+gbGUhMdBAhIf4E+CsaLTaqqpsoKm4gJMSPUcMiuHxaPzIu\n7UdoiPeeJUiTFx6t0QJrt8BFQ2TZA4CC4kYWP5fH974Vy5VTI80uB4CmJhur1p7m3Y/P8OA9Q7hs\nsus3zW602Cg720hZuYWGBiuWJk1wkB9hYQEkxgU75aKwp5AmLzxeWTW8vQ2umwQJfc2uxnx5RY08\n8UIe370uhlnT3OOMHmDPwUqW/2826UPC+OmdqURHBZpdkk+QJi+8Qk4R/OeAcSE2TCZdnGv0C+fE\ncI2bRDdgzH9fueY0H2YVc9f3BnDtzHinZvXiQtLkhdfYfgxyz8ANUyDAeyPWLisobmTJi3nMuawv\nN1zlHhdjWxzJqebZ13Kx2TQ//WEqo9Ltrw8vHCdNXngNrWHTbuP51eN9eyGzFiVnLSx5MY8p48K5\n9fpYt5oqaLNpNn1awqt/O8G4kZHc+d2BpCTKj2HOJk1eeJUmq5HPD4yDycPMrsY9VFZb+c2f8umf\nEMh930s4bx12d1BbZ2XdhgLWvVvAFdNiue3GFIfmvYvzSZMXXqe2Ad7cApcMgxEpZlfjHuobbPx+\nRQFWK/z8zkTC3HBpgIoqC6vfzmfDx0VcPiWWm65PZmCK7y0N7Gwua/JKqWjgH8Ag4DiQqbWusDPu\nOFAB2ACL1npyB+8pTV50ScuMm6vHwwDzbwJ1C1ar5rV1xew/Vsfj9yUTH+Oes1sqqiy89V4h73xY\nyOjhEXznuiQmjI50q6jJk3TU5B1dfPlR4COt9XBgM/BYO+NsQIbWemJHDb67srKynPVWvcLT6gX3\nrjkmHOZMNDL64spvXnfnmu1xZr3+/op7MuO4+tJIHv39KQ5ku2aHdEdrjooI5I6bBrD6jxcxeUJf\nnn0tlx88sJs1/8qnvNLinCLb8NWvC0eb/DxgVfPzVcD8dsYpJxzrAr76j9ab3L3m5BhjEbN3t0Nl\nrfGau9fclrPrVUrx7ZnR3H9rAstfLWDj5xf8cO0wZ9UcEuzPvGsSWfnMeH5+7xCO5tZwy/07efx3\nh8jaUkp9g9UpxwHf/bpw9JaweK11EYDWulApFd/OOA1sUkpZgVe01q86eFwhzhmaBHWNsP4rY7MR\nYbhoVBhPPdifZa8WcPREPXcvjOtwb1UzKaUYNyqScaMiqaltImtLKRs2FfE/f8zmojFRTJ8czZSJ\n0XJzVQ902uSVUpuA1ou9Koym/bid4e2F6dO11gVKqTiMZn9Qa/1Zt6sVoh1jBxmNfvsxsytxLykJ\nQfzu5wN48a9F/O1fpdy+oPPt9swW1ieA62YlcN2sBCqqLGz5+iyff3mWF1YcZ9LYKP7rF8PNLtGj\nOHrh9SBG1l6klEoEPtFad7jNg1LqCaBKa/10O78vV12FEKKb2rvw6mhcsx64HVgO/AB4p+0ApVQf\nwE9rXa2UCgNmA092t1AhhBDd5+iZfAywBhgAnMCYQlmulEoCXtVaf0splQr8EyPKCQD+qrVe5njp\nQgghOuN2N0MJIYRwHve81N4OpVS0UmqjUuqwUupDpZTdpfeUUn9WShUppfb0do3Nx5+jlDqklDqi\nlHqknTHPK6WOKqV2KaUm9HaNdurpsGal1HCl1BdKqXql1ENm1Nimns7q/Z5Sanfz4zOl1Fgz6mxT\nU2c1X99c706l1JdKqelm1Nmmpk6/lpvHXaKUsiilFvRmfXbq6OwzvkIpVa6U2tH8sDeBpFd1sV9k\nNH9d7FNKfdKtA2itPeaBkf3/svn5I8CydsbNACYAe0yo0Q84hnEXcCCwCxjRZsy1wLvNz6cAW03+\nXLtScz9gErAUeMgD6p0KRDU/n+Mhn3GfVs/HAgfdveZW4z4GNgAL3Lle4ApgvZmfaw9qjgL2AynN\nv+7XnWN41Jk8Xbz5ShvTM8/2VlFtTAaOaq1PaK0twGqMulubB7wOoLXeBkQppRIwT6c1a61LtNZf\nA01mFNhGV+rdqr9ZYmMrYPYKN12pubbVL8Mx7hQ3U1e+lgF+AqwDzvRmcXZ0tV53mtzRlZq/B7yp\ntc4D4//F7hzA05r8eTdfAe3dfGWmFOBUq1+f5sIG03ZMnp0xvakrNbuT7tZ7F/C+SyvqXJdqVkrN\nb56a/C/gh71UW3s6rVkplQzM11r/EfObZ1e/LqY1x6TvKqVG9U5p7epKzelAjFLqE6XUV0qp27pz\nALfbBNFJN18JAYBSaiZwB0aE5/a01m8DbyulZgD/DVxtckmdeRYjOm1hdqPvzNfAQK11rVLqWuBt\njCbqzgKAi4ArgTBgi1Jqi9a6S7f+uV2T11q3+0XdfDE1QX9z85XZPx7akwcMbPXr/s2vtR0zoJMx\nvakrNbuTLtWrlBoHvALM0VqbFd+16NZnrLX+TCk1RCkVo7Uuc3l19nWl5ouB1cpYPrIfcK1SyqK1\nXt9LNbbWab1a6+pWz99XSr3kAZ/xaaBEa10P1CulPgXGY2T5nfK0uKbl5ito5+arVhTmnFV8BQxV\nSg1SSgUBN2PU3dp64PsASqmpQHlLDGWSrtTcmtlna53Wq5QaCLwJ3Ka1zjahxra6UnNaq+cXAUEm\nNh/oQs1a6yHNj1SMXP5HJjV46NpnnNDq+WSMaeRu/Rlj9LkZSin/5ptLpwAHu3wEs68ud/NKdAzw\nEXAY2Aj0bX49CdjQatzfgHygATgJ3NHLdc5prvEo8Gjza/cC97Qa8yLGd+LdwEVu8Nl2WDNGhHYK\nKAfKmj/XcDeu91WgFNgB7AS+9IDP+JfAvuaaPwemuXvNbcauwMTZNV38jH/c/BnvBL4ApnjCZwz8\nHGOGzR7gJ915f7kZSgghvJinxTVCCCG6QZq8EEJ4MWnyQgjhxaTJCyGEF5MmL4QQXkyavBBCeDFp\n8kII4cWkyQshhBf7/9jzT0yj+OLKAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import numpy.linalg as la\n", "import matplotlib.pyplot as plt\n", "\n", "# The function graddesc performs the gradient descent algorithm with starting value x and tolerance tol\n", "def graddesc(A, b, x, tol):\n", " \"\"\"\n", " Input: matrix A, vector b, vector x, tolarance tol\n", " Output: trajectory of points x_0,x_1,... of the gradient descent method\n", " \"\"\"\n", " # Compute the negative gradient r = A^T(b-Ax)\n", " r = np.dot(A.transpose(),b-np.dot(A,x))\n", " # Start with an empty array\n", " xout = [x]\n", " while la.norm(r,2) > tol:\n", " # If the gradient is bigger than the tolerance\n", " Ar = np.dot(A,r)\n", " alpha = np.dot(r,r)/np.dot(Ar,Ar)\n", " x = x + alpha*r\n", " xout.append(x)\n", " r = r-alpha*np.dot(A.transpose(),Ar)\n", " return np.array(xout).transpose()\n", "\n", "# Define the matrix A and the vector b for the problem we consider, as well as the zero starting point\n", "A = np.array([[2,0], [1,3], [0,1]])\n", "b = np.array([1, -1, 0])\n", "tol = 1e-4\n", "x = np.zeros(2)\n", "\n", "# Call the gradient descent function with input A and with starting value x=0\n", "traj = graddesc(A, b, x, tol)\n", "\n", "# Plot\n", "%matplotlib inline\n", "\n", "# Define the function we aim to minimize\n", "def f(x):\n", " return np.dot(np.dot(A,x)-b,np.dot(A,x)-b)\n", "\n", "# Determine the values of the function for the first 7 steps\n", "# This is used to specify the location of the level sets\n", "fvals = []\n", "for i in range(7):\n", " fvals.append(f(traj[:,i]))\n", " \n", "# Create a mesh grid for plotting the contours / level sets \n", "xx = np.linspace(0.1,0.5,100)\n", "yy = np.linspace(-0.5,-0.1,100)\n", "X, Y = np.meshgrid(xx, yy)\n", "Z = np.zeros(X.shape)\n", "for i in range(Z.shape[0]):\n", " for j in range(Z.shape[1]):\n", " Z[i,j] = f(np.array([X[i,j], Y[i,j]]))\n", "\n", "# Get a nice monotone colormap\n", "cmap = plt.cm.get_cmap(\"coolwarm\")\n", "\n", "# Plot the contours and the trajectory\n", "plt.contour(X, Y, Z, cmap = cmap, levels = list(reversed(fvals)))\n", "plt.plot(traj[0,:], traj[1,:], 'o-k')\n", "#plt.axis('off')\n", "plt.xlim([0.1,0.5])\n", "plt.ylim([-0.5,-0.15])\n", "plt.axis('equal')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.12" } }, "nbformat": 4, "nbformat_minor": 0 }