1
2
3
4
5
6
7
8
9
10
11
12
13 """optimize.py
14
15 A collection of general-purpose optimization routines using Numeric
16
17
18 N-D Algorithms
19 ===============
20 fmin --- Nelder-Mead Simplex algorithm (uses only function calls).
21 fmin_powell --- Powell (direction set) method (uses only function calls).
22 fmin_bfgs --- Quasi-Newton method (uses function and gradient).
23 fmin_ncg --- Line-search Newton Conjugate Gradient (uses function,
24 gradient and hessian (if it's provided)).
25
26 1-D Algorithms
27 ===============
28 brent --- Use Brent's method (does not need inital guess)
29 fminbound --- Bounded minimization for scalar functions on an
30 interval using Brent's parabolic/golden_mean method.
31 golden -- Use Golden Section method (does not need initial guess)
32
33 bracket --- Find a bracket containing the minimum.
34
35 """
36
37
38 __all__ = ['fmin', 'fmin_powell','fmin_bfgs', 'fmin_ncg', 'fminbound', 'brent',
39 'golden','bracket','rosen','rosen_der', 'rosen_hess','rosen_hess_prod']
40
41 from __future__ import nested_scopes
42 import Numeric
43 import MLab
44 import sys
45
46 from Numeric import absolute, sqrt, asarray
47 from MLab import squeeze
48 Num = Numeric
49 max = MLab.max
50 min = MLab.min
51 abs = absolute
52 __version__="0.4.2"
53
55 return MLab.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
56
58 xm = x[1:-1]
59 xm_m1 = x[:-2]
60 xm_p1 = x[2:]
61 der = MLab.zeros(x.shape,x.typecode())
62 der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)
63 der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])
64 der[-1] = 200*(x[-1]-x[-2]**2)
65 return der
66
68 x = atleast_1d(x)
69 H = MLab.diag(-400*x[:-1],1) - MLab.diag(400*x[:-1],-1)
70 diagonal = Num.zeros(len(x),x.typecode())
71 diagonal[0] = 1200*x[0]-400*x[1]+2
72 diagonal[-1] = 200
73 diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
74 H = H + MLab.diag(diagonal)
75 return H
76
78 x = atleast_1d(x)
79 Hp = Num.zeros(len(x),x.typecode())
80 Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
81 Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] \
82 -400*x[1:-1]*p[2:]
83 Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
84 return Hp
85
86
87 -def fmin(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None,
88 full_output=0, disp=1):
89 """Minimize a function using the downhill simplex algorithm.
90
91 Description:
92
93 Uses a Nelder-Mead simplex algorithm to find the minimum of function
94 of one or more variables.
95
96 Inputs:
97
98 func -- the Python function or method to be minimized.
99 x0 -- the initial guess.
100 args -- extra arguments for func.
101
102 Outputs: (xopt, {fopt, iter, funcalls, warnflag})
103
104 xopt -- minimizer of function
105
106 fopt -- value of function at minimum: fopt = func(xopt)
107 iter -- number of iterations
108 funcalls -- number of function calls
109 warnflag -- Integer warning flag:
110 1 : 'Maximum number of function evaluations.'
111 2 : 'Maximum number of iterations.'
112
113 Additional Inputs:
114
115 xtol -- acceptable relative error in xopt for convergence.
116 ftol -- acceptable relative error in func(xopt) for convergence.
117 maxiter -- the maximum number of iterations to perform.
118 maxfun -- the maximum number of function evaluations.
119 full_output -- non-zero if fval and warnflag outputs are desired.
120 disp -- non-zero to print convergence messages.
121
122 """
123 x0 = asarray(x0)
124 N = len(x0)
125 print "n: " + `N`
126 print "m: " + `N+1`
127 rank = len(x0.shape)
128 if not -1 < rank < 2:
129 raise ValueError, "Initial guess must be a scalar or rank-1 sequence."
130 if maxiter is None:
131 maxiter = N * 200
132 if maxfun is None:
133 maxfun = N * 200
134
135 rho = 1; chi = 2; psi = 0.5; sigma = 0.5;
136 one2np1 = range(1,N+1)
137
138 if rank == 0:
139 sim = Num.zeros((N+1,),x0.typecode())
140 else:
141 sim = Num.zeros((N+1,N),x0.typecode())
142 fsim = Num.zeros((N+1,),'d')
143 sim[0] = x0
144 fsim[0] = apply(func,(x0,)+args)
145 nonzdelt = 0.05
146 zdelt = 0.00025
147 for k in range(0,N):
148 y = Num.array(x0,copy=1)
149 if y[k] != 0:
150 y[k] = (1+nonzdelt)*y[k]
151 else:
152 y[k] = zdelt
153
154 sim[k+1] = y
155 f = apply(func,(y,)+args)
156 fsim[k+1] = f
157
158 ind = Num.argsort(fsim)
159 fsim = Num.take(fsim,ind)
160 sim = Num.take(sim,ind,0)
161 print "\n\nMinimisation run number: " + `0`
162 print "Initial simplex:\n" + `sim`
163 print "Initial function vector:\n" + `fsim`
164
165 iterations = 1
166 funcalls = N+1
167
168 while (funcalls < maxfun and iterations < maxiter):
169
170
171
172 if max(abs(fsim[0]-fsim[1:])) <= ftol:
173 break
174
175 print "\n\n< Minimisation run number: " + `iterations` + " >"
176 print "Function diff: " + `max(abs(fsim[0]-fsim[1:]))`
177 print "Function diff limit:" + `ftol`
178 print "Simplex:\n" + `sim`
179 print "Function vector:\n" + `fsim`
180 xbar = Num.add.reduce(sim[:-1],0) / N
181 print "Pivot point: " + `xbar`
182 xr = (1+rho)*xbar - rho*sim[-1]
183 fxr = apply(func,(xr,)+args)
184 funcalls = funcalls + 1
185 doshrink = 0
186 print "\nReflecting."
187 print "\t%-29s%-40s" % ("Reflection vector:", `xr`)
188 print "\t%-29s%-40s" % ("Reflection function value:", `fxr`)
189
190 if fxr < fsim[0]:
191 xe = (1+rho*chi)*xbar - rho*chi*sim[-1]
192 fxe = apply(func,(xe,)+args)
193 funcalls = funcalls + 1
194 print "\nExtending."
195 print "\t%-29s%-40s" % ("Extension vector:", `xe`)
196 print "\t%-29s%-40s" % ("Extension function value:", `fxe`)
197
198 if fxe < fxr:
199 sim[-1] = xe
200 fsim[-1] = fxe
201 else:
202 sim[-1] = xr
203 fsim[-1] = fxr
204 else:
205 if fxr < fsim[-2]:
206 sim[-1] = xr
207 fsim[-1] = fxr
208 else:
209
210 if fxr < fsim[-1]:
211 xc = (1+psi*rho)*xbar - psi*rho*sim[-1]
212 fxc = apply(func,(xc,)+args)
213 funcalls = funcalls + 1
214 print "\nContracting."
215 print "\t%-29s%-40s" % ("Contraction vector:", `xc`)
216 print "\t%-29s%-40s" % ("Contraction function value:", `fxc`)
217
218 if fxc <= fxr:
219 sim[-1] = xc
220 fsim[-1] = fxc
221 else:
222 doshrink=1
223 else:
224
225 xcc = (1-psi)*xbar + psi*sim[-1]
226 fxcc = apply(func,(xcc,)+args)
227 funcalls = funcalls + 1
228 print "\nOriginal contracting."
229 print "\t%-29s%-40s" % ("Contraction orig vector:", `xcc`)
230 print "\t%-29s%-40s" % ("Contraction orig function value:", `fxcc`)
231
232 if fxcc < fsim[-1]:
233 sim[-1] = xcc
234 fsim[-1] = fxcc
235 else:
236 doshrink = 1
237
238 if doshrink:
239 for j in one2np1:
240 sim[j] = sim[0] + sigma*(sim[j] - sim[0])
241 fsim[j] = apply(func,(sim[j],)+args)
242 funcalls = funcalls + N
243 print "\nShrinking."
244
245 ind = Num.argsort(fsim)
246 sim = Num.take(sim,ind,0)
247 fsim = Num.take(fsim,ind)
248 print "Final simplex:\n" + `sim`
249 print "Final function vector:\n" + `fsim`
250 iterations = iterations + 1
251
252 x = sim[0]
253 fval = min(fsim)
254 warnflag = 0
255
256 if funcalls >= maxfun:
257 warnflag = 1
258 if disp:
259 print "Warning: Maximum number of function evaluations has "\
260 "been exceeded."
261 elif iterations >= maxiter:
262 warnflag = 2
263 if disp:
264 print "Warning: Maximum number of iterations has been exceeded"
265 else:
266 if disp:
267 print "Optimization terminated successfully."
268 print " Current function value: %f" % fval
269 print " Iterations: %d" % iterations
270 print " Function evaluations: %d" % funcalls
271
272 if full_output:
273 return x, fval, iterations, funcalls, warnflag
274 else:
275 return x
276
277
278 -def zoom(a_lo, a_hi):
280
281
282 -def line_search(f, fprime, xk, pk, gfk, args=(), c1=1e-4, c2=0.9, amax=50):
283 """Minimize the function f(xk+alpha pk)
284
285 Uses the line search algorithm of
286 Wright and Nocedal in 'Numerical Optimization', 1999, pg. 59-60
287
288 Outputs: (alpha0, gc, fc)
289 """
290
291 fc = 0
292 gc = 0
293 alpha0 = 1.0
294 phi0 = apply(f,(xk,)+args)
295 phi_a0 = apply(f,(xk+alpha0*pk,)+args)
296 fc = fc + 2
297 derphi0 = Num.dot(gfk,pk)
298 derphi_a0 = Num.dot(apply(fprime,(xk+alpha0*pk,)+args),pk)
299 gc = gc + 1
300
301
302 if (phi_a0 <= phi0 + c1*alpha0*derphi0) \
303 and (abs(derphi_a0) <= c2*abs(derphi0)):
304 return alpha0, fc, gc
305
306 alpha0 = 0
307 alpha1 = 1
308 phi_a1 = phi_a0
309 phi_a0 = phi0
310
311 i = 1
312 while 1:
313 if (phi_a1 > phi0 + c1*alpha1*derphi0) or \
314 ((phi_a1 >= phi_a0) and (i > 1)):
315 return zoom(alpha0, alpha1)
316
317 derphi_a1 = Num.dot(apply(fprime,(xk+alpha1*pk,)+args),pk)
318 gc = gc + 1
319 if (abs(derphi_a1) <= -c2*derphi0):
320 return alpha1
321
322 if (derphi_a1 >= 0):
323 return zoom(alpha1, alpha0)
324
325 alpha2 = (amax-alpha1)*0.25 + alpha1
326 i = i + 1
327 alpha0 = alpha1
328 alpha1 = alpha2
329 phi_a0 = phi_a1
330 phi_a1 = apply(f,(xk+alpha1*pk,)+args)
331
332
333
335 """Minimize over alpha, the function f(xk+alpha pk)
336
337 Uses the interpolation algorithm (Armiijo backtracking) as suggested by
338 Wright and Nocedal in 'Numerical Optimization', 1999, pg. 56-57
339
340 Outputs: (alpha, fc, gc)
341 """
342
343 fc = 0
344 phi0 = apply(f,(xk,)+args)
345 phi_a0 = apply(f,(xk+alpha0*pk,)+args)
346 fc = fc + 2
347 derphi0 = Num.dot(gfk,pk)
348
349 if (phi_a0 <= phi0 + c1*alpha0*derphi0):
350 return alpha0, fc, 0
351
352
353
354 alpha1 = -(derphi0) * alpha0**2 / 2.0 / (phi_a0 - phi0 - derphi0 * alpha0)
355 phi_a1 = apply(f,(xk+alpha1*pk,)+args)
356 fc = fc + 1
357
358 if (phi_a1 <= phi0 + c1*alpha1*derphi0):
359 return alpha1, fc, 0
360
361
362
363
364
365
366 while 1:
367 factor = alpha0**2 * alpha1**2 * (alpha1-alpha0)
368 a = alpha0**2 * (phi_a1 - phi0 - derphi0*alpha1) - \
369 alpha1**2 * (phi_a0 - phi0 - derphi0*alpha0)
370 a = a / factor
371 b = -alpha0**3 * (phi_a1 - phi0 - derphi0*alpha1) + \
372 alpha1**3 * (phi_a0 - phi0 - derphi0*alpha0)
373 b = b / factor
374
375 alpha2 = (-b + Num.sqrt(abs(b**2 - 3 * a * derphi0))) / (3.0*a)
376 phi_a2 = apply(f,(xk+alpha2*pk,)+args)
377 fc = fc + 1
378
379 if (phi_a2 <= phi0 + c1*alpha2*derphi0):
380 return alpha2, fc, 0
381
382 if (alpha1 - alpha2) > alpha1 / 2.0 or (1 - alpha2/alpha1) < 0.96:
383 alpha2 = alpha1 / 2.0
384
385 alpha0 = alpha1
386 alpha1 = alpha2
387 phi_a0 = phi_a1
388 phi_a1 = phi_a2
389
390 epsilon = 1e-8
391
393 f0 = apply(f,(xk,)+args)
394 grad = Num.zeros((len(xk),),'d')
395 ei = Num.zeros((len(xk),),'d')
396 for k in range(len(xk)):
397 ei[k] = 1.0
398 grad[k] = (apply(f,(xk+epsilon*ei,)+args) - f0)/epsilon
399 ei[k] = 0.0
400 return grad
401
403 f2 = apply(fprime,(x0+epsilon*p,)+args)
404 f1 = apply(fprime,(x0,)+args)
405 return (f2 - f1)/epsilon
406
407
408 -def fmin_bfgs(f, x0, fprime=None, args=(), avegtol=1e-5, maxiter=None,
409 full_output=0, disp=1):
410 """Minimize a function using the BFGS algorithm.
411
412 Description:
413
414 Optimize the function, f, whose gradient is given by fprime using the
415 quasi-Newton method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS)
416 See Wright, and Nocedal 'Numerical Optimization', 1999, pg. 198.
417
418 Inputs:
419
420 f -- the Python function or method to be minimized.
421 x0 -- the initial guess for the minimizer.
422 fprime -- a function to compute the gradient of f.
423 args -- extra arguments to f and fprime.
424
425 Outputs: (xopt, {fopt, func_calls, grad_calls, warnflag})
426
427 xopt -- the minimizer of f.
428
429 fopt -- the value of f(xopt).
430 func_calls -- the number of function_calls.
431 grad_calls -- the number of gradient calls.
432 warnflag -- an integer warning flag:
433 1 : 'Maximum number of iterations exceeded.'
434
435 Additional Inputs:
436
437 avegtol -- the minimum occurs when fprime(xopt)==0. This specifies how
438 close to zero the average magnitude of fprime(xopt) needs
439 to be.
440 maxiter -- the maximum number of iterations.
441 full_output -- if non-zero then return fopt, func_calls, grad_calls,
442 and warnflag in addition to xopt.
443 disp -- print convergence message if non-zero.
444 """
445 app_fprime = 0
446 if fprime is None:
447 app_fprime = 1
448
449 x0 = asarray(x0)
450 if maxiter is None:
451 maxiter = len(x0)*200
452 func_calls = 0
453 grad_calls = 0
454 k = 0
455 N = len(x0)
456 gtol = N*avegtol
457 I = MLab.eye(N)
458 Hk = I
459
460 if app_fprime:
461 gfk = apply(approx_fprime,(x0,f)+args)
462 func_calls = func_calls + len(x0) + 1
463 else:
464 gfk = apply(fprime,(x0,)+args)
465 grad_calls = grad_calls + 1
466 xk = x0
467 sk = [2*gtol]
468 warnflag = 0
469 while (Num.add.reduce(abs(gfk)) > gtol) and (k < maxiter):
470 pk = -Num.dot(Hk,gfk)
471 alpha_k, fc, gc = line_search_BFGS(f,xk,pk,gfk,args)
472 func_calls = func_calls + fc
473 xkp1 = xk + alpha_k * pk
474 sk = xkp1 - xk
475 xk = xkp1
476 if app_fprime:
477 gfkp1 = apply(approx_fprime,(xkp1,f)+args)
478 func_calls = func_calls + gc + len(x0) + 1
479 else:
480 gfkp1 = apply(fprime,(xkp1,)+args)
481 grad_calls = grad_calls + gc + 1
482
483 yk = gfkp1 - gfk
484 k = k + 1
485
486 try:
487 rhok = 1 / Num.dot(yk,sk)
488 except ZeroDivisionError:
489 warnflag = 2
490 break
491 A1 = I - sk[:,Num.NewAxis] * yk[Num.NewAxis,:] * rhok
492 A2 = I - yk[:,Num.NewAxis] * sk[Num.NewAxis,:] * rhok
493 Hk = Num.dot(A1,Num.dot(Hk,A2)) + rhok * sk[:,Num.NewAxis] \
494 * sk[Num.NewAxis,:]
495 gfk = gfkp1
496
497
498 if disp or full_output:
499 fval = apply(f,(xk,)+args)
500 if warnflag == 2:
501 if disp:
502 print "Warning: Desired error not necessarily achieved due to precision loss"
503 print " Current function value: %f" % fval
504 print " Iterations: %d" % k
505 print " Function evaluations: %d" % func_calls
506 print " Gradient evaluations: %d" % grad_calls
507
508 elif k >= maxiter:
509 warnflag = 1
510 if disp:
511 print "Warning: Maximum number of iterations has been exceeded"
512 print " Current function value: %f" % fval
513 print " Iterations: %d" % k
514 print " Function evaluations: %d" % func_calls
515 print " Gradient evaluations: %d" % grad_calls
516 else:
517 if disp:
518 print "Optimization terminated successfully."
519 print " Current function value: %f" % fval
520 print " Iterations: %d" % k
521 print " Function evaluations: %d" % func_calls
522 print " Gradient evaluations: %d" % grad_calls
523
524 if full_output:
525 return xk, fval, func_calls, grad_calls, warnflag
526 else:
527 return xk
528
529
530 -def fmin_ncg(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5,
531 maxiter=None, full_output=0, disp=1):
532 """Description:
533
534 Minimize the function, f, whose gradient is given by fprime using the
535 Newton-CG method. fhess_p must compute the hessian times an arbitrary
536 vector. If it is not given, finite-differences on fprime are used to
537 compute it. See Wright, and Nocedal 'Numerical Optimization', 1999,
538 pg. 140.
539
540 Inputs:
541
542 f -- the Python function or method to be minimized.
543 x0 -- the initial guess for the minimizer.
544 fprime -- a function to compute the gradient of f: fprime(x, *args)
545 fhess_p -- a function to compute the Hessian of f times an
546 arbitrary vector: fhess_p (x, p, *args)
547 fhess -- a function to compute the Hessian matrix of f.
548 args -- extra arguments for f, fprime, fhess_p, and fhess (the same
549 set of extra arguments is supplied to all of these functions).
550
551 Outputs: (xopt, {fopt, fcalls, gcalls, hcalls, warnflag})
552
553 xopt -- the minimizer of f
554
555 fopt -- the value of the function at xopt: fopt = f(xopt)
556 fcalls -- the number of function calls.
557 gcalls -- the number of gradient calls.
558 hcalls -- the number of hessian calls.
559 warnflag -- algorithm warnings:
560 1 : 'Maximum number of iterations exceeded.'
561
562 Additional Inputs:
563
564 avextol -- Convergence is assumed when the average relative error in
565 the minimizer falls below this amount.
566 maxiter -- Maximum number of iterations to allow.
567 full_output -- If non-zero return the optional outputs.
568 disp -- If non-zero print convergence message.
569
570 Remarks:
571
572 Only one of fhess_p or fhess need be given. If fhess is provided,
573 then fhess_p will be ignored. If neither fhess nor fhess_p is
574 provided, then the hessian product will be approximated using finite
575 differences on fprime.
576
577 """
578
579 x0 = asarray(x0)
580 fcalls = 0
581 gcalls = 0
582 hcalls = 0
583 if maxiter is None:
584 maxiter = len(x0)*200
585
586 xtol = len(x0)*avextol
587 update = [2*xtol]
588 xk = x0
589 k = 0
590 while (Num.add.reduce(abs(update)) > xtol) and (k < maxiter):
591
592
593 b = -apply(fprime,(xk,)+args)
594 gcalls = gcalls + 1
595 maggrad = Num.add.reduce(abs(b))
596 eta = min([0.5,Num.sqrt(maggrad)])
597 termcond = eta * maggrad
598 xsupi = 0
599 ri = -b
600 psupi = -ri
601 i = 0
602 dri0 = Num.dot(ri,ri)
603
604 if fhess is not None:
605 A = apply(fhess,(xk,)+args)
606 hcalls = hcalls + 1
607
608 while Num.add.reduce(abs(ri)) > termcond:
609 if fhess is None:
610 if fhess_p is None:
611 Ap = apply(approx_fhess_p,(xk,psupi,fprime)+args)
612 gcalls = gcalls + 2
613 else:
614 Ap = apply(fhess_p,(xk,psupi)+args)
615 hcalls = hcalls + 1
616 else:
617 Ap = Num.dot(A,psupi)
618
619 curv = Num.dot(psupi,Ap)
620 if (curv <= 0):
621 if (i > 0):
622 break
623 else:
624 xsupi = xsupi + dri0/curv * psupi
625 break
626 alphai = dri0 / curv
627 xsupi = xsupi + alphai * psupi
628 ri = ri + alphai * Ap
629 dri1 = Num.dot(ri,ri)
630 betai = dri1 / dri0
631 psupi = -ri + betai * psupi
632 i = i + 1
633 dri0 = dri1
634
635 pk = xsupi
636 gfk = -b
637 alphak, fc, gc = line_search_BFGS(f,xk,pk,gfk,args)
638 fcalls = fcalls + fc
639 gcalls = gcalls + gc
640
641 update = alphak * pk
642 xk = xk + update
643 k = k + 1
644
645 if disp or full_output:
646 fval = apply(f,(xk,)+args)
647 if k >= maxiter:
648 warnflag = 1
649 if disp:
650 print "Warning: Maximum number of iterations has been exceeded"
651 print " Current function value: %f" % fval
652 print " Iterations: %d" % k
653 print " Function evaluations: %d" % fcalls
654 print " Gradient evaluations: %d" % gcalls
655 print " Hessian evaluations: %d" % hcalls
656 else:
657 warnflag = 0
658 if disp:
659 print "Optimization terminated successfully."
660 print " Current function value: %f" % fval
661 print " Iterations: %d" % k
662 print " Function evaluations: %d" % fcalls
663 print " Gradient evaluations: %d" % gcalls
664 print " Hessian evaluations: %d" % hcalls
665
666 if full_output:
667 return xk, fval, fcalls, gcalls, hcalls, warnflag
668 else:
669 return xk
670
671 -def fminbound(func, x1, x2, args=(), xtol=1e-5, maxfun=500,
672 full_output=0, disp=1):
673 """Bounded minimization for scalar functions.
674
675 Description:
676
677 Finds a local minimizer of the scalar function func in the interval
678 x1 < xopt < x2 using Brent's method. (See brent for auto-bracketing).
679
680 Inputs:
681
682 func -- the function to be minimized (must accept scalar input and return
683 scalar output).
684 x1, x2 -- the optimization bounds.
685 args -- extra arguments to pass to function.
686 xtol -- the convergence tolerance.
687 maxfun -- maximum function evaluations.
688 full_output -- Non-zero to return optional outputs.
689 disp -- Non-zero to print messages.
690 0 : no message printing.
691 1 : non-convergence notification messages only.
692 2 : print a message on convergence too.
693 3 : print iteration results.
694
695
696 Outputs: (xopt, {fval, ierr, numfunc})
697
698 xopt -- The minimizer of the function over the interval.
699 fval -- The function value at the minimum point.
700 ierr -- An error flag (0 if converged, 1 if maximum number of
701 function calls reached).
702 numfunc -- The number of function calls.
703 """
704
705 if x1 > x2:
706 raise ValueError, "The lower bound exceeds the upper bound."
707
708 flag = 0
709 header = ' Func-count x f(x) Procedure'
710 step=' initial'
711
712 sqrt_eps = sqrt(2.2e-16)
713 golden_mean = 0.5*(3.0-sqrt(5.0))
714 a, b = x1, x2
715 fulc = a + golden_mean*(b-a)
716 nfc, xf = fulc, fulc
717 rat = e = 0.0
718 x = xf
719 fx = func(x,*args)
720 num = 1
721 fmin_data = (1, xf, fx)
722
723 ffulc = fnfc = fx
724 xm = 0.5*(a+b)
725 tol1 = sqrt_eps*abs(xf) + xtol / 3.0
726 tol2 = 2.0*tol1
727
728 if disp > 2:
729 print (" ")
730 print (header)
731 print "%5.0f %12.6g %12.6g %s" % (fmin_data + (step,))
732
733
734 while ( abs(xf-xm) > (tol2 - 0.5*(b-a)) ):
735 golden = 1
736
737 if abs(e) > tol1:
738 golden = 0
739 r = (xf-nfc)*(fx-ffulc)
740 q = (xf-fulc)*(fx-fnfc)
741 p = (xf-fulc)*q - (xf-nfc)*r
742 q = 2.0*(q-r)
743 if q > 0.0: p = -p
744 q = abs(q)
745 r = e
746 e = rat
747
748
749 if ( (abs(p) < abs(0.5*q*r)) and (p > q*(a-xf)) and (p < q*(b-xf))):
750 rat = (p+0.0) / q;
751 x = xf + rat
752 step = ' parabolic'
753
754 if ((x-a) < tol2) or ((b-x) < tol2):
755 si = Numeric.sign(xm-xf) + ((xm-xf)==0)
756 rat = tol1*si
757 else:
758 golden = 1
759
760 if golden:
761 if xf >= xm:
762 e=a-xf
763 else:
764 e=b-xf
765 rat = golden_mean*e
766 step = ' golden'
767
768 si = Numeric.sign(rat) + (rat == 0)
769 x = xf + si*max([abs(rat), tol1])
770 fu = func(x,*args)
771 num += 1
772 fmin_data = (num, x, fu)
773 if disp > 2:
774 print "%5.0f %12.6g %12.6g %s" % (fmin_data + (step,))
775
776 if fu <= fx:
777 if x >= xf:
778 a = xf
779 else:
780 b = xf
781 fulc, ffulc = nfc, fnfc
782 nfc, fnfc = xf, fx
783 xf, fx = x, fu
784 else:
785 if x < xf:
786 a = x
787 else:
788 b = x
789 if (fu <= fnfc) or (nfc == xf):
790 fulc, ffulc = nfc, fnfc
791 nfc, fnfc = x, fu
792 elif (fu <= ffulc) or (fulc == xf) or (fulc == nfc):
793 fulc, ffulc = x, fu
794
795 xm = 0.5*(a+b)
796 tol1 = sqrt_eps*abs(xf) + xtol/3.0
797 tol2 = 2.0*tol1
798
799 if num >= maxfun:
800 flag = 1
801 fval = fx
802 if disp > 0:
803 _endprint(x, flag, fval, maxfun, tol, disp)
804 if full_output:
805 return xf, fval, flag, num
806 else:
807 return xf
808
809 fval = fx
810 if disp > 0:
811 _endprint(x, flag, fval, maxfun, xtol, disp)
812
813 if full_output:
814 return xf, fval, flag, num
815 else:
816 return xf
817
818 -def brent(func, args=(), brack=None, tol=1.48e-8, full_output=0, maxiter=500):
819 """ Given a function of one-variable and a possible bracketing interval,
820 return the minimum of the function isolated to a fractional precision of
821 tol. A bracketing interval is a triple (a,b,c) where (a<b<c) and
822 func(b) < func(a),func(c). If bracket is two numbers then they are
823 assumed to be a starting interval for a downhill bracket search
824 (see bracket)
825
826 Uses inverse parabolic interpolation when possible to speed up convergence
827 of golden section method.
828
829 """
830 _mintol = 1.0e-11
831 _cg = 0.3819660
832 if brack is None:
833 xa,xb,xc,fa,fb,fc,funcalls = bracket(func, args=args)
834 elif len(brack) == 2:
835 xa,xb,xc,fa,fb,fc,funcalls = bracket(func, xa=brack[0], xb=brack[1], args=args)
836 elif len(brack) == 3:
837 xa,xb,xc = brack
838 if (xa > xc):
839 dum = xa; xa=xc; xc=dum
840 assert ((xa < xb) and (xb < xc)), "Not a bracketing interval."
841 fa = apply(func, (xa,)+args)
842 fb = apply(func, (xb,)+args)
843 fc = apply(func, (xc,)+args)
844 assert ((fb<fa) and (fb < fc)), "Not a bracketing interval."
845 funcalls = 3
846 else:
847 raise ValuError, "Bracketing interval must be length 2 or 3 sequence."
848
849 x=w=v=xb
850 fw=fv=fx=apply(func, (x,)+args)
851 if (xa < xc):
852 a = xa; b = xc
853 else:
854 a = xc; b = xa
855 deltax= 0.0
856 funcalls = 1
857 iter = 0
858 while (iter < maxiter):
859 tol1 = tol*abs(x) + _mintol
860 tol2 = 2.0*tol1
861 xmid = 0.5*(a+b)
862 if abs(x-xmid) < (tol2-0.5*(b-a)):
863 xmin=x; fval=fx
864 break
865 if (abs(deltax) <= tol1):
866 if (x>=xmid): deltax=a-x
867 else: deltax=b-x
868 rat = _cg*deltax
869 else:
870 tmp1 = (x-w)*(fx-fv)
871 tmp2 = (x-v)*(fx-fw)
872 p = (x-v)*tmp2 - (x-w)*tmp1;
873 tmp2 = 2.0*(tmp2-tmp1)
874 if (tmp2 > 0.0): p = -p
875 tmp2 = abs(tmp2)
876 dx_temp = deltax
877 deltax= rat
878
879 if ((p > tmp2*(a-x)) and (p < tmp2*(b-x)) and (abs(p) < abs(0.5*tmp2*dx_temp))):
880 rat = p*1.0/tmp2
881 u = x + rat
882 if ((u-a) < tol2 or (b-u) < tol2):
883 if xmid-x >= 0: rat = tol1
884 else: rat = -tol1
885 else:
886 if (x>=xmid): deltax=a-x
887 else: deltax=b-x
888 rat = _cg*deltax
889
890 if (abs(rat) < tol1):
891 if rat >= 0: u = x + tol1
892 else: u = x - tol1
893 else:
894 u = x + rat
895 fu = apply(func, (u,)+args)
896 funcalls += 1
897
898 if (fu > fx):
899 if (u<x): a=u
900 else: b=u
901 if (fu<=fw) or (w==x):
902 v=w; w=u; fv=fw; fw=fu
903 elif (fu<=fv) or (v==x) or (v==w):
904 v=u; fv=fu
905 else:
906 if (u >= x): a = x
907 else: b = x
908 v=w; w=x; x=u
909 fv=fw; fw=fx; fx=fu
910
911 xmin = x
912 fval = fx
913 if full_output:
914 return xmin, fval, iter, funcalls
915 else:
916 return xmin
917
918 -def golden(func, args=(), brack=None, tol=1.49e-8, full_output=0):
919 """ Given a function of one-variable and a possible bracketing interval,
920 return the minimum of the function isolated to a fractional precision of
921 tol. A bracketing interval is a triple (a,b,c) where (a<b<c) and
922 func(b) < func(a),func(c). If bracket is two numbers then they are
923 assumed to be a starting interval for a downhill bracket search
924 (see bracket)
925
926 Uses analog of bisection method to decrease the bracketed interval.
927 """
928 if brack is None:
929 xa,xb,xc,fa,fb,fc,funcalls = bracket(func, args=args)
930 elif len(brack) == 2:
931 xa,xb,xc,fa,fb,fc,funcalls = bracket(func, xa=brack[0], xb=brack[1], args=args)
932 elif len(brack) == 3:
933 xa,xb,xc = brack
934 if (xa > xc):
935 dum = xa; xa=xc; xc=dum
936 assert ((xa < xb) and (xb < xc)), "Not a bracketing interval."
937 fa = apply(func, (xa,)+args)
938 fb = apply(func, (xb,)+args)
939 fc = apply(func, (xc,)+args)
940 assert ((fb<fa) and (fb < fc)), "Not a bracketing interval."
941 funcalls = 3
942 else:
943 raise ValuError, "Bracketing interval must be length 2 or 3 sequence."
944
945 _gR = 0.61803399
946 _gC = 1.0-_gR
947 x3 = xc
948 x0 = xa
949 if (abs(xc-xb) > abs(xb-xa)):
950 x1 = xb
951 x2 = xb + _gC*(xc-xb)
952 else:
953 x2 = xb
954 x1 = xb - _gC*(xb-xa)
955 f1 = apply(func, (x1,)+args)
956 f2 = apply(func, (x2,)+args)
957 funcalls += 2
958 while (abs(x3-x0) > tol*(abs(x1)+abs(x2))):
959 if (f2 < f1):
960 x0 = x1; x1 = x2; x2 = _gR*x1 + _gC*x3
961 f1 = f2; f2 = apply(func, (x2,)+args)
962 else:
963 x3 = x2; x2 = x1; x1 = _gR*x2 + _gC*x0
964 f2 = f1; f1 = apply(func, (x1,)+args)
965 funcalls += 1
966 if (f1 < f2):
967 xmin = x1
968 fval = f1
969 else:
970 xmin = x2
971 fval = f2
972 if full_output:
973 return xmin, fval, funcalls
974 else:
975 return xmin
976
977
978 -def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0):
979 """Given a function and distinct initial points, search in the downhill
980 direction (as defined by the initital points) and return new points
981 xa, xb, xc that bracket the minimum of the function:
982 f(xa) > f(xb) < f(xc)
983 """
984 _gold = 1.618034
985 _verysmall_num = 1e-21
986 fa = apply(func, (xa,)+args)
987 fb = apply(func, (xb,)+args)
988 if (fa < fb):
989 dum = xa; xa = xb; xb = dum
990 dum = fa; fa = fb; fb = dum
991 xc = xb + _gold*(xb-xa)
992 fc = apply(func, (xc,)+args)
993 funcalls = 3
994 iter = 0
995 while (fc < fb):
996 tmp1 = (xb - xa)*(fb-fc)
997 tmp2 = (xb - xc)*(fb-fa)
998 val = tmp2-tmp1
999 if abs(val) < _verysmall_num:
1000 denom = 2.0*_verysmall_num
1001 else:
1002 denom = 2.0*val
1003 w = xb - ((xb-xc)*tmp2-(xb-xa)*tmp1)/denom
1004 wlim = xb + grow_limit*(xc-xb)
1005 if iter > 1000:
1006 raise RunTimeError, "Too many iterations."
1007 if (w-xc)*(xb-w) > 0.0:
1008 fw = apply(func, (w,)+args)
1009 funcalls += 1
1010 if (fw < fc):
1011 xa = xb; xb=w; fa=fb; fb=fw
1012 return xa, xb, xc, fa, fb, fc, funcalls
1013 elif (fw > fb):
1014 xc = w; fc=fw
1015 return xa, xb, xc, fa, fb, fc, funcalls
1016 w = xc + _gold*(xc-xb)
1017 fw = apply(func, (w,)+args)
1018 funcalls += 1
1019 elif (w-wlim)*(wlim-xc) >= 0.0:
1020 w = wlim
1021 fw = apply(func, (w,)+args)
1022 funcalls += 1
1023 elif (w-wlim)*(xc-w) > 0.0:
1024 fw = apply(func, (w,)+args)
1025 funcalls += 1
1026 if (fw < fc):
1027 xb=xc; xc=w; w=xc+_gold*(xc-xb)
1028 fb=fc; fc=fw; fw=apply(func, (w,)+args)
1029 funcalls += 1
1030 else:
1031 w = xc + _gold*(xc-xb)
1032 fw = apply(func, (w,)+args)
1033 funcalls += 1
1034 xa=xb; xb=xc; xc=w
1035 fa=fb; fb=fc; fc=fw
1036 return xa, xb, xc, fa, fb, fc, funcalls
1037
1038
1039 global _powell_funcalls
1040
1041 -def _myfunc(alpha, func, x0, direc, args=()):
1042 funcargs = (x0 + alpha * direc,)+args
1043 return func(*funcargs)
1044
1046
1047
1048
1049 global _powell_funcalls
1050 extra_args = (func, p, xi) + args
1051 alpha_min, fret, iter, num = brent(_myfunc, args=extra_args,
1052 full_output=1, tol=tol)
1053 xi = alpha_min*xi
1054 _powell_funcalls += num
1055 return squeeze(fret), p+xi, xi
1056
1057
1058 -def fmin_powell(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None,
1059 full_output=0, disp=1):
1060 """Minimize a function using modified Powell's method.
1061
1062 Description:
1063
1064 Uses a modification of Powell's method to find the minimum of a function
1065 of N
1066
1067 Inputs:
1068
1069 func -- the Python function or method to be minimized.
1070 x0 -- the initial guess.
1071 args -- extra arguments for func.
1072
1073 Outputs: (xopt, {fopt, xi, iter, funcalls, warnflag})
1074
1075 xopt -- minimizer of function
1076
1077 fopt -- value of function at minimum: fopt = func(xopt)
1078 direc -- current direction set
1079 iter -- number of iterations
1080 funcalls -- number of function calls
1081 warnflag -- Integer warning flag:
1082 1 : 'Maximum number of function evaluations.'
1083 2 : 'Maximum number of iterations.'
1084
1085 Additional Inputs:
1086
1087 xtol -- line-search error tolerance.
1088 ftol -- acceptable relative error in func(xopt) for convergence.
1089 maxiter -- the maximum number of iterations to perform.
1090 maxfun -- the maximum number of function evaluations.
1091 full_output -- non-zero if fval and warnflag outputs are desired.
1092 disp -- non-zero to print convergence messages.
1093
1094 """
1095 global _powell_funcalls
1096 x = asarray(x0)
1097 N = len(x)
1098 rank = len(x.shape)
1099 if not -1 < rank < 2:
1100 raise ValueError, "Initial guess must be a scalar or rank-1 sequence."
1101 if maxiter is None:
1102 maxiter = N * 1000
1103 if maxfun is None:
1104 maxfun = N * 1000
1105
1106 direc = eye(N,typecode='d')
1107 fval = squeeze(apply(func, (x,)+args))
1108 _powell_funcalls = 1
1109 x1 = x.copy()
1110 iter = 0;
1111 ilist = range(N)
1112 while 1:
1113 fx = fval
1114 bigind = 0
1115 delta = 0.0
1116 for i in ilist:
1117 direc1 = direc[i]
1118 fx2 = fval
1119 fval, x, direc1 = _linesearch_powell(func, x, direc1, args=args, tol=xtol)
1120 if (fx2 - fval) > delta:
1121 delta = fx2 - fval
1122 bigind = i
1123 iter += 1
1124 if (2.0*(fx - fval) <= ftol*(abs(fx)+abs(fval))+1e-20): break
1125 if _powell_funcalls >= maxfun: break
1126 if iter >= maxiter: break
1127
1128
1129 direc1 = x - x1
1130 x2 = 2*x - x1
1131 x1 = x.copy()
1132 fx2 = squeeze(apply(func, (x2,)+args))
1133 _powell_funcalls +=1
1134
1135 if (fx > fx2):
1136 t = 2.0*(fx+fx2-2.0*fval)
1137 temp = (fx-fval-delta)
1138 t *= temp*temp
1139 temp = fx-fx2
1140 t -= delta*temp*temp
1141 if t < 0.0:
1142 fval, x, direc1 = _linesearch_powell(func, x, direc1, args=args, tol=xtol)
1143 direc[bigind] = direc[-1]
1144 direc[-1] = direc1
1145
1146 warnflag = 0
1147 if _powell_funcalls >= maxfun:
1148 warnflag = 1
1149 if disp:
1150 print "Warning: Maximum number of function evaluations has "\
1151 "been exceeded."
1152 elif iter >= maxiter:
1153 warnflag = 2
1154 if disp:
1155 print "Warning: Maximum number of iterations has been exceeded"
1156 else:
1157 if disp:
1158 print "Optimization terminated successfully."
1159 print " Current function value: %f" % fval
1160 print " Iterations: %d" % iter
1161 print " Function evaluations: %d" % _powell_funcalls
1162
1163 x = squeeze(x)
1164 if full_output:
1165 return x, fval, direc, iter, _powell_funcalls, warnflag
1166 else:
1167 return x
1168
1169
1170
1171 -def _endprint(x, flag, fval, maxfun, xtol, disp):
1172 if flag == 0:
1173 if disp > 1:
1174 print "\nOptimization terminated successfully;\n the returned value" + \
1175 " satisfies the termination criteria\n (using xtol = ", xtol, ")"
1176 if flag == 1:
1177 print "\nMaximum number of function evaluations exceeded --- increase maxfun argument.\n"
1178 return
1179
1180
1181 if __name__ == "__main__":
1182 import string
1183 import time
1184
1185
1186 times = []
1187 algor = []
1188 x0 = [0.8,1.2,0.7]
1189 start = time.time()
1190 x = fmin(rosen,x0)
1191 print x
1192 times.append(time.time() - start)
1193 algor.append('Nelder-Mead Simplex\t')
1194
1195 start = time.time()
1196 x = fmin_powell(rosen,x0)
1197 print x
1198 times.append(time.time() - start)
1199 algor.append('Powell Direction Set Method.')
1200
1201 start = time.time()
1202 x = fmin_bfgs(rosen, x0, fprime=rosen_der, maxiter=80)
1203 print x
1204 times.append(time.time() - start)
1205 algor.append('BFGS Quasi-Newton\t')
1206
1207 start = time.time()
1208 x = fmin_bfgs(rosen, x0, avegtol=1e-4, maxiter=100)
1209 print x
1210 times.append(time.time() - start)
1211 algor.append('BFGS without gradient\t')
1212
1213
1214 start = time.time()
1215 x = fmin_ncg(rosen, x0, rosen_der, fhess_p=rosen_hess_prod, maxiter=80)
1216 print x
1217 times.append(time.time() - start)
1218 algor.append('Newton-CG with hessian product')
1219
1220
1221 start = time.time()
1222 x = fmin_ncg(rosen, x0, rosen_der, fhess=rosen_hess, maxiter=80)
1223 print x
1224 times.append(time.time() - start)
1225 algor.append('Newton-CG with full hessian')
1226
1227 print "\nMinimizing the Rosenbrock function of order 3\n"
1228 print " Algorithm \t\t\t Seconds"
1229 print "===========\t\t\t ========="
1230 for k in range(len(algor)):
1231 print algor[k], "\t -- ", times[k]
1232