{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def newton_gll_step( f,\n", " df,\n", " ddf,\n", " x,\n", " l,\n", " beta = 0.1,\n", " sigma = 0.1,\n", " eps = 1e-10,\n", " return_niter = False,\n", " maxiter = 100000,\n", " kappa = 1e-2,\n", " gamma = 1e-2,\n", " M = 0\n", " ):\n", " \n", " grad = df( x )\n", " f_x = f( x )\n", " \n", " val_list = [ f_x ]\n", " val_list = val_list + [ float( \"-Inf\" ) ] * M\n", " val_list_idx = 0\n", " \n", " grad_sq = np.sum( grad ** 2.0 )\n", " niter = 0\n", " nfeval = 1\n", " while ( grad_sq ** 0.5 > eps ) and ( niter < maxiter ):\n", " \n", " l_k = l\n", " try:\n", " d_k = -np.linalg.solve( ddf( x ), df( x ) )\n", " except np.linalg.LinAlgError:\n", " d_k = -grad\n", " nrm_d_k = np.linalg.norm( d_k )\n", " if d_k @ grad >= -gamma * ( grad_sq ** 0.5 ) * nrm_d_k:\n", " d_k = -grad\n", " if nrm_d_k <= kappa * ( grad_sq ** 0.5 ):\n", " d_k = -grad\n", " \n", " desc = sigma * np.sum( d_k * grad )\n", " \n", " x_next = x + l_k * d_k\n", " f_x_next = f( x_next )\n", " nfeval = nfeval + 1\n", " mx = max( val_list )\n", " while f_x_next >= mx + l_k * desc:\n", " l_k = l_k * beta\n", " x_next = x + l_k * d_k\n", " f_x_next = f( x_next )\n", " nfeval = nfeval + 1\n", " \n", " niter = niter + 1\n", "\n", " x = x_next\n", " grad = df( x )\n", " grad_sq = np.sum( grad ** 2.0 )\n", " f_x = f_x_next\n", " \n", " val_list_idx = ( val_list_idx + 1 ) % ( M + 1 )\n", " val_list[ val_list_idx ] = f_x\n", " \n", " if return_niter :\n", " return ( x, niter, nfeval )\n", " else:\n", " return x" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "l = [ 1, 2, 3, 4, 5 ]\n", "l_idx = 0\n", "val = 6" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[6, 2, 3, 4, 5]\n" ] } ], "source": [ "l[ l_idx ] = val\n", "l_idx = ( l_idx + 1 ) % 5\n", "val = val + 1\n", "print( l )" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[11.1]\n", "11.1\n" ] } ], "source": [ "val_list = [ 11.1 ] + [ float( \"-Inf\" ) ] * 0\n", "print( val_list )\n", "print( max( val_list ) )" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def rosenbrock( x, a = 1, b = 100 ):\n", " \n", " return ( a - x[ 0 ] ) ** 2.0 + b * ( x[ 1 ] - x[ 0 ] ** 2.0 ) ** 2.0\n", "\n", "def rosenbrock_grad( x, a = 1, b = 100 ):\n", " \n", " retval = np.empty( 2 )\n", " \n", " retval[ 0 ] = 2.0 * ( x[ 0 ] - a ) - 4.0 * x[ 0 ] * b * ( x[ 1 ] - x[ 0 ] ** 2.0 )\n", " retval[ 1 ] = 2.0 * b * ( x[ 1 ] - x[ 0 ] ** 2.0 )\n", " \n", " return retval\n", "\n", "def rosenbrock_hess( x, a = 1, b = 100 ):\n", " \n", " retval = np.empty( ( 2, 2 ) )\n", " \n", " retval[ 0, 0 ] = 2.0 * x[ 0 ] - 4.0 * b * ( x[ 1 ] - x[ 0 ] ** 2.0 ) + 8.0 * b * x[ 0 ]\n", " retval[ 0, 1 ] = -4.0 * x[ 0 ] * b\n", " retval[ 1, 0 ] = -4.0 * b * x[ 0 ]\n", " retval[ 1, 1 ] = 2.0 * b\n", " \n", " return retval" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.69203435 6.17892732]\n", "229\n", "114\n", "1.0087719298245614\n", "[1. 1.]\n" ] } ], "source": [ "x_0 = np.array( [ 1.0, 1.0 ] ) + np.random.randn( 2 ) * 3.0\n", "print( x_0 )\n", "x, niter, nfeval = newton_gll_step( \n", " rosenbrock,\n", " rosenbrock_grad,\n", " rosenbrock_hess,\n", " x_0,\n", " 1.0,\n", " return_niter = True,\n", " sigma = 1e-5,\n", " beta = 0.5,\n", " kappa = 1e-6,\n", " gamma = 1e-5\n", " )\n", "\n", "print( niter + nfeval )\n", "print( niter )\n", "print( nfeval / niter )\n", "print( x )" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.69203435 6.17892732]\n", "229\n", "114\n", "1.0087719298245614\n", "[1. 1.]\n" ] } ], "source": [ "print( x_0 )\n", "x, niter, nfeval = newton_gll_step( \n", " rosenbrock,\n", " rosenbrock_grad,\n", " rosenbrock_hess,\n", " x_0,\n", " 1.0,\n", " return_niter = True,\n", " sigma = 1e-5,\n", " beta = 0.5,\n", " kappa = 1e-6,\n", " gamma = 1e-5,\n", " M = 10\n", " )\n", "\n", "print( niter + nfeval )\n", "print( niter )\n", "print( nfeval / niter )\n", "print( x )" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def newton_zh_step( f,\n", " df,\n", " ddf,\n", " x,\n", " l,\n", " beta = 0.1,\n", " sigma = 0.1,\n", " eps = 1e-10,\n", " return_niter = False,\n", " maxiter = 100000,\n", " kappa = 1e-2,\n", " gamma = 1e-2,\n", " eta = 0\n", " ):\n", " \n", " grad = df( x )\n", " f_x = f( x )\n", "\n", " C = f_x\n", " Q = 1\n", " \n", " grad_sq = np.sum( grad ** 2.0 )\n", " niter = 0\n", " nfeval = 1\n", " while ( grad_sq ** 0.5 > eps ) and ( niter < maxiter ):\n", " \n", " l_k = l\n", " try:\n", " d_k = -np.linalg.solve( ddf( x ), df( x ) )\n", " except np.linalg.LinAlgError:\n", " d_k = -grad\n", " nrm_d_k = np.linalg.norm( d_k )\n", " if d_k @ grad >= -gamma * ( grad_sq ** 0.5 ) * nrm_d_k:\n", " d_k = -grad\n", " if nrm_d_k <= kappa * ( grad_sq ** 0.5 ):\n", " d_k = -grad\n", " \n", " desc = sigma * np.sum( d_k * grad )\n", " \n", " x_next = x + l_k * d_k\n", " f_x_next = f( x_next )\n", " nfeval = nfeval + 1\n", "\n", " while f_x_next >= C + l_k * desc:\n", " l_k = l_k * beta\n", " x_next = x + l_k * d_k\n", " f_x_next = f( x_next )\n", " nfeval = nfeval + 1\n", " \n", " niter = niter + 1\n", "\n", " x = x_next\n", " grad = df( x )\n", " grad_sq = np.sum( grad ** 2.0 )\n", " f_x = f_x_next\n", "\n", " Q_next = eta * Q + 1\n", " C = ( eta * Q * C + f_x ) / Q_next\n", " Q = Q_next\n", " \n", " if return_niter :\n", " return ( x, niter, nfeval )\n", " else:\n", " return x" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.69203435 6.17892732]\n", "229\n", "114\n", "1.0087719298245614\n", "[1. 1.]\n" ] } ], "source": [ "print( x_0 )\n", "x, niter, nfeval = newton_zh_step( \n", " rosenbrock,\n", " rosenbrock_grad,\n", " rosenbrock_hess,\n", " x_0,\n", " 1.0,\n", " return_niter = True,\n", " sigma = 1e-4,\n", " beta = 0.5,\n", " kappa = 1e-6,\n", " gamma = 1e-5,\n", " eta = 0.85\n", " )\n", "\n", "print( niter + nfeval )\n", "print( niter )\n", "print( nfeval / niter )\n", "print( x )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.6.9" } }, "nbformat": 4, "nbformat_minor": 2 }