{ "cells": [ { "cell_type": "markdown", "id": "c57b2ce5", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Introduction to Optimization\n", "\n", "\n", "Feng Li\n", "\n", "School of Statistics and Mathematics\n", "\n", "Central University of Finance and Economics\n", "\n", "[feng.li@cufe.edu.cn](mailto:feng.li@cufe.edu.cn)\n", "\n", "[https://feng.li/statcomp](https://feng.li/statcomp)\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "c085cebf", "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "options(jupyter.plot_scale=0.8)" ] }, { "cell_type": "markdown", "id": "b31a0a73", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Motivation: Why Optimize?\n", "\n", "### Optimization in Business\n", "\n", "- Many problems in business require something to be minimized or\n", " maximized\n", "\n", " - Maximizing Revenue\n", "\n", " - Minimizing Costs\n", "\n", " - Minimizing Delivery Time\n", "\n", " - Maximizing Financial Returns" ] }, { "cell_type": "markdown", "id": "2a3de9fc", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Input and output\n", "\n", "- For many of these problems there is some control over the input\n", "\n", " - Maximizing Revenue - Price\n", "\n", " - Minimizing Costs - Number of Workers\n", "\n", " - Minimizing Delivery Time - Driving Route\n", "\n", " - Maximizing Financial Returns - Portfolio weights" ] }, { "cell_type": "markdown", "id": "7f11faed", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Optimization in Statistics\n", "\n", "- In statistics, many estimators maximize or minimize a function\n", "\n", " - Maximum Likelihood\n", "\n", " - Least Squares\n", "\n", " - Method of Moments\n", "\n", " - Posterior Mode" ] }, { "cell_type": "markdown", "id": "93738c5b", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- Suppose we want to find an minimum or maximum of a function $f(x)$\n", "\n", "- Sometimes $f(x)$ will be very complicated\n", "\n", "- Are there computer algorithms that can help?\n", "\n", "- YES!\n", "\n", " - Newton’s Method\n", "\n", " - Quasi-Newton\n", "\n", " - Nelder Mead\n" ] }, { "cell_type": "markdown", "id": "09dd92f9", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Newton’s Method\n", "\n", "### Finding a root\n", "\n", "\n", "- Consider the problem of finding the **root** or **zero** a function.\n", "\n", "- For the function $g(x)$ the **root** is the point $x^*$ such that\n", " $g(x^*)=0$\n", "\n", "- An algorithm for solving this problem was proposed by Newton and\n", " Raphson nearly 500 years ago.\n", "\n", "- We will use this algorithm to find the root of\n", " $g(x)=3x^5-4x^4+6x^3+4x-4$\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "d08f0aea", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "gx<-function(x){\n", " y=(3*(x^5))-(4*(x^4))+(6*(x^3))+(4*x)-4\n", " return(y) \n", "}\n", "\n", "dgx<-function(x){\n", " y=(15*(x^4))-(16*(x^3))+(18*(x^2))+4\n", " return(y) \n", "}" ] }, { "cell_type": "code", "execution_count": 12, "id": "67b3430b", "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFoCAMAAAC46dgSAAAC/VBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKiqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4uLi5\nubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnKysrL\ny8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc3Nzd\n3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u7u7v\n7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7///9osWa/\nAAAACXBIWXMAABJ0AAASdAHeZh94AAAWRklEQVR4nO2dd3wVVfbATxJIIY2EUAKEEliJAVl+\ngoKEgGCkRIoiSHGVKAiRAK6iIvwQUWRFqkqCgIIs6AooLhYWAQXpCgpCQEmCKLDSQyf93c/O\nvJRXkjdvyp2Zl/PO94+ZuS93Zs6737zpcy4wAjVgdgCEvpBg5JBg5JBg5JBg5JBg5JBg5JBg\n5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg\n5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5JBg5FR3wWehnKC/vnRT4cw7\nwY4riwFCXNfdKFS5yJh0JQ8Ej2CBOwrd1C7p27fvdluRBHs+DoJhtpvaxUKd1bYiCfZ8RMFT\niouLiy4tEqa6uKldheBlxWUwizhwSZlg6UoeCAbBU0snuwFESFfOEo3O2Hu9vCwWV8hcUZng\nagciwY8DRIrjwnn9m4fd89Shshp2xUGlG+O95TM7Ci7b+k4HaFkyNz6w6UO/lH6ePTiq7oDP\nHDfRlSrdnHRHSNJR4a+ddP7CSkEkOBHgHmF0qG2pxpqvWJhTUa7gFsOs9QJ+Fj/eGmktDHYW\n7FDpdCtxul4aCeZNheCrS4WpDxi7FSuMm9wVIAyXM+eixEFWvp07AJ/oGsKwl1C8ECpMNAix\n1nEQbF+J9Rf1BosfkmC+OBxFjxU+mArgu4yxPzsAROU6F+UKfugCu3gXQJhQfBYg+HNWNKmy\nYLtKa4XiW5bC8SSYO/aCrb/kOIDh4viwD8A656JMwf7iYdiqUqPCD3iaWPf/nAXbVxoCkCCU\nSuJIMG8cfsFpjBUIW81PrH+5DWCmc1HiNMneXaz4x83CH8+y08LwoFic7SzYrhKLB/iHWJxK\ngnlTvg8uOd5TmPqV/SYM91j/0gNgpHOxCsFVHkWLxS1Wd98Iwwtica2zYLtKlpply1lKgnlj\nO4reK0x9av3JfmotC8e1052LygV/LwytJ1zLJQSzBgBzxeIMEswbm+ATwtT7VpF/E4tHfAHW\nOheVC84Vhq+KxYelBHcD6CkWO5Fg3tgEZwlTbzM2WThs/idjZ4Tj29rnnYui4GW2mWUIZn8R\nPt3ILG9XOoq2r/SGdUElr9BRNHccBc9m7GYzYdyic5AwfI9VKtYGaDrut/KZ5QheLJptHAGS\ngm82ECuFAwnmjk3weWFqtDA+0Lr0kLrmVOuVLMfiI+Kk+0uVYrHMXX630vn7SQkuu9zlP5AE\n88buUqVwpup7WBgXzO7bLLTjyJ/LajgUL4yIDoo7Uj6zHMGMffpI86g+//qPpGB2/IlWdR7Y\n+S4Jxs3LAIPNjsEJEsyDtJYtO95irChevJriWZBgHrwjbKkHfr0pCSD8hNmxOEGCeVAyuOxi\nafBnZofiDAnmw5b+t9dqcM/Ec2bHUQkSjBwSjJzKgt/faUIYhF5UFgypJoRB6EW54FNflgN9\nhIGpMREcKRe8AhwwNSaCI+Uqr6VAyJRZItBRGJgaE8ER2291bWTzHdZPaB+MCbuN8akevpML\nSTAy7Pe2ljn+7TJJMC4cD6cO3B64gASjwul4+VYakGBUVDoh2jJ3sxlxEDpBZ7zIIcHIIcHI\nIcHIIcHIIcHIcRKcP2aVOXEQOuEk+DqMNCcOQidIMHJIMHJIMHJIMHJIMHJIMHJIMHJIMHJI\nMHJIMHKcBN+o83dz4iB0gu4mIYcEI8cAwQf3E5w4qLz19Re8Dwhu7FPc/PoL3gUFuq/DSyiA\nXYrnIcHVCK2Cz1RMHeUQTTkkmBtaBUd+WLacV/21hnItt4KNJJgXWgW3gX5/CqNd8XC3xkhy\nfOyPDJT2+kq4QKvgwpmBESuvjvUJSy/RGkqm7dB+Clx3X5+Qg/aDrOweEAyD/ssrICuLSTAv\ntAs+NxzAfz7f/lNJMDc0C14e6ftsdn9o/xO3kBgJ5ohWwfdCm++F0YeRfs9zi4kEc0SrYP/p\nhdbx2Yd4Xv4gwZqwnLZNaxWcWTG1unJF1ZBgLVhGx9oKfC5V8s42S4I1YJkQsttW4iOYd5Yd\nEqweyzPB39kVtQjWL9ssCVaNZWzIdvuyFsErHO87ypz76uk/3V70IsFqKRoR7ihUi2AV2WYP\nPy722OfXaJj0TpsEqyRvQNR+x0+07YOVZpsd5wPRHZOTOzUGGCVVjwSr41KXJr84faTxIEtZ\nttkM6FV2wStzCMyTqEiCVXE8ru1p58+0HkUryjbbuVVRxXyJCRIVSbAadtTtebXSh9pPkxRk\nmw0bYZueEi5RkQSr4D3/tKLKn3I4D5afbbZznO2mU3f6BXMlf4z/4qo+53KhQ2622Qzoc6h0\n6thwsedtl5BgpfzWoVHVIo19qjIVIKZL/wFdmwOkWCTqkWCFrKmd5KILRIMfmz0wLEo8D44e\ntk2yGglWxLUna7zu6uKR8c9FXz55pspgrowfXUEiCVbA1uYt9rj8oxkPvpdkHanicO/io4Mr\naA/XNK7De7iW5ve0xM+Bg+Ab7/4hc8apy4RB0ZshAAGjr0hVpE20bNbFxH4j9XcOgk/BF3Jn\nvFcYjIeIQWM6QXy+REUSLJPs5JqTpB8hN1xwps/dF4XJZTBNoiIJlsXlFwK6Z7qpY7jgpVD6\nuEHCXRIVSbAM8hdExX7itpZGwesWLlw4A1KFoawZBcHTytylhkpUJMFuKVzaJGq+1G6uDI2C\nuym63S8KXgWlG5UH75CoSILdkPdus/BXKt9ZqAKNgouLiop+h/XCUNaMDWes/aHuUHHyh5pP\nSlQkwZJcmtkgcnquvLqG7oNjSl8g/Jaxl4LqnJSoSIIlODSmVtMFstvHUMHs1qFP3niyy3bG\n4mIkr1WSYFfcXJEAXdfI2l6WYqzgCo5IP3dHgqukZNvIsNrjDiuaxzNzdJDgylj2Tozx6/3R\nLYWzkeBqQf7XaY19E94+476mMyTY88la1D/E//4Mde/Yk2DPJmv5402g6VOfqr67RoI9lqtb\n33yoPjR45N0sLUvRIriYb94GO7xd8KkNs4a28q3RbswKTXJF1AkufKPXczen1g564obW9VeN\n1wouOPbF3FGda0PAnU+8s5NLKil1gl+IfPr+xNs2bGg7kUcMlfE6wYW/71j5+sgezfwguN2Q\nV9ccVXAlww3qBDdew4qi5zO2uxm3QBzwEsHXc3avXzw9tV+HaF/waXjPsCnvf3uK9zrUCQ45\nwljfPYydkHo9QSl/Hq/gdXvBl4+jIXPvlnX/XDhzUuqwPp3i6gcAQEjsXQ888dL8j7f9wnM9\nl22tp05wYmrp8dWs7lqt2shxeNnY7rQgVoccysjRnKNjf8OQHMbOtA91/bimck7a/gMdfsEX\neP5va+PgT99v27Zh/fqPVq7MWLhw5szXX3zxxTFjnho6dGBy8r0JCW3atIxpHB7uV97ONcLr\nxtzeNuG+5KEjnn5x+uyMles27T5sRJwXbK2n8jTp5nbhJ3Zxdg4Hr1Vh9j742snDOzesWTJ/\n5qS0UYP7JyW0vz02OiLE/jcSHBERExvbtn37bklJAwc/Mnp02qRJr82as2TJR2vWb/52/0/H\nf8uVdT9edzwzX7QZgm9lbV/99v+PGpAQV6+GqDCwXuydCfcPfmz0xElvvLlkyZo1Gzdv378/\n5/ip3FwZT8p4DB6UL9oOIwXnH/n3/An92kYC+NS/4/5Hn5mRsXrTvuwLhYYFoC+eky/aHmME\n5+9f/nzv5r4Q3HbAMws+2X2a38mnB+FB+aLt0F1w0Y8ZKW1rQKNezy39TsU9uGqEF+aLztv6\nSvdgaDF0zjeX9FuJx2B8vmhz82QdndcrqEbHF9af12n5HofB+aJNzZNVsmtiC/jL+C+86s1F\nY/NFm5kna/+zjXwT52VzX66HY2i+aPPyZJ2bE++TkI77cKpqDM0XbVaerB1D/Ju+cpznEqsP\nXB7ZOfejvBv/puTJKvrwTt9+/+F5Elet0Cb49xHvMravHYBvPzk3Mk3Ik5WX0SzkGS/98VrR\nJDi7DixgWUE+PVO7QbSMl6EMz5NVkNGozmsy39JCiibBg3zes7CHfcUsaKthnIw5jc2TZflX\nbORMrzonqgJNgut3FAaN+1ink1rLmdXIPFl7Oga9dNl9NeRoEhw6XBjUG2mdfipM5tyu8mTZ\nw0Hw+RTf4XKz/2BGk+DuDa8y1retuLEtadNV9vxV58myR7Ngy/LIO7a7reUNaBL8nX+nPexA\nyORilpcG893PaFierD/uD3oTyw1djWg7Tfq4BsQkxkJUhzBIkTOjVJ6s3LHcUhl+EJbodZck\nXaHxQscfzzW0Pt7S+2tZM0rlyeIm+PLggLlee12jEtqvZF0/9dtZme1pSJ6sH5rFH1Q9Mz4M\nfbvQiDxZiwJSqHt4OwwXrG+erPwnA95TNydWjBWsd56scwmNvlc1I14MTSesd56so83v4vts\nGAIw5cnaEfEg7X6dQZQn6/OgNN1SDlRfTBLsBjWCP6r5Ku8wMGBoOuEyvhzkroYKwctqvKN0\nFq/A0HTCZbzltrJywUv9liqcw0swNJ1wGToIXu63TNkMXoMZ+2D+glfXWKIsBO8BheANNWXc\nq/RSzBB886y7GsoE76k1VVkA3gSCVIbHokbrF0m1p/oLvtDiAZRvbnOi2gvOT7zTK5KmqaXa\nC05pdFrPSKo91V3wgqB9ukZS7dEieCPfdEKuEqFJsbXmKq4x4EOLYOC6cXSZylCC0/We4RkC\nRlQJntPOCsS3a8cxFFfJSF1TmNCFnn52gyrBh1tHvpGeng4z0tN1CInJ3ge/UJcOsNyhbhOd\nN6Hlbt6baHvkCd7gt1GvAPCgdh/8deOphSYLPlPvRb3WjwjVB1kXH+xgrmBLnw4Feq0fERqO\noj8Ypdu783IEpwf/qtfqMVFt0wln1VrEcYV4qa7phIs795LKAUGUU13TCc8Ll3pwnqigmqYT\nzq5FryDJo3qmE7b0uI820PIwPp2wHNwJXh6kV3cg6DA4nbBM3Ai+UGcWz7Whxth0wnJxIzil\nDd1jkIuh6YRlIy14py9lSJKNoemEZSMpuLjdCNd/JJzg88jO+9IZ+hUjKTgjzO1z1UQFfARD\nKo9YbEgJvhQplSuecEKL4FNflgN9hAHHqKQEj29FN5EUoEXwCsc+TTlGJSH415rcXzdHjRbB\n11IgZMosEegoDDRGcvHRwRW0d/3QXf/7NK7Hy9C2D14b2XyH9RMO++Ar4+WkMtzqe0D7qrwJ\njQdZp3r4Ti5UJlhTz2eWux+TvyaCaT+Ktszxb5cpX7DWns8+Cfxd5pqIUrSfJh24PXCBXMFa\nez4rjntW3oqIcjicB99KA5mCNfd8tjzUazqV5AWXCx1b5m6WNaPWns8Kmr4saz2EDUPfLtTa\n89miSMmOAIgqMFSwxp7P8hu/rnbN3ouhgjX2fJZRx9t7uVKBsS+Aa+r5rCCGfsDKMfgNfy09\nny2JoD2wcoxP4aC257Oi2Gma1uulcBCcP0ZRGgW1lypXhlxSshqiFA6Cr8NI2bOqv1RpaT1R\n9loIG8YK1nCp8vMA6o5BDYYK1nKpsov8zQRhh6GCNVyq/N7niMyVEA4YKljDpcqH+8lcB+GI\nsb9g1Zcqc/ykT5wJVxi8D1Z7qXJCB7nBEY4YexQteakyc38FU5wEXwn9SHGUhBWDz4MlLlXm\n+Ng/g+vYhdm8hvS2mUoMFswkLlVey61gIzg83F4SO0NhiEQ5Rgs+92vZmdIFqSRbuxwF/zuA\n3kZSC4/eR+v8Xe6cB9oCNFhhnewtdcvCSXASvU6oGkPvJuUE+iYlB0KGOK1AcJbPD2rXSBgq\neKjPBsbOtwwUU9RJCt4HBDeUp8RXLbh5L3F4LEi8KiUpmB20nTLtD5y4ihsTA/gta1XXrhwX\nFsDzWwbaNd9B5Z5UCw4tvYX0Mmx3J9ieYI4vpn4ZzG9ZLCWF48I86VuqFtwl3jq6EdO6gAQ7\n4UnfUrXgyTDO2o/HVzA0jwQ74knfskxN3gml2c/yEiG0rzjxMjSqS4Id8KRvCWxoNrv+mC8E\nTVbYqdzll+JKt9IrWsnPCOBJX90BxIJhL3u65VcnN8W9pnYZlhNb5Fb1pK/uAG7B0d8IU9ta\ncAlIGk/66g7gFhzxizB1qLaq9Q9SVN2TvroDmAWPXZf8HGMFTyWpmd99D+AOeNJXdwCx4BcG\nxAfADdYz+JCa+UlwFXjSt7T6KTlRwtaqS6uvUHDEJlVrqZJNEfyWxUbz7Hnck76l1l5XFAo+\nwbGzgJIT/JbFcnl2K+RJ31Jrrys36fa9Z6N/ryuEqejf6wphKvr3ukKYiv69rhCmon+vK4Sp\n6N/rCmEq+ve6QpiK/r2uEKZS+beaS+dImLAJzlv4xMws9llDCBlA50l4qBB8OR4A6v8YENa9\nDdTXrb93wmgqBD8Pzx3a3DK4ifDr/RgoyxEaKgTHd2LiM7D/EKfvbWdaPARnKgQHiYneT8Ma\ncfrpWqbFQ3CmQnCs2IfRrVTr2y8PR5kWj0KyF5odQVV4UlQVgofU/Lx8MicoWbf15U9JDIsd\nxq3P7wmqnhSszKKE8IRFfBbF+EXFo7kqBB+v5dPe2tHc4QnhPlu1BuaKK4kQP6qnTxCnHrE2\nBfBpylRo9fhtMI7LsvhFxaW5bOfB2QPrp4vjxVB/jebIXDEZ0oThV75/5bGwR1sBcGnKA9C7\niBX19DnMY2HcomJcmsvhSpb1GlbOLh2z4MSFWt9YS4JzHBb2UN++oVyachj8LAx/hMd5LIxb\nVIxLcxl8WyHe+r4aS4Zf+SyvDZemjGpsHUU34LEwxisqxqW5TLlvdD6wvsI33VzBpSkvQ2kq\nxo6u+0lVBjfBpWhqLjMEH2sJH3BaFJemPAn9reNkkEoHpQC+grU1l0GCb74lUPq8/41pQYHp\nnJbFpynPwADrOBn+5LA0xlew1uYySPBZMUWM9T21DU2gr7YdsG1ZnJqyxK+rddzJj9OdUo6C\nNTeX0ZvoadD6O46L49OU0bHWUUwjHgtjPAVrby6DBa+AoQXua8mHT1MOg2PCMBOG8VgY4yiY\nQ3MZK9jSqlEe1wXyacpt8DchtiGwg8fCGD/BPJrLWMEnoG7vUi7wWSCnpkyBHlO6Kku0KwUv\nwTyay1jB31Tk5POsExLLm53DOs/hsigRXoJ5NBc9IIscEowcEowcEowcEowcEowcEowcEowc\nEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEowc\nEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEowcEsxYpv+9wrCwTeQZsyPRARIsMB2W\nMzYTPjQ7Dj0gwQIFbSLPZwf2MzsMXSDBInt9h/eI4JRL2MMgwVaeBVhpdgz6QIKtZEPwVbNj\n0AcSbKW/P4w1OwZ9IMEiH8L8QT67zI5CF0iwwNk67Yv/GxbPNVW5p0CCBQb6/chYOrxqdhx6\nQIIZW23tELvk7oCjZkeiAyQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQY\nOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQYOSQY\nOf8DPLA9ERlYiSQAAAAASUVORK5CYII=", "text/plain": [ "Plot with title “Root Finding”" ] }, "metadata": { "image/png": { "height": 450, "width": 600 } }, "output_type": "display_data" } ], "source": [ "gridx<-seq(-2,2,0.01)\n", "gridy<-gx(gridx)\n", "plot(gridx,gridy,\"l\",main='Root Finding', ylim= c(-200, 100), \n", " xlab = 'x', ylab=expression(y=3*x^5-4*x^4+6*x^3+4*x-4))\n", "lines(c(-2,2),c(0,0))" ] }, { "cell_type": "markdown", "id": "b24f9dca", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Finding the Tangent\n", "\n", "- To find the tangent evaluate the first derivative of $g(x)$.\n", "\n", "- The function is $$g(x)=3x^5-4x^4+6x^3+4x-4$$\n", "\n", "- The first derivative is $$g'(x)=15x^4-16x^3+18x^2+4$$\n" ] }, { "cell_type": "markdown", "id": "e189bffd", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "![rootfinding](figures/rf7.png)" ] }, { "cell_type": "markdown", "id": "a6158085", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Find the crossing point\n", "\n", "- Find the crossing point From basic Geometry\n", "\n", "$$g'(x_0)=\\frac{g(x_0)}{x_0-x_1}$$\n", "\n", "\n", "- Rearrange\n", "\n", "$$\\begin{aligned}\n", " x_0-x_1&=\\frac{g(x_0)}{g'(x_0)}\\\\\n", " -x_1&=-x_0+\\frac{g(x_0)}{g'(x_0)}\\\\\n", " x_1&=x_0-\\frac{g(x_0)}{g'(x_0)}\n", " \\end{aligned}$$" ] }, { "cell_type": "markdown", "id": "8b5b2e43", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Stopping Rule\n", "\n", "- With each step the algorithm should get closer to the root.\n", "\n", "- However, it can run for a long time without reaching the **exact**\n", " root\n", "\n", "- There must be a **stopping rule** otherwise the program could run\n", " forever.\n", "\n", "- Let $\\epsilon$ be an extremely small number e.g. $1\\times 10^{-10}$\n", " called the **tolerance level**\n", "\n", "- If $|g(x^*)|<\\epsilon$ then the solution is close enough and there\n", " is a root at $x^*$" ] }, { "cell_type": "markdown", "id": "644826ac", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Newton-Raphson algorithm for root finding\n", "\n", "1. Select initial value $x_0$ and set $n=0$\n", "\n", "2. Set $x_{n+1}=x_{n}-\\frac{g(x_n)}{g'(x_n)}$\n", "\n", "3. Evaluate $|g(x_{n+1})|$\n", "\n", " - If $|g(x_{n+1})|\\leq\\epsilon$ then stop.\n", "\n", " - Otherwise set $n=n+1$ and go back to step 2." ] }, { "cell_type": "markdown", "id": "2e472147", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Your task: Write R/Python code to find the root\n", "\n", "$g(x)=3x^5-4x^4+6x^3+4x-4$\n", "\n", "Tips:\n", "\n", "- Write functions for $g(x)$ and $g'(x)$ first.\n", "\n", "- These can be inputs into a function that carries out the Newton\n", " Raphson method. Code should be flexible.\n", "\n", "- Use loops!" ] }, { "cell_type": "markdown", "id": "d22e51a2", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Your task: a different target function\n", "\n", "- Now use your Newton-Raphson code to find the root of\n", " $g(x)=\\sqrt{|x|}$\n", "\n", "- The derivative has two parts $$g'(x)=\\left\\{\n", " \\begin{array}{l}\n", " 1/(2\\sqrt{x})\\quad\\mbox{if $x>0$}\\\\\n", " -1/(2\\sqrt{-x})\\quad\\mbox{if $x<0$}\n", " \\end{array}\n", " \\right.$$\n", "\n", "- Use 0.25 as the starting value." ] }, { "cell_type": "markdown", "id": "aff0ae58", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Learn from mistakes\n", "\n", "- Newton-Raphson does not always converge\n", "\n", "- Be careful using *while*. Avoid infinite loops.\n", "\n", "- Don’t always assume the answer given by code is correct. Check\n", " carefully!\n", "\n", "- Print warning messages in code" ] }, { "cell_type": "markdown", "id": "ccde220c", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- For some functions, using some certain starting values leads to a\n", " series that **converges**, while other starting values lead to a\n", " series that **diverges**\n", "\n", "- For other functions different starting values converge to different\n", " roots.\n", "\n", "- Be careful when choosing the initial value.\n", "\n", "- Newton-Raphson doesn’t work if the first derivative is zero.\n", "\n", "- When can this happen?" ] }, { "cell_type": "markdown", "id": "b988bb7d", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Rough Proof of Quadratic Convergence\n", "\n", "- Can we prove anything about the rate of convergence for the Newton\n", " Raphson Method?\n", "\n", "- To do so requires the **Taylor Series**\n", "\n", "- Let $f(x)$ have a root at $\\alpha$. The Taylor approximation states\n", " that\n", "\n", " $$f(\\alpha)\\approx f(x_n)+f'(x_n)(\\alpha-x_n)+\\frac{1}{2}f''(x_n)(\\alpha-x_n)^2$$\n", "\n", "- The quality of the approximation depends on the function and how\n", " close $x_n$ is to $\\alpha$\n" ] }, { "cell_type": "markdown", "id": "d6db2a23", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- Since $\\alpha$ is a root, $f(\\alpha)=0$ This implies\n", "\n", " $$0\\approx f(x_n)+f'(x_n)(\\alpha-x_n)+\\frac{1}{2}f''(x_n)(\\alpha-x_n)^2$$\n", "\n", "- Dividing by $f'(x_n)$ and rearranging gives:\n", " $$\\frac{f(x_n)}{f'(x_n)}+(\\alpha-x_n)\\approx\\frac{-f''(x_n)}{2f'(x_n)}(\\alpha-x_n)^2$$\n", "\n", "- More rearranging\n", " $$\\alpha-\\left(x_n-\\frac{f(x_n)}{f'(x_n)}\\right)\\approx\\frac{-f''(x_n)}{2f'(x_n)}(\\alpha-x_n)^2$$\n", "\n", "- The term in brackets on the left hand side is the formula used to\n", " update $x$ in the Newton Raphson method\n", " $$(\\alpha-x_{n+1})\\approx\\frac{-f''(x_n)}{2f'(x_n)}(\\alpha-x_n)^2$$\n", "\n", "- This can be rewritten in terms of errors $e_{n+1}=\\alpha-x_{n+1}$\n", " and $e_n=\\alpha-x_n$\n", " $$e_{n+1}\\approx\\frac{-f''(x_n)}{2f'(x_n)}e_n^2$$" ] }, { "cell_type": "markdown", "id": "fe95d8ac", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Conclusion\n", "\n", "- Why did we spend so much time on finding roots of an equation?\n", "\n", "- Isn’t this topic meant to be about optimization?\n", "\n", "- Can we change this algorithm slightly so that it works for\n", " optimization?" ] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "R", "language": "R", "name": "ir" }, "language_info": { "codemirror_mode": "r", "file_extension": ".r", "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", "version": "4.3.2" } }, "nbformat": 4, "nbformat_minor": 5 }