# Nelder-Mead Simplex Optimization

The Nelder-Mead-Algorithm (also known as the “Simplex Algorithm” or even
as the “Amoeba Algorithm”) is an algorithm for the minimization of
non-linear functions in several variables. In contrast to other non-linear
minimization methods, it does *not* require gradient information. This
makes it less efficient, but also less prone to divergence problems. In
contrast to other methods, it is *not* necessary for the minimum to be
bracketed by the initial guess: the algorithm performs a limited “global”
search. (It may still converge to a local, rather than the global extremum,
of course.) Finally, the algorithm is fairly simple to implement as a
stand-alone routine, which makes it a natural choice for multi-dimensional
minimization *if function evaluations are not prohibitively expensive*.

Convention: The presentation in this post is phrased in terms of the
*minimization* of a function, for greater clarity. Everything is equally
applicable to *maximization* problems as well, with the objective function
*f* being replaced by *-f*.

### How It Works

The algorithm maintains a “simplex” of guesses for the location of the
minimum. A simplex is a convex shape that is spanned by *d+1* points,
where *d* is the dimension of the space. In other words, in 2 dimensions,
a simplex is a triangle, in 3 dimensions a simplex is a (generally
irregular) tetrahedron, and so on.

At each step, the algorithm selects the vertex with the greatest function value, and probes several alternative positions for that vertex. It selects the best of these positions, until the minimum has been sufficiently accurately determined.

The algorithm probes three candidate positions for the vertex in question.
All three lie on the line that passes through this vertex and through the
centroid of all the other vertices. (A “centroid” is the arithmetic mean
of vertices.) The first position is found by *reflection* of the initial
vertex through the centroid. The second and third positions are found by
*extension* and *contraction* along the same line. For instance, we may
choose a point that is twice as far from the centroid as the reflected
point (extension) or that is half-way between the centroid and the reflected
point (contraction). The algorithm chooses the best among those three
candidate points and replaces the initial vertex with that point. This
process is now repeated with the vertex that now corresponds to the
greatest function value, until the minimum has been determined.

It is because both the reflection and the extension step probe points
*outside* of the original simplex that the algorithm can find solutions
that are not bracketed by the original starting configuration. Instead,
the simplex “crawls down” the function, by extending and contracting
like an Amoeba (hence the name), until it reaches a (local) minimum.

### The Implementation

The implementation follows the description on Wikipedia quite closely.

```
package main
import (
"sort"
"fmt" )
// NelderMead finds the minimum of the function f for a given starting guess.
// - The first argument is the objective function f to minimize. The objective
// function takes a slice of length d, and returns a scalar float64.
// - The second argument is a slice of (d+1) slices, each of length d,
// representing the starting points of the iteration.
// - The last argument contains stopping conditions, as a slice of 3 values
// of type float64. The first is the maximum number of iteration steps (as
// a float, but interpreted as an int), the second is the maximum allowed
// relative deviation of the objective function at the centroid of the
// simplex from any of the centroid vertices (that is, f(p)/f(c)-1.0), the
// third is the maximum square distance of any of the vertices from the
// centroid.
// NelderMead returns the number of iteration steps taken (as an int), the
// value of the objective function at the centroid (as float64), and a slice
// of length d with the final centroid coordinates.
//
// NelderMead overwrites the values in the second argument (with the starting
// guess)!
func NelderMead( f func( []float64 ) float64, p [][]float64,
limits []float64 ) ( int, float64, []float64 ) {
a, b, c, d := 1.0, 2.0, 0.5, 0.5 // fixed parameters
maxitr, dplim, dflim := int(limits[0]), limits[1], limits[2]
n := len(p) - 1 // dimension - must have n+1 starting pts!
ctr := make( []float64, n, n ) // centroid
ref := make( []float64, n, n ) // reflection
ext := make( []float64, n, n ) // extension, expansion
con := make( []float64, n, n ) // contraction
for itr:=0; ; itr++ {
// sort points p by value of f
sort.Slice( p, func (i,j int) bool { return f(p[i]) < f(p[j]) } )
// centroid, without largest-valued pt
for i:=0; i<n; i++ { ctr[i] = 0.0 }
for j:=0; j<n-1; j++ {
for i:=0; i<n; i++ {
ctr[i] += p[j][i]
}
}
for i:=0; i<n; i++ { ctr[i] /= float64(n-1) }
fc := f(ctr)
// check termination
df := 0.0 // MAX (!) deviation of f(p) from fc
for j:=0; j<n+1; j++ {
tmp := f(p[j])/fc - 1.0 // ( f(p[j]) - fc )/fc
if tmp < 0.0 { tmp = -tmp } // abs value
if tmp > df { df = tmp }
}
dp := 0.0 // MAX of deviation of p from ctr
for j:=0; j<n+1; j++ {
tmp := 0.0
for i:=0; i<n; i++ {
tmp += ( p[j][i] - ctr[i] )*( p[j][i] - ctr[i] )
}
if tmp > dp { dp = tmp } // NOT relative, absolute bounding
}
if itr > maxitr || df < dflim || dp < dplim {
return itr, fc, ctr // terminate
}
// reflection
for i:=0; i<n; i++ {
ref[i] = ctr[i] + a * ( ctr[i] - p[n][i] )
}
fr := f(ref)
if f( p[0] ) <= fr && fr < f( p[n-1] ) {
for i:=0; i<n; i++ {
p[n][i] = ref[i]
}
continue
}
// extension
if fr < f(p[0]) {
for i:=0; i<n; i++ {
ext[i] = ctr[i] + b * ( ref[i] - ctr[i] )
}
if f(ext) < fr {
for i:=0; i<n; i++ {
p[n][i] = ext[i]
}
} else {
for i:=0; i<n; i++ {
p[n][i] = ref[i]
}
}
continue
}
// contraction
if fr < f( p[n] ) {
for i:=0; i<n; i++ {
con[i] = ctr[i] + c * ( ref[i] - ctr[i] ) // outside
}
if f(con) < fr {
for i:=0; i<n; i++ {
p[n][i] = con[i]
}
}
continue
} else {
for i:=0; i<n; i++ {
con[i] = ctr[i] + c * ( p[n][i] - ctr[i] ) // inside
}
if f(con) < f(p[n]) {
for i:=0; i<n; i++ {
p[n][i] = con[i]
}
}
continue
}
// shrink
for j:=1; j<n+1; j++ {
for i:=0; i<n; i++ {
p[j][i] = p[0][i] + d * ( p[j][i] - p[0][i] )
}
}
} // end for
panic( "Never get here!" )
}
```

The only somewhat delicate point concerns the question how to terminate the iteration. A number of different stopping criteria can be imagined, but some of them can be fooled into stopping the iteration prematurely.

In my opinion, the most robust condition is based on the *size* of the
simplex: an upper bound on the distances between the centroid and all
the vertices fixes the location of the minimum. If only the minimum
*value* of the function is needed, but not its location, then bounding
the variation of the function value across all vertices may lead to
faster termination, in particular for functions that are “flat” near
their minimum.

By contrast, my experience with conditions that focus on the “rate of improvement” (the reduction in the function value from one iteration step to the next, or the reduction in size of the simplex) is less good. The reason is that the present algorithm occasionally makes steps that predominantly reshape the simplex, without improving the objective function significantly. Stopping criteria that end the iteration when the “rate of improvement” is small are easily confused by such behavior and may end the iteration too early.

### Usage

The examples below show how to use the current implementation of the Nelder-Mead Simplex algorithm. Notice that the dimensionality of the space does not need to be specified explicitly; instead, it is inferred from the dimensionality of the starting points.

The examples below use a simple function, with a well-defined, isolated minimum, and the Rosenbrock function, which is often used as test function for optimization algorithms.

```
func simple() {
f := func( x []float64 ) float64 {
return (x[0]-.3)*(x[0]-.3) + 2.0*x[1]*x[1] + 0.5
}
start := [][]float64{ []float64{1,1}, []float64{2,1}, []float64{-1,2}, }
limits := []float64{50., 1e-9, 1e-9}
fmt.Println( "Simple, 2dim" )
fmt.Println( NelderMead( f, start, limits ) )
}
func rosenbrock( x []float64 ) float64 {
s := 0.0
for i:=0; i<len(x)-1; i++ {
s += 100.0*(x[i+1]-x[i])*(x[i+1]-x[i]) + (1.0 - x[i])*(1.0 - x[i])
}
return s
}
func rosenbrock_2d() {
start := [][]float64{ []float64{1,1}, []float64{2,1}, []float64{-1,2}, }
limits := []float64{50., 1.e-9, 1e-9}
fmt.Println( "Rosenbrock, 2dim" )
fmt.Println( NelderMead( rosenbrock, start, limits ) )
}
func rosenbrock_3d() {
start := [][]float64{ []float64{1,1,1}, []float64{2,1,2},
[]float64{-1,2,-1}, []float64{-1,-1,-1} }
limits := []float64{50., 1e-9, 1e-9}
fmt.Println( "Rosenbrock, 3dim" )
fmt.Println( NelderMead( rosenbrock, start, limits ) )
}
func main() {
simple()
rosenbrock_2d()
rosenbrock_3d()
}
```