{ "cells": [ { "cell_type": "markdown", "id": "1a523c57", "metadata": {}, "source": [ "\n", "*This notebook contains material from [CBE60499](https://ndcbe.github.io/CBE60499);\n", "content is available [on Github](git@github.com:ndcbe/CBE60499.git).*\n" ] }, { "cell_type": "markdown", "id": "a45391be", "metadata": {}, "source": [ "\n", "< [3.3 Unconstrained Optimality Conditions](https://ndcbe.github.io/CBE60499/03.03-Optimality.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [3.5 Quasi-Newton Methods for Unconstrained Optimization](https://ndcbe.github.io/CBE60499/03.05-Quasi-Newton-Method.html) >

\"Open

\"Download\"" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 1, "link": "[3.4 Newton-type Methods for Unconstrained Optimization](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4-Newton-type-Methods-for-Unconstrained-Optimization)", "section": "3.4 Newton-type Methods for Unconstrained Optimization" } }, "source": [ "# 3.4 Newton-type Methods for Unconstrained Optimization\n", "\n", "**Reference**: Section 2.4.2 in Biegler (2010)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.1 Test Problem: Example 2.19](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.1-Test-Problem:-Example-2.19)", "section": "3.4.1 Test Problem: Example 2.19" } }, "source": [ "## 3.4.1 Test Problem: Example 2.19\n", "\n", "![ex2-19](figures/ex2-19.png)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "nbpages": { "level": 2, "link": "[3.4.1 Test Problem: Example 2.19](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.1-Test-Problem:-Example-2.19)", "section": "3.4.1 Test Problem: Example 2.19" } }, "outputs": [], "source": [ "# Load required Python libraries.\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from scipy import linalg\n", "\n", "## Define Python function to calculate objective\n", "def my_f(x,verbose=False):\n", " ''' Evaluate function given above at point x\n", "\n", " Inputs:\n", " x - vector with 2 elements\n", " \n", " Outputs:\n", " f - function value (scalar)\n", " '''\n", " # Constants\n", " a = np.array([0.3, 0.6, 0.2])\n", " b = np.array([5, 26, 3])\n", " c = np.array([40, 1, 10])\n", " \n", " # Intermediates. Recall Python indicies start at 0\n", " u = x[0] - 0.8\n", " s = np.sqrt(1-u)\n", " s2 = np.sqrt(1+u)\n", " v = x[1] -(a[0] + a[1]*u**2*s-a[2]*u)\n", " alpha = -b[0] + b[1]*u**2*s2+ b[2]*u # September 5, 2018: changed 's' to 's2'\n", " beta = c[0]*v**2*(1-c[1]*v)/(1+c[2]*u**2)\n", " \n", " if verbose:\n", " print(\"##### my_f at x = \",x, \"#####\")\n", " print(\"u = \",u)\n", " print(\"sqrt(1-u) = \",s)\n", " print(\"sqrt(1+u) = \",s2)\n", " print(\"v = \",v)\n", " print(\"alpha = \",alpha)\n", " print(\"beta = \",beta)\n", " print(\"f(x) = \",)\n", " print(\"##### Done. #####\\n\")\n", " \n", " return alpha*np.exp(-beta)" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.2 Helper Functions](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.2-Helper-Functions)", "section": "3.4.2 Helper Functions" } }, "source": [ "## 3.4.2 Helper Functions" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "nbpages": { "level": 2, "link": "[3.4.2 Helper Functions](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.2-Helper-Functions)", "section": "3.4.2 Helper Functions" } }, "outputs": [], "source": [ "## Calculate gradient with central finite difference\n", "def my_grad_approx(x,f,eps1,verbose=False):\n", " '''\n", " Calculate gradient of function my_f using central difference formula\n", " \n", " Inputs:\n", " x - point for which to evaluate gradient\n", " f - function to consider\n", " eps1 - perturbation size\n", " \n", " Outputs:\n", " grad - gradient (vector)\n", " '''\n", " \n", " n = len(x)\n", " grad = np.zeros(n)\n", " \n", " if(verbose):\n", " print(\"***** my_grad_approx at x = \",x,\"*****\")\n", " \n", " for i in range(0,n):\n", " \n", " # Create vector of zeros except eps in position i\n", " e = np.zeros(n)\n", " e[i] = eps1\n", " \n", " # Finite difference formula\n", " my_f_plus = f(x + e)\n", " my_f_minus = f(x - e)\n", " \n", " # Diagnostics\n", " if(verbose):\n", " print(\"e[\",i,\"] = \",e)\n", " print(\"f(x + e[\",i,\"]) = \",my_f_plus)\n", " print(\"f(x - e[\",i,\"]) = \",my_f_minus)\n", " \n", " \n", " grad[i] = (my_f_plus - my_f_minus)/(2*eps1)\n", " \n", " if(verbose):\n", " print(\"***** Done. ***** \\n\")\n", " \n", " return grad\n", "\n", "## Calculate gradient using central finite difference and my_hes_approx\n", "def my_hes_approx(x,grad,eps2):\n", " '''\n", " Calculate gradient of function my_f using central difference formula and my_grad\n", " \n", " Inputs:\n", " x - point for which to evaluate gradient\n", " grad - function to calculate the gradient\n", " eps2 - perturbation size (for Hessian NOT gradient approximation)\n", " \n", " Outputs:\n", " H - Hessian (matrix)\n", " '''\n", " \n", " n = len(x)\n", " H = np.zeros([n,n])\n", " \n", " for i in range(0,n):\n", " # Create vector of zeros except eps in position i\n", " e = np.zeros(n)\n", " e[i] = eps2\n", " \n", " # Evaluate gradient twice\n", " grad_plus = grad(x + e)\n", " grad_minus = grad(x - e)\n", " \n", " # Notice we are building the Hessian by column (or row)\n", " H[:,i] = (grad_plus - grad_minus)/(2*eps2)\n", "\n", " return H\n", "\n", "def check_nan(A):\n", " return np.sum(np.isnan(A))\n", "\n", "## Analyze Hessian.\n", "def analyze_hes(B):\n", " print(B,\"\\n\")\n", " \n", " if not check_nan(B):\n", " # Calculate eigenvalues\n", " l = linalg.eigvals(B)\n", " print(\"Eigenvalues: \",l,\"\\n\")" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.3 Algorithm 2.1: Basic Newton Method](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.3-Algorithm-2.1:-Basic-Newton-Method)", "section": "3.4.3 Algorithm 2.1: Basic Newton Method" } }, "source": [ "## 3.4.3 Algorithm 2.1: Basic Newton Method\n", "\n", "![alg2-1](figures/alg2-1.png)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "nbpages": { "level": 2, "link": "[3.4.3 Algorithm 2.1: Basic Newton Method](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.3-Algorithm-2.1:-Basic-Newton-Method)", "section": "3.4.3 Algorithm 2.1: Basic Newton Method" } }, "outputs": [], "source": [ "## Implement Alg 1 in a function\n", "def alg1(x0,calc_f,calc_grad,calc_hes,eps1,eps2,verbose=False,max_iter=250):\n", " '''\n", " Arguments:\n", " x0 - starting point\n", " calc_f - funcation that calculates f(x)\n", " calc_grad - function that calculates gradient(x)\n", " calc_hes - function that calculates hessian(x)\n", " \n", " Outputs:\n", " x - iteration history of x (decision variables)\n", " f - iteration history of f(x) (objective value)\n", " p - iteration history of p (steps)\n", " '''\n", " \n", " # Allocate outputs as empty lists\n", " x = []\n", " f = []\n", " p = []\n", " \n", " # Store starting point\n", " x.append(x0)\n", " k = 0\n", " \n", " flag = True\n", " \n", " print(\"Iter. \\tf(x) \\t\\t||grad(x)|| \\t||p|| \\t\\tmin(lambda)\")\n", " \n", " while flag:\n", " # Evaluate f(x) at current iteration\n", " f.append(calc_f(x[k]))\n", " \n", " # Evaluate gradient\n", " grad = calc_grad(x[k])\n", " \n", " if(check_nan(grad)):\n", " print(\"WARNING: gradiant calculation returned NaN\")\n", " break\n", " \n", " # Evaluate Hessian\n", " hes = calc_hes(x[k])\n", " \n", " if(check_nan(hes)):\n", " print(\"WARNING: Hessian calculation returned NaN\")\n", " break\n", " \n", " if verbose:\n", " print(\"\\n\")\n", " print(\"k = \",k)\n", " print(\"x = \",x[k])\n", " print(\"grad = \",grad)\n", " print(\"hes = \\n\",hes)\n", " \n", " # Check if singular via condition number\n", " c = np.linalg.cond(hes)\n", " if c > 1E12:\n", " flag = False\n", " print(\"Warning: Hessian is near singular.\")\n", " \n", " else:\n", " # Calculate step\n", " \n", " p.append(linalg.solve(hes,-grad))\n", " \n", " if verbose:\n", " print(\"p = \",p[k])\n", " \n", " # Take step. x[k+1] = x[k] + p[k]\n", " x.append(x[k] + p[k])\n", " \n", " # Calculate norms\n", " norm_p = linalg.norm(p[k])\n", " norm_g = linalg.norm(grad)\n", " \n", " # Calculate eigenvalues (for display only)\n", " ev = np.real(linalg.eigvals(hes))\n", " \n", " # print(\"k = \",k,\"\\t\"f[k],\"\\t\",norm_g,\"\\t\",norm_p)\n", " print(k,' \\t{0: 1.4e} \\t{1:1.4e} \\t{2:1.4e} \\t{3: 1.4e}'.format(f[k],norm_g,norm_p,np.min(ev)))\n", " \n", " # Check convergence criteria\n", " flag = (norm_p > eps1) and (norm_g > eps2)\n", " \n", " # Update iteration counter\n", " k = k + 1\n", "\n", " if k > max_iter:\n", " flag = False\n", " print(\"Maximum number of iterations.\")\n", " print(\"Done.\")\n", " print(\"x* = \",x[-1])\n", " \n", " return x,f,p" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.4 Starting Point Near Optimal Solution](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.4-Starting-Point-Near-Optimal-Solution)", "section": "3.4.4 Starting Point Near Optimal Solution" } }, "source": [ "## 3.4.4 Starting Point Near Optimal Solution" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "nbpages": { "level": 2, "link": "[3.4.4 Starting Point Near Optimal Solution](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.4-Starting-Point-Near-Optimal-Solution)", "section": "3.4.4 Starting Point Near Optimal Solution" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iter. \tf(x) \t\t||grad(x)|| \t||p|| \t\tmin(lambda)\n", "0 \t-4.9246e+00 \t1.0874e+01 \t4.1946e-02 \t 4.9368e+01\n", "1 \t-5.0888e+00 \t6.2736e-01 \t1.9752e-03 \t 4.2541e+01\n", "2 \t-5.0893e+00 \t2.4210e-03 \t3.0032e-05 \t 4.3424e+01\n", "3 \t-5.0893e+00 \t2.9199e-07 \t3.1620e-09 \t 4.3418e+01\n", "Done.\n", "x* = [0.73950546 0.3143601 ]\n", "Hessian at x*= \n", "[[ 77.01173033 108.33423048]\n", " [108.33423048 392.76693009]] \n", "\n", "Eigenvalues: [ 43.41702924+0.j 426.36163117+0.j] \n", "\n" ] } ], "source": [ "## Test on example\n", "\n", "# Specify convergence criteria\n", "eps1 = 1E-8\n", "eps2 = 1E-4\n", "\n", "# Create a Lambda (anonymous) function for gradient calculation\n", "calc_grad = lambda x : my_grad_approx(x,my_f,1E-6)\n", "\n", "# Create a Lambda (anonymous) function for Hessian calculation\n", "calc_hes = lambda x : my_hes_approx(x,calc_grad,1E-6)\n", "\n", "# Specify starting point\n", "x0 = np.array([0.7, 0.3])\n", "\n", "# Call optimization routine\n", "x,f,p = alg1(x0,my_f,calc_grad,calc_hes,eps1,eps2);\n", "\n", "# Actual Hessian\n", "print(\"Hessian at x*= \")\n", "analyze_hes(my_hes_approx(x[-1],calc_grad,1E-6))" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.5 Activity: A Different Starting Point](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.5-Activity:-A-Different-Starting-Point)", "section": "3.4.5 Activity: A Different Starting Point" } }, "source": [ "## 3.4.5 Activity: A Different Starting Point\n", "\n", "Try with $x_0 = [0,0]^T$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "nbpages": { "level": 2, "link": "[3.4.5 Activity: A Different Starting Point](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.5-Activity:-A-Different-Starting-Point)", "section": "3.4.5 Activity: A Different Starting Point" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iter. \tf(x) \t\t||grad(x)|| \t||p|| \t\tmin(lambda)\n", "0 \t 1.6212e-06 \t1.5568e-04 \t3.8590e-02 \t-2.9496e-03\n", "1 \t 5.4102e-07 \t5.5405e-05 \t3.6202e-02 \t-1.0900e-03\n", "Done.\n", "x* = [-0.00124597 -0.07478209]\n", "Hessian at x* = \n", "[[0.0001531 0.00055527]\n", " [0.00055527 0.00014766]] \n", "\n", "Eigenvalues: [ 0.00070565+0.j -0.00040489+0.j] \n", "\n" ] } ], "source": [ "# Specify starting point\n", "x0 = np.array([0, 0])\n", "\n", "# Call optimization routine\n", "x,f,p = alg1(x0,my_f,calc_grad,calc_hes,eps1,eps2);\n", "\n", "# Actual Hessian\n", "print(\"Hessian at x* = \")\n", "analyze_hes(my_hes_approx(x[-1],calc_grad,1E-6))" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.6 Activity: Let's break it](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.6-Activity:-Let's-break-it)", "section": "3.4.6 Activity: Let's break it" } }, "source": [ "## 3.4.6 Activity: Let's break it\n", "\n", "**Activity**: Try $x_0 = [-0.2, -0.2]^T$. Why does the gradient or Hessian return NaN? *Hint*: redefine `calc_grad` and create `my_f_verbose` using `verbose=True`." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "nbpages": { "level": 2, "link": "[3.4.6 Activity: Let's break it](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.6-Activity:-Let's-break-it)", "section": "3.4.6 Activity: Let's break it" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iter. \tf(x) \t\t||grad(x)|| \t||p|| \t\tmin(lambda)\n", "***** my_grad_approx at x = [-0.2 -0.2] *****\n", "##### my_f at x = [-0.199999 -0.2 ] #####\n", "u = -0.9999990000000001\n", "sqrt(1-u) = 1.4142132088196604\n", "sqrt(1+u) = 0.0009999999999588667\n", "v = -1.5485260282367943\n", "alpha = -7.973997052001044\n", "beta = 22.222565160998283\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.200001 -0.2 ] #####\n", "u = -1.0000010000000001\n", "sqrt(1-u) = 1.4142139159264415\n", "sqrt(1+u) = nan\n", "v = -1.5485302466134128\n", "alpha = nan\n", "beta = 22.222642208927322\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 0 ] = [1.e-06 0.e+00]\n", "f(x + e[ 0 ]) = -1.780486346234917e-09\n", "f(x - e[ 0 ]) = nan\n", "##### my_f at x = [-0.2 -0.199999] #####\n", "u = -1.0\n", "sqrt(1-u) = 1.4142135623730951\n", "sqrt(1+u) = 0.0\n", "v = -1.548527137423857\n", "alpha = -8.0\n", "beta = 22.222566263573963\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.2 -0.200001] #####\n", "u = -1.0\n", "sqrt(1-u) = 1.4142135623730951\n", "sqrt(1+u) = 0.0\n", "v = -1.548529137423857\n", "alpha = -8.0\n", "beta = 22.22264110629725\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 1 ] = [0.e+00 1.e-06]\n", "f(x + e[ 1 ]) = -1.786290485440547e-09\n", "f(x - e[ 1 ]) = -1.7861567995988046e-09\n", "***** Done. ***** \n", "\n", "WARNING: gradiant calculation returned NaN\n", "Done.\n", "x* = [-0.2 -0.2]\n", "Hessian at x* = \n", "***** my_grad_approx at x = [-0.199999 -0.2 ] *****\n", "##### my_f at x = [-0.199998 -0.2 ] #####\n", "u = -0.999998\n", "sqrt(1-u) = 1.414212855266137\n", "sqrt(1+u) = 0.0014142135623541761\n", "v = -1.5485239190522238\n", "alpha = -7.963224594456855\n", "beta = 22.22252663717692\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.2 -0.2] #####\n", "u = -1.0\n", "sqrt(1-u) = 1.4142135623730951\n", "sqrt(1+u) = 0.0\n", "v = -1.5485281374238569\n", "alpha = -8.0\n", "beta = 22.22260368491507\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 0 ] = [1.e-06 0.e+00]\n", "f(x + e[ 0 ]) = -1.778149501075135e-09\n", "f(x - e[ 0 ]) = -1.7862236413056787e-09\n", "##### my_f at x = [-0.199999 -0.199999] #####\n", "u = -0.9999990000000001\n", "sqrt(1-u) = 1.4142132088196604\n", "sqrt(1+u) = 0.0009999999999588667\n", "v = -1.5485250282367944\n", "alpha = -7.973997052001044\n", "beta = 22.222527739675733\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.199999 -0.200001] #####\n", "u = -0.9999990000000001\n", "sqrt(1-u) = 1.4142132088196604\n", "sqrt(1+u) = 0.0009999999999588667\n", "v = -1.5485270282367944\n", "alpha = -7.973997052001044\n", "beta = 22.222602582361898\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 1 ] = [0.e+00 1.e-06]\n", "f(x + e[ 1 ]) = -1.7805529756354463e-09\n", "f(x - e[ 1 ]) = -1.7804197192545872e-09\n", "***** Done. ***** \n", "\n", "***** my_grad_approx at x = [-0.200001 -0.2 ] *****\n", "##### my_f at x = [-0.2 -0.2] #####\n", "u = -1.0\n", "sqrt(1-u) = 1.4142135623730951\n", "sqrt(1+u) = 0.0\n", "v = -1.5485281374238569\n", "alpha = -8.0\n", "beta = 22.22260368491507\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.200002 -0.2 ] #####\n", "u = -1.000002\n", "sqrt(1-u) = 1.4142142694796995\n", "sqrt(1+u) = nan\n", "v = -1.5485323558054604\n", "alpha = nan\n", "beta = 22.222680733035\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 0 ] = [1.e-06 0.e+00]\n", "f(x + e[ 0 ]) = -1.7862236413056787e-09\n", "f(x - e[ 0 ]) = nan\n", "##### my_f at x = [-0.200001 -0.199999] #####\n", "u = -1.0000010000000001\n", "sqrt(1-u) = 1.4142139159264415\n", "sqrt(1+u) = nan\n", "v = -1.5485292466134128\n", "alpha = nan\n", "beta = 22.22260478756765\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.200001 -0.200001] #####\n", "u = -1.0000010000000001\n", "sqrt(1-u) = 1.4142139159264415\n", "sqrt(1+u) = nan\n", "v = -1.548531246613413\n", "alpha = nan\n", "beta = 22.222679630328063\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 1 ] = [0.e+00 1.e-06]\n", "f(x + e[ 1 ]) = nan\n", "f(x - e[ 1 ]) = nan\n", "***** Done. ***** \n", "\n", "***** my_grad_approx at x = [-0.2 -0.199999] *****\n", "##### my_f at x = [-0.199999 -0.199999] #####\n", "u = -0.9999990000000001\n", "sqrt(1-u) = 1.4142132088196604\n", "sqrt(1+u) = 0.0009999999999588667\n", "v = -1.5485250282367944\n", "alpha = -7.973997052001044\n", "beta = 22.222527739675733\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.200001 -0.199999] #####\n", "u = -1.0000010000000001\n", "sqrt(1-u) = 1.4142139159264415\n", "sqrt(1+u) = nan\n", "v = -1.5485292466134128\n", "alpha = nan\n", "beta = 22.22260478756765\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 0 ] = [1.e-06 0.e+00]\n", "f(x + e[ 0 ]) = -1.7805529756354463e-09\n", "f(x - e[ 0 ]) = nan\n", "##### my_f at x = [-0.2 -0.199998] #####\n", "u = -1.0\n", "sqrt(1-u) = 1.4142135623730951\n", "sqrt(1+u) = 0.0\n", "v = -1.548526137423857\n", "alpha = -8.0\n", "beta = 22.222528842273906\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.2 -0.2] #####\n", "u = -1.0\n", "sqrt(1-u) = 1.4142135623730951\n", "sqrt(1+u) = 0.0\n", "v = -1.5485281374238569\n", "alpha = -8.0\n", "beta = 22.22260368491507\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 1 ] = [0.e+00 1.e-06]\n", "f(x + e[ 1 ]) = -1.7863573320035263e-09\n", "f(x - e[ 1 ]) = -1.7862236413056787e-09\n", "***** Done. ***** \n", "\n", "***** my_grad_approx at x = [-0.2 -0.200001] *****\n", "##### my_f at x = [-0.199999 -0.200001] #####\n", "u = -0.9999990000000001\n", "sqrt(1-u) = 1.4142132088196604\n", "sqrt(1+u) = 0.0009999999999588667\n", "v = -1.5485270282367944\n", "alpha = -7.973997052001044\n", "beta = 22.222602582361898\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.200001 -0.200001] #####\n", "u = -1.0000010000000001\n", "sqrt(1-u) = 1.4142139159264415\n", "sqrt(1+u) = nan\n", "v = -1.548531246613413\n", "alpha = nan\n", "beta = 22.222679630328063\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 0 ] = [1.e-06 0.e+00]\n", "f(x + e[ 0 ]) = -1.7804197192545872e-09\n", "f(x - e[ 0 ]) = nan\n", "##### my_f at x = [-0.2 -0.2] #####\n", "u = -1.0\n", "sqrt(1-u) = 1.4142135623730951\n", "sqrt(1+u) = 0.0\n", "v = -1.5485281374238569\n", "alpha = -8.0\n", "beta = 22.22260368491507\n", "f(x) = \n", "##### Done. #####\n", "\n", "##### my_f at x = [-0.2 -0.200002] #####\n", "u = -1.0\n", "sqrt(1-u) = 1.4142135623730951\n", "sqrt(1+u) = 0.0\n", "v = -1.548530137423857\n", "alpha = -8.0\n", "beta = 22.222678527720475\n", "f(x) = \n", "##### Done. #####\n", "\n", "e[ 1 ] = [0.e+00 1.e-06]\n", "f(x + e[ 1 ]) = -1.7862236413056787e-09\n", "f(x - e[ 1 ]) = -1.7860899603198772e-09\n", "***** Done. ***** \n", "\n", "[[ nan nan]\n", " [ nan -2.42801152e-06]] \n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:24: RuntimeWarning: invalid value encountered in sqrt\n" ] } ], "source": [ "# Specify starting point\n", "x0 = np.array([-0.2, -0.2])\n", "\n", "# Create a Lambda (anonymous) function for f(x) calculation with verbose output\n", "my_f_verbose = lambda x : my_f(x,verbose=True);\n", "\n", "# Create a Lambda (anonymous) function for gradient calculation\n", "calc_grad = lambda x : my_grad_approx(x,my_f_verbose,1E-6,verbose=True);\n", "\n", "# Call optimization routine\n", "x,f,p = alg1(x0,my_f,calc_grad,calc_hes,eps1,eps2,verbose=True);\n", "\n", "# Actual Hessian\n", "print(\"Hessian at x* = \")\n", "analyze_hes(my_hes_approx(x[-1],calc_grad,1E-6))" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.7 Activity: Use $I$ in place of Hessian](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.7-Activity:-Use-$I$-in-place-of-Hessian)", "section": "3.4.7 Activity: Use $I$ in place of Hessian" } }, "source": [ "## 3.4.7 Activity: Use $I$ in place of Hessian\n", "\n", "Create an alternative function that returns the identity matrix for the Hessian." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "nbpages": { "level": 2, "link": "[3.4.7 Activity: Use $I$ in place of Hessian](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.7-Activity:-Use-$I$-in-place-of-Hessian)", "section": "3.4.7 Activity: Use $I$ in place of Hessian" } }, "outputs": [], "source": [ "def eye_hes(x):\n", " return 50*np.eye(len(x))" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.7 Activity: Use $I$ in place of Hessian](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.7-Activity:-Use-$I$-in-place-of-Hessian)", "section": "3.4.7 Activity: Use $I$ in place of Hessian" } }, "source": [ "Test with starting point $x=[0.7,0.3]^T$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "nbpages": { "level": 2, "link": "[3.4.7 Activity: Use $I$ in place of Hessian](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.7-Activity:-Use-$I$-in-place-of-Hessian)", "section": "3.4.7 Activity: Use $I$ in place of Hessian" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iter. \tf(x) \t\t||grad(x)|| \t||p|| \t\tmin(lambda)\n", "0 \t-4.9246e+00 \t1.0874e+01 \t2.1749e-01 \t 5.0000e+01\n", "1 \t-1.4795e+00 \t1.6883e+01 \t3.3765e-01 \t 5.0000e+01\n", "2 \t-1.8403e+00 \t2.6363e+01 \t5.2725e-01 \t 5.0000e+01\n", "3 \t-1.2461e-01 \t1.6246e+00 \t3.2493e-02 \t 5.0000e+01\n", "4 \t-1.9169e-01 \t2.5753e+00 \t5.1507e-02 \t 5.0000e+01\n", "5 \t-3.8498e-01 \t5.2112e+00 \t1.0422e-01 \t 5.0000e+01\n", "6 \t-1.4594e+00 \t1.6836e+01 \t3.3673e-01 \t 5.0000e+01\n", "7 \t-1.9249e+00 \t2.6467e+01 \t5.2934e-01 \t 5.0000e+01\n", "8 \t-1.1531e-01 \t1.5097e+00 \t3.0195e-02 \t 5.0000e+01\n", "9 \t-1.7247e-01 \t2.3351e+00 \t4.6703e-02 \t 5.0000e+01\n", "10 \t-3.2703e-01 \t4.5030e+00 \t9.0059e-02 \t 5.0000e+01\n", "11 \t-1.0837e+00 \t1.3485e+01 \t2.6971e-01 \t 5.0000e+01\n", "12 \t-4.4497e+00 \t2.0946e+01 \t4.1891e-01 \t 5.0000e+01\n", "13 \t-1.6236e-01 \t2.1013e+00 \t4.2026e-02 \t 5.0000e+01\n", "14 \t-2.8114e-01 \t3.6906e+00 \t7.3812e-02 \t 5.0000e+01\n", "15 \t-7.2912e-01 \t9.1154e+00 \t1.8231e-01 \t 5.0000e+01\n", "16 \t-4.1506e+00 \t2.1547e+01 \t4.3094e-01 \t 5.0000e+01\n", "17 \t-5.4382e-03 \t2.2424e-01 \t4.4848e-03 \t 5.0000e+01\n", "18 \t-6.5332e-03 \t2.6508e-01 \t5.3017e-03 \t 5.0000e+01\n", "19 \t-8.0846e-03 \t3.2179e-01 \t6.4357e-03 \t 5.0000e+01\n", "20 \t-1.0414e-02 \t4.0482e-01 \t8.0965e-03 \t 5.0000e+01\n", "21 \t-1.4199e-02 \t5.3561e-01 \t1.0712e-02 \t 5.0000e+01\n", "22 \t-2.1098e-02 \t7.6415e-01 \t1.5283e-02 \t 5.0000e+01\n", "23 \t-3.6112e-02 \t1.2322e+00 \t2.4644e-02 \t 5.0000e+01\n", "24 \t-8.0386e-02 \t2.4796e+00 \t4.9592e-02 \t 5.0000e+01\n", "25 \t-3.1729e-01 \t7.8096e+00 \t1.5619e-01 \t 5.0000e+01\n", "Maximum number of iterations.\n", "Done.\n", "x* = [0.71010106 0.23551614]\n", "Hessian at x* = \n", "[[ 45.25146924 7.93809463]\n", " [ 7.93809463 147.99461656]] \n", "\n", "Eigenvalues: [ 44.64177775+0.j 148.60430805+0.j] \n", "\n" ] } ], "source": [ "# Specify starting point\n", "x0 = np.array([0.7, 0.3])\n", "\n", "# Create a Lambda (anonymous) function for gradient calculation\n", "calc_grad = lambda x : my_grad_approx(x,my_f,1E-6,verbose=False);\n", "\n", "# Call optimization routine\n", "x,f,p = alg1(x0,my_f,calc_grad,eye_hes,eps1,eps2,max_iter=25,verbose=False);\n", "\n", "# Actual Hessian\n", "print(\"Hessian at x* = \")\n", "analyze_hes(my_hes_approx(x[-1],calc_grad,1E-6))" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.7 Activity: Use $I$ in place of Hessian](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.7-Activity:-Use-$I$-in-place-of-Hessian)", "section": "3.4.7 Activity: Use $I$ in place of Hessian" } }, "source": [ "Test with starting point $x=[0.0,0.0]^T$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "nbpages": { "level": 2, "link": "[3.4.7 Activity: Use $I$ in place of Hessian](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.7-Activity:-Use-$I$-in-place-of-Hessian)", "section": "3.4.7 Activity: Use $I$ in place of Hessian" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iter. \tf(x) \t\t||grad(x)|| \t||p|| \t\tmin(lambda)\n", "0 \t 1.6212e-06 \t1.5568e-04 \t3.8590e-02 \t-2.9496e-03\n", "1 \t 5.4102e-07 \t5.5405e-05 \t3.6202e-02 \t-1.0900e-03\n", "Done.\n", "x* = [-0.00124597 -0.07478209]\n", "Hessian at x*= \n", "[[0.0001531 0.00055527]\n", " [0.00055527 0.00014766]] \n", "\n", "Eigenvalues: [ 0.00070565+0.j -0.00040489+0.j] \n", "\n" ] } ], "source": [ "# Specify starting point\n", "x0 = np.array([0.0, 0.0])\n", "\n", "# Call optimization routine\n", "x,f,p = alg1(x0,my_f,calc_grad,calc_hes,eps1,eps2);\n", "\n", "# Actual Hessian\n", "print(\"Hessian at x*= \")\n", "analyze_hes(my_hes_approx(x[-1],calc_grad,1E-6))" ] }, { "cell_type": "markdown", "metadata": { "nbpages": { "level": 2, "link": "[3.4.8 Adjust Hessian with Levenberg-Marquardt Correction](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.8-Adjust-Hessian-with-Levenberg-Marquardt-Correction)", "section": "3.4.8 Adjust Hessian with Levenberg-Marquardt Correction" } }, "source": [ "## 3.4.8 Adjust Hessian with Levenberg-Marquardt Correction" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "nbpages": { "level": 2, "link": "[3.4.8 Adjust Hessian with Levenberg-Marquardt Correction](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.8-Adjust-Hessian-with-Levenberg-Marquardt-Correction)", "section": "3.4.8 Adjust Hessian with Levenberg-Marquardt Correction" } }, "outputs": [], "source": [ "def adjusted_hes(x,grad,eps2,eps3,verbose=False):\n", " \n", " # Estimate Hessian with finite difference\n", " hes = my_hes_approx(x,grad,eps2)\n", " \n", " # Calculate eigenvalues\n", " l, V = linalg.eig(hes)\n", " smallest_ev = np.min(np.real(l))\n", " \n", " # Calculate modification\n", " delta = 0\n", " if(smallest_ev - eps3 < 0):\n", " delta = -smallest_ev + eps3\n", " \n", " if(verbose):\n", " print(\"Added ({0:1.4e})*I in LM correction.\".format(delta))\n", " \n", " # Adjust hessian with Levenberg-Marquardt Correction\n", " return V.dot(np.diag(np.real(l)) + delta*np.eye(len(x))).dot(V.T)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "nbpages": { "level": 2, "link": "[3.4.8 Adjust Hessian with Levenberg-Marquardt Correction](https://ndcbe.github.io/CBE60499/03.04-Netwon-Methods.html#3.4.8-Adjust-Hessian-with-Levenberg-Marquardt-Correction)", "section": "3.4.8 Adjust Hessian with Levenberg-Marquardt Correction" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Iter. \tf(x) \t\t||grad(x)|| \t||p|| \t\tmin(lambda)\n", "Added (1.2030e-03)*I in LM correction.\n", "0 \t-2.0284e-07 \t6.5946e-06 \t6.5923e-03 \t 1.0000e-03\n", "Done.\n", "x* = [-0.09849776 -0.0935812 ]\n", "Hessian at x*= \n", "[[ 4.20491243e-05 -5.31304857e-05]\n", " [-5.31304857e-05 -2.35437850e-04]] \n", "\n", "Eigenvalues: [ 5.18741524e-05+0.j -2.45262878e-04+0.j] \n", "\n" ] } ], "source": [ "# Specify starting point\n", "x0 = np.array([-0.1, -0.1])\n", "\n", "# Create a Lambda (anonymous) function for Hessian calculation\n", "calc_hes = lambda x : adjusted_hes(x,calc_grad,1E-6,1E-3,verbose=True)\n", "\n", "# Call optimization routine\n", "x,f,p = alg1(x0,my_f,calc_grad,calc_hes,eps1,eps2,verbose=False);\n", "\n", "# Actual Hessian\n", "print(\"Hessian at x*= \")\n", "analyze_hes(my_hes_approx(x[-1],calc_grad,1E-6))" ] }, { "cell_type": "markdown", "id": "d5894f6a", "metadata": {}, "source": [ "\n", "< [3.3 Unconstrained Optimality Conditions](https://ndcbe.github.io/CBE60499/03.03-Optimality.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [3.5 Quasi-Newton Methods for Unconstrained Optimization](https://ndcbe.github.io/CBE60499/03.05-Quasi-Newton-Method.html) >

\"Open

\"Download\"" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }