# Copyright (C) 2009-2010, University of Maryland
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/ or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# Author: Ismet Sahin
"""
BFGS quasi-newton optimizer.
All modules in this file are implemented from the book
"Numerical Methods for Unconstrained Optimization and Nonlinear Equations" by
J.E. Dennis and Robert B. Schnabel (Only a few minor modifications are done).
The interface is through the :func:`quasinewton` function. Here is an
example call::
n = 2
x0 = [-0.9 0.9]'
fn = lambda p: (1-p[0])**2 + 100*(p[1]-p[0]**2)**2
grad = lambda p: array([-2*(1-p[0]) - 400*(p[1]-p[0]**2)*p[0], 200*p[1]])
Sx = ones(n,1)
typf = 1 # todo. see what default value is the best
macheps = eps
eta = eps
maxstep = 100
gradtol = 1e-6
steptol = 1e-12 # do not let steptol larger than 1e-9
itnlimit = 1000
result = quasinewton(fn, x0, grad, Sx, typf,
macheps, eta, maxstep, gradtola, steptol, itnlimit)
print("status code %d"%result['status'])
print("x_min=%s, f(x_min)=%g"%(str(result['x']),result['fx']))
print("iterations, function calls, linesearch function calls",
result['iterations'],result['evals'],result['linesearch_evals'])
"""
from __future__ import print_function
__all__ = ["quasinewton"]
from numpy import inf, sqrt, isnan, isinf, finfo, diag, zeros, ones
from numpy import array, linalg, inner, outer, dot, amax, maximum
STATUS = {
1: "Gradient < tolerance",
2: "Step size < tolerance",
3: "Invalid point in line search",
4: "Iterations exceeded",
5: "Max step taken --- function unbounded?",
6: "User abort",
7: "Iterations exceeded in line search",
8: "Line search step size is too small",
9: "Singular Hessian",
}
[docs]
def quasinewton(fn, x0=None, grad=None, Sx=None, typf=1, macheps=None, eta=None,
maxstep=100, gradtol=1e-6, steptol=1e-12, itnlimit=2000,
abort_test=None, monitor=lambda **kw: True):
r"""
Run a quasinewton optimization on the problem.
*fn(x)* is the cost function, which takes a point x and returns a scalar fx.
*x0* is the initial point
*grad* is the analytic gradient (if available)
*Sx* is a scale vector indicating the typical values for parameters in
the fitted result. This is used for a variety of things such as setting
the step size in the finite difference approximation to the gradient, and
controlling numerical accuracy in calculating the Hessian matrix. If for
example some of your model parameters are in the order of 1e-6, then Sx
for those parameters should be set to 1e-6. Default: [1, ...]
*typf* is the typical value for f(x) near the minimum. This is used along
with gradtol to check the gradient stopping condition. Default: 1
*macheps* is the minimum value that can be added to 1 to produce a number
not equal to 1. Default: numpy.finfo(float).eps
*eta* adapts the numerical gradient calculations to machine precision.
Default: *macheps*
*maxstep* is the maximum step size in any gradient step, after normalizing
by *Sx*. Default: 100
*gradtol* is a stopping condition for the fit based on the amount of
improvement expected at the next step. Default: 1e-6
*steptol* is a stopping condition for the fit based on the size
of the step. Default: 1e-12
*itnlimit* is the maximum number of steps to take before stopping.
Default: 2000
*abort_test* is a function which tests whether the user has requested
abort. Default: None.
*monitor(x,fx,step)* is called every iteration so that a user interface
function can monitor the progress of the fit. Default: lambda \*\*kw: True
Returns the fit result as a dictionary:
*status* is a status code indicating why the fit terminated. Turn the
status code into a string with *STATUS[result.status]*. Status values
vary from 1 to 9, with 1 and 2 indicating convergence and the remaining
codes indicating some form of premature termination.
*x* is the minimum point
*fx* is the value fn(x) at the minimum
*H* is the approximate Hessian matrix, which is the inverse of the
covariance matrix
*L* is the cholesky decomposition of H+D, where D is a small correction
to force H+D to be positive definite. To compute parameter uncertainty
*iterations* is the number of iterations
*evals* is the number of function evaluations
*linesearch_evals* is the number of function evaluations for line search
"""
# print("starting QN")
# If some input parameters are not specified, define default values for them
# here. First and second parameters fn and x0 must be defined, others may be
# passed. If you want to set a value to a parameter, say to typf, make
# sure all the parameters before this parameter are specified, in this
# case fn, x0, grad, and Sx if you want to have default values for grad
# and Sx, for each enter [].
# important for also computing fcount (function count)
n = len(x0)
if x0 is None:
x0 = zeros(n)
if grad is None:
analgrad = 0
else:
analgrad = 1
if Sx is None:
Sx = ones(n)
#Sx = x0 + (x0==0.)
elif len(Sx) != n:
raise ValueError("sizes of x0 and Sx must be the same")
if macheps is None:
# PAK: use finfo rather than macheps
macheps = finfo('d').eps
if eta is None:
eta = macheps
fcount = 0 # total function count
fcount_ls = 0 # funciton count due to line search
# If analytic gradient is available then fn will return both function
# value and analytic gradient. Otherwise, use finite difference method
# for estimating the gradient
if analgrad == 1:
fc = fn(x0)
gc = grad(x0)
fcount = fcount + 1
else:
fc = fn(x0)
gc = fdgrad(n, x0, fc, fn, Sx, eta)
fcount = fcount + n + 1
# Check if the initial guess is a local minimizer
termcode = umstop0(n, x0, fc, gc, Sx, typf, gradtol)
consecmax = 0
# Value to return if we fail early
# Approximately x0 is a critical point
xf = x0
ff = fc
H = L = None
if termcode == 0:
H = inithessunfac(n, fc, typf, Sx)
# STEP 9.
xc = x0
# Iterate until convergence in the following loop
itncount = 0
while termcode == 0: # todo. increase itncount
# print("update",itncount)
itncount = itncount + 1
# disp(['Iteration = ' num2str(itncount)])
# Find Newton step sN
H, L = modelhess(n, Sx, macheps, H)
# the vector obtained in the middle
middle_step_v = linalg.solve(L, -gc)
sN = linalg.solve(L.transpose(), middle_step_v) # the last step
if isnan(sN).any():
# print("H",H)
# print("L",L)
# print("v",middle_step_v)
# print("Sx",Sx)
# print("gc",gc)
termcode = 9
break
# Perform line search (Alg.6.3.1). todo. put param order as in the book
# print("calling linesearch",xc,fc,gc,sN,Sx,H,L,middle_step_v)
# print("linesearch",xc,fc)
retcode, xp, fp, maxtaken, fcnt \
= linesearch(fn, n, xc, fc, gc, sN, Sx, maxstep, steptol)
fcount += fcnt
fcount_ls += fcnt
#plot(xp(1), xp(2), 'g.')
# Evaluate gradient at new point xp
if analgrad == 1:
gp = grad(xp)
else:
gp = fdgrad(n, xp, fp, fn, Sx, eta)
fcount = fcount + n
# Check stopping criteria (alg.7.2.1)
consecmax = consecmax + 1 if maxtaken else 0
termcode = umstop(n, xc, xp, fp, gp, Sx, typf, retcode, gradtol,
steptol, itncount, itnlimit, consecmax)
if abort_test():
termcode = 6
# STEP 10.6
# If termcode is larger than zero, we found a point satisfying one
# of the termination criteria, return from here. Otherwise evaluate
# the next Hessian approximation (Alg. 9.4.1).
if termcode > 0:
xf = xp # x final
ff = fp # f final
elif not monitor(x=xp, fx=fp, step=itncount):
termcode = 6
else:
H = bfgsunfac(n, xc, xp, gc, gp, macheps, eta, analgrad, H)
xc = xp
fc = fp
gc = gp
# STOPHERE
result = dict(status=termcode, x=xf, fx=ff, H=H, L=L,
iterations=itncount, evals=fcount, linesearch_evals=fcount_ls)
#print("result",result, steptol, macheps)
return result
#------------------------------------------------------------------------------
#@author: Ismet Sahin
# Alg. 9.4.1
# NOTE:
# BFCG Hessian update is performed unless the following two conditions hold
# (i) y'*s < sqrt(macheps)*norm(s)*norm(y)
# (ii)
def bfgsunfac(n, xc, xp, gc, gp, macheps, eta, analgrad, H):
s = xp - xc
y = gp - gc
temp1 = inner(y, s)
# ISMET : I added condition of having temp1 != 0
if temp1 >= sqrt(macheps) * linalg.norm(s) * linalg.norm(y) and temp1 != 0:
if analgrad == 1:
tol = eta
else:
tol = sqrt(eta)
# deal with noise levels in y
skipupdate = 1
t = dot(H, s)
temp_logicals = (abs(y - t) >= tol * maximum(abs(gc), abs(gp)))
if sum(temp_logicals):
skipupdate = 0
# do the BFGS update if skipdate is false
if skipupdate == 0:
temp2 = dot(s, t)
H = H + outer(y, y) / temp1 - outer(t, t) / temp2
return H
#------------------------------------------------------------------------------
'''
@author: Ismet Sahin
'''
def choldecomp(n, H, maxoffl, macheps):
minl = (macheps) ** (0.25) * maxoffl
if maxoffl == 0:
# H is known to be a positive definite matrix
maxoffl = sqrt(max(abs(diag(H))))
minl2 = sqrt(macheps) * maxoffl
# 3. maxadd is the number (R) specifying the maximum amount added to any
# diagonal entry of Hessian matrix H
maxadd = 0
# 4. form column j of L
L = zeros((n, n))
for j in range(1, n + 1):
L[j - 1, j - 1] = H[j - 1, j - 1] - sum(L[j - 1, 0:j - 1] ** 2)
minljj = 0
for i in range(j + 1, n + 1):
L[i - 1, j - 1] = H[j - 1, i - 1] - \
sum(L[i - 1, 0:j - 1] * L[j - 1, 0:j - 1])
minljj = max(abs(L[i - 1, j - 1]), minljj)
# 4.4
minljj = max(minljj / maxoffl, minl)
# 4.5
if L[j - 1, j - 1] > minljj ** 2:
# normal Cholesky iteration
L[j - 1, j - 1] = sqrt(L[j - 1, j - 1])
else:
# augment H[j-1,j-1]
if minljj < minl2:
minljj = minl2 # occurs only if maxoffl = 0
maxadd = max(maxadd, minljj ** 2 - L[j - 1, j - 1])
L[j - 1, j - 1] = minljj
# 4.6
L[j:n, j - 1] = L[j:n, j - 1] / L[j - 1, j - 1]
return L, maxadd
#------------------------------------------------------------------------------
# ALGORITHM 5.6.3
# Ismet Sahin
# function g = fdgrad(n, xc, fc, objfunc, sx, eta)
# g = fdgrad(@obj_function1, 2, [1 -1]', 10, [1 1], eps)
# NOTATION:
# N : Natural number
# R : Real number
# Rn: nx1 real vector
# Rnxm : nxm real matrix
# INPUTS:
# n : the dimension of the gradient vector (N)
# xc : the current point at which the value of gradient is computed (Rn)
# fc : function value at xc (R)
# objfunc : a function handle which is used to compute function values
# Sx : a n-dim vector, jth entry specifies the typical value of jth param.
# (Rn)
# eta: equals to 1e-DIGITS where DIGITS is an integer specifying the
# number of reliable digits (R)
# OUTPUT:
# g : the n-dim finite difference gradient vector (Rn)
# NOTES :
# hj : is the constant specifying the step size in the direction of jth
# coordinate (R)
# ej : the unit vector, jth column of the identity matrix (Rn)
# COMMENTS:
#--- FIND STEP SIZE hj
# 1.a : sign(x) does not work for us when x = 0 since this makes the step
# size hj zero which is not allowed. (Step size = 0 => gj = inf.)
# 1.b : evaluation of the step size
# 1.c : a trick to reduce error due to finite precision. The line xc(j) =
# xc(j) + hj is equivalent to xc = xc + hj * ej where ej is the jth column
# of identity matrix.
#
#--- EVALUATE APPR. GRADIENT
# First evaluate function at xc + hj * ej and then estimate jth entry of
# the gradient.
def fdgrad(n, xc, fc, fn, Sx, eta):
# create memory for gradient
g = zeros(n)
sqrteta = sqrt(eta)
for j in range(1, n + 1):
#--- FIND STEP SIZE hj
if xc[j - 1] >= 0:
signxcj = 1
else:
signxcj = -1 # 1.a
# 1.b
hj = sqrteta * max(abs(xc[j - 1]), 1 / Sx[j - 1]) * signxcj
# 1.c
tempj = xc[j - 1]
xc[j - 1] = xc[j - 1] + hj
hj = xc[j - 1] - tempj
#--- EVALUATE APPR. GRADIENT
fj = fn(xc)
# PAK: hack for infeasible region: point the other way
if isinf(fj):
fj = fc + hj
g[j - 1] = (fj - fc) / hj
# if isinf(g[j-1]):
# print("fc,fj,hj,Sx,xc",fc,fj,hj,Sx[j-1],xc[j-1])
# now reset the current
xc[j - 1] = tempj
#print("gradient", g)
return g
#------------------------------------------------------------------------------
# @author: Ismet Sahin
# Example call:
# H = inithessunfac(2, f, 1, [1 0.1]')
def inithessunfac(n, f, typf, Sx):
temp = max(abs(f), typf)
H = diag(temp * Sx ** 2)
return H
#------------------------------------------------------------------------------
def linesearch(cost_func, n, xc, fc, g, p, Sx, maxstep, steptol):
"""
ALGORITHM 6.3.1
Ismet Sahin
THE PURPOSE
is to find a step size which yields the new function value smaller than the
current function value, i.e. f(xc + alfa*p) <= f(xc) + alfa * lambda * g'p
CONDITIONS
g'p < 0
alfa < 0.5
NOTATION:
N : Natural number
R : Real number
Rn: nx1 real vector
Rnxm : nxm real matrix
Str: a string
INPUTS
n : dimensionality (N)
xc : the current point ( Rn)
fc : the function value at xc (R)
obj_func : the function handle to evaluate function values (str like :
'@costfunction1')
g : gradient (Rn)
p : the descent direction (Rn)
Sx : scale factors (Rn)
maxstep : maximum step size allowed (R)
steptol : step tolerance in order to break infinite loop in line search (R)
OUTPUTS
retcode : boolean indicating a new point xp found (0) or not (1) (N).
xp : the new point (Rn)
fp : function value at xp (R)
maxtaken : boolean (N)
NOTES:
alfa : is used to prevent function value reductions which are too small.
Here we'll use a very small number in order to accept very small
reductions but not too small.
"""
maxtaken = 0
# alfa specifies how much function value reduction is allowable. The
# smaller the alfa, the smaller the function value reduction we allow.
alfa = 1e-4
# the magnitude of the Newton step
Newtlen = linalg.norm(Sx * p)
if Newtlen > maxstep:
# Newton step is larger than the max acceptable step size (maxstep).
# Make it equal or smaller than maxstep
p = p * (maxstep / Newtlen)
Newtlen = maxstep
initslope = inner(g, p)
# "Relative length of p as calculated in the stopping routine"
# rellength = amax(abs(p) / maximum(abs(xc), Sx)) # this was a bug
rellength = amax(abs(p) / maximum(abs(xc), 1 / Sx))
minlambda = steptol / rellength
lambdaM = 1.0
# In this loop, we try to find an acceptable next point
# xp = xc + lambda * p by finding an optimal lambda based on one
# dimensional quadratic and cubic models
fcount = 0
while True: # 10 starts.
# next point candidate
xp = xc + lambdaM * p
if isnan(xp).any():
#print("nan xp")
retcode = 1
xp, fp = xc, fc
break
if fcount > 20:
#print("too many cycles in linesearch",xp)
retcode = 2
xp, fp = xc, fc
break
# function value at xp
fp = cost_func(xp)
#print("linesearch",fcount,xp,xc,lambdaM,p,fp,fc)
if isinf(fp):
fp = 2 * fc # PAK: infeasible region hack
fcount = fcount + 1
if fp <= fc + alfa * lambdaM * initslope:
# satisfactory xp is found
retcode = 0
if lambdaM == 1.0 and Newtlen > 0.99 * maxstep:
maxtaken = 1
# return from here
break
elif lambdaM < minlambda:
# step length is too small, so a satisfactory xp cannot be found
#print("step",lambdaM,minlambda,steptol,rellength)
retcode = 3
xp, fp = xc, fc
break
else: # 10.3c starts
# reduce lambda by a factor between 0.1 and 0.5
if lambdaM == 1.0:
# first backtrack with one dimensional quadratic fit
lambda_temp = -initslope / (2.0 * (fp - fc - initslope))
#print("L1",lambda_temp)
else:
# perform second and following backtracks with cubic fit
Mt = array([[1.0/lambdaM**2, -1.0/lambda_prev**2],
[-lambda_prev/lambdaM**2, lambdaM/lambda_prev**2]])
vt = array([[fp - fc - lambdaM * initslope],
[fp_prev - fc - lambda_prev * initslope]])
ab = (1.0 / (lambdaM - lambda_prev)) * dot(Mt, vt)
# a = ab(1) and b = ab(2)
disc = ab[1, 0] ** 2 - 3.0 * ab[0, 0] * initslope
#print("Mt,vt,ab,disc",Mt,vt,ab,disc)
if ab[0, 0] == 0.0:
# cubic model turn out to be a quadratic
lambda_temp = -initslope / (2.0 * ab[1, 0])
#print("L2",lambda_temp)
else:
# the model is a legitimate cubic
lambda_temp = (-ab[1, 0] + sqrt(disc)) / (3.0 * ab[0, 0])
#print("L3",lambda_temp)
if lambda_temp > 0.5 * lambdaM:
# larger than half of previous lambda is not allowed.
lambda_temp = 0.5 * lambdaM
#print("L4",lambda_temp)
lambda_prev = lambdaM
fp_prev = fp
if lambda_temp <= 0.1 * lambdaM:
# smaller than 1/10 th of previous lambda is not allowed.
lambdaM = 0.1 * lambdaM
else:
lambdaM = lambda_temp
#print('lambda = ', lambdaM)
# return xp, fp, retcode
return retcode, xp, fp, maxtaken, fcount
#------------------------------------------------------------------------------
# @author: Ismet Sahin
# ALGORITHM 1.3.1
def machineeps():
macheps = 1.0
while (macheps + 1) != 1:
macheps = macheps / 2
macheps = 2 * macheps
return macheps
#------------------------------------------------------------------------------
def modelhess(n, Sx, macheps, H):
"""
@author: Ismet Sahin.
Thanks to Christopher Meeting for his help in converting this module from
Matlab to Python
ALGORITHM 5.5.1
NOTES:
Currently we are not implementing steps 1, 14, and 15 (TODO)
This function performs perturbed Cholesky decomposition (CD) as if the input
Hessian matrix is positive definite. The code for perturbed CD resides in
choldecomp.m file which returns the factored lower triangle matrix L and a
number, maxadd, specifying the largest number added to a diagonal element of
H during the CD decomposition. This function checks if the decomposition is
completed without adding any positive number to the diagonal elements of H,
i.e. maxadd <= 0. Otherwise, this function adds the least number to the
diagonals of H which makes it positive definite based on maxadd and other
entries in H.
EXAMPLE CALLS::
A1 =[2 0 2.4
0 2 0
2.4 0 3]
A2 =[2 0 2.5
0 2 0
2.5 0 3]
A3 =[2 0 10
0 2 0
10 0 3]
"""
# SCALING
scale_needed = 0 # ISMET uses this parameter
if sum(Sx - ones(n)) != 0:
# scaling is requested by the user
scale_needed = 1
Dx = diag(Sx)
Dx_inv = diag(1.0 / Sx)
H = dot(Dx_inv, dot(H, Dx_inv))
# STEP I.
sqrteps = sqrt(macheps)
# 2-4.
H_diag = diag(H)
maxdiag = max(H_diag)
mindiag = min(H_diag)
# 5.
maxposdiag = max(0, maxdiag)
# 6. mu is the amount to be added to diagonal of H before the
# Cholesky decomp. If the minimum diagonal is much much smaller than
# the maximum diagonal element then adjust mu accordingly otherwise mu = 0.
if mindiag <= sqrteps * maxposdiag:
mu = 2 * (maxposdiag - mindiag) * sqrteps - mindiag
maxdiag = maxdiag + mu
else:
mu = 0
# 7. maximum of off-diagonal elements of H
diag_infinite = diag(inf * ones(n))
maxoff = (H - diag_infinite).max()
# 8. if maximum off diagonal element is much much larger than the maximum
# diagonal element of the Hessian H
if maxoff * (1 + 2 * sqrteps) > maxdiag:
mu = mu + (maxoff - maxdiag) + 2 * sqrteps * maxoff
maxdiag = maxoff * (1 + 2 * sqrteps)
# 9.
if maxdiag == 0: # if H == 0
mu = 1
maxdiag = 1
# 10. mu>0 => need to add mu amount to the diagonal elements: H = H + mu*I
if mu > 0:
diag_mu = diag(mu * ones(n))
H = H + diag_mu
# 11.
maxoffl = sqrt(max(maxdiag, maxoff / n))
# STEP II. Perform perturbed Cholesky decomposition H + D = LL' where D is
# a diagonal matrix which is implicitly added to H if H is not positive
# definite. Matrix D has only positive elements. The output variable maxadd
# indicates the maximum number added to a diagonal entry of the Hesian,
# i.e. the maximum of D. If maxadd is returned 0, then H was indeed pd
# and L is the resulting factor.
# 12.
L, maxadd = choldecomp(n, H, maxoffl, macheps)
# STEP III.
# 13. If maxadd <= 0, we are done H was positive definite.
if maxadd > 0:
# H was not positive definite
# print('WARNING: Hessian is not pd. Max number added to H is ',maxadd)
maxev = H[0, 0]
minev = H[0, 0]
for i in range(1, n + 1):
offrow = sum(abs(H[0:i - 1, i - 1])) + sum(abs(H[i - 1, i:n]))
maxev = max(maxev, H[i - 1, i - 1] + offrow)
minev = min(minev, H[i - 1, i - 1] - offrow)
sdd = (maxev - minev) * sqrteps - minev
sdd = max(sdd, 0)
mu = min(maxadd, sdd)
H = H + diag(mu * ones(n))
L, maxadd = choldecomp(n, H, 0, macheps)
if scale_needed: # todo. this calculation can be done faster
H = dot(Dx, dot(H, Dx))
L = dot(Dx, L)
return H, L
#------------------------------------------------------------------------------
def umstop(n, xc, xp, f, g, Sx, typf, retcode, gradtol, steptol,
itncount, itnlimit, consecmax):
"""
#@author: Ismet Sahin
ALGORITHM 7.2.1
Return codes:
Note that return codes are nonnegative integers. When it is not zero, there is
a termination condition which is satisfied.
0 : None of the termination conditions is satisfied
1 : Magnitude of scaled grad is less than gradtol; this is the primary
condition. The new point xp is most likely a local minimizer. If gradtol
is too large, then this condition can be satisfied easier and therefore
xp may not be a local minimizer
2 : Scaled distance between last two points is less than steptol; xp might be
a local minimizer. This condition may also be satisfied if step is
chosen too large or the algorithm is far from the minimizer and making
small progress
3 : The algorithm cannot find a new point giving smaller function value than
the current point. The current may be a local minimizer, or analytic
gradient implementation has some mistakes, or finite difference gradient
estimation is not accurate, or steptol is too large.
4 : Maximum number of iterations are completed
5 : The maximum step length maxstep is taken for last ten consecutive
iterations. This may happen if the function is not bounded from below,
or the function has a finite asymptote in some direction, or maxstep is
too small.
"""
termcode = 0
if retcode == 1:
termcode = 3
elif retcode == 2:
termcode = 7
elif retcode == 3:
termcode = 8
elif retcode > 0:
raise ValueError("Unknown linesearch return code")
elif max(abs(g) * maximum(abs(xp), 1 / Sx) / max(abs(f), typf)) <= gradtol:
# maximum component of scaled gradient is smaller than gradtol.
# TODO: make sure not to use a too large typf value which leads to the
# satisfaction of this algorithm easily.
termcode = 1
elif max(abs(xp - xc) / maximum(abs(xp), 1 / Sx)) <= steptol:
# maximum component of scaled step is smaller than steptol
termcode = 2
elif itncount >= itnlimit:
# maximum number of iterations are performed
termcode = 4
elif consecmax == 10:
# not more than 10 steps will be taken consecutively.
termcode = 5
return termcode
#------------------------------------------------------------------------------
#@author: Ismet Sahin
# This function checks whether initial conditions are acceptable for
# continuing unconstrained optimization
# f : the function value at x0, i.e. f = f(x0), (R)
# g : the gradient at x0, (Rn)
# termcode = 0 : x0 is not a critical point of f(x), (Z)
# termcode = 1 : x0 is a critical point of f(x), (Z)
# Note that x0 may be a critical point of the function; in this case, it is
# either a local minimizer or a saddle point of the function. If the Hessian
# at x0 is positive definite than it is indeed a local minimizer. Instead of
# checking Hessian, we can also restart the driver program umexample from
# another point which is close to x0. If x0 is the local minimizer, the
# algorithm will approach it.
def umstop0(n, x0, f, g, Sx, typf, gradtol):
#consecmax = 0
if max(abs(g) * maximum(abs(x0), 1./Sx)/max(abs(f), typf)) <= 1e-3*gradtol:
termcode = 1
else:
termcode = 0
return termcode
#------------------------------------------------------------------------------
def example_call():
print('***********************************')
# Rosenbrock function
fn = lambda p: (1 - p[0])**2 + 100*(p[1] - p[0]**2)**2
grad = lambda p: array([-2*(1 - p[0]) - 400*(p[1] - p[0]**2)*p[0],
200*(p[1] - p[0]**2)])
x0 = array([2.320894, -0.534223])
# x0 = array([2.0,1.0])
result = quasinewton(fn=fn, x0=x0, grad=grad)
#result = quasinewton(fn=fn, x0=x0)
print('\n\nInitial point x0 = ', x0, ', f(x0) = ', fn(x0))
for k in sorted(result.keys()):
print(k, "=", result[k])
if __name__ == "__main__":
example_call()