sesh6.ipynb (7849B)
1 { 2 "cells": [ 3 { 4 "cell_type": "markdown", 5 "id": "bcc85c40", 6 "metadata": {}, 7 "source": [ 8 "# MILUTIN POPOVIC: EXERCISE SHEET 6" 9 ] 10 }, 11 { 12 "cell_type": "code", 13 "execution_count": 1, 14 "id": "77f12d86", 15 "metadata": {}, 16 "outputs": [], 17 "source": [ 18 "import numpy as np\n", 19 "import numdifftools as nd" 20 ] 21 }, 22 { 23 "cell_type": "markdown", 24 "id": "921657fc", 25 "metadata": {}, 26 "source": [ 27 "Exercise 41:\n", 28 "-----------------\n", 29 "Implement the local Newton algorithm. Use as input data the starting vector x0,\n", 30 "the parameter for the stopping criterion ε and the parameter for the maximal\n", 31 "number of allowed iterations kmax. The sequence $x_0$, $x_1$, $x_2$, ...\n", 32 "containing the iteration history and the number of performed iterations should be returned.\n", 33 "The implemented algorithm should be tested for ε = 10^(−6) and kmax = 200, and the following\n", 34 "\n", 35 "Exercise 42:\n", 36 "-----------------\n", 37 "Implement the globalized Newton algorithm. Use as input data the starting vector $x_0$ the\n", 38 "parameter for the stopping criterion ε, the parameter for the maximal number of allowed\n", 39 "iterations kmax, the parameters for the determination of the Armijo step size σ and β, and\n", 40 "the parameters ρ and p. The sequence $x_0$, $x_1$, $x_2$, ... containing the iteration history and the\n", 41 "number of performed iterations should be returned.\n", 42 "The implemented algorithm should be tested for ε=10−6, kmax=200, ρ=10−8, p=2.1, σ=10−4" 43 ] 44 }, 45 { 46 "cell_type": "code", 47 "execution_count": 2, 48 "id": "bc388246", 49 "metadata": {}, 50 "outputs": [], 51 "source": [ 52 "def local_newton(f, x_0, eps=10e-6, kmax=200):\n", 53 " f_grad = nd.Gradient(f)\n", 54 " f_hess = nd.Hessian(f)\n", 55 " x_k = x_0\n", 56 " \n", 57 " for k in range(kmax):\n", 58 " if np.linalg.norm(f_grad(x_k)) <= eps:\n", 59 " return x_k\n", 60 " break\n", 61 " \n", 62 " if type(x_0) == float:\n", 63 " d_k = -float(f_hess(x_k))/float(f_grad(x_k))\n", 64 " else:\n", 65 " d_k = np.linalg.solve(f_hess(x_k), -f_grad(x_k))\n", 66 " \n", 67 " x_k = x_k + np.array(d_k)\n", 68 " return x_k\n", 69 "\n", 70 "def global_newton(f, x_0, eps=10e-6, kmax=200, rho=10e-8, p=2.1, sig=10e-4, beta=0.5):\n", 71 " f_grad = nd.Gradient(f)\n", 72 " f_hess = nd.Hessian(f)\n", 73 " x_k = x_0\n", 74 " for k in range(kmax):\n", 75 " if np.linalg.norm(f_grad(x_k)) <= eps:\n", 76 " return x_k\n", 77 " break\n", 78 " \n", 79 " try:\n", 80 " if type(x_0) == float:\n", 81 " d_k = -float(f_hess(x_k))/float(f_grad(x_k))\n", 82 " else:\n", 83 " d_k = np.linalg.solve(f_hess(x_k), -f_grad(x_k))\n", 84 " if (np.dot(f_grad(x_k), d_k) > -rho*np.linalg.norm(d_k)**p):\n", 85 " d_k = -f_grad(x_k)\n", 86 " except np.linalg.LinAlgError:\n", 87 " d_k = -f_grad(x_k)\n", 88 " \n", 89 " # Find step size (Arnijno Rule)\n", 90 " t_k = 1\n", 91 " while f(x_k + t_k*d_k) > f(x_k) + sig*t_k * np.dot(f_grad(x_k), d_k):\n", 92 " t_k = t_k * beta\n", 93 " \n", 94 " x_k = x_k + t_k * d_k\n", 95 " \n", 96 " return x_k" 97 ] 98 }, 99 { 100 "cell_type": "code", 101 "execution_count": 3, 102 "id": "58828079", 103 "metadata": {}, 104 "outputs": [ 105 { 106 "name": "stdout", 107 "output_type": "stream", 108 "text": [ 109 "a)\n", 110 "Local Newton Algorithm for x_0 = [-1.2, 1] is: [0.9999957 0.99999139]\n", 111 "Global Newton Algorithm for x_0 = [-1.2, 1] is: [1. 1.]\n", 112 "\n", 113 "b)\n", 114 "Local Newton Algorithm for x_0 = [0.25, 0.25, 0.25, 0.25] is: [1.50383723 0.7951105 0.47235406 0.32296666]\n", 115 "Global Newton Algorithm for x_0 = [0.25, 0.25, 0.25, 0.25] is: [1.50383723 0.7951105 0.47235406 0.32296666]\n", 116 "\n", 117 "c)\n", 118 "Local Newton Algorithm for x_0 = [1, 1] is: [158.75270666 79.36690168]\n", 119 "Global Newton Algorithm for x_0 = [1, 1] is: [1.24996948e-06 2.00000000e+06]\n", 120 "\n", 121 "d)\n", 122 "Local Newton Algorithm for x_0 = [-3, -1, -3, -1] is: [-1.41934291 2.02305705 0.36974654 0.12724277]\n", 123 "Global Newton Algorithm for x_0 = [-3, -1, -3, -1] is: [-1.41934305 2.02305745 0.36974706 0.12724316]\n", 124 "\n", 125 "e)\n", 126 "Local Newton Algorithm for x_0 = 2.0 is: -10.69705641570642\n", 127 "Local Newton Algorithm for x_0 = 1.0 is: 21.27903082235975\n", 128 "Local Newton Algorithm for x_0 = 0.5 is: 21.26959143622489\n", 129 "Global Newton Algorithm for x_0 = 2.0 is: -4.437802149414097e-11\n", 130 "Global Newton Algorithm for x_0 = 1.0 is: 8.446339466597688e-14\n", 131 "Global Newton Algorithm for x_0 = 0.5 is: 7.344043315106116e-14\n", 132 "\n" 133 ] 134 } 135 ], 136 "source": [ 137 "f_a = lambda x: (1-x[0])**2 + 100*(x[1] - x[0]**2)**2\n", 138 "x_a_0 = [-1.2, 1]\n", 139 "print('a)')\n", 140 "print(f'Local Newton Algorithm for x_0 = {x_a_0} is: {local_newton(f_a, x_a_0)}')\n", 141 "print(f'Global Newton Algorithm for x_0 = {x_a_0} is: {global_newton(f_a, x_a_0)}')\n", 142 "print('')\n", 143 "\n", 144 "f_b = lambda x: sum( (sum(np.cos(x[j]) + (i+1)*(1 - np.cos(x[i])) - np.sin(x[i]) for j in range(4) ) )**2 for i in range(4) )\n", 145 "x_b_0 = [1/4, 1/4, 1/4, 1/4]\n", 146 "print('b)')\n", 147 "print(f'Local Newton Algorithm for x_0 = {x_b_0} is: {local_newton(f_b, x_b_0)}')\n", 148 "print(f'Global Newton Algorithm for x_0 = {x_b_0} is: {global_newton(f_b, x_b_0)}')\n", 149 "print('')\n", 150 "\n", 151 "f_c = lambda x: (x[0] - 10**6)**2 + (x[1] - 2*10**6)**2 + (x[0]*x[1] - 2)**2\n", 152 "x_c_0 = [1, 1]\n", 153 "print('c)')\n", 154 "print(f'Local Newton Algorithm for x_0 = {x_c_0} is: {local_newton(f_c, x_c_0)}')\n", 155 "print(f'Global Newton Algorithm for x_0 = {x_c_0} is: {global_newton(f_c, x_c_0)}')\n", 156 "print('')\n", 157 "\n", 158 "f_d = lambda x: 100*(x[1] - x[0]**2)**2 + (1-x[0])**2 + 90*(x[3] - x[2]**2)**2 + (1 - x[2])**2 + 10*(x[1] - x[3] - 2)**2 + 1/10*(x[1] - x[3])**2\n", 159 "x_d_0 = [-3, -1, -3, -1]\n", 160 "print('d)')\n", 161 "print(f'Local Newton Algorithm for x_0 = {x_d_0} is: {local_newton(f_d, x_d_0)}')\n", 162 "print(f'Global Newton Algorithm for x_0 = {x_d_0} is: {global_newton(f_d, x_d_0)}')\n", 163 "print('')\n", 164 "\n", 165 "f_e = lambda x: np.sqrt(1 + x**2)\n", 166 "x_e_00 = 2.0\n", 167 "x_e_01 = 1.0\n", 168 "x_e_02 = 1/2\n", 169 "print('e)')\n", 170 "print(f'Local Newton Algorithm for x_0 = {x_e_00} is: {local_newton(f_e, x_e_00)}')\n", 171 "print(f'Local Newton Algorithm for x_0 = {x_e_01} is: {local_newton(f_e, x_e_01)}')\n", 172 "print(f'Local Newton Algorithm for x_0 = {x_e_02} is: {local_newton(f_e, x_e_02)}')\n", 173 "print(f'Global Newton Algorithm for x_0 = {x_e_00} is: {global_newton(f_e, x_e_00)}')\n", 174 "print(f'Global Newton Algorithm for x_0 = {x_e_01} is: {global_newton(f_e, x_e_01)}')\n", 175 "print(f'Global Newton Algorithm for x_0 = {x_e_02} is: {global_newton(f_e, x_e_02)}')\n", 176 "print('')" 177 ] 178 }, 179 { 180 "cell_type": "code", 181 "execution_count": null, 182 "id": "b1f45770", 183 "metadata": {}, 184 "outputs": [], 185 "source": [] 186 } 187 ], 188 "metadata": { 189 "kernelspec": { 190 "display_name": "Python 3", 191 "language": "python", 192 "name": "python3" 193 }, 194 "language_info": { 195 "codemirror_mode": { 196 "name": "ipython", 197 "version": 3 198 }, 199 "file_extension": ".py", 200 "mimetype": "text/x-python", 201 "name": "python", 202 "nbconvert_exporter": "python", 203 "pygments_lexer": "ipython3", 204 "version": "3.10.9" 205 } 206 }, 207 "nbformat": 4, 208 "nbformat_minor": 5 209 }