Source code for pyopus.optimizer.nm

.. inheritance-diagram:: pyopus.optimizer.nm
    :parts: 1

**Unconstrained Nelder-Mead optimizer (PyOPUS subsystem name: NMOPT)**

A very popular unconstrained optimization algorithm first published in [nm]_,

Unfortunately no convergence theory is available. There is even a 
counterexample available showing how the algorithm can fail. See [mck]_.

.. [nm] Nelder, J. A.,Mead, R.: A simplex method for function minimization. 
        Computer Journal, vol. 7, pp. 308-313, 1965. 

.. [mck] McKinnon, K. I. M.: Convergence of the Nelder-Mead Simplex Method to a 
        Nonstationary Point. SIAM Journal on Optimization, vol. 9, pp. 148-158, 1998. 

from ..misc.debug import DbgMsgOut, DbgMsg

from .base import Optimizer

from numpy import array, abs, lexsort, zeros, where

__all__ = [ 'NelderMead' ]

[docs]class NelderMead(Optimizer): """ Nelder-Mead optimizer class *reflect*, *expand*, *outerContract*, *innerContract*, and *shrink* are step size factors for the reflection, expansion, outer contraction, inner contraction, and shrink step, respectively. *expansion* must be above 1. *reflection* must be greater than 0 and smaller than *expansion*. *outerContraction* must be between 0 and 1, while *innerContraction* must be between -1 and 0. *shrink* must be between 0 and 1. *reltol* is the relative stopping tolerance. *ftol* and *xtol* are the absolute stopping tolerances for cost function values at simplex points and simlex side lengths. See the :meth:`checkFtol` and :meth:`checkXtol` methods. *simplex* is the initial simplex given by a (*ndim*+1) times *ndim* array where every row corresponds to one simplex point. If *simplex is not given an initial simplex is constructed around the initial point *x0*. See the :meth:`buildSimplex` method for details. If *looseContraction* is ``True`` the acceptance condition for contraction steps requres that the new point is not worse than the worst point. This is the behavior of the original algorithm. If this parameter is ``False`` (which is also the default) the new point is accepted if it is better than the worst point. See the :class:`~pyopus.optimizer.base.Optimizer` class for more information. """ # Note: shrink coefficient should be <0.5, because larger values may cause stagnation # due to roundoff errors (succesive shrinks do not result in a zero-diameter # simplex after infinite number of steps). If the value of the coefficient is 0.5 # or larger, roundoff errors may keep the simplex size at floating point precision # limit (relative tolerance 2**-52 = 2.22e-16) and it never reaches 0. def __init__(self, function, debug=0, fstop=None, maxiter=None, reflect=1.0, expand=2.0, outerContract=0.5, innerContract=-0.5, shrink=0.5, reltol=1e-15, ftol=1e-15, xtol=1e-9, simplex=None, looseContraction=False): Optimizer.__init__(self, function, debug, fstop, maxiter) # Coefficients self.reflect=reflect self.expand=expand self.outerContract=outerContract self.innerContract=innerContract self.shrink=shrink # Stopping condition self.reltol=reltol self.ftol=ftol self.xtol=xtol # Simplex self.simplex=simplex # Modifications self.looseContraction=looseContraction
[docs] def check(self): """ Checks the optimization algorithm's settings and raises an exception if something is wrong. """ Optimizer.check(self) if self.expand<=1.0: raise Exception(DbgMsg("NMOPT", "Expansion coefficient should be gerater than 1.")) if self.reflect>self.expand: raise Exception(DbgMsg("NMOPT", "Reflection coefficient should be smaller than expansion coefficient.")) if self.reflect<=0.0: raise Exception(DbgMsg("NMOPT", "Reflection coefficient should be greater than 0.")) if (self.outerContract<=0.0) or (self.outerContract>=self.reflect): raise Exception(DbgMsg("NMOPT", "Outer contraction coefficient should be between 0 and reflection coefficient.")) if (self.innerContract>=0.0) or (self.innerContract<=-1.0): raise Exception(DbgMsg("NMOPT", "Inner contraction coefficient must be from (-1,0).")) if (self.shrink<=0.0) or (self.shrink>=1.0): raise Exception(DbgMsg("NMOPT", "Shrink coefficient must be from (0,1).")) if self.reltol<0: raise Exception(DbgMsg("NMOPT", "Negative relative tolerance.")) if self.ftol<0: raise Exception(DbgMsg("NMOPT", "Negative f tolerance.")) if (self.xtol<0).any(): raise Exception(DbgMsg("NMOPT", "Negative x tolerance."))
def _setSimplex(self, sim): """ Sets the initial simplex to the array given by *sim* and checks it. """ self.npts=sim.shape[0] if sim.ndim!=2: raise Exception(DbgMsg("NMOPT", "Simplex must be a 2-dimensional array.")) if sim.shape[0]!=sim.shape[1]+1: raise Exception(DbgMsg("NMOPT", "Simplex must have dimension+1 points.")) self.simplexf=None self.simplex=sim
[docs] def buildSimplex(self, x0, rel=0.05, abs=0.00025): """ Builds an initial simplex around point given by a 1-dimensional array *x0*. *rel* and *abs* are used for the relative and absolute simplex size. The initial simplex has its first point positioned at *x0*. The *i*-th point among the remaining *ndim* points is obtained by movin along the *i*-th coordinate direction by :math:`x_0^i \\cdot rel` or *abs* if :math:`x_0^i` is zero. """ ndim=x0.shape[0] sim=zeros([ndim+1, ndim]) sim[0,:]=x0 for i in range(1,ndim+1): x=x0.copy() c=x[i-1] if c==0.0: x[i-1]+=abs else: x[i-1]+=c*rel sim[i,:]=x return sim
[docs] def orderSimplex(self): """ Reorders the points and the corresponding cost function values of the simplex in such way that the point with the lowest cost function value is the first point in the simplex. """ # Order simplex i=lexsort((-self.simplexmoves, self.simplexf), 0) self.simplexf=self.simplexf[i] self.simplex=self.simplex[i,:] self.simplexmoves=self.simplexmoves[i]
[docs] def checkFtol(self): """ Checks the function value tolerance and returns ``True`` if the function values are within :math:`\\max(ftol, reltol \\cdot |f_{best}|)` of the point with the lowest cost function value (:math:`f_{best}`). """ tol=max(self.ftol, self.reltol*abs(self.simplexf[0])) if abs(self.simplexf[1:]-self.simplexf[0]).max()<tol: return True else: return False
[docs] def checkXtol(self): """ Returns ``True`` if the components of points in the simplex are within :math:`max(reltol \\cdot |x_{best}^i|, xtol)` of the corresponding components of the point with the lowest cost function value (:math:`x_{best}`). """ tolr=self.xtol*abs(self.simplex[0,:]) tol=where(tolr>=self.xtol, tolr, self.xtol) if (abs(self.simplex[1:,:]-self.simplex[0,:]).max(0)<tol).all(): return True else: return False
[docs] def reset(self, x0): """ Puts the optimizer in its initial state and sets the initial point to be the 1-dimensional array *x0*. The length of the array becomes the dimension of the optimization problem (:attr:`ndim` member). The initial simplex is built around *x0* by calling the :meth:`buildSimplex` method with default values for the *rel* and *abs* arguments. If *x0* is a 2-dimensional array of size (*ndim*+1) times *ndim* it specifies the initial simplex. """ # Debug message if self.debug: DbgMsgOut("NM", "Resetting.") # Make it an array x0=array(x0) # Is x0 a point or a simplex? if x0.ndim==1: # Point # Set x now Optimizer.reset(self, x0) if self.debug: DbgMsgOut("NM", "Generating initial simplex from initial point.") sim=self.buildSimplex(x0) self._setSimplex(sim) else: # Simplex or error (handled in _setSimplex()) self._setSimplex(x0) if self.debug: DbgMsgOut("NM", "Using specified initial simplex.") # Set x to first point in simplex after it was checked in _setSimplex() Optimizer.reset(self, x0[0,:]) # Reset point moves counter self.simplexmoves=zeros(self.ndim+1) # Make x tolerance an array self.xtol=array(self.xtol) # Reset counters self.noc=0 self.nic=0 self.ns=0 self.nrok=0 self.neok=0 self.nocok=0 self.nicok=0 self.icconv=0 self.occonv=0
[docs] def run(self): """ Runs the optimization algorithm. """ # Debug message if self.debug: DbgMsgOut("NM", "Starting a run at i="+str(self.niter)) # Checks self.check() # Reset stop flag self.stop=False # Evaluate initial simplex if needed if self.simplexf is None: self.simplexf=zeros(self.npts) for i in range(0, self.ndim+1): self.simplexf[i][i,:]) if self.debug: DbgMsgOut("NM", "Initial simplex point i="+str(self.niter)+": f="+str(self.simplexf[i])) # Loop while not self.stop: # Order simplex (best point first) self.orderSimplex() # Centroid xc=self.simplex[:-1,:].sum(0)/self.ndim # Worst point xw=self.simplex[-1,:] fw=self.simplexf[-1] # Second worst point xsw=self.simplex[-2,:] fsw=self.simplexf[-2] # Best point xb=self.simplex[0,:] fb=self.simplexf[0] # No shrink shrink=False # Reflect xr=xc+(xc-xw)*self.reflect if self.debug: DbgMsgOut("NM", "Iteration i="+str(self.niter)+": reflect : f="+str(fr)) if fr<fb: # Try expansion xe=xc+(xc-xw)*self.expand if self.debug: DbgMsgOut("NM", "Iteration i="+str(self.niter)+": expand : f="+str(fe)) if fe<fr: # Accept expansion self.simplex[-1,:]=xe self.simplexf[-1]=fe self.simplexmoves[-1]+=1 self.neok+=1 else: # Accept reflection self.simplex[-1,:]=xr self.simplexf[-1]=fr self.simplexmoves[-1]+=1 self.nrok+=1 elif fb<=fr and fr<fsw: # Accept reflection self.simplex[-1,:]=xr self.simplexf[-1]=fr self.simplexmoves[-1]+=1 self.nrok+=1 elif fsw<=fr and fr<fw: # Try outer contraction xo=xc+(xc-xw)*self.outerContract self.noc+=1 if fo<((1+self.outerContract)*fw+(self.reflect-self.outerContract)*fr): self.occonv+=1 if self.debug: DbgMsgOut("NM", "Iteration i="+str(self.niter)+": outer con : f="+str(fo)) if fo<fw or (self.looseContraction and fo==fw): # Accept self.simplex[-1,:]=xo self.simplexf[-1]=fo self.simplexmoves[-1]+=1 self.nocok+=1 else: # Shrink shrink=True elif fw<=fr: # Try inner contraction xi=xc+(xc-xw)*self.innerContract self.nic+=1 if fi<((1+self.innerContract)*fw+(self.reflect-self.innerContract)*fr): self.icconv+=1 if self.debug: DbgMsgOut("NM", "Iteration i="+str(self.niter)+": inner con : f="+str(fi)) if fi<fw or (self.looseContraction and fi==fw): # Accept self.simplex[-1,:]=xi self.simplexf[-1]=fi self.simplexmoves[-1]+=1 self.nicok+=1 else: # Shrink shrink=True # Shrink if shrink: for i in range(1, self.ndim+1): xs=xb+(self.simplex[i,:]-xb)*self.shrink self.ns+=1 self.simplex[i,:]=xs self.simplexf[i]=fs self.simplexmoves[i]+=1 if self.debug: DbgMsgOut("NM", "Iteration i="+str(self.niter)+": shrink : f="+str(fs)) # Check stopping condition if self.checkFtol() and self.checkXtol(): if self.debug: DbgMsgOut("NM", "Iteration i="+str(self.niter)+": simplex x and f tolerance reached, stopping.") break # Debug message if self.debug: DbgMsgOut("NM", "Finished.")