{"id":1774,"date":"2023-01-16T01:01:44","date_gmt":"2023-01-16T09:01:44","guid":{"rendered":"https:\/\/gantovnik.com\/bio-tips\/?p=1774"},"modified":"2023-01-16T01:03:17","modified_gmt":"2023-01-16T09:03:17","slug":"338-minimization-using-newtons-method-on-the-non-quadratic-problem-in-python","status":"publish","type":"post","link":"https:\/\/gantovnik.com\/bio-tips\/2023\/01\/338-minimization-using-newtons-method-on-the-non-quadratic-problem-in-python\/","title":{"rendered":"#338 Minimization using Newton\u2019s method on the non-quadratic problem in python"},"content":{"rendered":"<p><a href=\"https:\/\/gantovnik.com\/bio-tips\/2023\/01\/338-minimization-using-newtons-method-on-the-non-quadratic-problem-in-python\/ex338\/\" rel=\"attachment wp-att-1776\"><img data-recalc-dims=\"1\" decoding=\"async\" src=\"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/01\/ex338.png?resize=1080%2C540&#038;ssl=1\" alt=\"\" width=\"1080\" height=\"540\" class=\"alignnone size-full wp-image-1776\" srcset=\"https:\/\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/01\/ex338.png 1600w, https:\/\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/01\/ex338-1280x640.png 1280w, https:\/\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/01\/ex338-980x490.png 980w, https:\/\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/01\/ex338-480x240.png 480w\" sizes=\"(min-width: 0px) and (max-width: 480px) 480px, (min-width: 481px) and (max-width: 980px) 980px, (min-width: 981px) and (max-width: 1280px) 1280px, (min-width: 1281px) 1600px, 100vw\" \/><\/a><\/p>\n<p>ex338.py<\/p>\n<pre class=\"brush: python; title: ; notranslate\" title=\"\">\r\nimport numpy as np\r\nimport math\r\nimport matplotlib.pyplot as plt\r\nimport warnings\r\nwarnings.filterwarnings('ignore')\r\n\r\ndef my_f(x):\r\n    #$f(x_1, x_2) = e^{x_1+3x_2-0.1}+e^{x_1-3x_2-0.1}+e^{-x_1-0.1}$\r\n    x1 = x&#x5B;0, 0]\r\n    x2 = x&#x5B;1, 0]\r\n    return math.exp(x1+3*x2-0.1)+math.exp(x1-3*x2-0.1)+math.exp(-x1-0.1)\r\n\r\ndef my_gradient_f(x):\r\n    #$\\nabla f(x_1, x_2)$\r\n    x1 = x&#x5B;0, 0]\r\n    x2 = x&#x5B;1, 0]\r\n    gradient_1=1*math.exp(x1+3*x2-0.1)+1*math.exp(x1-3*x2-0.1)-math.exp(-x1-0.1)\r\n    gradient_2=3*math.exp(x1+3*x2-0.1)-3*math.exp(x1-3*x2-0.1)\r\n    return np.array(&#x5B;&#x5B;gradient_1], &#x5B;gradient_2]])\r\n\r\ndef my_hessian_f(x):\r\n    #$\\nabla^2 f(x_1, x_2)$\r\n    x1 = x&#x5B;0, 0]\r\n    x2 = x&#x5B;1, 0]\r\n    x11 = math.exp(x1 + 3 * x2 - 0.1) + math.exp(x1 - 3 * x2 - 0.1) + math.exp(-x1 - 0.1)\r\n    x12 = 3 * math.exp(x1 + 3 * x2 - 0.1) - 3 * math.exp(x1 - 3 * x2 - 0.1)\r\n    x21 = 3 * math.exp(x1 + 3 * x2 - 0.1) - 3 * math.exp(x1 - 3 * x2 - 0.1)\r\n    x22 = 9 * math.exp(x1 + 3 * x2 - 0.1) + 9 * math.exp(x1 - 3 * x2 - 0.1)\r\n    return np.array(&#x5B;&#x5B;x11, x12], &#x5B;x21, x22]])\r\n\r\ndef my_newton_step_decrement(x):\r\n    # calculate newton step $\\Delta x_{nt} = -&#x5B;\\nabla^2 f(x)]^{-1}\\nabla f(x) $\r\n    gradient = my_gradient_f(x)\r\n    hessian = my_hessian_f(x)\r\n    hessian_inv = np.linalg.inv(hessian)\r\n    x_nt = -1 * hessian_inv @ gradient\r\n    decrement = -1 * gradient.T @ x_nt\r\n    return x_nt, decrement.item()\r\n\r\ndef backtracking_line_search(f, df, x, delta_x, alpha, beta):\r\n    # backtracking line search\r\n    t = 1\r\n    while f(x + delta_x * t) &gt; f(x) + alpha * t * ((df(x).T @ delta_x).item()):\r\n        t = beta * t\r\n    return t\r\n\r\n\r\ndef plot_ellipse_by_matrix(P, N):\r\n    # Plot a 2D ellipse according to the $X^TPX=1$\r\n    W, V = np.linalg.eig(P)\r\n    a = 1 \/ np.sqrt(W&#x5B;0])\r\n    b = 1 \/ np.sqrt(W&#x5B;1])\r\n    n_p = N\r\n    y = np.zeros((n_p + 1, 2, 1), dtype='complex_')\r\n    for i in range(n_p + 1):\r\n        omega = i * 2 * np.pi \/ n_p\r\n        y&#x5B;i, 0, 0] = a * np.cos(omega)\r\n        y&#x5B;i, 1, 0] = b * np.sin(omega)\r\n\r\n    x = np.zeros((n_p + 1, 2, 1), dtype='complex_')\r\n    for i in range(n_p+1):\r\n        x&#x5B;i] = V @ y&#x5B;i]\r\n\r\n    return x\r\n        \r\ndef plot_convergence_Newton_method_non_quadratic_problem(ax1, ax2):\r\n    # Plot the convergence of the Newton's method on the non-quadratic problem\r\n    min_x1 = -6\r\n    max_x1 = 5\r\n    min_x2 = -4\r\n    max_x2 = 4\r\n    delta = 0.1\r\n    x = np.arange(min_x1, max_x1, delta)\r\n    y = np.arange(min_x2, max_x2, delta)\r\n    X, Y = np.meshgrid(x, y)\r\n    Z = np.zeros(X.shape)\r\n    f_min = my_f(np.array(&#x5B;&#x5B;0.5 * math.log(0.5)], &#x5B;0]]))\r\n    for x1_i in range(len(x)):\r\n        for x2_i in range(len(y)):\r\n            x1 = x&#x5B;x1_i]\r\n            x2 = y&#x5B;x2_i]\r\n            Z&#x5B;x2_i]&#x5B;x1_i] = my_f(np.array(&#x5B;&#x5B;x1], &#x5B;x2]]))\r\n\r\n    levels = np.arange(0, 200, 1)\r\n    CS = ax1.contour(X, Y, Z, linestyles='dashed', levels=levels, vmin=0,vmax=100)\r\n    ax1.clabel(CS, fontsize=4, inline=True)\r\n    alpha = 0.1\r\n    beta = 0.7\r\n    x0 = np.array(&#x5B;&#x5B;-2], &#x5B;2]])\r\n    x = x0\r\n    epsilon = 1e-10\r\n    x_nt, decrement = my_newton_step_decrement(x)\r\n    x_points = &#x5B;x]\r\n    fx = &#x5B;my_f(x)]\r\n    iter = 0\r\n    while decrement\/2 &gt; epsilon:\r\n        if iter &lt; 4:\r\n            hessian = my_hessian_f(x)\r\n            v = plot_ellipse_by_matrix(hessian, 100)\r\n            if iter == 0:\r\n                ax1.plot(v&#x5B;:,0, 0] + x&#x5B;0, 0], v&#x5B;:, 1, 0] + x&#x5B;1, 0],color='red',\r\n                         label=r'ellipsoid {$x^{(k)}$+x||$x^T\\nabla^2f(x) x$|=1}')\r\n            else:\r\n                ax1.plot(v&#x5B;:, 0, 0] + x&#x5B;0, 0], v&#x5B;:, 1, 0] + x&#x5B;1, 0],color='red')\r\n        print('Info: decrement = {}, fx={}'.format(decrement, fx&#x5B;-1]))\r\n        t = backtracking_line_search(my_f, my_gradient_f, x, x_nt, alpha, beta)\r\n        x_next = x + t * x_nt\r\n        x = x_next\r\n        x_nt, decrement = my_newton_step_decrement(x)\r\n        x_points.append(x)\r\n        fx.append(my_f(x))\r\n        iter = iter + 1\r\n    x_points_arr = np.array(x_points)\r\n    ax1.plot(x_points_arr&#x5B;:, 0], x_points_arr&#x5B;:, 1], color='black',marker='o',\r\n             label='Newton descent')\r\n    ax1.set_aspect('auto', adjustable='box')\r\n    ax1.legend()\r\n    ax1.set_xlabel(r'$x_1$')\r\n    ax1.set_ylabel(r'$x_2$')\r\n    ax1.set_title('f(x)')\r\n    ax1.grid()\r\n    ax2.semilogy(np.array(fx) - f_min, marker='o', label=r&quot;Newton's method&quot;)\r\n    ax2.grid()\r\n    ax2.legend()\r\n    ax2.set_xlabel('iteration')\r\n    ax2.set_title(r'$f(x)-p^*$')\r\n\r\nif __name__ == &quot;__main__&quot;:\r\n    plt.figure(figsize=(16, 8))\r\n    ax1 = plt.subplot(1, 2, 1)\r\n    ax2 = plt.subplot(1, 2, 2)\r\n    plot_convergence_Newton_method_non_quadratic_problem(ax1, ax2)\r\n    plt.savefig('ex338.png', dpi=100)\r\n    plt.show()\r\n<\/pre>\n","protected":false},"excerpt":{"rendered":"<p>ex338.py import numpy as np import math import matplotlib.pyplot as plt import warnings warnings.filterwarnings(&#8216;ignore&#8217;) def my_f(x): #$f(x_1, x_2) = e^{x_1+3x_2-0.1}+e^{x_1-3x_2-0.1}+e^{-x_1-0.1}$ x1 = x&#x5B;0, 0] x2 = x&#x5B;1, 0] return math.exp(x1+3*x2-0.1)+math.exp(x1-3*x2-0.1)+math.exp(-x1-0.1) def my_gradient_f(x): #$\\nabla f(x_1, x_2)$ x1 = x&#x5B;0, 0] x2 = x&#x5B;1, 0] gradient_1=1*math.exp(x1+3*x2-0.1)+1*math.exp(x1-3*x2-0.1)-math.exp(-x1-0.1) gradient_2=3*math.exp(x1+3*x2-0.1)-3*math.exp(x1-3*x2-0.1) return np.array(&#x5B;&#x5B;gradient_1], &#x5B;gradient_2]]) def my_hessian_f(x): #$\\nabla^2 f(x_1, x_2)$ x1 = [&hellip;]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"nf_dc_page":"","_et_pb_use_builder":"","_et_pb_old_content":"","_et_gb_content_width":"","_lmt_disableupdate":"yes","_lmt_disable":"","jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[9,2],"tags":[75,25,3],"class_list":["post-1774","post","type-post","status-publish","format-standard","hentry","category-optimization","category-python","tag-newtons-method","tag-optimization","tag-python"],"modified_by":"gantovnik","jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"jetpack_shortlink":"https:\/\/wp.me\/p8bH0k-sC","jetpack_likes_enabled":true,"jetpack-related-posts":[{"id":104,"url":"https:\/\/gantovnik.com\/bio-tips\/2018\/12\/unconstrained-multivariate-optimization\/","url_meta":{"origin":1774,"position":0},"title":"Unconstrained multivariate optimization","author":"gantovnik","date":"2018-12-31","format":false,"excerpt":"import os import matplotlib.pyplot as plt import numpy as np import scipy import sympy os.chdir('\/home\/vg\/Downloads\/projects\/ex18') os.getcwd() x1,x2=sympy.symbols(\"x_1,x_2\") f_sym=(x1-1)**4 + 5*(x2-1)**2 - 2*x1*x2 fprime_sym=[f_sym.diff(x_) for x_ in (x1,x2)] sympy.Matrix(fprime_sym) fhess_sym=[[f_sym.diff(x1_,x2_) for x1_ in (x1,x2)] for x2_ in (x1,x2)] sympy.Matrix(fhess_sym) f_lmbda=sympy.lambdify((x1,x2),f_sym,'numpy') fprime_lmbda=sympy.lambdify((x1,x2),fprime_sym,'numpy') fhess_lmbda=sympy.lambdify((x1,x2),fhess_sym,'numpy') def func_XY_to_xy(f): return lambda X: np.array(f(X[0],X[1])) f=func_XY_to_xy(f_lmbda) fprime=func_XY_to_xy(fprime_lmbda) fhess=func_XY_to_xy(fhess_lmbda)\u2026","rel":"","context":"In &quot;python&quot;","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"example18","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2018\/12\/example18.png?resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2018\/12\/example18.png?resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2018\/12\/example18.png?resize=525%2C300 1.5x"},"classes":[]},{"id":85,"url":"https:\/\/gantovnik.com\/bio-tips\/2018\/12\/subplots-in-matplotlib\/","url_meta":{"origin":1774,"position":1},"title":"#13 Subplots in matplotlib","author":"gantovnik","date":"2018-12-29","format":false,"excerpt":"[code language=\"python\"] import os import matplotlib.pyplot as plt import numpy as np os.chdir('\/home\/vg\/Downloads\/projects\/ex13') os.getcwd() fig,axes = plt.subplots(2,2,figsize=(6,6),sharex=True,sharey=True,squeeze=False) x1=np.random.randn(100) x2=np.random.randn(100) axes[0,0].set_title(\"Uncorrelated\") axes[0,0].scatter(x1,x2) axes[0,1].set_title(\"Weakly positively correlated\") axes[0,1].scatter(x1,x1+x2) axes[1,0].set_title(\"Weakly negatively correlated\") axes[1,0].scatter(x1,-x1+x2) axes[1,1].set_title(\"Strongly correlated\") axes[1,1].scatter(x1,x1+0.15*x2) axes[1,1].set_xlabel(\"x\") axes[1,0].set_xlabel(\"x\") axes[0,0].set_ylabel(\"y\") axes[1,0].set_ylabel(\"y\") plt.subplots_adjust(left=0.1,right=0.95,bottom=0.1,top=0.95,wspace=0.1,hspace=0.2) plt.savefig(\"example13.png\", dpi=100) plt.show() plt.close() [\/code]","rel":"","context":"In &quot;python&quot;","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"example13","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2018\/12\/example13.png?resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2018\/12\/example13.png?resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2018\/12\/example13.png?resize=525%2C300 1.5x"},"classes":[]},{"id":1738,"url":"https:\/\/gantovnik.com\/bio-tips\/2023\/01\/204-mandelbrot-fractal-using-python-2-2-2-2-2-2-2-2\/","url_meta":{"origin":1774,"position":2},"title":"#329 Golden section method using python","author":"gantovnik","date":"2023-01-03","format":false,"excerpt":"golden.py [code language=\"python\"] import math as mt def golden(f,a,b,tol=1.0e-10): # Golden section method for determining x # that minimizes the scalar function f(x). # The minimum must be bracketed in (a,b). c1 = (mt.sqrt(5.0)-1.0)\/2.0 c2 = 1.0 - c1 nIt = int(mt.ceil(mt.log(tol\/abs(a-b))\/mt.log(c1))) # first step x1 = c1*a + c2*b\u2026","rel":"","context":"In &quot;optimization&quot;","block_context":{"text":"optimization","link":"https:\/\/gantovnik.com\/bio-tips\/category\/optimization\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/01\/ex329.png?resize=350%2C200&ssl=1","width":350,"height":200},"classes":[]},{"id":1742,"url":"https:\/\/gantovnik.com\/bio-tips\/2023\/01\/204-mandelbrot-fractal-using-python-2-2-2-2-2-2-2-2-2\/","url_meta":{"origin":1774,"position":3},"title":"#330 Gradient method using python","author":"gantovnik","date":"2023-01-03","format":false,"excerpt":"golden.py [code language=\"python\"] import math as mt def golden(f,a,b,tol=1.0e-10): # Golden section method for determining x # that minimizes the scalar function f(x). # The minimum must be bracketed in (a,b). c1 = (mt.sqrt(5.0)-1.0)\/2.0 c2 = 1.0 - c1 nIt = int(mt.ceil(mt.log(tol\/abs(a-b))\/mt.log(c1))) # first step x1 = c1*a + c2*b\u2026","rel":"","context":"In &quot;optimization&quot;","block_context":{"text":"optimization","link":"https:\/\/gantovnik.com\/bio-tips\/category\/optimization\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/01\/ex330.png?resize=350%2C200&ssl=1","width":350,"height":200},"classes":[]},{"id":1984,"url":"https:\/\/gantovnik.com\/bio-tips\/2023\/12\/398-find-points-of-intersection-of-two-circles\/","url_meta":{"origin":1774,"position":4},"title":"#398 Find points of intersection of two circles","author":"gantovnik","date":"2023-12-26","format":false,"excerpt":"- how to find an equation of the common chord of two intersected circles. [code language=\"python\"] import numpy as np import matplotlib.pyplot as plt import math def get_intersections(x0, y0, r0, x1, y1, r1): # circle 1: (x0, y0), radius r0 # circle 2: (x1, y1), radius r1 d=math.sqrt((x1-x0)**2 + (y1-y0)**2)\u2026","rel":"","context":"In &quot;matplotlib&quot;","block_context":{"text":"matplotlib","link":"https:\/\/gantovnik.com\/bio-tips\/category\/matplotlib\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/12\/ex398.png?fit=640%2C480&ssl=1&resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/12\/ex398.png?fit=640%2C480&ssl=1&resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2023\/12\/ex398.png?fit=640%2C480&ssl=1&resize=525%2C300 1.5x"},"classes":[]},{"id":1695,"url":"https:\/\/gantovnik.com\/bio-tips\/2022\/11\/210-parametric-curve-in-3d-2-2-2-2-2-2-2-2-2-2-2-2-2-3-3-2-2-3-2-2-6\/","url_meta":{"origin":1774,"position":5},"title":"#321 unittest &#8211; unit testing framework in python","author":"gantovnik","date":"2022-11-30","format":false,"excerpt":"calculator.py [code language=\"python\"] #A simple calculator class Calculator: #empty constructor def __init__(self): pass #add method - given two numbers, return the addition def add(self, x1, x2): return x1 + x2 #multiply method - given two numbers, return the #multiplication of the two def multiply(self, x1, x2): return x1 * x2\u2026","rel":"","context":"In &quot;python&quot;","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"","src":"","width":0,"height":0},"classes":[]}],"_links":{"self":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts\/1774","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/comments?post=1774"}],"version-history":[{"count":0,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts\/1774\/revisions"}],"wp:attachment":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/media?parent=1774"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/categories?post=1774"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/tags?post=1774"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}