optimize: Change initialization, remove Needser, and update Problem f… (#779)

* optimize: Change initialization, remove Needser, and update Problem function calls We need a better way to express the Hessian function call so that sparse Hessians can be provided. This change updates the Problem function definitions to allow an arbitrary Symmetric matrix. With this change, we need to change how Location is used, so that we do not allocate a SymDense. Once this location is changed, we no longer need Needser to allocate the appropriate memory, and can shift that to initialization, further simplifying the interfaces. A 'fake' Problem is passed to Method to continue to make it impossible for the Method to call the functions directly. Fixes #727, #593.
2025-09-26 19:21:17 +08:00 · 2019-02-01 15:26:26 +00:00
parent 199b7405a3
commit c07f678f3f
20 changed files with 427 additions and 192 deletions
--- a/optimize/bfgs.go
+++ b/optimize/bfgs.go
@@ -10,6 +10,11 @@ import (
 	"gonum.org/v1/gonum/mat"
 )

+var (
+	_ Method      = (*BFGS)(nil)
+	_ localMethod = (*BFGS)(nil)
+)
+
 // BFGS implements the Broyden–Fletcher–Goldfarb–Shanno optimization method. It
 // is a quasi-Newton method that performs successive rank-one updates to an
 // estimate of the inverse Hessian of the objective function. It exhibits
@@ -46,6 +51,10 @@ func (b *BFGS) Status() (Status, error) {
 	return b.status, b.err
 }

+func (*BFGS) Uses(has Available) (uses Available, err error) {
+	return has.gradient()
+}
+
 func (b *BFGS) Init(dim, tasks int) int {
 	b.status = NotTerminated
 	b.err = nil
@@ -172,7 +181,7 @@ func (b *BFGS) NextDirection(loc *Location, dir []float64) (stepSize float64) {
 	return 1
 }

-func (*BFGS) Needs() struct {
+func (*BFGS) needs() struct {
 	Gradient bool
 	Hessian  bool
 } {
--- a/optimize/cg.go
+++ b/optimize/cg.go
@@ -15,6 +15,11 @@ const (
 	angleRestartThreshold  = -0.9
 )

+var (
+	_ Method      = (*CG)(nil)
+	_ localMethod = (*CG)(nil)
+)
+
 // CGVariant calculates the scaling parameter, β, used for updating the
 // conjugate direction in the nonlinear conjugate gradient (CG) method.
 type CGVariant interface {
@@ -109,6 +114,10 @@ func (cg *CG) Status() (Status, error) {
 	return cg.status, cg.err
 }

+func (*CG) Uses(has Available) (uses Available, err error) {
+	return has.gradient()
+}
+
 func (cg *CG) Init(dim, tasks int) int {
 	cg.status = NotTerminated
 	cg.err = nil
@@ -234,7 +243,7 @@ func (cg *CG) NextDirection(loc *Location, dir []float64) (stepSize float64) {
 	return stepSize
 }

-func (*CG) Needs() struct {
+func (*CG) needs() struct {
 	Gradient bool
 	Hessian  bool
 } {
--- a/optimize/cmaes.go
+++ b/optimize/cmaes.go
@@ -116,10 +116,6 @@ var (
 	_ Method   = (*CmaEsChol)(nil)
 )

-func (cma *CmaEsChol) Needs() struct{ Gradient, Hessian bool } {
-	return struct{ Gradient, Hessian bool }{false, false}
-}
-
 func (cma *CmaEsChol) methodConverged() Status {
 	sd := cma.StopLogDet
 	switch {
@@ -142,6 +138,10 @@ func (cma *CmaEsChol) Status() (Status, error) {
 	return cma.methodConverged(), nil
 }

+func (*CmaEsChol) Uses(has Available) (uses Available, err error) {
+	return has.function()
+}
+
 func (cma *CmaEsChol) Init(dim, tasks int) int {
 	if dim <= 0 {
 		panic(nonpositiveDimension)
--- a/optimize/errors.go
+++ b/optimize/errors.go
@@ -31,6 +31,14 @@ var (
 	// ErrLinesearcherBound signifies that a Linesearcher reached a step that
 	// lies out of allowed bounds.
 	ErrLinesearcherBound = errors.New("linesearch: step out of bounds")
+
+	// ErrMissingGrad signifies that a Method requires a Gradient function that
+	// is not supplied by Problem.
+	ErrMissingGrad = errors.New("optimize: problem does not provide needed Grad function")
+
+	// ErrMissingHess signifies that a Method requires a Hessian function that
+	// is not supplied by Problem.
+	ErrMissingHess = errors.New("optimize: problem does not provide needed Hess function")
 )

 // ErrFunc is returned when an initial function value is invalid. The error
--- a/optimize/functions/functions.go
+++ b/optimize/functions/functions.go
@@ -36,10 +36,13 @@ func (Beale) Func(x []float64) float64 {
 	return f1*f1 + f2*f2 + f3*f3
 }

-func (Beale) Grad(grad, x []float64) {
+func (Beale) Grad(grad, x []float64) []float64 {
 	if len(x) != 2 {
 		panic("dimension of the problem must be 2")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -54,12 +57,16 @@ func (Beale) Grad(grad, x []float64) {

 	grad[0] = -2 * (f1*t1 + f2*t2 + f3*t3)
 	grad[1] = 2 * x[0] * (f1 + 2*f2*x[1] + 3*f3*x[1]*x[1])
+	return grad
 }

-func (Beale) Hess(hess mat.MutableSymmetric, x []float64) {
+func (Beale) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
 	if len(x) != 2 {
 		panic("dimension of the problem must be 2")
 	}
+	if hess == nil {
+		hess = mat.NewSymDense(len(x), nil)
+	}
 	if len(x) != hess.Symmetric() {
 		panic("incorrect size of the Hessian")
 	}
@@ -74,9 +81,11 @@ func (Beale) Hess(hess mat.MutableSymmetric, x []float64) {
 	h00 := 2 * (t1*t1 + t2*t2 + t3*t3)
 	h01 := 2 * (f1 + x[1]*(2*f2+3*x[1]*f3) - x[0]*(t1+x[1]*(2*t2+3*x[1]*t3)))
 	h11 := 2 * x[0] * (x[0] + 2*f2 + x[1]*(6*f3+x[0]*x[1]*(4+9*x[1]*x[1])))
-	hess.SetSym(0, 0, h00)
-	hess.SetSym(0, 1, h01)
-	hess.SetSym(1, 1, h11)
+	h := hess.(*mat.SymDense)
+	h.SetSym(0, 0, h00)
+	h.SetSym(0, 1, h01)
+	h.SetSym(1, 1, h11)
+	return h
 }

 func (Beale) Minima() []Minimum {
@@ -113,10 +122,13 @@ func (BiggsEXP2) Func(x []float64) (sum float64) {
 	return sum
 }

-func (BiggsEXP2) Grad(grad, x []float64) {
+func (BiggsEXP2) Grad(grad, x []float64) []float64 {
 	if len(x) != 2 {
 		panic("dimension of the problem must be 2")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -135,6 +147,7 @@ func (BiggsEXP2) Grad(grad, x []float64) {
 		grad[0] += 2 * f * dfdx0
 		grad[1] += 2 * f * dfdx1
 	}
+	return grad
 }

 func (BiggsEXP2) Minima() []Minimum {
@@ -171,10 +184,13 @@ func (BiggsEXP3) Func(x []float64) (sum float64) {
 	return sum
 }

-func (BiggsEXP3) Grad(grad, x []float64) {
+func (BiggsEXP3) Grad(grad, x []float64) []float64 {
 	if len(x) != 3 {
 		panic("dimension of the problem must be 3")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -195,6 +211,7 @@ func (BiggsEXP3) Grad(grad, x []float64) {
 		grad[1] += 2 * f * dfdx1
 		grad[2] += 2 * f * dfdx2
 	}
+	return grad
 }

 func (BiggsEXP3) Minima() []Minimum {
@@ -231,10 +248,13 @@ func (BiggsEXP4) Func(x []float64) (sum float64) {
 	return sum
 }

-func (BiggsEXP4) Grad(grad, x []float64) {
+func (BiggsEXP4) Grad(grad, x []float64) []float64 {
 	if len(x) != 4 {
 		panic("dimension of the problem must be 4")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -257,6 +277,7 @@ func (BiggsEXP4) Grad(grad, x []float64) {
 		grad[2] += 2 * f * dfdx2
 		grad[3] += 2 * f * dfdx3
 	}
+	return grad
 }

 func (BiggsEXP4) Minima() []Minimum {
@@ -293,10 +314,13 @@ func (BiggsEXP5) Func(x []float64) (sum float64) {
 	return sum
 }

-func (BiggsEXP5) Grad(grad, x []float64) {
+func (BiggsEXP5) Grad(grad, x []float64) []float64 {
 	if len(x) != 5 {
 		panic("dimension of the problem must be 5")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -321,6 +345,7 @@ func (BiggsEXP5) Grad(grad, x []float64) {
 		grad[3] += 2 * f * dfdx3
 		grad[4] += 2 * f * dfdx4
 	}
+	return grad
 }

 func (BiggsEXP5) Minima() []Minimum {
@@ -360,10 +385,13 @@ func (BiggsEXP6) Func(x []float64) (sum float64) {
 	return sum
 }

-func (BiggsEXP6) Grad(grad, x []float64) {
+func (BiggsEXP6) Grad(grad, x []float64) []float64 {
 	if len(x) != 6 {
 		panic("dimension of the problem must be 6")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -390,6 +418,7 @@ func (BiggsEXP6) Grad(grad, x []float64) {
 		grad[4] += 2 * f * dfdx4
 		grad[5] += 2 * f * dfdx5
 	}
+	return grad
 }

 func (BiggsEXP6) Minima() []Minimum {
@@ -440,10 +469,13 @@ func (Box3D) Func(x []float64) (sum float64) {
 	return sum
 }

-func (Box3D) Grad(grad, x []float64) {
+func (Box3D) Grad(grad, x []float64) []float64 {
 	if len(x) != 3 {
 		panic("dimension of the problem must be 3")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -459,6 +491,7 @@ func (Box3D) Grad(grad, x []float64) {
 		grad[1] += -2 * f * c * math.Exp(c*x[1])
 		grad[2] += -2 * f * y
 	}
+	return grad
 }

 func (Box3D) Minima() []Minimum {
@@ -543,10 +576,13 @@ func (BrownBadlyScaled) Func(x []float64) float64 {
 	return f1*f1 + f2*f2 + f3*f3
 }

-func (BrownBadlyScaled) Grad(grad, x []float64) {
+func (BrownBadlyScaled) Grad(grad, x []float64) []float64 {
 	if len(x) != 2 {
 		panic("dimension of the problem must be 2")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -556,12 +592,16 @@ func (BrownBadlyScaled) Grad(grad, x []float64) {
 	f3 := x[0]*x[1] - 2
 	grad[0] = 2*f1 + 2*f3*x[1]
 	grad[1] = 2*f2 + 2*f3*x[0]
+	return grad
 }

-func (BrownBadlyScaled) Hess(hess mat.MutableSymmetric, x []float64) {
+func (BrownBadlyScaled) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
 	if len(x) != 2 {
 		panic("dimension of the problem must be 2")
 	}
+	if hess == nil {
+		hess = mat.NewSymDense(len(x), nil)
+	}
 	if len(x) != hess.Symmetric() {
 		panic("incorrect size of the Hessian")
 	}
@@ -569,9 +609,11 @@ func (BrownBadlyScaled) Hess(hess mat.MutableSymmetric, x []float64) {
 	h00 := 2 + 2*x[1]*x[1]
 	h01 := 4*x[0]*x[1] - 4
 	h11 := 2 + 2*x[0]*x[0]
-	hess.SetSym(0, 0, h00)
-	hess.SetSym(0, 1, h01)
-	hess.SetSym(1, 1, h11)
+	h := hess.(*mat.SymDense)
+	h.SetSym(0, 0, h00)
+	h.SetSym(0, 1, h01)
+	h.SetSym(1, 1, h11)
+	return h
 }

 func (BrownBadlyScaled) Minima() []Minimum {
@@ -612,10 +654,13 @@ func (BrownAndDennis) Func(x []float64) (sum float64) {
 	return sum
 }

-func (BrownAndDennis) Grad(grad, x []float64) {
+func (BrownAndDennis) Grad(grad, x []float64) []float64 {
 	if len(x) != 4 {
 		panic("dimension of the problem must be 4")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -633,19 +678,24 @@ func (BrownAndDennis) Grad(grad, x []float64) {
 		grad[2] += 4 * f * f2
 		grad[3] += 4 * f * f2 * math.Sin(c)
 	}
+	return grad
 }

-func (BrownAndDennis) Hess(hess mat.MutableSymmetric, x []float64) {
+func (BrownAndDennis) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
 	if len(x) != 4 {
 		panic("dimension of the problem must be 4")
 	}
+	if hess == nil {
+		hess = mat.NewSymDense(len(x), nil)
+	}
 	if len(x) != hess.Symmetric() {
 		panic("incorrect size of the Hessian")
 	}

+	h := hess.(*mat.SymDense)
 	for i := 0; i < 4; i++ {
 		for j := i; j < 4; j++ {
-			hess.SetSym(i, j, 0)
+			h.SetSym(i, j, 0)
 		}
 	}
 	for i := 1; i <= 20; i++ {
@@ -657,22 +707,23 @@ func (BrownAndDennis) Hess(hess mat.MutableSymmetric, x []float64) {
 		s3 := 2 * t1 * t2
 		r1 := t + 2*t1*t1
 		r2 := t + 2*t2*t2
-		hess.SetSym(0, 0, hess.At(0, 0)+r1)
-		hess.SetSym(0, 1, hess.At(0, 1)+d1*r1)
-		hess.SetSym(1, 1, hess.At(1, 1)+d1*d1*r1)
-		hess.SetSym(0, 2, hess.At(0, 2)+s3)
-		hess.SetSym(1, 2, hess.At(1, 2)+d1*s3)
-		hess.SetSym(2, 2, hess.At(2, 2)+r2)
-		hess.SetSym(0, 3, hess.At(0, 3)+d2*s3)
-		hess.SetSym(1, 3, hess.At(1, 3)+d1*d2*s3)
-		hess.SetSym(2, 3, hess.At(2, 3)+d2*r2)
-		hess.SetSym(3, 3, hess.At(3, 3)+d2*d2*r2)
+		h.SetSym(0, 0, h.At(0, 0)+r1)
+		h.SetSym(0, 1, h.At(0, 1)+d1*r1)
+		h.SetSym(1, 1, h.At(1, 1)+d1*d1*r1)
+		h.SetSym(0, 2, h.At(0, 2)+s3)
+		h.SetSym(1, 2, h.At(1, 2)+d1*s3)
+		h.SetSym(2, 2, h.At(2, 2)+r2)
+		h.SetSym(0, 3, h.At(0, 3)+d2*s3)
+		h.SetSym(1, 3, h.At(1, 3)+d1*d2*s3)
+		h.SetSym(2, 3, h.At(2, 3)+d2*r2)
+		h.SetSym(3, 3, h.At(3, 3)+d2*d2*r2)
 	}
 	for i := 0; i < 4; i++ {
 		for j := i; j < 4; j++ {
-			hess.SetSym(i, j, 4*hess.At(i, j))
+			h.SetSym(i, j, 4*h.At(i, j))
 		}
 	}
+	return h
 }

 func (BrownAndDennis) Minima() []Minimum {
@@ -716,10 +767,13 @@ func (ExtendedPowellSingular) Func(x []float64) (sum float64) {
 	return sum
 }

-func (ExtendedPowellSingular) Grad(grad, x []float64) {
+func (ExtendedPowellSingular) Grad(grad, x []float64) []float64 {
 	if len(x)%4 != 0 {
 		panic("dimension of the problem must be a multiple of 4")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -736,6 +790,7 @@ func (ExtendedPowellSingular) Grad(grad, x []float64) {
 		grad[i+2] = 10*f2 - 8*f3*t1
 		grad[i+3] = -10*f2 - 40*f4*t2
 	}
+	return grad
 }

 func (ExtendedPowellSingular) Minima() []Minimum {
@@ -780,7 +835,10 @@ func (ExtendedRosenbrock) Func(x []float64) (sum float64) {
 	return sum
 }

-func (ExtendedRosenbrock) Grad(grad, x []float64) {
+func (ExtendedRosenbrock) Grad(grad, x []float64) []float64 {
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -796,6 +854,7 @@ func (ExtendedRosenbrock) Grad(grad, x []float64) {
 	for i := 1; i < dim; i++ {
 		grad[i] += 200 * (x[i] - x[i-1]*x[i-1])
 	}
+	return grad
 }

 func (ExtendedRosenbrock) Minima() []Minimum {
@@ -902,10 +961,13 @@ func (g Gaussian) Func(x []float64) (sum float64) {
 	return sum
 }

-func (g Gaussian) Grad(grad, x []float64) {
+func (g Gaussian) Grad(grad, x []float64) []float64 {
 	if len(x) != 3 {
 		panic("dimension of the problem must be 3")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -923,6 +985,7 @@ func (g Gaussian) Grad(grad, x []float64) {
 		grad[1] -= f * e * d * x[0]
 		grad[2] += 2 * f * e * x[0] * x[1] * b
 	}
+	return grad
 }

 func (Gaussian) Minima() []Minimum {
@@ -964,10 +1027,13 @@ func (GulfResearchAndDevelopment) Func(x []float64) (sum float64) {
 	return sum
 }

-func (GulfResearchAndDevelopment) Grad(grad, x []float64) {
+func (GulfResearchAndDevelopment) Grad(grad, x []float64) []float64 {
 	if len(x) != 3 {
 		panic("dimension of the problem must be 3")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -989,6 +1055,7 @@ func (GulfResearchAndDevelopment) Grad(grad, x []float64) {
 	grad[0] *= 2 / x[0]
 	grad[1] *= 2 * x[2]
 	grad[2] *= 2
+	return grad
 }

 func (GulfResearchAndDevelopment) Minima() []Minimum {
@@ -1042,10 +1109,13 @@ func (HelicalValley) Func(x []float64) float64 {
 	return f1*f1 + f2*f2 + f3*f3
 }

-func (HelicalValley) Grad(grad, x []float64) {
+func (HelicalValley) Grad(grad, x []float64) []float64 {
 	if len(x) != 3 {
 		panic("dimension of the problem must be 3")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1064,6 +1134,7 @@ func (HelicalValley) Grad(grad, x []float64) {
 	grad[0] = 200 * (5*s*q*x[1] + (h-1)*r*x[0])
 	grad[1] = 200 * (-5*s*q*x[0] + (h-1)*r*x[1])
 	grad[2] = 2 * (100*s + x[2])
+	return grad
 }

 func (HelicalValley) Minima() []Minimum {
@@ -1083,7 +1154,7 @@ func (Linear) Func(x []float64) float64 {
 	return floats.Sum(x)
 }

-func (Linear) Grad(grad, x []float64) {
+func (Linear) Grad(grad, x []float64) []float64 {
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1091,6 +1162,7 @@ func (Linear) Grad(grad, x []float64) {
 	for i := range grad {
 		grad[i] = 1
 	}
+	return grad
 }

 // PenaltyI implements the first penalty function by Gill, Murray and Pitfield.
@@ -1120,7 +1192,10 @@ func (PenaltyI) Func(x []float64) (sum float64) {
 	return sum
 }

-func (PenaltyI) Grad(grad, x []float64) {
+func (PenaltyI) Grad(grad, x []float64) []float64 {
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1132,6 +1207,7 @@ func (PenaltyI) Grad(grad, x []float64) {
 	for i, v := range x {
 		grad[i] = 2 * (2*s*v + 1e-5*(v-1))
 	}
+	return grad
 }

 func (PenaltyI) Minima() []Minimum {
@@ -1185,7 +1261,10 @@ func (PenaltyII) Func(x []float64) (sum float64) {
 	return sum
 }

-func (PenaltyII) Grad(grad, x []float64) {
+func (PenaltyII) Grad(grad, x []float64) []float64 {
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1209,6 +1288,7 @@ func (PenaltyII) Grad(grad, x []float64) {
 		grad[i] += 1e-5 * f * math.Exp(x[i]/10) / 5
 	}
 	grad[0] += 2 * (x[0] - 0.2)
+	return grad
 }

 func (PenaltyII) Minima() []Minimum {
@@ -1254,10 +1334,13 @@ func (PowellBadlyScaled) Func(x []float64) float64 {
 	return f1*f1 + f2*f2
 }

-func (PowellBadlyScaled) Grad(grad, x []float64) {
+func (PowellBadlyScaled) Grad(grad, x []float64) []float64 {
 	if len(x) != 2 {
 		panic("dimension of the problem must be 2")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1266,12 +1349,16 @@ func (PowellBadlyScaled) Grad(grad, x []float64) {
 	f2 := math.Exp(-x[0]) + math.Exp(-x[1]) - 1.0001
 	grad[0] = 2 * (1e4*f1*x[1] - f2*math.Exp(-x[0]))
 	grad[1] = 2 * (1e4*f1*x[0] - f2*math.Exp(-x[1]))
+	return grad
 }

-func (PowellBadlyScaled) Hess(hess mat.MutableSymmetric, x []float64) {
+func (PowellBadlyScaled) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
 	if len(x) != 2 {
 		panic("dimension of the problem must be 2")
 	}
+	if hess == nil {
+		hess = mat.NewSymDense(len(x), nil)
+	}
 	if len(x) != hess.Symmetric() {
 		panic("incorrect size of the Hessian")
 	}
@@ -1281,12 +1368,14 @@ func (PowellBadlyScaled) Hess(hess mat.MutableSymmetric, x []float64) {
 	s2 := math.Exp(-x[1])
 	t2 := s1 + s2 - 1.0001

+	h := hess.(*mat.SymDense)
 	h00 := 2 * (1e8*x[1]*x[1] + s1*(s1+t2))
 	h01 := 2 * (1e4*(1+2*t1) + s1*s2)
 	h11 := 2 * (1e8*x[0]*x[0] + s2*(s2+t2))
-	hess.SetSym(0, 0, h00)
-	hess.SetSym(0, 1, h01)
-	hess.SetSym(1, 1, h11)
+	h.SetSym(0, 0, h00)
+	h.SetSym(0, 1, h01)
+	h.SetSym(1, 1, h11)
+	return h
 }

 func (PowellBadlyScaled) Minima() []Minimum {
@@ -1324,7 +1413,10 @@ func (Trigonometric) Func(x []float64) (sum float64) {
 	return sum
 }

-func (Trigonometric) Grad(grad, x []float64) {
+func (Trigonometric) Grad(grad, x []float64) []float64 {
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1344,6 +1436,7 @@ func (Trigonometric) Grad(grad, x []float64) {
 	for i, v := range x {
 		grad[i] += 2 * s2 * math.Sin(v)
 	}
+	return grad
 }

 func (Trigonometric) Minima() []Minimum {
@@ -1396,7 +1489,10 @@ func (VariablyDimensioned) Func(x []float64) (sum float64) {
 	return sum
 }

-func (VariablyDimensioned) Grad(grad, x []float64) {
+func (VariablyDimensioned) Grad(grad, x []float64) []float64 {
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1408,6 +1504,7 @@ func (VariablyDimensioned) Grad(grad, x []float64) {
 	for i, v := range x {
 		grad[i] = 2 * (v - 1 + s*float64(i+1)*(1+2*s*s))
 	}
+	return grad
 }

 func (VariablyDimensioned) Minima() []Minimum {
@@ -1480,7 +1577,10 @@ func (Watson) Func(x []float64) (sum float64) {
 	return sum
 }

-func (Watson) Grad(grad, x []float64) {
+func (Watson) Grad(grad, x []float64) []float64 {
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1516,17 +1616,21 @@ func (Watson) Grad(grad, x []float64) {
 	t := x[1] - x[0]*x[0] - 1
 	grad[0] += x[0] * (2 - 4*t)
 	grad[1] += 2 * t
+	return grad
 }

-func (Watson) Hess(hess mat.MutableSymmetric, x []float64) {
+func (Watson) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
 	dim := len(x)
+	if hess == nil {
+		hess = mat.NewSymDense(len(x), nil)
+	}
 	if dim != hess.Symmetric() {
 		panic("incorrect size of the Hessian")
 	}
-
+	h := hess.(*mat.SymDense)
 	for j := 0; j < dim; j++ {
 		for k := j; k < dim; k++ {
-			hess.SetSym(j, k, 0)
+			h.SetSym(j, k, 0)
 		}
 	}
 	for i := 1; i <= 29; i++ {
@@ -1553,16 +1657,17 @@ func (Watson) Hess(hess mat.MutableSymmetric, x []float64) {
 			v := float64(j) - s3
 			d3 := 1 / d1
 			for k := 0; k <= j; k++ {
-				hess.SetSym(k, j, hess.At(k, j)+d2*d3*(v*(float64(k)-s3)-th))
+				h.SetSym(k, j, h.At(k, j)+d2*d3*(v*(float64(k)-s3)-th))
 				d3 *= d1
 			}
 			d2 *= d1
 		}
 	}
 	t1 := x[1] - x[0]*x[0] - 1
-	hess.SetSym(0, 0, hess.At(0, 0)+8*x[0]*x[0]+2-4*t1)
-	hess.SetSym(0, 1, hess.At(0, 1)-4*x[0])
-	hess.SetSym(1, 1, hess.At(1, 1)+2)
+	h.SetSym(0, 0, h.At(0, 0)+8*x[0]*x[0]+2-4*t1)
+	h.SetSym(0, 1, h.At(0, 1)-4*x[0])
+	h.SetSym(1, 1, h.At(1, 1)+2)
+	return h
 }

 func (Watson) Minima() []Minimum {
@@ -1618,10 +1723,13 @@ func (Wood) Func(x []float64) (sum float64) {
 	return 100*f1*f1 + f2*f2 + 90*f3*f3 + f4*f4 + 10*f5*f5 + 0.1*f6*f6
 }

-func (Wood) Grad(grad, x []float64) {
+func (Wood) Grad(grad, x []float64) []float64 {
 	if len(x) != 4 {
 		panic("dimension of the problem must be 4")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1636,26 +1744,32 @@ func (Wood) Grad(grad, x []float64) {
 	grad[1] = 2 * (100*f1 + 10*f5 + 0.1*f6)
 	grad[2] = -2 * (180*f3*x[2] + f4)
 	grad[3] = 2 * (90*f3 + 10*f5 - 0.1*f6)
+	return grad
 }

-func (Wood) Hess(hess mat.MutableSymmetric, x []float64) {
+func (Wood) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
 	if len(x) != 4 {
 		panic("dimension of the problem must be 4")
 	}
+	if hess == nil {
+		hess = mat.NewSymDense(len(x), nil)
+	}
 	if len(x) != hess.Symmetric() {
 		panic("incorrect size of the Hessian")
 	}
+	h := hess.(*mat.SymDense)

-	hess.SetSym(0, 0, 400*(3*x[0]*x[0]-x[1])+2)
-	hess.SetSym(0, 1, -400*x[0])
-	hess.SetSym(1, 1, 220.2)
-	hess.SetSym(0, 2, 0)
-	hess.SetSym(1, 2, 0)
-	hess.SetSym(2, 2, 360*(3*x[2]*x[2]-x[3])+2)
-	hess.SetSym(0, 3, 0)
-	hess.SetSym(1, 3, 19.8)
-	hess.SetSym(2, 3, -360*x[2])
-	hess.SetSym(3, 3, 200.2)
+	h.SetSym(0, 0, 400*(3*x[0]*x[0]-x[1])+2)
+	h.SetSym(0, 1, -400*x[0])
+	h.SetSym(1, 1, 220.2)
+	h.SetSym(0, 2, 0)
+	h.SetSym(1, 2, 0)
+	h.SetSym(2, 2, 360*(3*x[2]*x[2]-x[3])+2)
+	h.SetSym(0, 3, 0)
+	h.SetSym(1, 3, 19.8)
+	h.SetSym(2, 3, -360*x[2])
+	h.SetSym(3, 3, 200.2)
+	return h
 }

 func (Wood) Minima() []Minimum {
@@ -1683,7 +1797,7 @@ func (ConcaveRight) Func(x []float64) float64 {
 	return -x[0] / (x[0]*x[0] + 2)
 }

-func (ConcaveRight) Grad(grad, x []float64) {
+func (ConcaveRight) Grad(grad, x []float64) []float64 {
 	if len(x) != 1 {
 		panic("dimension of the problem must be 1")
 	}
@@ -1692,6 +1806,7 @@ func (ConcaveRight) Grad(grad, x []float64) {
 	}
 	xSqr := x[0] * x[0]
 	grad[0] = (xSqr - 2) / (xSqr + 2) / (xSqr + 2)
+	return grad
 }

 // ConcaveLeft implements an univariate function that is concave to the left of
@@ -1709,14 +1824,18 @@ func (ConcaveLeft) Func(x []float64) float64 {
 	return math.Pow(x[0]+0.004, 4) * (x[0] - 1.996)
 }

-func (ConcaveLeft) Grad(grad, x []float64) {
+func (ConcaveLeft) Grad(grad, x []float64) []float64 {
 	if len(x) != 1 {
 		panic("dimension of the problem must be 1")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
 	grad[0] = math.Pow(x[0]+0.004, 3) * (5*x[0] - 7.98)
+	return grad
 }

 // Plassmann implements an univariate oscillatory function where the value of L
@@ -1753,10 +1872,13 @@ func (f Plassmann) Func(x []float64) float64 {
 	return r
 }

-func (f Plassmann) Grad(grad, x []float64) {
+func (f Plassmann) Grad(grad, x []float64) []float64 {
 	if len(x) != 1 {
 		panic("dimension of the problem must be 1")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1772,6 +1894,7 @@ func (f Plassmann) Grad(grad, x []float64) {
 	default: // a > 1+b
 		grad[0]++
 	}
+	return grad
 }

 // YanaiOzawaKaneko is an univariate convex function where the values of Beta1
@@ -1802,10 +1925,13 @@ func (f YanaiOzawaKaneko) Func(x []float64) float64 {
 	return g1*math.Sqrt((a-1)*(a-1)+b2*b2) + g2*math.Sqrt(a*a+b1*b1)
 }

-func (f YanaiOzawaKaneko) Grad(grad, x []float64) {
+func (f YanaiOzawaKaneko) Grad(grad, x []float64) []float64 {
 	if len(x) != 1 {
 		panic("dimension of the problem must be 1")
 	}
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
 	if len(x) != len(grad) {
 		panic("incorrect size of the gradient")
 	}
@@ -1815,4 +1941,5 @@ func (f YanaiOzawaKaneko) Grad(grad, x []float64) {
 	g1 := math.Sqrt(1+b1*b1) - b1
 	g2 := math.Sqrt(1+b2*b2) - b2
 	grad[0] = g1*(a-1)/math.Sqrt(b2*b2+(a-1)*(a-1)) + g2*a/math.Sqrt(b1*b1+a*a)
+	return grad
 }
--- a/optimize/functions/minsurf.go
+++ b/optimize/functions/minsurf.go
@@ -75,12 +75,15 @@ func (ms *MinimalSurface) Func(x []float64) (area float64) {
 }

 // Grad evaluates the area gradient of the surface represented by the vector.
-func (ms *MinimalSurface) Grad(grad, x []float64) {
+func (ms *MinimalSurface) Grad(grad, x []float64) []float64 {
 	nx, ny := ms.Dims()
 	if len(x) != (nx-2)*(ny-2) {
 		panic("functions: problem size mismatch")
 	}
-	if grad != nil && len(x) != len(grad) {
+	if grad == nil {
+		grad = make([]float64, len(x))
+	}
+	if len(x) != len(grad) {
 		panic("functions: unexpected size mismatch")
 	}

@@ -128,6 +131,7 @@ func (ms *MinimalSurface) Grad(grad, x []float64) {
 	for i := range grad {
 		grad[i] *= cellSize
 	}
+	return grad
 }

 // InitX returns a starting location for the minimization problem. Length of
--- a/optimize/functions/validate.go
+++ b/optimize/functions/validate.go
@@ -18,7 +18,7 @@ type function interface {
 }

 type gradient interface {
-	Grad(grad, x []float64)
+	Grad(grad, x []float64) []float64
 }

 // minimumer is an objective function that can also provide information about
--- a/optimize/gradientdescent.go
+++ b/optimize/gradientdescent.go
@@ -6,6 +6,11 @@ package optimize

 import "gonum.org/v1/gonum/floats"

+var (
+	_ Method      = (*GradientDescent)(nil)
+	_ localMethod = (*GradientDescent)(nil)
+)
+
 // GradientDescent implements the steepest descent optimization method that
 // performs successive steps along the direction of the negative gradient.
 type GradientDescent struct {
@@ -30,6 +35,10 @@ func (g *GradientDescent) Status() (Status, error) {
 	return g.status, g.err
 }

+func (*GradientDescent) Uses(has Available) (uses Available, err error) {
+	return has.gradient()
+}
+
 func (g *GradientDescent) Init(dim, tasks int) int {
 	g.status = NotTerminated
 	g.err = nil
@@ -75,7 +84,7 @@ func (g *GradientDescent) NextDirection(loc *Location, dir []float64) (stepSize
 	return g.StepSizer.StepSize(loc, dir)
 }

-func (*GradientDescent) Needs() struct {
+func (*GradientDescent) needs() struct {
 	Gradient bool
 	Hessian  bool
 } {
--- a/optimize/guessandcheck.go
+++ b/optimize/guessandcheck.go
@@ -10,6 +10,8 @@ import (
 	"gonum.org/v1/gonum/stat/distmv"
 )

+var _ Method = (*GuessAndCheck)(nil)
+
 // GuessAndCheck is a global optimizer that evaluates the function at random
 // locations. Not a good optimizer, but useful for comparison and debugging.
 type GuessAndCheck struct {
@@ -19,8 +21,8 @@ type GuessAndCheck struct {
 	bestX []float64
 }

-func (g *GuessAndCheck) Needs() struct{ Gradient, Hessian bool } {
-	return struct{ Gradient, Hessian bool }{false, false}
+func (*GuessAndCheck) Uses(has Available) (uses Available, err error) {
+	return has.function()
 }

 func (g *GuessAndCheck) Init(dim, tasks int) int {
--- a/optimize/interfaces.go
+++ b/optimize/interfaces.go
@@ -44,15 +44,15 @@ type localMethod interface {
 	// updates loc and returns the next operation.
 	iterateLocal(loc *Location) (Operation, error)

-	Needser
+	needser
 }

-type Needser interface {
-	// Needs specifies information about the objective function needed by the
+type needser interface {
+	// needs specifies information about the objective function needed by the
 	// optimizer beyond just the function value. The information is used
 	// internally for initialization and must match evaluation types returned
 	// by Init and Iterate during the optimization process.
-	Needs() struct {
+	needs() struct {
 		Gradient bool
 		Hessian  bool
 	}
--- a/optimize/lbfgs.go
+++ b/optimize/lbfgs.go
@@ -8,6 +8,11 @@ import (
 	"gonum.org/v1/gonum/floats"
 )

+var (
+	_ Method      = (*LBFGS)(nil)
+	_ localMethod = (*LBFGS)(nil)
+)
+
 // LBFGS implements the limited-memory BFGS method for gradient-based
 // unconstrained minimization.
 //
@@ -52,6 +57,10 @@ func (l *LBFGS) Status() (Status, error) {
 	return l.status, l.err
 }

+func (*LBFGS) Uses(has Available) (uses Available, err error) {
+	return has.gradient()
+}
+
 func (l *LBFGS) Init(dim, tasks int) int {
 	l.status = NotTerminated
 	l.err = nil
@@ -179,7 +188,7 @@ func (l *LBFGS) NextDirection(loc *Location, dir []float64) (stepSize float64) {
 	return 1
 }

-func (*LBFGS) Needs() struct {
+func (*LBFGS) needs() struct {
 	Gradient bool
 	Hessian  bool
 } {
--- a/optimize/linesearcher_test.go
+++ b/optimize/linesearcher_test.go
@@ -41,7 +41,7 @@ func TestBacktracking(t *testing.T) {

 type funcGrader interface {
 	Func([]float64) float64
-	Grad([]float64, []float64)
+	Grad([]float64, []float64) []float64
 }

 type linesearcherTest struct {
--- a/optimize/listsearch.go
+++ b/optimize/listsearch.go
@@ -10,6 +10,8 @@ import (
 	"gonum.org/v1/gonum/mat"
 )

+var _ Method = (*ListSearch)(nil)
+
 // ListSearch finds the optimum location from a specified list of possible
 // optimum locations.
 type ListSearch struct {
@@ -24,8 +26,8 @@ type ListSearch struct {
 	bestIdx int
 }

-func (*ListSearch) Needs() struct{ Gradient, Hessian bool } {
-	return struct{ Gradient, Hessian bool }{false, false}
+func (*ListSearch) Uses(has Available) (uses Available, err error) {
+	return has.function()
 }

 // InitGlobal initializes the method for optimization. The input dimension
--- a/optimize/local.go
+++ b/optimize/local.go
@@ -39,7 +39,6 @@ func (l localOptimizer) run(method localMethod, gradThresh float64, operation ch
 		l.finish(operation, result)
 		return NotTerminated, nil
 	}
-
 	op, err := method.initLocal(task.Location)
 	if err != nil {
 		l.finishMethodDone(operation, result, task)
@@ -77,13 +76,13 @@ Loop:

 // initialOperation returns the Operation needed to fill the initial location
 // based on the needs of the method and the values already supplied.
-func (localOptimizer) initialOperation(task Task, needser Needser) Operation {
+func (localOptimizer) initialOperation(task Task, n needser) Operation {
 	var newOp Operation
 	op := task.Op
 	if op&FuncEvaluation == 0 {
 		newOp |= FuncEvaluation
 	}
-	needs := needser.Needs()
+	needs := n.needs()
 	if needs.Gradient && op&GradEvaluation == 0 {
 		newOp |= GradEvaluation
 	}
@@ -95,8 +94,8 @@ func (localOptimizer) initialOperation(task Task, needser Needser) Operation {

 // initialLocation fills the initial location based on the needs of the method.
 // The task passed to initialLocation should be the first task sent in RunGlobal.
-func (l localOptimizer) initialLocation(operation chan<- Task, result <-chan Task, task Task, needser Needser) Task {
-	task.Op = l.initialOperation(task, needser)
+func (l localOptimizer) initialLocation(operation chan<- Task, result <-chan Task, task Task, needs needser) Task {
+	task.Op = l.initialOperation(task, needs)
 	operation <- task
 	return <-result
 }
--- a/optimize/local_example_test.go
+++ b/optimize/local_example_test.go
@@ -27,12 +27,12 @@ func ExampleMinimize() {
 		log.Fatal(err)
 	}
 	fmt.Printf("result.Status: %v\n", result.Status)
-	fmt.Printf("result.X: %v\n", result.X)
-	fmt.Printf("result.F: %v\n", result.F)
+	fmt.Printf("result.X: %0.4g\n", result.X)
+	fmt.Printf("result.F: %0.4g\n", result.F)
 	fmt.Printf("result.Stats.FuncEvaluations: %d\n", result.Stats.FuncEvaluations)
 	// Output:
 	// result.Status: GradientThreshold
 	// result.X: [1 1 1 1 1]
-	// result.F: 0
-	// result.Stats.FuncEvaluations: 35
+	// result.F: 4.98e-30
+	// result.Stats.FuncEvaluations: 31
 }
--- a/optimize/minimize.go
+++ b/optimize/minimize.go
@@ -33,17 +33,26 @@ type Task struct {
 	*Location
 }

+// Location represents a location in the optimization procedure.
+type Location struct {
+	X        []float64
+	F        float64
+	Gradient []float64
+	Hessian  mat.Symmetric
+}
+
 // Method is a type which can search for an optimum of an objective function.
 type Method interface {
-	Needser
-	// Init takes as input the problem dimension and number of available
-	// concurrent tasks, and returns the number of concurrent processes to be used.
-	// The returned value must be less than or equal to tasks.
-	Init(dim, tasks int) int
+	// Init initializes the method for optimization. The inputs are
+	// the problem dimension and number of available concurrent tasks.
+	//
+	// Init returns the number of concurrent processes to use, which must be
+	// less than or equal to tasks.
+	Init(dim, tasks int) (concurrent int)
 	// Run runs an optimization. The method sends Tasks on
 	// the operation channel (for performing function evaluations, major
 	// iterations, etc.). The result of the tasks will be returned on Result.
-	// See the documentation for Operation types for the possible tasks.
+	// See the documentation for Operation types for the possible operations.
 	//
 	// The caller of Run will signal the termination of the optimization
 	// (i.e. convergence from user settings) by sending a task with a PostIteration
@@ -70,6 +79,13 @@ type Method interface {
 	// composed to specify the valid fields. Methods are free to use or
 	// ignore these values.
 	//
+	// Successful execution of an Operation may require the Method to modify
+	// fields a Location. MajorIteration calls will not modify the values in
+	// the Location, but Evaluation operations will. Methods are encouraged to
+	// leave Location fields untouched to allow memory re-use. If data needs to
+	// be stored, the respective field should be set to nil -- Methods should
+	// not allocate Location memory themselves.
+	//
 	// Method may have its own specific convergence criteria, which can
 	// be communicated using a MethodDone operation. This will trigger a
 	// PostIteration to be sent on result, and the MethodDone task will not be
@@ -79,6 +95,11 @@ type Method interface {
 	// The operation and result tasks are guaranteed to have a buffer length
 	// equal to the return from Init.
 	Run(operation chan<- Task, result <-chan Task, tasks []Task)
+	// Uses checks if the Method is suited to the optimization problem. The
+	// input is the available functions in Problem to call, and the returns are
+	// the functions which may be used and an error if there is a mismatch
+	// between the Problem and the Method's capabilities.
+	Uses(has Available) (uses Available, err error)
 }

 // Minimize uses an optimizer to search for a minimum of a function. A
@@ -110,7 +131,9 @@ type Method interface {
 //
 // The final argument is the optimization method to use. If method == nil, then
 // an appropriate default is chosen based on the properties of the other arguments
-// (dimension, gradient-free or gradient-based, etc.).
+// (dimension, gradient-free or gradient-based, etc.). If method is not nil,
+// Minimize panics if the Problem is not consistent with the Method (Uses
+// returns an error).
 //
 // Minimize returns a Result struct and any error that occurred. See the
 // documentation of Result for more information.
@@ -131,15 +154,15 @@ func Minimize(p Problem, initX []float64, settings *Settings, method Method) (*R
 	}
 	stats := &Stats{}
 	dim := len(initX)
-	err := checkOptimization(p, dim, method, settings.Recorder)
+	err := checkOptimization(p, dim, settings.Recorder)
 	if err != nil {
 		return nil, err
 	}

-	optLoc := newLocation(dim, method)
+	optLoc := newLocation(dim) // This must have an allocated X field.
 	optLoc.F = math.Inf(1)

-	initOp, initLoc := getInitLocation(dim, initX, settings.InitValues, method)
+	initOp, initLoc := getInitLocation(dim, initX, settings.InitValues)

 	converger := settings.Converger
 	if converger == nil {
@@ -175,7 +198,7 @@ func Minimize(p Problem, initX []float64, settings *Settings, method Method) (*R

 func getDefaultMethod(p *Problem) Method {
 	if p.Grad != nil {
-		return &BFGS{}
+		return &LBFGS{}
 	}
 	return &NelderMead{}
 }
@@ -188,6 +211,11 @@ func minimize(prob *Problem, method Method, settings *Settings, converger Conver
 	if nTasks == 0 {
 		nTasks = 1
 	}
+	has := availFromProblem(*prob)
+	_, initErr := method.Uses(has)
+	if initErr != nil {
+		panic(fmt.Sprintf("optimize: specified method inconsistent with Problem: %v", initErr))
+	}
 	newNTasks := method.Init(dim, nTasks)
 	if newNTasks > nTasks {
 		panic("optimize: too many tasks returned by Method")
@@ -203,7 +231,7 @@ func minimize(prob *Problem, method Method, settings *Settings, converger Conver
 		tasks[0].Location = initLoc
 		tasks[0].Op = initOp
 		for i := 1; i < len(tasks); i++ {
-			tasks[i].Location = newLocation(dim, method)
+			tasks[i].Location = newLocation(dim)
 		}
 		method.Run(operations, results, tasks)
 	}()
@@ -370,44 +398,18 @@ func defaultFunctionConverge() *FunctionConverge {
 	}
 }

-// newLocation allocates a new locatian structure of the appropriate size. It
-// allocates memory based on the dimension and the values in Needs.
-func newLocation(dim int, method Needser) *Location {
-	// TODO(btracey): combine this with Local.
-	loc := &Location{
+// newLocation allocates a new locatian structure with an X field of the
+// appropriate size.
+func newLocation(dim int) *Location {
+	return &Location{
 		X: make([]float64, dim),
 	}
-	if method.Needs().Gradient {
-		loc.Gradient = make([]float64, dim)
-	}
-	if method.Needs().Hessian {
-		loc.Hessian = mat.NewSymDense(dim, nil)
-	}
-	return loc
-}
-
-func copyLocation(dst, src *Location) {
-	dst.X = resize(dst.X, len(src.X))
-	copy(dst.X, src.X)
-
-	dst.F = src.F
-
-	dst.Gradient = resize(dst.Gradient, len(src.Gradient))
-	copy(dst.Gradient, src.Gradient)
-
-	if src.Hessian != nil {
-		if dst.Hessian == nil || dst.Hessian.Symmetric() != len(src.X) {
-			dst.Hessian = mat.NewSymDense(len(src.X), nil)
-		}
-		dst.Hessian.CopySym(src.Hessian)
-	}
 }

 // getInitLocation checks the validity of initLocation and initOperation and
 // returns the initial values as a *Location.
-func getInitLocation(dim int, initX []float64, initValues *Location, method Needser) (Operation, *Location) {
-	needs := method.Needs()
-	loc := newLocation(dim, method)
+func getInitLocation(dim int, initX []float64, initValues *Location) (Operation, *Location) {
+	loc := newLocation(dim)
 	if initX == nil {
 		if initValues != nil {
 			panic("optimize: initValues is non-nil but no initial location specified")
@@ -428,33 +430,26 @@ func getInitLocation(dim int, initX []float64, initValues *Location, method Need
 		if len(initValues.Gradient) != dim {
 			panic("optimize: initial gradient does not match problem dimension")
 		}
-		if needs.Gradient {
-			copy(loc.Gradient, initValues.Gradient)
-			op |= GradEvaluation
-		}
+		loc.Gradient = initValues.Gradient
+		op |= GradEvaluation
 	}
 	if initValues.Hessian != nil {
 		if initValues.Hessian.Symmetric() != dim {
 			panic("optimize: initial Hessian does not match problem dimension")
 		}
-		if needs.Hessian {
-			loc.Hessian.CopySym(initValues.Hessian)
-			op |= HessEvaluation
-		}
+		loc.Hessian = initValues.Hessian
+		op |= HessEvaluation
 	}
 	return op, loc
 }

-func checkOptimization(p Problem, dim int, method Needser, recorder Recorder) error {
+func checkOptimization(p Problem, dim int, recorder Recorder) error {
 	if p.Func == nil {
 		panic(badProblem)
 	}
 	if dim <= 0 {
 		panic("optimize: impossible problem dimension")
 	}
-	if err := p.satisfies(method); err != nil {
-		return err
-	}
 	if p.Status != nil {
 		_, err := p.Status()
 		if err != nil {
@@ -482,10 +477,10 @@ func evaluate(p *Problem, loc *Location, op Operation, x []float64) {
 		loc.F = p.Func(x)
 	}
 	if op&GradEvaluation != 0 {
-		p.Grad(loc.Gradient, x)
+		loc.Gradient = p.Grad(loc.Gradient, x)
 	}
 	if op&HessEvaluation != 0 {
-		p.Hess(loc.Hessian, x)
+		loc.Hessian = p.Hess(loc.Hessian, x)
 	}
 }

@@ -559,7 +554,16 @@ func checkIterationLimits(loc *Location, stats *Stats, settings *Settings) Statu
 // It increments the iteration count, updates the optimal location, and checks
 // the necessary convergence criteria.
 func performMajorIteration(optLoc, loc *Location, stats *Stats, converger Converger, startTime time.Time, settings *Settings) Status {
-	copyLocation(optLoc, loc)
+	optLoc.F = loc.F
+	copy(optLoc.X, loc.X)
+	if loc.Gradient == nil {
+		optLoc.Gradient = nil
+	} else {
+		if optLoc.Gradient == nil {
+			optLoc.Gradient = make([]float64, len(loc.Gradient))
+		}
+		copy(optLoc.Gradient, loc.Gradient)
+	}
 	stats.MajorIterations++
 	stats.Runtime = time.Since(startTime)
 	status := checkLocationConvergence(optLoc, settings, converger)
--- a/optimize/neldermead.go
+++ b/optimize/neldermead.go
@@ -42,6 +42,8 @@ func (n nmVertexSorter) Swap(i, j int) {
 	n.vertices[i], n.vertices[j] = n.vertices[j], n.vertices[i]
 }

+var _ Method = (*NelderMead)(nil)
+
 // NelderMead is an implementation of the Nelder-Mead simplex algorithm for
 // gradient-free nonlinear optimization (not to be confused with Danzig's
 // simplex algorithm for linear programming). The implementation follows the
@@ -90,6 +92,10 @@ func (n *NelderMead) Status() (Status, error) {
 	return n.status, n.err
 }

+func (*NelderMead) Uses(has Available) (uses Available, err error) {
+	return has.function()
+}
+
 func (n *NelderMead) Init(dim, tasks int) int {
 	n.status = NotTerminated
 	n.err = nil
@@ -332,7 +338,7 @@ func (n *NelderMead) replaceWorst(x []float64, f float64) {
 	floats.AddScaled(n.centroid, 1/float64(dim), x)
 }

-func (*NelderMead) Needs() struct {
+func (*NelderMead) needs() struct {
 	Gradient bool
 	Hessian  bool
 } {
--- a/optimize/newton.go
+++ b/optimize/newton.go
@@ -12,6 +12,11 @@ import (

 const maxNewtonModifications = 20

+var (
+	_ Method      = (*Newton)(nil)
+	_ localMethod = (*Newton)(nil)
+)
+
 // Newton implements a modified Newton's method for Hessian-based unconstrained
 // minimization. It applies regularization when the Hessian is not positive
 // definite, and it can converge to a local minimum from any starting point.
@@ -64,6 +69,10 @@ func (n *Newton) Status() (Status, error) {
 	return n.status, n.err
 }

+func (*Newton) Uses(has Available) (uses Available, err error) {
+	return has.hessian()
+}
+
 func (n *Newton) Init(dim, tasks int) int {
 	n.status = NotTerminated
 	n.err = nil
@@ -91,7 +100,6 @@ func (n *Newton) initLocal(loc *Location) (Operation, error) {
 	}
 	n.ls.Linesearcher = n.Linesearcher
 	n.ls.NextDirectioner = n
-
 	return n.ls.Init(loc)
 }

@@ -159,7 +167,7 @@ func (n *Newton) NextDirection(loc *Location, dir []float64) (stepSize float64)
 	return 1
 }

-func (n *Newton) Needs() struct {
+func (n *Newton) needs() struct {
 	Gradient bool
 	Hessian  bool
 } {
--- a/optimize/types.go
+++ b/optimize/types.go
@@ -5,7 +5,6 @@
 package optimize

 import (
-	"errors"
 	"fmt"
 	"time"

@@ -85,17 +84,9 @@ var operationNames = map[Operation]string{
 	signalDone:     "signalDone",
 }

-// Location represents a location in the optimization procedure.
-type Location struct {
-	X        []float64
-	F        float64
-	Gradient []float64
-	Hessian  *mat.SymDense
-}
-
 // Result represents the answer of an optimization run. It contains the optimum
-// location as well as the Status at convergence and Statistics taken during the
-// run.
+// function value, X location, and gradient as well as the Status at convergence
+// and Statistics taken during the run.
 type Result struct {
 	Location
 	Stats
@@ -132,13 +123,17 @@ type Problem struct {
 	// must not modify x.
 	Func func(x []float64) float64

-	// Grad evaluates the gradient at x and stores the result in-place in grad.
-	// Grad must not modify x.
-	Grad func(grad []float64, x []float64)
+	// Grad evaluates the gradient at x and returns the result. Grad may use
+	// (and return) the provided slice if it is non-nil, or must allocate a new
+	// slice otherwise. Grad must not modify x.
+	Grad func(grad []float64, x []float64) []float64

 	// Hess evaluates the Hessian at x and stores the result in-place in hess.
-	// Hess must not modify x.
-	Hess func(hess mat.MutableSymmetric, x []float64)
+	// Hess must not modify x. Hess may use (and return) the provided Symmetric
+	// if it is non-nil, or must allocate a new Symmetric otherwise. Minimize
+	// will 'give back' the returned Symmetric where possible, allowing Hess
+	// to use a type assertion on the provided matrix.
+	Hess func(hess mat.Symmetric, x []float64) mat.Symmetric

 	// Status reports the status of the objective function being optimized and any
 	// error. This can be used to terminate early, for example when the function is
@@ -147,16 +142,47 @@ type Problem struct {
 	Status func() (Status, error)
 }

-// TODO(btracey): Think about making this an exported function when the
-// constraint interface is designed.
-func (p Problem) satisfies(method Needser) error {
-	if method.Needs().Gradient && p.Grad == nil {
-		return errors.New("optimize: problem does not provide needed Grad function")
+// Available describes the functions available to call in Problem.
+type Available struct {
+	Grad bool
+	Hess bool
+}
+
+func availFromProblem(prob Problem) Available {
+	return Available{Grad: prob.Grad != nil, Hess: prob.Hess != nil}
+}
+
+// function tests if the Problem described by the receiver is suitable for an
+// unconstrained Method that only calls the function, and returns the result.
+func (has Available) function() (uses Available, err error) {
+	// TODO(btracey): This needs to be modified when optimize supports
+	// constrained optimization.
+	return Available{}, nil
+}
+
+// gradient tests if the Problem described by the receiver is suitable for an
+// unconstrained gradient-based Method, and returns the result.
+func (has Available) gradient() (uses Available, err error) {
+	// TODO(btracey): This needs to be modified when optimize supports
+	// constrained optimization.
+	if !has.Grad {
+		return Available{}, ErrMissingGrad
 	}
-	if method.Needs().Hessian && p.Hess == nil {
-		return errors.New("optimize: problem does not provide needed Hess function")
+	return Available{Grad: true}, nil
+}
+
+// hessian tests if the Problem described by the receiver is suitable for an
+// unconstrained Hessian-based Method, and returns the result.
+func (has Available) hessian() (uses Available, err error) {
+	// TODO(btracey): This needs to be modified when optimize supports
+	// constrained optimization.
+	if !has.Grad {
+		return Available{}, ErrMissingGrad
 	}
-	return nil
+	if !has.Hess {
+		return Available{}, ErrMissingHess
+	}
+	return Available{Grad: true, Hess: true}, nil
 }

 // Settings represents settings of the optimization run. It contains initial
@@ -167,7 +193,8 @@ type Settings struct {
 	// InitValues specifies properties (function value, gradient, etc.) known
 	// at the initial location passed to Minimize. If InitValues is non-nil, then
 	// the function value F must be provided, the location X must not be specified
-	// and other fields may be specified.
+	// and other fields may be specified. The values in Location may be modified
+	// during the call to Minimize.
 	InitValues *Location

 	// GradientThreshold stops optimization with GradientThreshold status if the
--- a/optimize/unconstrained_test.go
+++ b/optimize/unconstrained_test.go
@@ -1162,9 +1162,21 @@ func testLocal(t *testing.T, tests []unconstrainedTest, method Method) {

 		settings := &Settings{}
 		settings.Converger = defaultFunctionConverge()
-		if method != nil && method.Needs().Gradient {
+		var uses Available
+		if method != nil {
+			var err error
+			has := availFromProblem(test.p)
+			uses, err = method.Uses(has)
+			if err != nil {
+				t.Errorf("problem and method mismatch: %v", err)
+				continue
+			}
+		}
+		if method != nil {
 			// Turn off function convergence checks for gradient-based methods.
-			settings.Converger = NeverTerminate{}
+			if uses.Grad {
+				settings.Converger = NeverTerminate{}
+			}
 		} else {
 			if test.fIter == 0 {
 				test.fIter = 20
@@ -1222,7 +1234,7 @@ func testLocal(t *testing.T, tests []unconstrainedTest, method Method) {
 			continue
 		}

-		if !method.Needs().Gradient && !method.Needs().Hessian {
+		if !uses.Grad && !uses.Hess {
 			// Gradient-free tests can correctly terminate only with
 			// FunctionConvergence status.
 			if result.Status != FunctionConvergence {
@@ -1234,11 +1246,11 @@ func testLocal(t *testing.T, tests []unconstrainedTest, method Method) {
 		// evaluate them.
 		settings.InitValues = &Location{}
 		settings.InitValues.F = test.p.Func(test.x)
-		if method.Needs().Gradient {
+		if uses.Grad {
 			settings.InitValues.Gradient = resize(settings.InitValues.Gradient, len(test.x))
 			test.p.Grad(settings.InitValues.Gradient, test.x)
 		}
-		if method.Needs().Hessian {
+		if uses.Hess {
 			settings.InitValues.Hessian = mat.NewSymDense(len(test.x), nil)
 			test.p.Hess(settings.InitValues.Hessian, test.x)
 		}
@@ -1266,13 +1278,13 @@ func testLocal(t *testing.T, tests []unconstrainedTest, method Method) {
 			t.Errorf("Providing initial data does not reduce the number of Func calls for:\n%v", test)
 			continue
 		}
-		if method.Needs().Gradient {
+		if uses.Grad {
 			if result.GradEvaluations != result2.GradEvaluations+1 {
 				t.Errorf("Providing initial data does not reduce the number of Grad calls for:\n%v", test)
 				continue
 			}
 		}
-		if method.Needs().Hessian {
+		if uses.Hess {
 			if result.HessEvaluations != result2.HessEvaluations+1 {
 				t.Errorf("Providing initial data does not reduce the number of Hess calls for:\n%v", test)
 				continue