#__docformat__ = "restructuredtext en" # ******NOTICE*************** # optimize.py module by Travis E. Oliphant # # You may copy and use this module as you see fit with no # guarantee implied provided you keep this notice in all copies. # *****END NOTICE************ # A collection of optimization algorithms. Version 0.5 # CHANGES # Added fminbound (July 2001) # Added brute (Aug. 2002) # Finished line search satisfying strong Wolfe conditions (Mar. 2004) # Updated strong Wolfe conditions line search to use cubic-interpolation (Mar. 2004) # Minimization routines __all__ = ['fmin', 'fmin_powell','fmin_bfgs', 'fmin_ncg', 'fmin_cg', 'fminbound','brent', 'golden','bracket','rosen','rosen_der', 'rosen_hess', 'rosen_hess_prod', 'brute', 'approx_fprime', 'line_search', 'check_grad'] __docformat__ = "restructuredtext en" import numpy from numpy import atleast_1d, eye, mgrid, argmin, zeros, shape, \ squeeze, vectorize, asarray, absolute, sqrt, Inf, asfarray, isinf from linesearch import \ line_search_BFGS, line_search_wolfe1, line_search_wolfe2, \ line_search_wolfe2 as line_search # These have been copied from Numeric's MLab.py # I don't think they made the transition to scipy_core def max(m,axis=0): """max(m,axis=0) returns the maximum of m along dimension axis. """ m = asarray(m) return numpy.maximum.reduce(m,axis) def min(m,axis=0): """min(m,axis=0) returns the minimum of m along dimension axis. """ m = asarray(m) return numpy.minimum.reduce(m,axis) def is_array_scalar(x): """Test whether `x` is either a scalar or an array scalar. """ return len(atleast_1d(x) == 1) abs = absolute import __builtin__ pymin = __builtin__.min pymax = __builtin__.max __version__="0.7" _epsilon = sqrt(numpy.finfo(float).eps) def vecnorm(x, ord=2): if ord == Inf: return numpy.amax(abs(x)) elif ord == -Inf: return numpy.amin(abs(x)) else: return numpy.sum(abs(x)**ord,axis=0)**(1.0/ord) def rosen(x): """The Rosenbrock function. The function computed is sum(100.0*(x[1:] - x[:-1]**2.0)**2.0 + (1 - x[:-1])**2.0 Parameters ---------- x : array_like, 1D The point at which the Rosenbrock function is to be computed. Returns ------- f : float The value of the Rosenbrock function See Also -------- rosen_der, rosen_hess, rosen_hess_prod """ x = asarray(x) return numpy.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0,axis=0) def rosen_der(x): """The derivative (i.e. gradient) of the Rosenbrock function. Parameters ---------- x : array_like, 1D The point at which the derivative is to be computed. Returns ------- der : 1D numpy array The gradient of the Rosenbrock function at `x`. See Also -------- rosen, rosen_hess, rosen_hess_prod """ x = asarray(x) xm = x[1:-1] xm_m1 = x[:-2] xm_p1 = x[2:] der = numpy.zeros_like(x) der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm) der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0]) der[-1] = 200*(x[-1]-x[-2]**2) return der def rosen_hess(x): """The Hessian matrix of the Rosenbrock function. Parameters ---------- x : array_like, 1D The point at which the Hessian matrix is to be computed. Returns ------- hess : 2D numpy array The Hessian matrix of the Rosenbrock function at `x`. See Also -------- rosen, rosen_der, rosen_hess_prod """ x = atleast_1d(x) H = numpy.diag(-400*x[:-1],1) - numpy.diag(400*x[:-1],-1) diagonal = numpy.zeros(len(x), dtype=x.dtype) diagonal[0] = 1200*x[0]**2 - 400*x[1] + 2 diagonal[-1] = 200 diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:] H = H + numpy.diag(diagonal) return H def rosen_hess_prod(x,p): """Product of the Hessian matrix of the Rosenbrock function with a vector. Parameters ---------- x : array_like, 1D The point at which the Hessian matrix is to be computed. p : array_like, 1D, same size as `x`. The vector to be multiplied by the Hessian matrix. Returns ------- v : 1D numpy array The Hessian matrix of the Rosenbrock function at `x` multiplied by the vector `p`. See Also -------- rosen, rosen_der, rosen_hess """ x = atleast_1d(x) Hp = numpy.zeros(len(x), dtype=x.dtype) Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1] Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] \ -400*x[1:-1]*p[2:] Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1] return Hp def wrap_function(function, args): ncalls = [0] def function_wrapper(x): ncalls[0] += 1 return function(x, *args) return ncalls, function_wrapper def fmin(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None): """ Minimize a function using the downhill simplex algorithm. This algorithm only uses function values, not derivatives or second derivatives. Parameters ---------- func : callable func(x,*args) The objective function to be minimized. x0 : ndarray Initial guess. args : tuple Extra arguments passed to func, i.e. ``f(x,*args)``. callback : callable Called after each iteration, as callback(xk), where xk is the current parameter vector. Returns ------- xopt : ndarray Parameter that minimizes function. fopt : float Value of function at minimum: ``fopt = func(xopt)``. iter : int Number of iterations performed. funcalls : int Number of function calls made. warnflag : int 1 : Maximum number of function evaluations made. 2 : Maximum number of iterations reached. allvecs : list Solution at each iteration. Other parameters ---------------- xtol : float Relative error in xopt acceptable for convergence. ftol : number Relative error in func(xopt) acceptable for convergence. maxiter : int Maximum number of iterations to perform. maxfun : number Maximum number of function evaluations to make. full_output : bool Set to True if fopt and warnflag outputs are desired. disp : bool Set to True to print convergence messages. retall : bool Set to True to return list of solutions at each iteration. Notes ----- Uses a Nelder-Mead simplex algorithm to find the minimum of function of one or more variables. This algorithm has a long history of successful use in applications. But it will usually be slower than an algorithm that uses first or second derivative information. In practice it can have poor performance in high-dimensional problems and is not robust to minimizing complicated functions. Additionally, there currently is no complete theory describing when the algorithm will successfully converge to the minimum, or how fast it will if it does. References ---------- Nelder, J.A. and Mead, R. (1965), "A simplex method for function minimization", The Computer Journal, 7, pp. 308-313 Wright, M.H. (1996), "Direct Search Methods: Once Scorned, Now Respectable", in Numerical Analysis 1995, Proceedings of the 1995 Dundee Biennial Conference in Numerical Analysis, D.F. Griffiths and G.A. Watson (Eds.), Addison Wesley Longman, Harlow, UK, pp. 191-208. """ fcalls, func = wrap_function(func, args) x0 = asfarray(x0).flatten() N = len(x0) rank = len(x0.shape) if not -1 < rank < 2: raise ValueError("Initial guess must be a scalar or rank-1 sequence.") if maxiter is None: maxiter = N * 200 if maxfun is None: maxfun = N * 200 rho = 1; chi = 2; psi = 0.5; sigma = 0.5; one2np1 = range(1,N+1) if rank == 0: sim = numpy.zeros((N+1,), dtype=x0.dtype) else: sim = numpy.zeros((N+1,N), dtype=x0.dtype) fsim = numpy.zeros((N+1,), float) sim[0] = x0 if retall: allvecs = [sim[0]] fsim[0] = func(x0) nonzdelt = 0.05 zdelt = 0.00025 for k in range(0,N): y = numpy.array(x0,copy=True) if y[k] != 0: y[k] = (1+nonzdelt)*y[k] else: y[k] = zdelt sim[k+1] = y f = func(y) fsim[k+1] = f ind = numpy.argsort(fsim) fsim = numpy.take(fsim,ind,0) # sort so sim[0,:] has the lowest function value sim = numpy.take(sim,ind,0) iterations = 1 while (fcalls[0] < maxfun and iterations < maxiter): if (max(numpy.ravel(abs(sim[1:]-sim[0]))) <= xtol \ and max(abs(fsim[0]-fsim[1:])) <= ftol): break xbar = numpy.add.reduce(sim[:-1],0) / N xr = (1+rho)*xbar - rho*sim[-1] fxr = func(xr) doshrink = 0 if fxr < fsim[0]: xe = (1+rho*chi)*xbar - rho*chi*sim[-1] fxe = func(xe) if fxe < fxr: sim[-1] = xe fsim[-1] = fxe else: sim[-1] = xr fsim[-1] = fxr else: # fsim[0] <= fxr if fxr < fsim[-2]: sim[-1] = xr fsim[-1] = fxr else: # fxr >= fsim[-2] # Perform contraction if fxr < fsim[-1]: xc = (1+psi*rho)*xbar - psi*rho*sim[-1] fxc = func(xc) if fxc <= fxr: sim[-1] = xc fsim[-1] = fxc else: doshrink=1 else: # Perform an inside contraction xcc = (1-psi)*xbar + psi*sim[-1] fxcc = func(xcc) if fxcc < fsim[-1]: sim[-1] = xcc fsim[-1] = fxcc else: doshrink = 1 if doshrink: for j in one2np1: sim[j] = sim[0] + sigma*(sim[j] - sim[0]) fsim[j] = func(sim[j]) ind = numpy.argsort(fsim) sim = numpy.take(sim,ind,0) fsim = numpy.take(fsim,ind,0) if callback is not None: callback(sim[0]) iterations += 1 if retall: allvecs.append(sim[0]) x = sim[0] fval = min(fsim) warnflag = 0 if fcalls[0] >= maxfun: warnflag = 1 if disp: print "Warning: Maximum number of function evaluations has "\ "been exceeded." elif iterations >= maxiter: warnflag = 2 if disp: print "Warning: Maximum number of iterations has been exceeded" else: if disp: print "Optimization terminated successfully." print " Current function value: %f" % fval print " Iterations: %d" % iterations print " Function evaluations: %d" % fcalls[0] if full_output: retlist = x, fval, iterations, fcalls[0], warnflag if retall: retlist += (allvecs,) else: retlist = x if retall: retlist = (x, allvecs) return retlist def approx_fprime(xk,f,epsilon,*args): f0 = f(*((xk,)+args)) grad = numpy.zeros((len(xk),), float) ei = numpy.zeros((len(xk),), float) for k in range(len(xk)): ei[k] = epsilon grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon ei[k] = 0.0 return grad def check_grad(func, grad, x0, *args): """Check the correctness of a gradient function by comparing it against a finite-difference approximation of the gradient. Parameters ---------- func: callable func(x0,*args) Function whose derivative is to be checked grad: callable grad(x0, *args) Gradient of func x0: ndarray Points to check grad against finite difference approximation of grad using func. args: optional Extra arguments passed to func and grad Returns ------- err: float The square root of the sum of squares (i.e. the 2-norm) of the difference between grad(x0, *args) and the finite difference approximation of grad using func at the points x0. """ return sqrt(sum((grad(x0,*args)-approx_fprime(x0,func,_epsilon,*args))**2)) def approx_fhess_p(x0,p,fprime,epsilon,*args): f2 = fprime(*((x0+epsilon*p,)+args)) f1 = fprime(*((x0,)+args)) return (f2 - f1)/epsilon def fmin_bfgs(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf, epsilon=_epsilon, maxiter=None, full_output=0, disp=1, retall=0, callback=None): """Minimize a function using the BFGS algorithm. Parameters ---------- f : callable f(x,*args) Objective function to be minimized. x0 : ndarray Initial guess. fprime : callable f'(x,*args) Gradient of f. args : tuple Extra arguments passed to f and fprime. gtol : float Gradient norm must be less than gtol before succesful termination. norm : float Order of norm (Inf is max, -Inf is min) epsilon : int or ndarray If fprime is approximated, use this value for the step size. callback : callable An optional user-supplied function to call after each iteration. Called as callback(xk), where xk is the current parameter vector. Returns ------- xopt : ndarray Parameters which minimize f, i.e. f(xopt) == fopt. fopt : float Minimum value. gopt : ndarray Value of gradient at minimum, f'(xopt), which should be near 0. Bopt : ndarray Value of 1/f''(xopt), i.e. the inverse hessian matrix. func_calls : int Number of function_calls made. grad_calls : int Number of gradient calls made. warnflag : integer 1 : Maximum number of iterations exceeded. 2 : Gradient and/or function calls not changing. allvecs : list Results at each iteration. Only returned if retall is True. Other Parameters ---------------- maxiter : int Maximum number of iterations to perform. full_output : bool If True,return fopt, func_calls, grad_calls, and warnflag in addition to xopt. disp : bool Print convergence message if True. retall : bool Return a list of results at each iteration if True. Notes ----- Optimize the function, f, whose gradient is given by fprime using the quasi-Newton method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS) References ---------- Wright, and Nocedal 'Numerical Optimization', 1999, pg. 198. """ x0 = asarray(x0).flatten() if x0.ndim == 0: x0.shape = (1,) if maxiter is None: maxiter = len(x0)*200 func_calls, f = wrap_function(f, args) if fprime is None: grad_calls, myfprime = wrap_function(approx_fprime, (f, epsilon)) else: grad_calls, myfprime = wrap_function(fprime, args) gfk = myfprime(x0) k = 0 N = len(x0) I = numpy.eye(N,dtype=int) Hk = I old_fval = f(x0) old_old_fval = old_fval + 5000 xk = x0 if retall: allvecs = [x0] sk = [2*gtol] warnflag = 0 gnorm = vecnorm(gfk,ord=norm) while (gnorm > gtol) and (k < maxiter): pk = -numpy.dot(Hk,gfk) alpha_k, fc, gc, old_fval2, old_old_fval2, gfkp1 = \ line_search_wolfe1(f,myfprime,xk,pk,gfk, old_fval,old_old_fval) if alpha_k is not None: old_fval = old_fval2 old_old_fval = old_old_fval2 else: # line search failed: try different one. alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \ line_search_wolfe2(f,myfprime,xk,pk,gfk, old_fval,old_old_fval) if alpha_k is None: # This line search also failed to find a better solution. warnflag = 2 break xkp1 = xk + alpha_k * pk if retall: allvecs.append(xkp1) sk = xkp1 - xk xk = xkp1 if gfkp1 is None: gfkp1 = myfprime(xkp1) yk = gfkp1 - gfk gfk = gfkp1 if callback is not None: callback(xk) k += 1 gnorm = vecnorm(gfk,ord=norm) if (gnorm <= gtol): break if not numpy.isfinite(old_fval): # We correctly found +-Inf as optimal value, or something went # wrong. warnflag = 2 break try: # this was handled in numeric, let it remaines for more safety rhok = 1.0 / (numpy.dot(yk,sk)) except ZeroDivisionError: rhok = 1000.0 print "Divide-by-zero encountered: rhok assumed large" if isinf(rhok): # this is patch for numpy rhok = 1000.0 print "Divide-by-zero encountered: rhok assumed large" A1 = I - sk[:,numpy.newaxis] * yk[numpy.newaxis,:] * rhok A2 = I - yk[:,numpy.newaxis] * sk[numpy.newaxis,:] * rhok Hk = numpy.dot(A1,numpy.dot(Hk,A2)) + rhok * sk[:,numpy.newaxis] \ * sk[numpy.newaxis,:] if disp or full_output: fval = old_fval if warnflag == 2: if disp: print "Warning: Desired error not necessarily achieved" \ "due to precision loss" print " Current function value: %f" % fval print " Iterations: %d" % k print " Function evaluations: %d" % func_calls[0] print " Gradient evaluations: %d" % grad_calls[0] elif k >= maxiter: warnflag = 1 if disp: print "Warning: Maximum number of iterations has been exceeded" print " Current function value: %f" % fval print " Iterations: %d" % k print " Function evaluations: %d" % func_calls[0] print " Gradient evaluations: %d" % grad_calls[0] else: if disp: print "Optimization terminated successfully." print " Current function value: %f" % fval print " Iterations: %d" % k print " Function evaluations: %d" % func_calls[0] print " Gradient evaluations: %d" % grad_calls[0] if full_output: retlist = xk, fval, gfk, Hk, func_calls[0], grad_calls[0], warnflag if retall: retlist += (allvecs,) else: retlist = xk if retall: retlist = (xk, allvecs) return retlist def fmin_cg(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf, epsilon=_epsilon, maxiter=None, full_output=0, disp=1, retall=0, callback=None): """Minimize a function using a nonlinear conjugate gradient algorithm. Parameters ---------- f : callable f(x,*args) Objective function to be minimized. x0 : ndarray Initial guess. fprime : callable f'(x,*args) Function which computes the gradient of f. args : tuple Extra arguments passed to f and fprime. gtol : float Stop when norm of gradient is less than gtol. norm : float Order of vector norm to use. -Inf is min, Inf is max. epsilon : float or ndarray If fprime is approximated, use this value for the step size (can be scalar or vector). callback : callable An optional user-supplied function, called after each iteration. Called as callback(xk), where xk is the current parameter vector. Returns ------- xopt : ndarray Parameters which minimize f, i.e. f(xopt) == fopt. fopt : float Minimum value found, f(xopt). func_calls : int The number of function_calls made. grad_calls : int The number of gradient calls made. warnflag : int 1 : Maximum number of iterations exceeded. 2 : Gradient and/or function calls not changing. allvecs : ndarray If retall is True (see other parameters below), then this vector containing the result at each iteration is returned. Other Parameters ---------------- maxiter : int Maximum number of iterations to perform. full_output : bool If True then return fopt, func_calls, grad_calls, and warnflag in addition to xopt. disp : bool Print convergence message if True. retall : bool Return a list of results at each iteration if True. Notes ----- Optimize the function, f, whose gradient is given by fprime using the nonlinear conjugate gradient algorithm of Polak and Ribiere. See Wright & Nocedal, 'Numerical Optimization', 1999, pg. 120-122. """ x0 = asarray(x0).flatten() if maxiter is None: maxiter = len(x0)*200 func_calls, f = wrap_function(f, args) if fprime is None: grad_calls, myfprime = wrap_function(approx_fprime, (f, epsilon)) else: grad_calls, myfprime = wrap_function(fprime, args) gfk = myfprime(x0) k = 0 N = len(x0) xk = x0 old_fval = f(xk) old_old_fval = old_fval + 5000 if retall: allvecs = [xk] sk = [2*gtol] warnflag = 0 pk = -gfk gnorm = vecnorm(gfk,ord=norm) while (gnorm > gtol) and (k < maxiter): deltak = numpy.dot(gfk,gfk) # These values are modified by the line search, even if it fails old_fval_backup = old_fval old_old_fval_backup = old_old_fval alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \ line_search_wolfe1(f,myfprime,xk,pk,gfk,old_fval, old_old_fval,c2=0.4) if alpha_k is None: # line search failed -- use different one. alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \ line_search_wolfe2(f,myfprime,xk,pk,gfk, old_fval_backup,old_old_fval_backup) if alpha_k is None or alpha_k == 0: # This line search also failed to find a better solution. warnflag = 2 break xk = xk + alpha_k*pk if retall: allvecs.append(xk) if gfkp1 is None: gfkp1 = myfprime(xk) yk = gfkp1 - gfk beta_k = pymax(0,numpy.dot(yk,gfkp1)/deltak) pk = -gfkp1 + beta_k * pk gfk = gfkp1 gnorm = vecnorm(gfk,ord=norm) if callback is not None: callback(xk) k += 1 if disp or full_output: fval = old_fval if warnflag == 2: if disp: print "Warning: Desired error not necessarily achieved due to precision loss" print " Current function value: %f" % fval print " Iterations: %d" % k print " Function evaluations: %d" % func_calls[0] print " Gradient evaluations: %d" % grad_calls[0] elif k >= maxiter: warnflag = 1 if disp: print "Warning: Maximum number of iterations has been exceeded" print " Current function value: %f" % fval print " Iterations: %d" % k print " Function evaluations: %d" % func_calls[0] print " Gradient evaluations: %d" % grad_calls[0] else: if disp: print "Optimization terminated successfully." print " Current function value: %f" % fval print " Iterations: %d" % k print " Function evaluations: %d" % func_calls[0] print " Gradient evaluations: %d" % grad_calls[0] if full_output: retlist = xk, fval, func_calls[0], grad_calls[0], warnflag if retall: retlist += (allvecs,) else: retlist = xk if retall: retlist = (xk, allvecs) return retlist def fmin_ncg(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5, epsilon=_epsilon, maxiter=None, full_output=0, disp=1, retall=0, callback=None): """Unconstrained minimization of a function using the Newton-CG method. Parameters ---------- f : callable f(x,*args) Objective function to be minimized. x0 : ndarray Initial guess. fprime : callable f'(x,*args) Gradient of f. fhess_p : callable fhess_p(x,p,*args) Function which computes the Hessian of f times an arbitrary vector, p. fhess : callable fhess(x,*args) Function to compute the Hessian matrix of f. args : tuple Extra arguments passed to f, fprime, fhess_p, and fhess (the same set of extra arguments is supplied to all of these functions). epsilon : float or ndarray If fhess is approximated, use this value for the step size. callback : callable An optional user-supplied function which is called after each iteration. Called as callback(xk), where xk is the current parameter vector. Returns ------- xopt : ndarray Parameters which minimizer f, i.e. ``f(xopt) == fopt``. fopt : float Value of the function at xopt, i.e. ``fopt = f(xopt)``. fcalls : int Number of function calls made. gcalls : int Number of gradient calls made. hcalls : int Number of hessian calls made. warnflag : int Warnings generated by the algorithm. 1 : Maximum number of iterations exceeded. allvecs : list The result at each iteration, if retall is True (see below). Other Parameters ---------------- avextol : float Convergence is assumed when the average relative error in the minimizer falls below this amount. maxiter : int Maximum number of iterations to perform. full_output : bool If True, return the optional outputs. disp : bool If True, print convergence message. retall : bool If True, return a list of results at each iteration. Notes ----- Only one of `fhess_p` or `fhess` need to be given. If `fhess` is provided, then `fhess_p` will be ignored. If neither `fhess` nor `fhess_p` is provided, then the hessian product will be approximated using finite differences on `fprime`. `fhess_p` must compute the hessian times an arbitrary vector. If it is not given, finite-differences on `fprime` are used to compute it. Newton-CG methods are also called truncated Newton methods. This function differs from scipy.optimize.fmin_tnc because 1. scipy.optimize.fmin_ncg is written purely in python using numpy and scipy while scipy.optimize.fmin_tnc calls a C function. 2. scipy.optimize.fmin_ncg is only for unconstrained minimization while scipy.optimize.fmin_tnc is for unconstrained minimization or box constrained minimization. (Box constraints give lower and upper bounds for each variable seperately.) References ---------- Wright & Nocedal, 'Numerical Optimization', 1999, pg. 140. """ x0 = asarray(x0).flatten() fcalls, f = wrap_function(f, args) gcalls, fprime = wrap_function(fprime, args) hcalls = 0 if maxiter is None: maxiter = len(x0)*200 xtol = len(x0)*avextol update = [2*xtol] xk = x0 if retall: allvecs = [xk] k = 0 old_fval = f(x0) while (numpy.add.reduce(abs(update)) > xtol) and (k < maxiter): # Compute a search direction pk by applying the CG method to # del2 f(xk) p = - grad f(xk) starting from 0. b = -fprime(xk) maggrad = numpy.add.reduce(abs(b)) eta = min([0.5,numpy.sqrt(maggrad)]) termcond = eta * maggrad xsupi = zeros(len(x0), dtype=x0.dtype) ri = -b psupi = -ri i = 0 dri0 = numpy.dot(ri,ri) if fhess is not None: # you want to compute hessian once. A = fhess(*(xk,)+args) hcalls = hcalls + 1 while numpy.add.reduce(abs(ri)) > termcond: if fhess is None: if fhess_p is None: Ap = approx_fhess_p(xk,psupi,fprime,epsilon) else: Ap = fhess_p(xk,psupi, *args) hcalls = hcalls + 1 else: Ap = numpy.dot(A,psupi) # check curvature Ap = asarray(Ap).squeeze() # get rid of matrices... curv = numpy.dot(psupi,Ap) if 0 <= curv <= 3*numpy.finfo(numpy.float64).eps: break elif curv < 0: if (i > 0): break else: xsupi = xsupi + dri0/curv * psupi break alphai = dri0 / curv xsupi = xsupi + alphai * psupi ri = ri + alphai * Ap dri1 = numpy.dot(ri,ri) betai = dri1 / dri0 psupi = -ri + betai * psupi i = i + 1 dri0 = dri1 # update numpy.dot(ri,ri) for next time. pk = xsupi # search direction is solution to system. gfk = -b # gradient at xk alphak, fc, gc, old_fval = line_search_BFGS(f,xk,pk,gfk,old_fval) update = alphak * pk xk = xk + update # upcast if necessary if callback is not None: callback(xk) if retall: allvecs.append(xk) k += 1 if disp or full_output: fval = old_fval if k >= maxiter: warnflag = 1 if disp: print "Warning: Maximum number of iterations has been exceeded" print " Current function value: %f" % fval print " Iterations: %d" % k print " Function evaluations: %d" % fcalls[0] print " Gradient evaluations: %d" % gcalls[0] print " Hessian evaluations: %d" % hcalls else: warnflag = 0 if disp: print "Optimization terminated successfully." print " Current function value: %f" % fval print " Iterations: %d" % k print " Function evaluations: %d" % fcalls[0] print " Gradient evaluations: %d" % gcalls[0] print " Hessian evaluations: %d" % hcalls if full_output: retlist = xk, fval, fcalls[0], gcalls[0], hcalls, warnflag if retall: retlist += (allvecs,) else: retlist = xk if retall: retlist = (xk, allvecs) return retlist def fminbound(func, x1, x2, args=(), xtol=1e-5, maxfun=500, full_output=0, disp=1): """Bounded minimization for scalar functions. Parameters ---------- func : callable f(x,*args) Objective function to be minimized (must accept and return scalars). x1, x2 : float or array scalar The optimization bounds. args : tuple Extra arguments passed to function. xtol : float The convergence tolerance. maxfun : int Maximum number of function evaluations allowed. full_output : bool If True, return optional outputs. disp : int If non-zero, print messages. 0 : no message printing. 1 : non-convergence notification messages only. 2 : print a message on convergence too. 3 : print iteration results. Returns ------- xopt : ndarray Parameters (over given interval) which minimize the objective function. fval : number The function value at the minimum point. ierr : int An error flag (0 if converged, 1 if maximum number of function calls reached). numfunc : int The number of function calls made. Notes ----- Finds a local minimizer of the scalar function `func` in the interval x1 < xopt < x2 using Brent's method. (See `brent` for auto-bracketing). """ # Test bounds are of correct form if not (is_array_scalar(x1) and is_array_scalar(x2)): raise ValueError("Optimisation bounds must be scalars" " or array scalars.") if x1 > x2: raise ValueError("The lower bound exceeds the upper bound.") flag = 0 header = ' Func-count x f(x) Procedure' step=' initial' sqrt_eps = sqrt(2.2e-16) golden_mean = 0.5*(3.0-sqrt(5.0)) a, b = x1, x2 fulc = a + golden_mean*(b-a) nfc, xf = fulc, fulc rat = e = 0.0 x = xf fx = func(x,*args) num = 1 fmin_data = (1, xf, fx) ffulc = fnfc = fx xm = 0.5*(a+b) tol1 = sqrt_eps*abs(xf) + xtol / 3.0 tol2 = 2.0*tol1 if disp > 2: print (" ") print (header) print "%5.0f %12.6g %12.6g %s" % (fmin_data + (step,)) while ( abs(xf-xm) > (tol2 - 0.5*(b-a)) ): golden = 1 # Check for parabolic fit if abs(e) > tol1: golden = 0 r = (xf-nfc)*(fx-ffulc) q = (xf-fulc)*(fx-fnfc) p = (xf-fulc)*q - (xf-nfc)*r q = 2.0*(q-r) if q > 0.0: p = -p q = abs(q) r = e e = rat # Check for acceptability of parabola if ( (abs(p) < abs(0.5*q*r)) and (p > q*(a-xf)) and \ (p < q*(b-xf))): rat = (p+0.0) / q; x = xf + rat step = ' parabolic' if ((x-a) < tol2) or ((b-x) < tol2): si = numpy.sign(xm-xf) + ((xm-xf)==0) rat = tol1*si else: # do a golden section step golden = 1 if golden: # Do a golden-section step if xf >= xm: e=a-xf else: e=b-xf rat = golden_mean*e step = ' golden' si = numpy.sign(rat) + (rat == 0) x = xf + si*max([abs(rat), tol1]) fu = func(x,*args) num += 1 fmin_data = (num, x, fu) if disp > 2: print "%5.0f %12.6g %12.6g %s" % (fmin_data + (step,)) if fu <= fx: if x >= xf: a = xf else: b = xf fulc, ffulc = nfc, fnfc nfc, fnfc = xf, fx xf, fx = x, fu else: if x < xf: a = x else: b = x if (fu <= fnfc) or (nfc == xf): fulc, ffulc = nfc, fnfc nfc, fnfc = x, fu elif (fu <= ffulc) or (fulc == xf) or (fulc == nfc): fulc, ffulc = x, fu xm = 0.5*(a+b) tol1 = sqrt_eps*abs(xf) + xtol/3.0 tol2 = 2.0*tol1 if num >= maxfun: flag = 1 fval = fx if disp > 0: _endprint(x, flag, fval, maxfun, xtol, disp) if full_output: return xf, fval, flag, num else: return xf fval = fx if disp > 0: _endprint(x, flag, fval, maxfun, xtol, disp) if full_output: return xf, fval, flag, num else: return xf class Brent: #need to rethink design of __init__ def __init__(self, func, args=(), tol=1.48e-8, maxiter=500, full_output=0): self.func = func self.args = args self.tol = tol self.maxiter = maxiter self._mintol = 1.0e-11 self._cg = 0.3819660 self.xmin = None self.fval = None self.iter = 0 self.funcalls = 0 #need to rethink design of set_bracket (new options, etc) def set_bracket(self, brack = None): self.brack = brack def get_bracket_info(self): #set up func = self.func args = self.args brack = self.brack ### BEGIN core bracket_info code ### ### carefully DOCUMENT any CHANGES in core ## if brack is None: xa,xb,xc,fa,fb,fc,funcalls = bracket(func, args=args) elif len(brack) == 2: xa,xb,xc,fa,fb,fc,funcalls = bracket(func, xa=brack[0], xb=brack[1], args=args) elif len(brack) == 3: xa,xb,xc = brack if (xa > xc): # swap so xa < xc can be assumed dum = xa; xa=xc; xc=dum if not ((xa < xb) and (xb < xc)): raise ValueError("Not a bracketing interval.") fa = func(*((xa,)+args)) fb = func(*((xb,)+args)) fc = func(*((xc,)+args)) if not ((fb=xmid): deltax=a-x # do a golden section step else: deltax=b-x rat = _cg*deltax else: # do a parabolic step tmp1 = (x-w)*(fx-fv) tmp2 = (x-v)*(fx-fw) p = (x-v)*tmp2 - (x-w)*tmp1; tmp2 = 2.0*(tmp2-tmp1) if (tmp2 > 0.0): p = -p tmp2 = abs(tmp2) dx_temp = deltax deltax= rat # check parabolic fit if ((p > tmp2*(a-x)) and (p < tmp2*(b-x)) and (abs(p) < abs(0.5*tmp2*dx_temp))): rat = p*1.0/tmp2 # if parabolic step is useful. u = x + rat if ((u-a) < tol2 or (b-u) < tol2): if xmid-x >= 0: rat = tol1 else: rat = -tol1 else: if (x>=xmid): deltax=a-x # if it's not do a golden section step else: deltax=b-x rat = _cg*deltax if (abs(rat) < tol1): # update by at least tol1 if rat >= 0: u = x + tol1 else: u = x - tol1 else: u = x + rat fu = func(*((u,)+self.args)) # calculate new output value funcalls += 1 if (fu > fx): # if it's bigger than current if (u= x): a = x else: b = x v=w; w=x; x=u fv=fw; fw=fx; fx=fu iter += 1 ################################# #END CORE ALGORITHM ################################# self.xmin = x self.fval = fx self.iter = iter self.funcalls = funcalls def get_result(self, full_output=False): if full_output: return self.xmin, self.fval, self.iter, self.funcalls else: return self.xmin def brent(func, args=(), brack=None, tol=1.48e-8, full_output=0, maxiter=500): """Given a function of one-variable and a possible bracketing interval, return the minimum of the function isolated to a fractional precision of tol. Parameters ---------- func : callable f(x,*args) Objective function. args Additional arguments (if present). brack : tuple Triple (a,b,c) where (a xc): # swap so xa < xc can be assumed dum = xa; xa=xc; xc=dum if not ((xa < xb) and (xb < xc)): raise ValueError("Not a bracketing interval.") fa = func(*((xa,)+args)) fb = func(*((xb,)+args)) fc = func(*((xc,)+args)) if not ((fb abs(xb-xa)): x1 = xb x2 = xb + _gC*(xc-xb) else: x2 = xb x1 = xb - _gC*(xb-xa) f1 = func(*((x1,)+args)) f2 = func(*((x2,)+args)) funcalls += 2 while (abs(x3-x0) > tol*(abs(x1)+abs(x2))): if (f2 < f1): x0 = x1; x1 = x2; x2 = _gR*x1 + _gC*x3 f1 = f2; f2 = func(*((x2,)+args)) else: x3 = x2; x2 = x1; x1 = _gR*x2 + _gC*x0 f2 = f1; f1 = func(*((x1,)+args)) funcalls += 1 if (f1 < f2): xmin = x1 fval = f1 else: xmin = x2 fval = f2 if full_output: return xmin, fval, funcalls else: return xmin def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000): """Given a function and distinct initial points, search in the downhill direction (as defined by the initital points) and return new points xa, xb, xc that bracket the minimum of the function f(xa) > f(xb) < f(xc). It doesn't always mean that obtained solution will satisfy xa<=x<=xb Parameters ---------- func : callable f(x,*args) Objective function to minimize. xa, xb : float Bracketing interval. args : tuple Additional arguments (if present), passed to `func`. grow_limit : float Maximum grow limit. maxiter : int Maximum number of iterations to perform. Returns ------- xa, xb, xc : float Bracket. fa, fb, fc : float Objective function values in bracket. funcalls : int Number of function evaluations made. """ _gold = 1.618034 _verysmall_num = 1e-21 fa = func(*(xa,)+args) fb = func(*(xb,)+args) if (fa < fb): # Switch so fa > fb dum = xa; xa = xb; xb = dum dum = fa; fa = fb; fb = dum xc = xb + _gold*(xb-xa) fc = func(*((xc,)+args)) funcalls = 3 iter = 0 while (fc < fb): tmp1 = (xb - xa)*(fb-fc) tmp2 = (xb - xc)*(fb-fa) val = tmp2-tmp1 if abs(val) < _verysmall_num: denom = 2.0*_verysmall_num else: denom = 2.0*val w = xb - ((xb-xc)*tmp2-(xb-xa)*tmp1)/denom wlim = xb + grow_limit*(xc-xb) if iter > maxiter: raise RuntimeError("Too many iterations.") iter += 1 if (w-xc)*(xb-w) > 0.0: fw = func(*((w,)+args)) funcalls += 1 if (fw < fc): xa = xb; xb=w; fa=fb; fb=fw return xa, xb, xc, fa, fb, fc, funcalls elif (fw > fb): xc = w; fc=fw return xa, xb, xc, fa, fb, fc, funcalls w = xc + _gold*(xc-xb) fw = func(*((w,)+args)) funcalls += 1 elif (w-wlim)*(wlim-xc) >= 0.0: w = wlim fw = func(*((w,)+args)) funcalls += 1 elif (w-wlim)*(xc-w) > 0.0: fw = func(*((w,)+args)) funcalls += 1 if (fw < fc): xb=xc; xc=w; w=xc+_gold*(xc-xb) fb=fc; fc=fw; fw=func(*((w,)+args)) funcalls += 1 else: w = xc + _gold*(xc-xb) fw = func(*((w,)+args)) funcalls += 1 xa=xb; xb=xc; xc=w fa=fb; fb=fc; fc=fw return xa, xb, xc, fa, fb, fc, funcalls def _linesearch_powell(func, p, xi, tol=1e-3): """Line-search algorithm using fminbound. Find the minimium of the function ``func(x0+ alpha*direc)``. """ def myfunc(alpha): return func(p + alpha * xi) alpha_min, fret, iter, num = brent(myfunc, full_output=1, tol=tol) xi = alpha_min*xi return squeeze(fret), p+xi, xi def fmin_powell(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None, direc=None): """ Minimize a function using modified Powell's method. This method only uses function values, not derivatives. Parameters ---------- func : callable f(x,*args) Objective function to be minimized. x0 : ndarray Initial guess. args : tuple Extra arguments passed to func. callback : callable An optional user-supplied function, called after each iteration. Called as ``callback(xk)``, where ``xk`` is the current parameter vector. direc : ndarray Initial direction set. Returns ------- xopt : ndarray Parameter which minimizes `func`. fopt : number Value of function at minimum: ``fopt = func(xopt)``. direc : ndarray Current direction set. iter : int Number of iterations. funcalls : int Number of function calls made. warnflag : int Integer warning flag: 1 : Maximum number of function evaluations. 2 : Maximum number of iterations. allvecs : list List of solutions at each iteration. Other Parameters ---------------- xtol : float Line-search error tolerance. ftol : float Relative error in ``func(xopt)`` acceptable for convergence. maxiter : int Maximum number of iterations to perform. maxfun : int Maximum number of function evaluations to make. full_output : bool If True, fopt, xi, direc, iter, funcalls, and warnflag are returned. disp : bool If True, print convergence messages. retall : bool If True, return a list of the solution at each iteration. Notes ----- Uses a modification of Powell's method to find the minimum of a function of N variables. Powell's method is a conjugate direction method. The algorithm has two loops. The outer loop merely iterates over the inner loop. The inner loop minimizes over each current direction in the direction set. At the end of the inner loop, if certain conditions are met, the direction that gave the largest decrease is dropped and replaced with the difference between the current estiamted x and the estimated x from the beginning of the inner-loop. The technical conditions for replacing the direction of greatest increase amount to checking that 1. No further gain can be made along the direction of greatest increase from that iteration. 2. The direction of greatest increase accounted for a large sufficient fraction of the decrease in the function value from that iteration of the inner loop. References ---------- Powell M.J.D. (1964) An efficient method for finding the minimum of a function of several variables without calculating derivatives, Computer Journal, 7 (2):155-162. Press W., Teukolsky S.A., Vetterling W.T., and Flannery B.P.: Numerical Recipes (any edition), Cambridge University Press """ # we need to use a mutable object here that we can update in the # wrapper function fcalls, func = wrap_function(func, args) x = asarray(x0).flatten() if retall: allvecs = [x] N = len(x) rank = len(x.shape) if not -1 < rank < 2: raise ValueError("Initial guess must be a scalar or rank-1 sequence.") if maxiter is None: maxiter = N * 1000 if maxfun is None: maxfun = N * 1000 if direc is None: direc = eye(N, dtype=float) else: direc = asarray(direc, dtype=float) fval = squeeze(func(x)) x1 = x.copy() iter = 0; ilist = range(N) while True: fx = fval bigind = 0 delta = 0.0 for i in ilist: direc1 = direc[i] fx2 = fval fval, x, direc1 = _linesearch_powell(func, x, direc1, tol=xtol*100) if (fx2 - fval) > delta: delta = fx2 - fval bigind = i iter += 1 if callback is not None: callback(x) if retall: allvecs.append(x) if (2.0*(fx - fval) <= ftol*(abs(fx)+abs(fval))+1e-20): break if fcalls[0] >= maxfun: break if iter >= maxiter: break # Construct the extrapolated point direc1 = x - x1 x2 = 2*x - x1 x1 = x.copy() fx2 = squeeze(func(x2)) if (fx > fx2): t = 2.0*(fx+fx2-2.0*fval) temp = (fx-fval-delta) t *= temp*temp temp = fx-fx2 t -= delta*temp*temp if t < 0.0: fval, x, direc1 = _linesearch_powell(func, x, direc1, tol=xtol*100) direc[bigind] = direc[-1] direc[-1] = direc1 warnflag = 0 if fcalls[0] >= maxfun: warnflag = 1 if disp: print "Warning: Maximum number of function evaluations has "\ "been exceeded." elif iter >= maxiter: warnflag = 2 if disp: print "Warning: Maximum number of iterations has been exceeded" else: if disp: print "Optimization terminated successfully." print " Current function value: %f" % fval print " Iterations: %d" % iter print " Function evaluations: %d" % fcalls[0] x = squeeze(x) if full_output: retlist = x, fval, direc, iter, fcalls[0], warnflag if retall: retlist += (allvecs,) else: retlist = x if retall: retlist = (x, allvecs) return retlist def _endprint(x, flag, fval, maxfun, xtol, disp): if flag == 0: if disp > 1: print "\nOptimization terminated successfully;\n" \ "The returned value satisfies the termination criteria\n" \ "(using xtol = ", xtol, ")" if flag == 1: print "\nMaximum number of function evaluations exceeded --- " \ "increase maxfun argument.\n" return def brute(func, ranges, args=(), Ns=20, full_output=0, finish=fmin): """Minimize a function over a given range by brute force. Parameters ---------- func : callable ``f(x,*args)`` Objective function to be minimized. ranges : tuple Each element is a tuple of parameters or a slice object to be handed to ``numpy.mgrid``. args : tuple Extra arguments passed to function. Ns : int Default number of samples, if those are not provided. full_output : bool If True, return the evaluation grid. finish : callable, optional An optimization function that is called with the result of brute force minimization as initial guess. `finish` should take the initial guess as positional argument, and take take `args`, `full_output` and `disp` as keyword arguments. See Notes for more details. Returns ------- x0 : ndarray Value of arguments to `func`, giving minimum over the grid. fval : int Function value at minimum. grid : tuple Representation of the evaluation grid. It has the same length as x0. Jout : ndarray Function values over grid: ``Jout = func(*grid)``. Notes ----- The range is respected by the brute force minimization, but if the `finish` keyword specifies another optimization function (including the default `fmin`), the returned value may still be (just) outside the range. In order to ensure the range is specified, use ``finish=None``. """ N = len(ranges) if N > 40: raise ValueError("Brute Force not possible with more " \ "than 40 variables.") lrange = list(ranges) for k in range(N): if type(lrange[k]) is not type(slice(None)): if len(lrange[k]) < 3: lrange[k] = tuple(lrange[k]) + (complex(Ns),) lrange[k] = slice(*lrange[k]) if (N==1): lrange = lrange[0] def _scalarfunc(*params): params = squeeze(asarray(params)) return func(params,*args) vecfunc = vectorize(_scalarfunc) grid = mgrid[lrange] if (N==1): grid = (grid,) Jout = vecfunc(*grid) Nshape = shape(Jout) indx = argmin(Jout.ravel(),axis=-1) Nindx = zeros(N,int) xmin = zeros(N,float) for k in range(N-1,-1,-1): thisN = Nshape[k] Nindx[k] = indx % Nshape[k] indx = indx // thisN for k in range(N): xmin[k] = grid[k][tuple(Nindx)] Jmin = Jout[tuple(Nindx)] if (N==1): grid = grid[0] xmin = xmin[0] if callable(finish): vals = finish(func,xmin,args=args,full_output=1, disp=0) xmin = vals[0] Jmin = vals[1] if vals[-1] > 0: print "Warning: Final optimization did not succeed" if full_output: return xmin, Jmin, grid, Jout else: return xmin def main(): import time times = [] algor = [] x0 = [0.8,1.2,0.7] print "Nelder-Mead Simplex" print "===================" start = time.time() x = fmin(rosen,x0) print x times.append(time.time() - start) algor.append('Nelder-Mead Simplex\t') print print "Powell Direction Set Method" print "===========================" start = time.time() x = fmin_powell(rosen,x0) print x times.append(time.time() - start) algor.append('Powell Direction Set Method.') print print "Nonlinear CG" print "============" start = time.time() x = fmin_cg(rosen, x0, fprime=rosen_der, maxiter=200) print x times.append(time.time() - start) algor.append('Nonlinear CG \t') print print "BFGS Quasi-Newton" print "=================" start = time.time() x = fmin_bfgs(rosen, x0, fprime=rosen_der, maxiter=80) print x times.append(time.time() - start) algor.append('BFGS Quasi-Newton\t') print print "BFGS approximate gradient" print "=========================" start = time.time() x = fmin_bfgs(rosen, x0, gtol=1e-4, maxiter=100) print x times.append(time.time() - start) algor.append('BFGS without gradient\t') print print "Newton-CG with Hessian product" print "==============================" start = time.time() x = fmin_ncg(rosen, x0, rosen_der, fhess_p=rosen_hess_prod, maxiter=80) print x times.append(time.time() - start) algor.append('Newton-CG with hessian product') print print "Newton-CG with full Hessian" print "===========================" start = time.time() x = fmin_ncg(rosen, x0, rosen_der, fhess=rosen_hess, maxiter=80) print x times.append(time.time() - start) algor.append('Newton-CG with full hessian') print print "\nMinimizing the Rosenbrock function of order 3\n" print " Algorithm \t\t\t Seconds" print "===========\t\t\t =========" for k in range(len(algor)): print algor[k], "\t -- ", times[k] if __name__ == "__main__": main()