Added NelderMead method and settings for derivative-free optimization

Author: btracey

functionconvergence.go

@@ -0,0 +1,40 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package optimize
+
+import "math"
+
+// FunctionConverge tests for the convergence of function values. See comment
+// in Settings.
+type FunctionConverge struct {
+	Absolute   float64
+	Relative   float64
+	Iterations int
+
+	best float64
+	iter int
+}
+
+func (fc *FunctionConverge) Init(f float64) {
+	fc.best = f
+	fc.iter = 0
+}
+
+func (fc *FunctionConverge) FunctionConverged(f float64) Status {
+	if fc.Iterations == 0 {
+		return NotTerminated
+	}
+	maxAbs := math.Max(math.Abs(f), math.Abs(fc.best))
+	if f < fc.best && fc.best-f > fc.Relative*maxAbs+fc.Absolute {
+		fc.best = f
+		fc.iter = 0
+		return NotTerminated
+	}
+	fc.iter++
+	if fc.iter < fc.Iterations {
+		return NotTerminated
+	}
+	return FunctionConvergence
+}

gradfree_test.go

@@ -0,0 +1,100 @@
+// Copyright ©2015 The gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package optimize
+
+import (
+	"testing"
+
+	"github.com/gonum/floats"
+	"github.com/gonum/optimize/functions"
+)
+
+type gradFreeTest struct {
+	// f is the function that is being minimized.
+	f Function
+	// x is the initial guess.
+	x []float64
+	// absTol is the absolute function convergence for the test. If absTol == 0,
+	// the default value of 1e-6 will be used
+	absTol float64
+	// absIter is the number of iterations for function convergence. If iter == 0,
+	// the default value of 5 will be used
+	absIter int
+	// long indicates that the test takes long time to finish and will be
+	// excluded if testing.Short() is true.
+	long bool
+}
+
+// gradFree ensures that the function is gradient free
+type gradFree struct {
+	f Function
+}
+
+func (g gradFree) Func(x []float64) float64 {
+	return g.f.Func(x)
+}
+
+// makeGradFree ensures that a function contains no gradient method
+func makeGradFree(f Function) gradFree {
+	return gradFree{f}
+}
+
+// TODO(btracey): The gradient is still evaluated and tested if available
+// even if a gradient-free method is being used. This should be fixed. When that
+// is fixed, this should include Functions that also have gradients.
+var gradFreeTests = []gradFreeTest{
+	{
+		f: makeGradFree(functions.ExtendedRosenbrock{}),
+		x: []float64{-10, 10},
+	},
+	{
+		f: makeGradFree(functions.ExtendedRosenbrock{}),
+		x: []float64{-5, 4, 16, 3},
+	},
+}
+
+func TestLocalGradFree(t *testing.T) {
+	testLocalGradFree(t, gradFreeTests, nil)
+}
+
+func TestNelderMead(t *testing.T) {
+	testLocalGradFree(t, gradFreeTests, &NelderMead{})
+}
+
+func testLocalGradFree(t *testing.T, tests []gradFreeTest, method Method) {
+	for _, test := range tests {
+		if test.long && testing.Short() {
+			continue
+		}
+		settings := DefaultSettings()
+		settings.Recorder = nil
+		if test.absIter == 0 {
+			test.absIter = 5
+		}
+		if test.absTol == 0 {
+			test.absTol = 1e-6
+		}
+		result, err := Local(test.f, test.x, settings, method)
+		if err != nil {
+			t.Errorf("error finding minimum (%v) for \n%v", err, test)
+		}
+		if result == nil {
+			t.Errorf("nil result without error for:\n%v", test)
+			continue
+		}
+		if result.Status != FunctionConvergence {
+			t.Errorf("Status not %v, %v instead", FunctionConvergence, result.Status)
+		}
+
+		result2, err := Local(test.f, test.x, settings, method)
+		if err != nil {
+			t.Errorf("error finding minimum (%v) when reusing Method for \n%v", err, test)
+		}
+		if result.FuncEvaluations != result2.FuncEvaluations ||
+			result.F != result2.F || !floats.Equal(result.X, result2.X) {
+			t.Errorf("Different result when reuse method")
+		}
+	}
+}

local.go

@@ -89,6 +89,10 @@ func Local(f Function, initX []float64, settings *Settings, method Method) (*Res
 		return nil, err
 	}
 
+	if settings.FunctionConverge != nil {
+		settings.FunctionConverge.Init(optLoc.F)
+	}
+
 	// Runtime is the only Stats field that needs to be updated here.
 	stats.Runtime = time.Since(startTime)
 	// Send optLoc to Recorder before checking it for convergence.
@@ -198,8 +202,7 @@ func getDefaultMethod(funcInfo *functionInfo) Method {
 	if funcInfo.IsGradient || funcInfo.IsFunctionGradient {
 		return &BFGS{}
 	}
-	// TODO: Implement a gradient-free method
-	panic("optimize: gradient-free methods not yet coded")
+	return &NelderMead{}
 }
 
 // getStartingLocation allocates and initializes the starting location for the minimization.
@@ -256,12 +259,19 @@ func checkConvergence(loc *Location, iterType IterationType, stats *Stats, setti
 	if iterType == MajorIteration || iterType == InitIteration {
 		if loc.Gradient != nil {
 			norm := floats.Norm(loc.Gradient, math.Inf(1))
-			if norm < settings.GradientAbsTol {
-				return GradientThreshhold
+			if norm < settings.GradientThreshold {
+				return GradientThreshold
 			}
 		}
-		if loc.F < settings.FunctionAbsTol {
-			return FunctionThreshhold
+		if loc.F < settings.FunctionThreshold {
+			return FunctionThreshold
+		}
+	}
+
+	if iterType == MajorIteration && settings.FunctionConverge != nil {
+		status := settings.FunctionConverge.FunctionConverged(loc.F)
+		if status != NotTerminated {
+			return status
 		}
 	}
 

neldermead.go

@@ -1,4 +1,4 @@
-// Copyright ©2014 The gonum Authors. All rights reserved.
+// Copyright ©2015 The gonum Authors. All rights reserved.
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
 
@@ -49,7 +49,10 @@ func (n nmVertexSorter) Swap(i, j int) {
 //
 // If an initial simplex is provided, it is used and initLoc is ignored. If
 // InitialVertices and InitialValues are both nil, an initial simplex will be
-// generated automatically. If the simplex update parameters (Reflection, etc.)
+// generated automatically using the initial location as one vertex, and each
+// additional vertex as SimplexSize away in one dimension.
+//
+// If the simplex update parameters (Reflection, etc.)
 // are zero, they will be set automatically based on the dimension according to
 // the recommendations in
 //
@@ -61,6 +64,12 @@ type NelderMead struct {
 	Expansion       float64 // Expansion parameter (>1)
 	Contraction     float64 // Contraction parameter (>0, <1)
 	Shrink          float64 // Shrink parameter (>0, <1)
+	SimplexSize     float64 // size of auto-constructed initial simplex
+
+	reflection  float64
+	expansion   float64
+	contraction float64
+	shrink      float64
 
 	vertices [][]float64 // location of the vertices sorted in ascending f
 	values   []float64   // function values at the vertices sorted in ascending f
@@ -85,19 +94,27 @@ func (n *NelderMead) Init(loc *Location, f *FunctionInfo, xNext []float64) (Eval
 	n.centroid = resize(n.centroid, dim)
 	n.reflectedPoint = resize(n.reflectedPoint, dim)
 
+	if n.SimplexSize == 0 {
+		n.SimplexSize = 0.05
+	}
+
 	// Default parameter choices are chosen in a dimension-dependent way
 	// from http://www.webpages.uidaho.edu/~fuchang/res/ANMS.pdf
-	if n.Reflection == 0 {
-		n.Reflection = 1
+	n.reflection = n.Reflection
+	if n.reflection == 0 {
+		n.reflection = 1
 	}
-	if n.Expansion == 0 {
-		n.Expansion = 1 + 2.0/float64(dim)
+	n.expansion = n.Expansion
+	if n.expansion == 0 {
+		n.expansion = 1 + 2/float64(dim)
 	}
-	if n.Contraction == 0 {
-		n.Contraction = 0.75 - 1.0/(2.0*float64(dim))
+	n.contraction = n.Contraction
+	if n.contraction == 0 {
+		n.contraction = 0.75 - 1/(2*float64(dim))
 	}
-	if n.Shrink == 0 {
-		n.Shrink = 1 - 1.0/float64(dim)
+	n.shrink = n.Shrink
+	if n.shrink == 0 {
+		n.shrink = 1 - 1/float64(dim)
 	}
 
 	if n.InitialVertices != nil {
@@ -126,8 +143,7 @@ func (n *NelderMead) Init(loc *Location, f *FunctionInfo, xNext []float64) (Eval
 	n.values[dim] = loc.F
 	n.fillIdx = 0
 	copy(xNext, loc.X)
-	xNext[0] += 1
-	copy(n.vertices[0], xNext)
+	xNext[0] += n.SimplexSize
 	n.lastIter = nmInitialize
 	return FuncEvaluation, InitIteration, nil
 }
@@ -155,17 +171,16 @@ func (n *NelderMead) Iterate(loc *Location, xNext []float64) (EvaluationType, It
 	switch n.lastIter {
 	case nmInitialize:
 		n.values[n.fillIdx] = loc.F
+		copy(n.vertices[n.fillIdx], loc.X)
 		n.fillIdx++
 		if n.fillIdx == dim {
 			// Successfully finished building initial simplex.
 			sort.Sort(nmVertexSorter{n.vertices, n.values})
 			computeCentroid(n.vertices, n.centroid)
 			return n.returnNext(nmReflected, xNext)
 		}
-		copy(xNext, loc.X)
-		xNext[n.fillIdx-1] -= 1
-		xNext[n.fillIdx] += 1
-		copy(n.vertices[n.fillIdx], xNext)
+		copy(xNext, n.vertices[dim])
+		xNext[n.fillIdx] += n.SimplexSize
 		return FuncEvaluation, InitIteration, nil
 	case nmReflected:
 		n.reflectedValue = loc.F
@@ -228,13 +243,13 @@ func (n *NelderMead) returnNext(iter nmIterType, xNext []float64) (EvaluationTyp
 		var scale float64
 		switch iter {
 		case nmReflected:
-			scale = n.Reflection
+			scale = n.reflection
 		case nmExpanded:
-			scale = n.Reflection * n.Expansion
+			scale = n.reflection * n.expansion
 		case nmContractedOutside:
-			scale = n.Reflection * n.Contraction
+			scale = n.reflection * n.contraction
 		case nmContractedInside:
-			scale = -n.Contraction
+			scale = -n.contraction
 		}
 		floats.SubTo(xNext, n.centroid, n.vertices[dim])
 		floats.Scale(scale, xNext)
@@ -248,7 +263,7 @@ func (n *NelderMead) returnNext(iter nmIterType, xNext []float64) (EvaluationTyp
 	case nmShrink:
 		// x_shrink = x_best + delta * (x_i + x_best)
 		floats.SubTo(xNext, n.vertices[n.fillIdx], n.vertices[0])
-		floats.Scale(n.Shrink, xNext)
+		floats.Scale(n.shrink, xNext)
 		floats.Add(xNext, n.vertices[0])
 		return FuncEvaluation, SubIteration, nil
 	default:
@@ -280,3 +295,13 @@ func (n *NelderMead) replaceWorst(x []float64, f float64) {
 	floats.AddScaled(n.centroid, -1/float64(dim), n.vertices[dim])
 	floats.AddScaled(n.centroid, 1/float64(dim), x)
 }
+
+func (*NelderMead) Needs() struct {
+	Gradient bool
+	Hessian  bool
+} {
+	return struct {
+		Gradient bool
+		Hessian  bool
+	}{false, false}
+}

termination.go

@@ -13,10 +13,9 @@ type Status int
 const (
 	NotTerminated Status = iota
 	Success
-	FunctionAbsoluteConvergence
-	FunctionRelativeConvergence
-	FunctionThreshhold
-	GradientThreshhold
+	FunctionThreshold
+	FunctionConvergence
+	GradientThreshold
 	StepConvergence
 	FunctionNegativeInfinity
 	Failure
@@ -57,10 +56,13 @@ var statuses = []struct {
 		name: "Success",
 	},
 	{
-		name: "FunctionAbsoluteConvergence",
+		name: "FunctionThreshold",
 	},
 	{
-		name: "GradientAbsoluteConvergence",
+		name: "FunctionConvergence",
+	},
+	{
+		name: "GradientThreshold",
 	},
 	{
 		name: "StepConvergence",