14
14
# Updated strong Wolfe conditions line search to use cubic-interpolation (Mar. 2004)
16
16
# Minimization routines
19
A collection of general-purpose optimization routines using Numeric
24
fmin --- Nelder-Mead Simplex algorithm (uses only function calls).
25
fmin_powell --- Powell (direction set) method (uses only function calls).
26
fmin_cg --- Non-linear (Polak-Rubiere) conjugate gradient algorithm
27
(can use function and gradient).
28
fmin_bfgs --- Quasi-Newton method (Broyden-Fletcher-Goldfarb-Shanno
29
(can use function and gradient).
30
fmin_ncg --- Line-search Newton Conjugate Gradient (can use function,
31
gradient and hessian).
34
brute --- Perform a brute force search for the minimum
35
with final optimization if desired.
39
brent --- Use Brent's method (does not need inital guess)
40
fminbound --- Bounded minimization for scalar functions on an
41
interval using Brent's parabolic/golden_mean method.
42
golden -- Use Golden Section method (does not need initial guess)
44
bracket --- Find a bracket containing the minimum.
51
18
__all__ = ['fmin', 'fmin_powell','fmin_bfgs', 'fmin_ncg', 'fmin_cg',
52
19
'fminbound','brent', 'golden','bracket','rosen','rosen_der',
53
20
'rosen_hess', 'rosen_hess_prod', 'brute', 'approx_fprime',
54
21
'line_search', 'check_grad']
58
from scipy_base import atleast_1d, eye, mgrid, argmin, zeros, shape, \
59
squeeze, isscalar, vectorize, asarray, absolute, sqrt, Inf
24
from numpy import atleast_1d, eye, mgrid, argmin, zeros, shape, \
25
squeeze, isscalar, vectorize, asarray, absolute, sqrt, Inf, asfarray
28
# These have been copied from Numeric's MLab.py
29
# I don't think they made the transition to scipy_core
31
"""max(m,axis=0) returns the maximum of m along dimension axis.
34
return numpy.maximum.reduce(m,axis)
37
"""min(m,axis=0) returns the minimum of m along dimension axis.
40
return numpy.minimum.reduce(m,axis)
67
44
pymin = __builtin__.min
68
45
pymax = __builtin__.max
70
_epsilon = sqrt(scipy_base.limits.double_epsilon)
47
_epsilon = sqrt(numpy.finfo(float).eps)
72
49
def vecnorm(x, ord=2):
74
return scipy_base.amax(abs(x))
51
return numpy.amax(abs(x))
76
return scipy_base.amin(abs(x))
53
return numpy.amin(abs(x))
78
return scipy_base.sum(abs(x)**ord)**(1.0/ord)
55
return numpy.sum(abs(x)**ord,axis=0)**(1.0/ord)
80
57
def rosen(x): # The Rosenbrock function
82
return MLab.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
59
return numpy.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0,axis=0)
89
der = MLab.zeros(x.shape,x.typecode())
66
der = numpy.zeros_like(x)
90
67
der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)
91
68
der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])
92
69
der[-1] = 200*(x[-1]-x[-2]**2)
97
H = MLab.diag(-400*x[:-1],1) - MLab.diag(400*x[:-1],-1)
98
diagonal = Num.zeros(len(x),x.typecode())
74
H = numpy.diag(-400*x[:-1],1) - numpy.diag(400*x[:-1],-1)
75
diagonal = numpy.zeros(len(x), dtype=x.dtype)
99
76
diagonal[0] = 1200*x[0]-400*x[1]+2
100
77
diagonal[-1] = 200
101
78
diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
102
H = H + MLab.diag(diagonal)
79
H = H + numpy.diag(diagonal)
105
82
def rosen_hess_prod(x,p):
107
Hp = Num.zeros(len(x),x.typecode())
84
Hp = numpy.zeros(len(x), dtype=x.dtype)
108
85
Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
109
86
Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] \
110
87
-400*x[1:-1]*p[2:]
111
88
Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
114
def fmin(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None,
115
full_output=0, disp=1, retall=0):
91
def wrap_function(function, args):
93
def function_wrapper(x):
95
return function(x, *args)
96
return ncalls, function_wrapper
98
def fmin(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None,
99
full_output=0, disp=1, retall=0, callback=None):
116
100
"""Minimize a function using the downhill simplex algorithm.
120
104
Uses a Nelder-Mead simplex algorithm to find the minimum of function
121
105
of one or more variables.
147
134
full_output -- non-zero if fval and warnflag outputs are desired.
148
135
disp -- non-zero to print convergence messages.
149
136
retall -- non-zero to return list of solutions at each iteration
140
fmin, fmin_powell, fmin_cg,
141
fmin_bfgs, fmin_ncg -- multivariate local optimizers
142
leastsq -- nonlinear least squares minimizer
144
fmin_l_bfgs_b, fmin_tnc,
145
fmin_cobyla -- constrained multivariate optimizers
147
anneal, brute -- global optimizers
149
fminbound, brent, golden, bracket -- local scalar minimizers
151
fsolve -- n-dimenstional root-finding
153
brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
155
fixed_point -- scalar fixed-point finder
158
fcalls, func = wrap_function(func, args)
154
161
rank = len(x0.shape)
155
162
if not -1 < rank < 2:
163
170
one2np1 = range(1,N+1)
166
sim = Num.zeros((N+1,),x0.typecode())
168
sim = Num.zeros((N+1,N),x0.typecode())
169
fsim = Num.zeros((N+1,),'d')
173
sim = numpy.zeros((N+1,), dtype=x0.dtype)
175
sim = numpy.zeros((N+1,N), dtype=x0.dtype)
176
fsim = numpy.zeros((N+1,), float)
172
179
allvecs = [sim[0]]
173
fsim[0] = apply(func,(x0,)+args)
176
183
for k in range(0,N):
177
y = Num.array(x0,copy=1)
184
y = numpy.array(x0,copy=True)
179
186
y[k] = (1+nonzdelt)*y[k]
184
f = apply(func,(y,)+args)
187
ind = Num.argsort(fsim)
188
fsim = Num.take(fsim,ind) # sort so sim[0,:] has the lowest function value
189
sim = Num.take(sim,ind,0)
194
ind = numpy.argsort(fsim)
195
fsim = numpy.take(fsim,ind,0)
196
# sort so sim[0,:] has the lowest function value
197
sim = numpy.take(sim,ind,0)
194
while (funcalls < maxfun and iterations < maxiter):
195
if (max(Num.ravel(abs(sim[1:]-sim[0]))) <= xtol \
201
while (fcalls[0] < maxfun and iterations < maxiter):
202
if (max(numpy.ravel(abs(sim[1:]-sim[0]))) <= xtol \
196
203
and max(abs(fsim[0]-fsim[1:])) <= ftol):
199
xbar = Num.add.reduce(sim[:-1],0) / N
206
xbar = numpy.add.reduce(sim[:-1],0) / N
200
207
xr = (1+rho)*xbar - rho*sim[-1]
201
fxr = apply(func,(xr,)+args)
202
funcalls = funcalls + 1
205
211
if fxr < fsim[0]:
206
212
xe = (1+rho*chi)*xbar - rho*chi*sim[-1]
207
fxe = apply(func,(xe,)+args)
208
funcalls = funcalls + 1
245
248
for j in one2np1:
246
249
sim[j] = sim[0] + sigma*(sim[j] - sim[0])
247
fsim[j] = apply(func,(sim[j],)+args)
248
funcalls = funcalls + N
250
fsim[j] = func(sim[j])
250
ind = Num.argsort(fsim)
251
sim = Num.take(sim,ind,0)
252
fsim = Num.take(fsim,ind)
253
iterations = iterations + 1
252
ind = numpy.argsort(fsim)
253
sim = numpy.take(sim,ind,0)
254
fsim = numpy.take(fsim,ind,0)
255
if callback is not None:
255
allvecs.append(sim[0])
259
allvecs.append(sim[0])
261
if funcalls >= maxfun:
265
if fcalls[0] >= maxfun:
264
268
print "Warning: Maximum number of function evaluations has "\
272
276
print "Optimization terminated successfully."
273
277
print " Current function value: %f" % fval
274
278
print " Iterations: %d" % iterations
275
print " Function evaluations: %d" % funcalls
279
print " Function evaluations: %d" % fcalls[0]
279
retlist = x, fval, iterations, funcalls, warnflag
283
retlist = x, fval, iterations, fcalls[0], warnflag
281
285
retlist += (allvecs,)
285
289
retlist = (x, allvecs)
394
398
val_star = phi_aj
395
399
valprime_star = None
397
401
return a_star, val_star, valprime_star
399
403
def line_search(f, myfprime, xk, pk, gfk, old_fval, old_old_fval,
400
404
args=(), c1=1e-4, c2=0.9, amax=50):
401
"""Find alpha that satisfies strong Wolfe conditions.
403
Uses the line search algorithm to enforce strong Wolfe conditions
405
"""Find alpha that satisfies strong Wolfe conditions.
407
Uses the line search algorithm to enforce strong Wolfe conditions
404
408
Wright and Nocedal, 'Numerical Optimization', 1999, pg. 59-60
406
For the zoom phase it uses an algorithm by
410
For the zoom phase it uses an algorithm by
407
411
Outputs: (alpha0, gc, fc)
424
428
fprime = myfprime[0]
425
429
newargs = (f,eps) + args
426
430
_ls_ingfk = fprime(xk+alpha*pk,*newargs) # store for later use
427
return Num.dot(_ls_ingfk,pk)
431
return numpy.dot(_ls_ingfk,pk)
429
433
fprime = myfprime
430
434
def phiprime(alpha):
431
435
global _ls_gc, _ls_ingfk
433
437
_ls_ingfk = fprime(xk+alpha*pk,*args) # store for later use
434
return Num.dot(_ls_ingfk,pk)
438
return numpy.dot(_ls_ingfk,pk)
438
derphi0 = Num.dot(gfk,pk)
442
derphi0 = numpy.dot(gfk,pk)
440
444
alpha1 = pymin(1.0,1.01*2*(phi0-old_old_fval)/derphi0)
447
# This shouldn't happen. Perhaps the increment has slipped below
448
# machine precision? For now, set the return variables skip the
449
# useless while loop, and raise warnflag=2 due to possible imprecision.
452
old_fval = old_old_fval
441
455
phi_a1 = phi(alpha1)
442
456
#derphi_a1 = phiprime(alpha1) evaluated below
490
506
# this is the gradient at the next step
491
507
# no need to compute it again in the outer loop.
492
508
fprime_star = _ls_ingfk
494
510
return alpha_star, _ls_fc, _ls_gc, fval_star, old_fval, fprime_star
497
513
def line_search_BFGS(f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=1):
498
514
"""Minimize over alpha, the function f(xk+alpha pk)
500
516
Uses the interpolation algorithm (Armiijo backtracking) as suggested by
501
517
Wright and Nocedal in 'Numerical Optimization', 1999, pg. 56-57
503
519
Outputs: (alpha, fc, gc)
554
570
def approx_fprime(xk,f,epsilon,*args):
555
571
f0 = apply(f,(xk,)+args)
556
grad = Num.zeros((len(xk),),'d')
557
ei = Num.zeros((len(xk),),'d')
572
grad = numpy.zeros((len(xk),), float)
573
ei = numpy.zeros((len(xk),), float)
558
574
for k in range(len(xk)):
560
576
grad[k] = (apply(f,(xk+ei,)+args) - f0)/epsilon
592
608
norm -- order of norm (Inf is max, -Inf is min)
593
609
epsilon -- if fprime is approximated use this value for
594
610
the step size (can be scalar or vector)
611
callback -- an optional user-supplied function to call after each
612
iteration. It is called as callback(xk), where xk is the
613
current parameter vector.
596
615
Outputs: (xopt, {fopt, gopt, Hopt, func_calls, grad_calls, warnflag}, <allvecs>)
614
633
and warnflag in addition to xopt.
615
634
disp -- print convergence message if non-zero.
616
635
retall -- return a list of results at each iteration if non-zero
639
fmin, fmin_powell, fmin_cg,
640
fmin_bfgs, fmin_ncg -- multivariate local optimizers
641
leastsq -- nonlinear least squares minimizer
643
fmin_l_bfgs_b, fmin_tnc,
644
fmin_cobyla -- constrained multivariate optimizers
646
anneal, brute -- global optimizers
648
fminbound, brent, golden, bracket -- local scalar minimizers
650
fsolve -- n-dimenstional root-finding
652
brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
654
fixed_point -- scalar fixed-point finder
623
658
if maxiter is None:
624
659
maxiter = len(x0)*200
660
func_calls, f = wrap_function(f, args)
662
grad_calls, myfprime = wrap_function(approx_fprime, (f, epsilon))
664
grad_calls, myfprime = wrap_function(fprime, args)
668
I = numpy.eye(N,dtype=int)
631
old_fval = f(x0,*args)
632
671
old_old_fval = old_fval + 5000
635
gfk = apply(approx_fprime,(x0,f,epsilon)+args)
636
myfprime = (approx_fprime,epsilon)
637
func_calls = func_calls + len(x0) + 1
639
gfk = apply(fprime,(x0,)+args)
641
grad_calls = grad_calls + 1
647
677
gnorm = vecnorm(gfk,ord=norm)
648
678
while (gnorm > gtol) and (k < maxiter):
649
pk = -Num.dot(Hk,gfk)
679
pk = -numpy.dot(Hk,gfk)
650
680
alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
651
681
linesearch.line_search(f,myfprime,xk,pk,gfk,
652
old_fval,old_old_fval,args=args)
682
old_fval,old_old_fval)
653
683
if alpha_k is None: # line search failed try different one.
654
func_calls = func_calls + fc
655
grad_calls = grad_calls + gc
656
684
alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
657
685
line_search(f,myfprime,xk,pk,gfk,
658
old_fval,old_old_fval,args=args)
660
func_calls = func_calls + fc
661
grad_calls = grad_calls + gc
686
old_fval,old_old_fval)
688
# This line search also failed to find a better solution.
662
691
xkp1 = xk + alpha_k * pk
664
693
allvecs.append(xkp1)
667
696
if gfkp1 is None:
669
gfkp1 = apply(approx_fprime,(xkp1,f,epsilon)+args)
670
func_calls = func_calls + len(x0) + 1
672
gfkp1 = apply(fprime,(xkp1,)+args)
673
grad_calls = grad_calls + 1
697
gfkp1 = myfprime(xkp1)
701
if callback is not None:
678
704
gnorm = vecnorm(gfk,ord=norm)
679
705
if (gnorm <= gtol):
683
rhok = 1 / Num.dot(yk,sk)
709
rhok = 1 / (numpy.dot(yk,sk))
684
710
except ZeroDivisionError:
687
#print "Divide by zero encountered: Hessian calculation reset."
690
A1 = I - sk[:,Num.NewAxis] * yk[Num.NewAxis,:] * rhok
691
A2 = I - yk[:,Num.NewAxis] * sk[Num.NewAxis,:] * rhok
692
Hk = Num.dot(A1,Num.dot(Hk,A2)) + rhok * sk[:,Num.NewAxis] \
712
print "Divide-by-zero encountered: rhok assumed large"
713
A1 = I - sk[:,numpy.newaxis] * yk[numpy.newaxis,:] * rhok
714
A2 = I - yk[:,numpy.newaxis] * sk[numpy.newaxis,:] * rhok
715
Hk = numpy.dot(A1,numpy.dot(Hk,A2)) + rhok * sk[:,numpy.newaxis] \
716
* sk[numpy.newaxis,:]
695
718
if disp or full_output:
699
722
print "Warning: Desired error not necessarily achieved due to precision loss"
700
723
print " Current function value: %f" % fval
701
724
print " Iterations: %d" % k
702
print " Function evaluations: %d" % func_calls
703
print " Gradient evaluations: %d" % grad_calls
725
print " Function evaluations: %d" % func_calls[0]
726
print " Gradient evaluations: %d" % grad_calls[0]
705
728
elif k >= maxiter:
708
731
print "Warning: Maximum number of iterations has been exceeded"
709
732
print " Current function value: %f" % fval
710
733
print " Iterations: %d" % k
711
print " Function evaluations: %d" % func_calls
712
print " Gradient evaluations: %d" % grad_calls
734
print " Function evaluations: %d" % func_calls[0]
735
print " Gradient evaluations: %d" % grad_calls[0]
715
738
print "Optimization terminated successfully."
716
739
print " Current function value: %f" % fval
717
740
print " Iterations: %d" % k
718
print " Function evaluations: %d" % func_calls
719
print " Gradient evaluations: %d" % grad_calls
741
print " Function evaluations: %d" % func_calls[0]
742
print " Gradient evaluations: %d" % grad_calls[0]
722
retlist = xk, fval, gfk, Hk, func_calls, grad_calls, warnflag
745
retlist = xk, fval, gfk, Hk, func_calls[0], grad_calls[0], warnflag
724
747
retlist += (allvecs,)
728
751
retlist = (xk, allvecs)
771
797
and warnflag in addition to xopt.
772
798
disp -- print convergence message if non-zero.
773
799
retall -- return a list of results at each iteration if True
803
fmin, fmin_powell, fmin_cg,
804
fmin_bfgs, fmin_ncg -- multivariate local optimizers
805
leastsq -- nonlinear least squares minimizer
807
fmin_l_bfgs_b, fmin_tnc,
808
fmin_cobyla -- constrained multivariate optimizers
810
anneal, brute -- global optimizers
812
fminbound, brent, golden, bracket -- local scalar minimizers
814
fsolve -- n-dimenstional root-finding
816
brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
818
fixed_point -- scalar fixed-point finder
780
822
if maxiter is None:
781
823
maxiter = len(x0)*200
824
func_calls, f = wrap_function(f, args)
826
grad_calls, myfprime = wrap_function(approx_fprime, (f, epsilon))
828
grad_calls, myfprime = wrap_function(fprime, args)
787
old_fval = f(xk,*args)
788
834
old_old_fval = old_fval + 5000
792
gfk = apply(approx_fprime,(x0,f,epsilon)+args)
793
myfprime = (approx_fprime,epsilon)
794
func_calls = func_calls + len(x0) + 1
796
gfk = apply(fprime,(x0,)+args)
798
grad_calls = grad_calls + 1
804
841
gnorm = vecnorm(gfk,ord=norm)
805
842
while (gnorm > gtol) and (k < maxiter):
806
deltak = Num.dot(gfk,gfk)
843
deltak = numpy.dot(gfk,gfk)
845
# These values are modified by the line search, even if it fails
846
old_fval_backup = old_fval
847
old_old_fval_backup = old_old_fval
807
849
alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
808
850
linesearch.line_search(f,myfprime,xk,pk,gfk,old_fval,
809
old_old_fval,args=args,c2=0.4)
810
852
if alpha_k is None: # line search failed -- use different one.
813
853
alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
814
854
line_search(f,myfprime,xk,pk,gfk,
815
old_fval,old_old_fval,args=args)
855
old_fval_backup,old_old_fval_backup)
856
if alpha_k is None or alpha_k == 0:
857
# This line search also failed to find a better solution.
818
860
xk = xk + alpha_k*pk
820
862
allvecs.append(xk)
821
863
if gfkp1 is None:
823
gfkp1 = apply(approx_fprime,(xk,f,epsilon)+args)
824
func_calls = func_calls + len(x0) + 1
826
gfkp1 = apply(fprime,(xk,)+args)
827
grad_calls = grad_calls + 1
830
beta_k = pymax(0,Num.dot(yk,gfkp1)/deltak)
866
beta_k = pymax(0,numpy.dot(yk,gfkp1)/deltak)
831
867
pk = -gfkp1 + beta_k * pk
833
869
gnorm = vecnorm(gfk,ord=norm)
870
if callback is not None:
837
875
if disp or full_output:
839
877
if warnflag == 2:
841
879
print "Warning: Desired error not necessarily achieved due to precision loss"
842
880
print " Current function value: %f" % fval
843
881
print " Iterations: %d" % k
844
print " Function evaluations: %d" % func_calls
845
print " Gradient evaluations: %d" % grad_calls
882
print " Function evaluations: %d" % func_calls[0]
883
print " Gradient evaluations: %d" % grad_calls[0]
847
885
elif k >= maxiter:
850
888
print "Warning: Maximum number of iterations has been exceeded"
851
889
print " Current function value: %f" % fval
852
890
print " Iterations: %d" % k
853
print " Function evaluations: %d" % func_calls
854
print " Gradient evaluations: %d" % grad_calls
891
print " Function evaluations: %d" % func_calls[0]
892
print " Gradient evaluations: %d" % grad_calls[0]
857
895
print "Optimization terminated successfully."
858
896
print " Current function value: %f" % fval
859
897
print " Iterations: %d" % k
860
print " Function evaluations: %d" % func_calls
861
print " Gradient evaluations: %d" % grad_calls
898
print " Function evaluations: %d" % func_calls[0]
899
print " Gradient evaluations: %d" % grad_calls[0]
865
retlist = xk, fval, func_calls, grad_calls, warnflag
903
retlist = xk, fval, func_calls[0], grad_calls[0], warnflag
867
905
retlist += (allvecs,)
871
909
retlist = (xk, allvecs)
875
913
def fmin_ncg(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5,
876
epsilon=_epsilon, maxiter=None, full_output=0, disp=1, retall=0):
914
epsilon=_epsilon, maxiter=None, full_output=0, disp=1, retall=0,
879
918
Minimize the function, f, whose gradient is given by fprime using the
896
935
epsilon -- if fhess is approximated use this value for
897
936
the step size (can be scalar or vector)
937
callback -- an optional user-supplied function to call after each
938
iteration. It is called as callback(xk), where xk is the
939
current parameter vector.
899
941
Outputs: (xopt, {fopt, fcalls, gcalls, hcalls, warnflag},{allvecs})
901
943
xopt -- the minimizer of f
903
945
fopt -- the value of the function at xopt: fopt = f(xopt)
904
946
fcalls -- the number of function calls.
905
947
gcalls -- the number of gradient calls.
924
966
provided, then the hessian product will be approximated using finite
925
967
differences on fprime.
971
fmin, fmin_powell, fmin_cg,
972
fmin_bfgs, fmin_ncg -- multivariate local optimizers
973
leastsq -- nonlinear least squares minimizer
975
fmin_l_bfgs_b, fmin_tnc,
976
fmin_cobyla -- constrained multivariate optimizers
978
anneal, brute -- global optimizers
980
fminbound, brent, golden, bracket -- local scalar minimizers
982
fsolve -- n-dimenstional root-finding
984
brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
986
fixed_point -- scalar fixed-point finder
990
fcalls, f = wrap_function(f, args)
991
gcalls, fprime = wrap_function(fprime, args)
933
993
if maxiter is None:
934
994
maxiter = len(x0)*200
936
996
xtol = len(x0)*avextol
937
997
update = [2*xtol]
942
old_fval = f(x0,*args)
944
while (Num.add.reduce(abs(update)) > xtol) and (k < maxiter):
1003
while (numpy.add.reduce(abs(update)) > xtol) and (k < maxiter):
945
1004
# Compute a search direction pk by applying the CG method to
946
1005
# del2 f(xk) p = - grad f(xk) starting from 0.
947
b = -apply(fprime,(xk,)+args)
949
maggrad = Num.add.reduce(abs(b))
950
eta = min([0.5,Num.sqrt(maggrad)])
1007
maggrad = numpy.add.reduce(abs(b))
1008
eta = min([0.5,numpy.sqrt(maggrad)])
951
1009
termcond = eta * maggrad
1010
xsupi = zeros(len(x0), dtype=x0.dtype)
956
dri0 = Num.dot(ri,ri)
1014
dri0 = numpy.dot(ri,ri)
958
1016
if fhess is not None: # you want to compute hessian once.
959
1017
A = apply(fhess,(xk,)+args)
960
1018
hcalls = hcalls + 1
962
while Num.add.reduce(abs(ri)) > termcond:
1020
while numpy.add.reduce(abs(ri)) > termcond:
963
1021
if fhess is None:
964
1022
if fhess_p is None:
965
Ap = apply(approx_fhess_p,(xk,psupi,fprime,epsilon)+args)
1023
Ap = approx_fhess_p(xk,psupi,fprime,epsilon)
968
Ap = apply(fhess_p,(xk,psupi)+args)
1025
Ap = fhess_p(xk,psupi, *args)
969
1026
hcalls = hcalls + 1
971
Ap = Num.dot(A,psupi)
1028
Ap = numpy.dot(A,psupi)
972
1029
# check curvature
973
curv = Num.dot(psupi,Ap)
1030
Ap = asarray(Ap).squeeze() # get rid of matrices...
1031
curv = numpy.dot(psupi,Ap)
980
1040
alphai = dri0 / curv
981
1041
xsupi = xsupi + alphai * psupi
982
1042
ri = ri + alphai * Ap
983
dri1 = Num.dot(ri,ri)
1043
dri1 = numpy.dot(ri,ri)
984
1044
betai = dri1 / dri0
985
1045
psupi = -ri + betai * psupi
987
dri0 = dri1 # update Num.dot(ri,ri) for next time.
1047
dri0 = dri1 # update numpy.dot(ri,ri) for next time.
989
1049
pk = xsupi # search direction is solution to system.
990
1050
gfk = -b # gradient at xk
991
alphak, fc, gc, old_fval = line_search_BFGS(f,xk,pk,gfk,old_fval,args)
1051
alphak, fc, gc, old_fval = line_search_BFGS(f,xk,pk,gfk,old_fval)
995
1053
update = alphak * pk
1054
xk = xk + update # upcast if necessary
1055
if callback is not None:
998
1058
allvecs.append(xk)
1001
1061
if disp or full_output:
1002
1062
fval = old_fval
1015
1075
print "Optimization terminated successfully."
1016
1076
print " Current function value: %f" % fval
1017
1077
print " Iterations: %d" % k
1018
print " Function evaluations: %d" % fcalls
1019
print " Gradient evaluations: %d" % gcalls
1078
print " Function evaluations: %d" % fcalls[0]
1079
print " Gradient evaluations: %d" % gcalls[0]
1020
1080
print " Hessian evaluations: %d" % hcalls
1022
1082
if full_output:
1023
retlist = xk, fval, fcalls, gcalls, hcalls, warnflag
1083
retlist = xk, fval, fcalls[0], gcalls[0], hcalls, warnflag
1025
1085
retlist += (allvecs,)
1029
1089
retlist = (xk, allvecs)
1034
def fminbound(func, x1, x2, args=(), xtol=1e-5, maxfun=500,
1094
def fminbound(func, x1, x2, args=(), xtol=1e-5, maxfun=500,
1035
1095
full_output=0, disp=1):
1036
1096
"""Bounded minimization for scalar functions.
1048
1108
args -- extra arguments to pass to function.
1049
1109
xtol -- the convergence tolerance.
1050
1110
maxfun -- maximum function evaluations.
1051
full_output -- Non-zero to return optional outputs.
1111
full_output -- Non-zero to return optional outputs.
1052
1112
disp -- Non-zero to print messages.
1053
1113
0 : no message printing.
1054
1114
1 : non-convergence notification messages only.
1055
1115
2 : print a message on convergence too.
1056
3 : print iteration results.
1116
3 : print iteration results.
1059
1119
Outputs: (xopt, {fval, ierr, numfunc})
1063
1123
ierr -- An error flag (0 if converged, 1 if maximum number of
1064
1124
function calls reached).
1065
1125
numfunc -- The number of function calls.
1129
fmin, fmin_powell, fmin_cg,
1130
fmin_bfgs, fmin_ncg -- multivariate local optimizers
1131
leastsq -- nonlinear least squares minimizer
1133
fmin_l_bfgs_b, fmin_tnc,
1134
fmin_cobyla -- constrained multivariate optimizers
1136
anneal, brute -- global optimizers
1138
fminbound, brent, golden, bracket -- local scalar minimizers
1140
fsolve -- n-dimenstional root-finding
1142
brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
1144
fixed_point -- scalar fixed-point finder
1189
1269
Uses inverse parabolic interpolation when possible to speed up convergence
1190
1270
of golden section method.
1275
fmin, fmin_powell, fmin_cg,
1276
fmin_bfgs, fmin_ncg -- multivariate local optimizers
1277
leastsq -- nonlinear least squares minimizer
1279
fmin_l_bfgs_b, fmin_tnc,
1280
fmin_cobyla -- constrained multivariate optimizers
1282
anneal, brute -- global optimizers
1284
fminbound, brent, golden, bracket -- local scalar minimizers
1286
fsolve -- n-dimenstional root-finding
1288
brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
1290
fixed_point -- scalar fixed-point finder
1193
1293
_mintol = 1.0e-11
1194
1294
_cg = 0.3819660
1289
1389
Uses analog of bisection method to decrease the bracketed interval.
1393
fmin, fmin_powell, fmin_cg,
1394
fmin_bfgs, fmin_ncg -- multivariate local optimizers
1395
leastsq -- nonlinear least squares minimizer
1397
fmin_l_bfgs_b, fmin_tnc,
1398
fmin_cobyla -- constrained multivariate optimizers
1400
anneal, brute -- global optimizers
1402
fminbound, brent, golden, bracket -- local scalar minimizers
1404
fsolve -- n-dimenstional root-finding
1406
brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
1408
fixed_point -- scalar fixed-point finder
1291
1411
if brack is None:
1292
1412
xa,xb,xc,fa,fb,fc,funcalls = bracket(func, args=args)
1335
1455
if full_output:
1336
1456
return xmin, fval, funcalls
1341
def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0):
1461
def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
1342
1462
"""Given a function and distinct initial points, search in the downhill
1343
1463
direction (as defined by the initital points) and return new points
1344
1464
xa, xb, xc that bracket the minimum of the function:
1397
1518
xa=xb; xb=xc; xc=w
1398
1519
fa=fb; fb=fc; fc=fw
1399
1520
return xa, xb, xc, fa, fb, fc, funcalls
1402
global _powell_funcalls
1404
def _myfunc(alpha, func, x0, direc, args=()):
1405
funcargs = (x0 + alpha * direc,)+args
1406
return func(*funcargs)
1408
def _linesearch_powell(func, p, xi, args=(), tol=1e-3):
1524
def _linesearch_powell(func, p, xi, tol=1e-3):
1409
1525
# line-search algorithm using fminbound
1410
1526
# find the minimium of the function
1411
1527
# func(x0+ alpha*direc)
1412
global _powell_funcalls
1413
extra_args = (func, p, xi) + args
1414
alpha_min, fret, iter, num = brent(_myfunc, args=extra_args,
1415
full_output=1, tol=tol)
1529
return func(p + alpha * xi)
1530
alpha_min, fret, iter, num = brent(myfunc, full_output=1, tol=tol)
1416
1531
xi = alpha_min*xi
1417
_powell_funcalls += num
1418
1532
return squeeze(fret), p+xi, xi
1421
1535
def fmin_powell(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None,
1422
maxfun=None, full_output=0, disp=1, retall=0):
1536
maxfun=None, full_output=0, disp=1, retall=0, callback=None):
1423
1537
"""Minimize a function using modified Powell's method.
1427
1541
Uses a modification of Powell's method to find the minimum of a function
1456
1573
disp -- non-zero to print convergence messages.
1457
1574
retall -- non-zero to return a list of the solution at each iteration
1578
fmin, fmin_powell, fmin_cg,
1579
fmin_bfgs, fmin_ncg -- multivariate local optimizers
1580
leastsq -- nonlinear least squares minimizer
1582
fmin_l_bfgs_b, fmin_tnc,
1583
fmin_cobyla -- constrained multivariate optimizers
1585
anneal, brute -- global optimizers
1587
fminbound, brent, golden, bracket -- local scalar minimizers
1589
fsolve -- n-dimenstional root-finding
1591
brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
1593
fixed_point -- scalar fixed-point finder
1460
global _powell_funcalls
1596
# we need to use a mutable object here that we can update in the
1598
fcalls, func = wrap_function(func, args)
1461
1599
x = asarray(x0)
1470
1608
if maxfun is None:
1471
1609
maxfun = N * 1000
1473
direc = eye(N,typecode='d')
1474
fval = squeeze(apply(func, (x,)+args))
1475
_powell_funcalls = 1
1611
direc = eye(N,dtype=float)
1612
fval = squeeze(func(x))
1478
1615
ilist = range(N)
1483
1620
for i in ilist:
1484
1621
direc1 = direc[i]
1486
fval, x, direc1 = _linesearch_powell(func, x, direc1, args=args, tol=xtol*100)
1623
fval, x, direc1 = _linesearch_powell(func, x, direc1, tol=xtol*100)
1487
1624
if (fx2 - fval) > delta:
1488
1625
delta = fx2 - fval
1628
if callback is not None:
1492
1631
allvecs.append(x)
1493
1632
if (2.0*(fx - fval) <= ftol*(abs(fx)+abs(fval))+1e-20): break
1494
if _powell_funcalls >= maxfun: break
1633
if fcalls[0] >= maxfun: break
1495
1634
if iter >= maxiter: break
1497
1636
# Construct the extrapolated point
1498
1637
direc1 = x - x1
1501
fx2 = squeeze(apply(func, (x2,)+args))
1502
_powell_funcalls +=1
1640
fx2 = squeeze(func(x2))
1505
1643
t = 2.0*(fx+fx2-2.0*fval)
1527
1665
print "Optimization terminated successfully."
1528
1666
print " Current function value: %f" % fval
1529
1667
print " Iterations: %d" % iter
1530
print " Function evaluations: %d" % _powell_funcalls
1668
print " Function evaluations: %d" % fcalls[0]
1534
1672
if full_output:
1535
retlist = x, fval, direc, iter, _powell_funcalls, warnflag
1673
retlist = x, fval, direc, iter, fcalls[0], warnflag
1537
1675
retlist += (allvecs,)
1541
1679
retlist = (x, allvecs)
1548
1686
def _endprint(x, flag, fval, maxfun, xtol, disp):
1578
1716
grid -- tuple with same length as x0 representing the
1579
1717
evaluation grid
1580
1718
Jout -- Function values over grid: Jout = func(*grid)
1722
fmin, fmin_powell, fmin_cg,
1723
fmin_bfgs, fmin_ncg -- multivariate local optimizers
1724
leastsq -- nonlinear least squares minimizer
1726
fmin_l_bfgs_b, fmin_tnc,
1727
fmin_cobyla -- constrained multivariate optimizers
1729
anneal, brute -- global optimizers
1731
fminbound, brent, golden, bracket -- local scalar minimizers
1733
fsolve -- n-dimenstional root-finding
1735
brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
1737
fixed_point -- scalar fixed-point finder
1582
1740
N = len(ranges)
added
removed
Lib/optimize/optimize.py
