# 4.6. pyopus.optimizer.nm — Unconstrained Melder-Mead simplex optimizer¶

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.

class pyopus.optimizer.nm.NelderMead(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-09, simplex=None, looseContraction=False)[source]

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 checkFtol() and 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 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 Optimizer class for more information.

buildSimplex(x0, rel=0.05, abs=0.00025)[source]

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 $$x_0^i \cdot rel$$ or abs if $$x_0^i$$ is zero.

check()[source]

Checks the optimization algorithm’s settings and raises an exception if something is wrong.

checkFtol()[source]

Checks the function value tolerance and returns True if the function values are within $$\max(ftol, reltol \cdot |f_{best}|)$$ of the point with the lowest cost function value ($$f_{best}$$).

checkXtol()[source]

Returns True if the components of points in the simplex are within $$max(reltol \cdot |x_{best}^i|, xtol)$$ of the corresponding components of the point with the lowest cost function value ($$x_{best}$$).

orderSimplex()[source]

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.

reset(x0)[source]

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 (ndim member).

The initial simplex is built around x0 by calling the 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.

run()[source]

Runs the optimization algorithm.

Example file nm.py in folder demo/optimizer/

# Optimize MGH unconstrained test suite with Nelder-Mead optimizer.

from pyopus.optimizer.base import Reporter
from pyopus.problems.mgh import *
from numpy import array, sqrt
from platform import system
if system()=='Windows':
# perf_counter() is the most precise timer in Windows
from time import perf_counter as timer
else:
# time() is the most precise timer in Linux
from time import time as timer
from sys import stdout

# Custom reporter that prints dots
class MyReporter(Reporter):
def __init__(self, name, iterSpacing=1, concise=False, printImprovement=True):
Reporter.__init__(self)
self.name=name
self.iterSpacing=iterSpacing
self.concise=concise
self.printImprovement=printImprovement
self.fbest=None
self.lastPrintout=None

def reset(self):
self.fbest=None

def __call__(self, x, f, opt):
report=False

if self.fbest is None:
self.lastPrintout=opt.niter
self.fbest=f
report=True

if self.fbest>f and self.printImprovement:
self.fbest=f
report=True

if opt.niter-self.lastPrintout>=self.iterSpacing:
report=True

if report:
if self.concise:
stdout.write(".")
stdout.flush()
else:
print("%30s (%2d) iter=%-6d f=%12.4e fbest=%12.4e" % (
self.name[:30], x.size, opt.niter, f, self.fbest
))
self.lastPrintout=opt.niter

if __name__=='__main__':
suite=[
[ Rosenbrock() ],
[ FreudensteinAndRoth() ],
[ Beale() ],
[ JennrichAndSampson() ],
[ McKinnon() ],
[ McKinnon(), array([[0.0, 0.0], [1.0, 1.0], [(1.0+sqrt(33.0))/8.0, (1.0-sqrt(33.0))/8.0]]) ],
[ HelicalValley() ],
[ Bard() ],
[ Gaussian() ],
[ Meyer() ],
[ GulfResearchAndDevelopement() ],
[ Box3D() ],
[ PowellSingular() ],
[ Wood() ],
[ KowalikAndOsborne() ],
[ BrownAndDennis() ],
[ PenaltyI(4) ],
[ PenaltyII(4) ],
[ Osborne1() ],
[ BrownAlmostLinear(5) ],
[ BiggsEXP6() ],
[ ExtendedRosenbrock(6) ],
[ BrownAlmostLinear(7) ],
[ ExtendedRosenbrock(8) ],
[ VariablyDimensioned(8) ],
[ ExtendedPowellSingular(8) ],
[ Watson() ],
[ ExtendedRosenbrock(10) ],
[ PenaltyI(10) ],
[ PenaltyII(10) ],
[ Trigonometric() ],
[ Osborne2() ],
[ ExtendedPowellSingular(12) ],
]

t1=timer()

results=[]

for probdef in suite:
# Take problem function.
prob=probdef[0]

# Write a message.
print("\nProcessing: "+prob.name+" ("+str(prob.n)+") ")

# Create optimizer object.

# Install custom reporter plugin. Print dot for 1000 evaluations.
opt.installPlugin(MyReporter(prob.name, 1000, concise=True, printImprovement=False))

# If problem has a custom initial simplex, set it.
if len(probdef)==1:
opt.reset(prob.initial)
else:
# Custom simplex (for McKinnon)
opt.reset(probdef[1])

# Start timing, run, measure time.
dt=timer()
opt.run()
dt=timer()-dt

# Write number of function evaluations.
print(" %d evaluations" % opt.niter)

# Calculate initial and final gradient
gini=prob.g(prob.initial)
gend=prob.g(opt.x)

# Store results for this problem.
result={
'i': opt.niter,
'x': opt.x,
'f': opt.f,
'gi': sqrt((gini**2).sum()),
'ge': sqrt((gend**2).sum()),
't': dt
}
results.append(result)

# Print summary. Last column is iterations per second.
print("\n")
for i in range(0, len(suite)):
prob=suite[i][0]
print("%30s (%2d) : ni=%6d f=%16.8e gi=%11.3e ge=%11.3e  it/s=%11.3e" % (
prob.name[:30],
prob.initial.size,
results[i]['i'],
results[i]['f'],
results[i]['gi'],
results[i]['ge'],
results[i]['i']/results[i]['t']
)
)

t2=timer()

print("\ntime=%f" % (t2-t1))