diff/fd: implement Hessian finite difference, and code cleanups. (#109)

* diff/fd: implement Hessian finite difference, and code cleanups.

This commit primarily adds the Hessian function for finding a finite difference approximation to the Hessian. At the same time, it combines duplicated functionality across the difference routines so that the preludes to all the difference routines look similar

diff/fd/derivative.go

@@ -0,0 +1,71 @@
+// Copyright ©2017 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 fd
+
+import (
+	"math"
+	"runtime"
+	"sync"
+)
+
+// Derivative estimates the derivative of the function f at the given location.
+// The finite difference formula, the step size, and other options are
+// specified by settings. If settings is nil, the first derivative will be
+// estimated using the Forward formula and a default step size.
+func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
+	// Default settings.
+	formula := Forward
+	step := formula.Step
+	var originValue float64
+	var originKnown, concurrent bool
+
+	// Use user settings if provided.
+	if settings != nil {
+		if !settings.Formula.isZero() {
+			formula = settings.Formula
+			step = formula.Step
+			checkFormula(formula)
+		}
+		if settings.Step != 0 {
+			step = settings.Step
+		}
+		originKnown = settings.OriginKnown
+		originValue = settings.OriginValue
+		concurrent = settings.Concurrent
+	}
+
+	var deriv float64
+	if !concurrent || runtime.GOMAXPROCS(0) == 1 {
+		for _, pt := range formula.Stencil {
+			if originKnown && pt.Loc == 0 {
+				deriv += pt.Coeff * originValue
+				continue
+			}
+			deriv += pt.Coeff * f(x+step*pt.Loc)
+		}
+		return deriv / math.Pow(step, float64(formula.Derivative))
+	}
+
+	wg := &sync.WaitGroup{}
+	mux := &sync.Mutex{}
+	for _, pt := range formula.Stencil {
+		if originKnown && pt.Loc == 0 {
+			mux.Lock()
+			deriv += pt.Coeff * originValue
+			mux.Unlock()
+			continue
+		}
+		wg.Add(1)
+		go func(pt Point) {
+			defer wg.Done()
+			fofx := f(x + step*pt.Loc)
+			mux.Lock()
+			defer mux.Unlock()
+			deriv += pt.Coeff * fofx
+		}(pt)
+	}
+	wg.Wait()
+	return deriv / math.Pow(step, float64(formula.Derivative))
+}

diff/fd/diff.go

@@ -8,9 +8,6 @@ package fd // import "gonum.org/v1/gonum/diff/fd"
 import (
 	"math"
 	"runtime"
-	"sync"
-
-	"gonum.org/v1/gonum/floats"
 )
 
 // A Point is a stencil location in a finite difference formula.
@@ -22,6 +19,7 @@ type Point struct {
 // Formula represents a finite difference formula on a regularly spaced grid
 // that approximates the derivative of order k of a function f at x as
 //  d^k f(x) ≈ (1 / Step^k) * \sum_i Coeff_i * f(x + Step * Loc_i).
+// Step must be positive, or the finite difference formula will panic.
 type Formula struct {
 	// Stencil is the set of sampling Points which are used to estimate the
 	// derivative. The locations will be scaled by Step and are relative to x.
@@ -50,199 +48,6 @@ type Settings struct {
 	Concurrent bool // Should the function calls be executed concurrently.
 }
 
-// Derivative estimates the derivative of the function f at the given location.
-// The finite difference formula, the step size, and other options are
-// specified by settings. If settings is nil, the first derivative will be
-// estimated using the Forward formula and a default step size.
-func Derivative(f func(float64) float64, x float64, settings *Settings) float64 {
-	if settings == nil {
-		settings = &Settings{}
-	}
-	formula := settings.Formula
-	if formula.isZero() {
-		formula = Forward
-	}
-	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
-		panic("fd: bad formula")
-	}
-	step := settings.Step
-	if step == 0 {
-		step = formula.Step
-	}
-
-	var deriv float64
-	if !settings.Concurrent || runtime.GOMAXPROCS(0) == 1 {
-		for _, pt := range formula.Stencil {
-			if settings.OriginKnown && pt.Loc == 0 {
-				deriv += pt.Coeff * settings.OriginValue
-				continue
-			}
-			deriv += pt.Coeff * f(x+step*pt.Loc)
-		}
-		return deriv / math.Pow(step, float64(formula.Derivative))
-	}
-
-	wg := &sync.WaitGroup{}
-	mux := &sync.Mutex{}
-	for _, pt := range formula.Stencil {
-		if settings.OriginKnown && pt.Loc == 0 {
-			mux.Lock()
-			deriv += pt.Coeff * settings.OriginValue
-			mux.Unlock()
-			continue
-		}
-		wg.Add(1)
-		go func(pt Point) {
-			defer wg.Done()
-			fofx := f(x + step*pt.Loc)
-			mux.Lock()
-			defer mux.Unlock()
-			deriv += pt.Coeff * fofx
-		}(pt)
-	}
-	wg.Wait()
-	return deriv / math.Pow(step, float64(formula.Derivative))
-}
-
-// Gradient estimates the gradient of the multivariate function f at the
-// location x. If dst is not nil, the result will be stored in-place into dst
-// and returned, otherwise a new slice will be allocated first. Finite
-// difference kernel and other options are specified by settings. If settings is
-// nil, the gradient will be estimated using the Forward formula and a default
-// step size.
-//
-// Gradient panics if the length of dst and x is not equal, or if the derivative
-// order of the formula is not 1.
-func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *Settings) []float64 {
-	if dst == nil {
-		dst = make([]float64, len(x))
-	}
-	if len(dst) != len(x) {
-		panic("fd: slice length mismatch")
-	}
-	if settings == nil {
-		settings = &Settings{}
-	}
-
-	formula := settings.Formula
-	if formula.isZero() {
-		formula = Forward
-	}
-	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
-		panic("fd: bad formula")
-	}
-	if formula.Derivative != 1 {
-		panic("fd: invalid derivative order")
-	}
-
-	step := settings.Step
-	if step == 0 {
-		step = formula.Step
-	}
-
-	expect := len(formula.Stencil) * len(x)
-	nWorkers := 1
-	if settings.Concurrent {
-		nWorkers = runtime.GOMAXPROCS(0)
-		if nWorkers > expect {
-			nWorkers = expect
-		}
-	}
-
-	var hasOrigin bool
-	for _, pt := range formula.Stencil {
-		if pt.Loc == 0 {
-			hasOrigin = true
-			break
-		}
-	}
-	xcopy := make([]float64, len(x)) // So that x is not modified during the call.
-	originValue := settings.OriginValue
-	if hasOrigin && !settings.OriginKnown {
-		copy(xcopy, x)
-		originValue = f(xcopy)
-	}
-
-	if nWorkers == 1 {
-		for i := range xcopy {
-			var deriv float64
-			for _, pt := range formula.Stencil {
-				if pt.Loc == 0 {
-					deriv += pt.Coeff * originValue
-					continue
-				}
-				copy(xcopy, x)
-				xcopy[i] += pt.Loc * step
-				deriv += pt.Coeff * f(xcopy)
-			}
-			dst[i] = deriv / step
-		}
-		return dst
-	}
-
-	sendChan := make(chan fdrun, expect)
-	ansChan := make(chan fdrun, expect)
-	quit := make(chan struct{})
-	defer close(quit)
-
-	// Launch workers. Workers receive an index and a step, and compute the answer.
-	for i := 0; i < nWorkers; i++ {
-		go func(sendChan <-chan fdrun, ansChan chan<- fdrun, quit <-chan struct{}) {
-			xcopy := make([]float64, len(x))
-			for {
-				select {
-				case <-quit:
-					return
-				case run := <-sendChan:
-					copy(xcopy, x)
-					xcopy[run.idx] += run.pt.Loc * step
-					run.result = f(xcopy)
-					ansChan <- run
-				}
-			}
-		}(sendChan, ansChan, quit)
-	}
-
-	// Launch the distributor. Distributor sends the cases to be computed.
-	go func(sendChan chan<- fdrun, ansChan chan<- fdrun) {
-		for i := range x {
-			for _, pt := range formula.Stencil {
-				if pt.Loc == 0 {
-					// Answer already known. Send the answer on the answer channel.
-					ansChan <- fdrun{
-						idx:    i,
-						pt:     pt,
-						result: originValue,
-					}
-					continue
-				}
-				// Answer not known, send the answer to be computed.
-				sendChan <- fdrun{
-					idx: i,
-					pt:  pt,
-				}
-			}
-		}
-	}(sendChan, ansChan)
-
-	for i := range dst {
-		dst[i] = 0
-	}
-	// Read in all of the results.
-	for i := 0; i < expect; i++ {
-		run := <-ansChan
-		dst[run.idx] += run.pt.Coeff * run.result
-	}
-	floats.Scale(1/step, dst)
-	return dst
-}
-
-type fdrun struct {
-	idx    int
-	pt     Point
-	result float64
-}
-
 // Forward represents a first-order accurate forward approximation
 // to the first derivative.
 var Forward = Formula{
@@ -251,6 +56,14 @@ var Forward = Formula{
 	Step:       2e-8,
 }
 
+// Forward2nd represents a first-order accurate forward approximation
+// to the second derivative.
+var Forward2nd = Formula{
+	Stencil:    []Point{{Loc: 0, Coeff: 1}, {Loc: 1, Coeff: -2}, {Loc: 2, Coeff: 1}},
+	Derivative: 2,
+	Step:       1e-4,
+}
+
 // Backward represents a first-order accurate backward approximation
 // to the first derivative.
 var Backward = Formula{
@@ -259,6 +72,14 @@ var Backward = Formula{
 	Step:       2e-8,
 }
 
+// Backward2nd represents a first-order accurate forward approximation
+// to the second derivative.
+var Backward2nd = Formula{
+	Stencil:    []Point{{Loc: 0, Coeff: 1}, {Loc: -1, Coeff: -2}, {Loc: -2, Coeff: 1}},
+	Derivative: 2,
+	Step:       1e-4,
+}
+
 // Central represents a second-order accurate centered approximation
 // to the first derivative.
 var Central = Formula{
@@ -274,3 +95,55 @@ var Central2nd = Formula{
 	Derivative: 2,
 	Step:       1e-4,
 }
+
+var negativeStep = "fd: negative step"
+
+// checkFormula checks if the formula is valid, and panics otherwise.
+func checkFormula(formula Formula) {
+	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step <= 0 {
+		panic("fd: bad formula")
+	}
+}
+
+// computeWorkers returns the desired number of workers given the concurrency
+// level and number of evaluations.
+func computeWorkers(concurrent bool, evals int) int {
+	if !concurrent {
+		return 1
+	}
+	nWorkers := runtime.GOMAXPROCS(0)
+	if nWorkers > evals {
+		nWorkers = evals
+	}
+	return nWorkers
+}
+
+// usesOrigin returns whether the stencil uses the origin, which is true iff
+// one of the locations in the stencil equals 0.
+func usesOrigin(stencil []Point) bool {
+	for _, pt := range stencil {
+		if pt.Loc == 0 {
+			return true
+		}
+	}
+	return false
+}
+
+// getOrigin returns the value at the origin. It returns originValue if originKnown
+// is true. It returns the value returned by f if stencil contains a point with
+// zero location, and NaN otherwise.
+func getOrigin(originKnown bool, originValue float64, f func() float64, stencil []Point) float64 {
+	if originKnown {
+		return originValue
+	}
+	for _, pt := range stencil {
+		if pt.Loc == 0 {
+			return f()
+		}
+	}
+	return math.NaN()
+}
+
+const (
+	badDerivOrder = "fd: invalid derivative order"
+)

diff/fd/gradient.go

@@ -0,0 +1,144 @@
+// Copyright ©2017 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 fd
+
+import "gonum.org/v1/gonum/floats"
+
+// Gradient estimates the gradient of the multivariate function f at the
+// location x. If dst is not nil, the result will be stored in-place into dst
+// and returned, otherwise a new slice will be allocated first. Finite
+// difference formula and other options are specified by settings. If settings is
+// nil, the gradient will be estimated using the Forward formula and a default
+// step size.
+//
+// Gradient panics if the length of dst and x is not equal, or if the derivative
+// order of the formula is not 1.
+func Gradient(dst []float64, f func([]float64) float64, x []float64, settings *Settings) []float64 {
+	if dst == nil {
+		dst = make([]float64, len(x))
+	}
+	if len(dst) != len(x) {
+		panic("fd: slice length mismatch")
+	}
+
+	// Default settings.
+	formula := Forward
+	step := formula.Step
+	var originValue float64
+	var originKnown, concurrent bool
+
+	// Use user settings if provided.
+	if settings != nil {
+		if !settings.Formula.isZero() {
+			formula = settings.Formula
+			step = formula.Step
+			checkFormula(formula)
+			if formula.Derivative != 1 {
+				panic(badDerivOrder)
+			}
+		}
+		if settings.Step != 0 {
+			step = settings.Step
+		}
+		originKnown = settings.OriginKnown
+		originValue = settings.OriginValue
+		concurrent = settings.Concurrent
+	}
+
+	evals := len(formula.Stencil) * len(x)
+	nWorkers := computeWorkers(concurrent, evals)
+
+	hasOrigin := usesOrigin(formula.Stencil)
+	xcopy := make([]float64, len(x)) // So that x is not modified during the call.
+	if hasOrigin && !originKnown {
+		copy(xcopy, x)
+		originValue = f(xcopy)
+	}
+
+	if nWorkers == 1 {
+		for i := range xcopy {
+			var deriv float64
+			for _, pt := range formula.Stencil {
+				if pt.Loc == 0 {
+					deriv += pt.Coeff * originValue
+					continue
+				}
+				// Copying the code anew has two benefits. First, it
+				// avoids floating point issues where adding and then
+				// subtracting the step don't return to the exact same
+				// location. Secondly, it protects against the function
+				// modifying the input data.
+				copy(xcopy, x)
+				xcopy[i] += pt.Loc * step
+				deriv += pt.Coeff * f(xcopy)
+			}
+			dst[i] = deriv / step
+		}
+		return dst
+	}
+
+	sendChan := make(chan fdrun, evals)
+	ansChan := make(chan fdrun, evals)
+	quit := make(chan struct{})
+	defer close(quit)
+
+	// Launch workers. Workers receive an index and a step, and compute the answer.
+	for i := 0; i < nWorkers; i++ {
+		go func(sendChan <-chan fdrun, ansChan chan<- fdrun, quit <-chan struct{}) {
+			xcopy := make([]float64, len(x))
+			for {
+				select {
+				case <-quit:
+					return
+				case run := <-sendChan:
+					// See above comment on the copy.
+					copy(xcopy, x)
+					xcopy[run.idx] += run.pt.Loc * step
+					run.result = f(xcopy)
+					ansChan <- run
+				}
+			}
+		}(sendChan, ansChan, quit)
+	}
+
+	// Launch the distributor. Distributor sends the cases to be computed.
+	go func(sendChan chan<- fdrun, ansChan chan<- fdrun) {
+		for i := range x {
+			for _, pt := range formula.Stencil {
+				if pt.Loc == 0 {
+					// Answer already known. Send the answer on the answer channel.
+					ansChan <- fdrun{
+						idx:    i,
+						pt:     pt,
+						result: originValue,
+					}
+					continue
+				}
+				// Answer not known, send the answer to be computed.
+				sendChan <- fdrun{
+					idx: i,
+					pt:  pt,
+				}
+			}
+		}
+	}(sendChan, ansChan)
+
+	for i := range dst {
+		dst[i] = 0
+	}
+	// Read in all of the results.
+	for i := 0; i < evals; i++ {
+		run := <-ansChan
+		dst[run.idx] += run.pt.Coeff * run.result
+	}
+	floats.Scale(1/step, dst)
+	return dst
+}
+
+type fdrun struct {
+	idx    int
+	pt     Point
+	result float64
+}

diff/fd/hessian.go

@@ -0,0 +1,192 @@
+// Copyright ©2017 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 fd
+
+import (
+	"math"
+	"sync"
+
+	"gonum.org/v1/gonum/mat"
+)
+
+// Hessian approximates the Hessian matrix of the multivariate function f
+// at the location x. That is
+//  H_{i,j} = ∂^2 f(x)/∂x_i ∂x_j
+// If dst is not nil, the resulting H will be stored in-place into dst
+// and returned, otherwise a new matrix will be allocated first. Finite difference
+// formula and other options are specified by settings. If settings is nil,
+// the Hessian will be estimated using the Forward formula and a default step size.
+//
+// Hessian panics if the size of dst and x is not equal, or if the derivative
+// order of the formula is not 1.
+func Hessian(dst *mat.SymDense, f func(x []float64) float64, x []float64, settings *Settings) *mat.SymDense {
+	n := len(x)
+	if dst == nil {
+		dst = mat.NewSymDense(n, nil)
+	} else {
+		if n2 := dst.Symmetric(); n2 != n {
+			panic("hessian: dst size mismatch")
+		}
+		for i := 0; i < n; i++ {
+			for j := i; j < n; j++ {
+				dst.SetSym(i, j, 0)
+			}
+		}
+	}
+
+	// Default settings.
+	formula := Forward
+	step := math.Sqrt(formula.Step) // Use the sqrt because taking derivatives of derivatives.
+	var originValue float64
+	var originKnown, concurrent bool
+
+	// Use user settings if provided.
+	if settings != nil {
+		if !settings.Formula.isZero() {
+			formula = settings.Formula
+			step = math.Sqrt(formula.Step)
+			checkFormula(formula)
+			if formula.Derivative != 1 {
+				panic(badDerivOrder)
+			}
+		}
+		if settings.Step != 0 {
+			if settings.Step < 0 {
+				panic(negativeStep)
+			}
+			step = settings.Step
+		}
+		originKnown = settings.OriginKnown
+		originValue = settings.OriginValue
+		concurrent = settings.Concurrent
+	}
+
+	evals := n * (n + 1) / 2 * len(formula.Stencil) * len(formula.Stencil)
+	for _, pt := range formula.Stencil {
+		if pt.Loc == 0 {
+			evals -= n * (n + 1) / 2
+			break
+		}
+	}
+
+	nWorkers := computeWorkers(concurrent, evals)
+	if nWorkers == 1 {
+		hessianSerial(dst, f, x, formula.Stencil, step, originKnown, originValue)
+		return dst
+	}
+	hessianConcurrent(dst, nWorkers, evals, f, x, formula.Stencil, step, originKnown, originValue)
+	return dst
+}
+
+func hessianSerial(dst *mat.SymDense, f func(x []float64) float64, x []float64, stencil []Point, step float64, originKnown bool, originValue float64) {
+	n := len(x)
+	xCopy := make([]float64, n)
+	fo := func() float64 {
+		copy(xCopy, x)
+		return f(x)
+	}
+	is2 := 1 / (step * step)
+	origin := getOrigin(originKnown, originValue, fo, stencil)
+	for i := 0; i < n; i++ {
+		for j := i; j < n; j++ {
+			var hess float64
+			for _, pti := range stencil {
+				for _, ptj := range stencil {
+					var v float64
+					if pti.Loc == 0 && ptj.Loc == 0 {
+						v = origin
+					} else {
+						// Copying the code anew has two benefits. First, it
+						// avoids floating point issues where adding and then
+						// subtracting the step don't return to the exact same
+						// location. Secondly, it protects against the function
+						// modifying the input data.
+						copy(xCopy, x)
+						xCopy[i] += pti.Loc * step
+						xCopy[j] += ptj.Loc * step
+						v = f(xCopy)
+					}
+					hess += v * pti.Coeff * ptj.Coeff * is2
+				}
+			}
+			dst.SetSym(i, j, hess)
+		}
+	}
+}
+
+func hessianConcurrent(dst *mat.SymDense, nWorkers, evals int, f func(x []float64) float64, x []float64, stencil []Point, step float64, originKnown bool, originValue float64) {
+	n := dst.Symmetric()
+	type run struct {
+		i, j       int
+		iIdx, jIdx int
+		result     float64
+	}
+
+	send := make(chan run, evals)
+	ans := make(chan run)
+
+	var originWG sync.WaitGroup
+	hasOrigin := usesOrigin(stencil)
+	if hasOrigin {
+		originWG.Add(1)
+		// Launch worker to compute the origin.
+		go func() {
+			defer originWG.Done()
+			xCopy := make([]float64, len(x))
+			copy(xCopy, x)
+			originValue = f(xCopy)
+		}()
+	}
+
+	var workerWG sync.WaitGroup
+	// Launch workers.
+	for i := 0; i < nWorkers; i++ {
+		workerWG.Add(1)
+		go func(send <-chan run, ans chan<- run) {
+			defer workerWG.Done()
+			xCopy := make([]float64, len(x))
+			for r := range send {
+				if stencil[r.iIdx].Loc == 0 && stencil[r.jIdx].Loc == 0 {
+					originWG.Wait()
+					r.result = originValue
+				} else {
+					// See hessianSerial for comment on the copy.
+					copy(xCopy, x)
+					xCopy[r.i] += stencil[r.iIdx].Loc * step
+					xCopy[r.j] += stencil[r.jIdx].Loc * step
+					r.result = f(xCopy)
+				}
+				ans <- r
+			}
+		}(send, ans)
+	}
+
+	// Launch the distributor, which sends all of runs.
+	go func(send chan<- run) {
+		for i := 0; i < n; i++ {
+			for j := i; j < n; j++ {
+				for iIdx := range stencil {
+					for jIdx := range stencil {
+						send <- run{
+							i: i, j: j, iIdx: iIdx, jIdx: jIdx,
+						}
+					}
+				}
+			}
+		}
+		close(send)
+		// Wait for all the workers to quit, then close the ans channel.
+		workerWG.Wait()
+		close(ans)
+	}(send)
+
+	is2 := 1 / (step * step)
+	// Read in the results.
+	for r := range ans {
+		v := r.result * stencil[r.iIdx].Coeff * stencil[r.jIdx].Coeff * is2
+		v += dst.At(r.i, r.j)
+		dst.SetSym(r.i, r.j, v)
+	}
+}

diff/fd/hessian_test.go

@@ -0,0 +1,92 @@
+// Copyright ©2017 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 fd
+
+import (
+	"testing"
+
+	"gonum.org/v1/gonum/mat"
+)
+
+type HessianTester interface {
+	Func(x []float64) float64
+	Grad(grad, x []float64)
+	Hess(dst mat.MutableSymmetric, x []float64)
+}
+
+func TestHessian(t *testing.T) {
+	for cas, test := range []struct {
+		h        HessianTester
+		x        []float64
+		settings *Settings
+		tol      float64
+	}{
+		{
+			h:   Watson{},
+			x:   []float64{0.2, 0.3, 0.1, 0.4},
+			tol: 1e-3,
+		},
+		{
+			h:   Watson{},
+			x:   []float64{2, 3, 1, 4},
+			tol: 1e-3,
+			settings: &Settings{
+				Step:    1e-5,
+				Formula: Central,
+			},
+		},
+		{
+			h:   Watson{},
+			x:   []float64{2, 3, 1},
+			tol: 1e-3,
+			settings: &Settings{
+				OriginKnown: true,
+				OriginValue: 7606.529501201192,
+			},
+		},
+		{
+			h:   ConstFunc(5),
+			x:   []float64{1, 9},
+			tol: 1e-16,
+		},
+		{
+			h:   LinearFunc{w: []float64{10, 6, -1}, c: 5},
+			x:   []float64{3, 1, 8},
+			tol: 1e-6,
+		},
+		{
+			h: QuadFunc{
+				a: mat.NewSymDense(3, []float64{
+					10, 2, 1,
+					2, 5, -3,
+					1, -3, 6,
+				}),
+				b: mat.NewVector(3, []float64{3, -2, -1}),
+				c: 5,
+			},
+			x:   []float64{-1.6, -3, 2},
+			tol: 1e-6,
+		},
+	} {
+		n := len(test.x)
+		got := Hessian(nil, test.h.Func, test.x, test.settings)
+		want := mat.NewSymDense(n, nil)
+		test.h.Hess(want, test.x)
+		if !mat.EqualApprox(got, want, test.tol) {
+			t.Errorf("Cas %d: Hessian mismatch\ngot=\n%0.4v\nwant=\n%0.4v\n", cas, mat.Formatted(got), mat.Formatted(want))
+		}
+
+		// Test that concurrency works.
+		settings := test.settings
+		if settings == nil {
+			settings = &Settings{}
+		}
+		settings.Concurrent = true
+		got2 := Hessian(nil, test.h.Func, test.x, settings)
+		if !mat.EqualApprox(got, got2, 1e-5) {
+			t.Errorf("Cas %d: Hessian mismatch concurrent\ngot=\n%0.6v\nwant=\n%0.6v\n", cas, mat.Formatted(got2), mat.Formatted(got))
+		}
+	}
+}

diff/fd/jacobian.go

@@ -5,7 +5,6 @@
 package fd
 
 import (
-	"runtime"
 	"sync"
 
 	"gonum.org/v1/gonum/floats"
@@ -49,27 +48,30 @@ func Jacobian(dst *mat.Dense, f func(y, x []float64), x []float64, settings *Jac
 		panic("jacobian: mismatched matrix size")
 	}
 
-	if settings == nil {
-		settings = &JacobianSettings{}
+	// Default settings.
+	formula := Forward
+	step := formula.Step
+	var originValue []float64
+	var concurrent bool
+
+	// Use user settings if provided.
+	if settings != nil {
+		if !settings.Formula.isZero() {
+			formula = settings.Formula
+			step = formula.Step
+			checkFormula(formula)
+			if formula.Derivative != 1 {
+				panic(badDerivOrder)
 			}
-	if settings.OriginValue != nil && len(settings.OriginValue) != m {
+		}
+		if settings.Step != 0 {
+			step = settings.Step
+		}
+		originValue = settings.OriginValue
+		if originValue != nil && len(originValue) != m {
 			panic("jacobian: mismatched OriginValue slice length")
 		}
-
-	formula := settings.Formula
-	if formula.isZero() {
-		formula = Forward
-	}
-	if formula.Derivative == 0 || formula.Stencil == nil || formula.Step == 0 {
-		panic("jacobian: bad formula")
-	}
-	if formula.Derivative != 1 {
-		panic("jacobian: invalid derivative order")
-	}
-
-	step := settings.Step
-	if step == 0 {
-		step = formula.Step
+		concurrent = settings.Concurrent
 	}
 
 	evals := n * len(formula.Stencil)
@@ -79,18 +81,13 @@ func Jacobian(dst *mat.Dense, f func(y, x []float64), x []float64, settings *Jac
 			break
 		}
 	}
-	nWorkers := 1
-	if settings.Concurrent {
-		nWorkers = runtime.GOMAXPROCS(0)
-		if nWorkers > evals {
-			nWorkers = evals
-		}
-	}
+
+	nWorkers := computeWorkers(concurrent, evals)
 	if nWorkers == 1 {
-		jacobianSerial(dst, f, x, settings.OriginValue, formula, step)
-	} else {
-		jacobianConcurrent(dst, f, x, settings.OriginValue, formula, step, nWorkers)
+		jacobianSerial(dst, f, x, originValue, formula, step)
+		return
 	}
+	jacobianConcurrent(dst, f, x, originValue, formula, step, nWorkers)
 }
 
 func jacobianSerial(dst *mat.Dense, f func([]float64, []float64), x, origin []float64, formula Formula, step float64) {

diff/fd/simplefunctions_test.go

@@ -0,0 +1,87 @@
+// Copyright ©2017 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 fd
+
+import (
+	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/mat"
+)
+
+// ConstFunc is a constant function returning the value held by the type.
+type ConstFunc float64
+
+func (c ConstFunc) Func(x []float64) float64 {
+	return float64(c)
+}
+
+func (c ConstFunc) Grad(grad, x []float64) {
+	for i := range grad {
+		grad[i] = 0
+	}
+}
+
+func (c ConstFunc) Hess(dst mat.MutableSymmetric, x []float64) {
+	n := len(x)
+	for i := 0; i < n; i++ {
+		for j := i; j < n; j++ {
+			dst.SetSym(i, j, 0)
+		}
+	}
+}
+
+// LinearFunc is a linear function returning w*x+c.
+type LinearFunc struct {
+	w []float64
+	c float64
+}
+
+func (l LinearFunc) Func(x []float64) float64 {
+	return floats.Dot(l.w, x) + l.c
+}
+
+func (l LinearFunc) Grad(grad, x []float64) {
+	copy(grad, l.w)
+}
+
+func (l LinearFunc) Hess(dst mat.MutableSymmetric, x []float64) {
+	n := len(x)
+	for i := 0; i < n; i++ {
+		for j := i; j < n; j++ {
+			dst.SetSym(i, j, 0)
+		}
+	}
+}
+
+// QuadFunc is a quadratic function returning 0.5*x'*a*x + b*x + c.
+type QuadFunc struct {
+	a *mat.SymDense
+	b *mat.Vector
+	c float64
+}
+
+func (q QuadFunc) Func(x []float64) float64 {
+	v := mat.NewVector(len(x), x)
+	var tmp mat.Vector
+	tmp.MulVec(q.a, v)
+	return 0.5*mat.Dot(&tmp, v) + mat.Dot(q.b, v) + q.c
+}
+
+func (q QuadFunc) Grad(grad, x []float64) {
+	var tmp mat.Vector
+	v := mat.NewVector(len(x), x)
+	tmp.MulVec(q.a, v)
+	for i := range grad {
+		grad[i] = tmp.At(i, 0) + q.b.At(i, 0)
+	}
+}
+
+func (q QuadFunc) Hess(dst mat.MutableSymmetric, x []float64) {
+	n := len(x)
+	for i := 0; i < n; i++ {
+		for j := i; j < n; j++ {
+			dst.SetSym(i, j, q.a.At(i, j))
+		}
+	}
+}

diff/fd/watson_test.go

@@ -0,0 +1,132 @@
+// Copyright ©2017 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 fd
+
+import "gonum.org/v1/gonum/mat"
+
+// Watson implements the Watson's function.
+// Dimension of the problem should be 2 <= dim <= 31. For dim == 9, the problem
+// of minimizing the function is very ill conditioned.
+//
+// This is copied from gonum.org/v1/optimize/functions for testing Hessian-like
+// derivative methods.
+//
+// References:
+//  - Kowalik, J.S., Osborne, M.R.: Methods for Unconstrained Optimization
+//    Problems. Elsevier North-Holland, New York, 1968
+//  - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained
+//    optimization software. ACM Trans Math Softw 7 (1981), 17-41
+type Watson struct{}
+
+func (Watson) Func(x []float64) (sum float64) {
+	for i := 1; i <= 29; i++ {
+		d1 := float64(i) / 29
+
+		d2 := 1.0
+		var s1 float64
+		for j := 1; j < len(x); j++ {
+			s1 += float64(j) * d2 * x[j]
+			d2 *= d1
+		}
+
+		d2 = 1.0
+		var s2 float64
+		for _, v := range x {
+			s2 += d2 * v
+			d2 *= d1
+		}
+
+		t := s1 - s2*s2 - 1
+		sum += t * t
+	}
+	t := x[1] - x[0]*x[0] - 1
+	sum += x[0]*x[0] + t*t
+	return sum
+}
+
+func (Watson) Grad(grad, x []float64) {
+	if len(x) != len(grad) {
+		panic("incorrect size of the gradient")
+	}
+
+	for i := range grad {
+		grad[i] = 0
+	}
+	for i := 1; i <= 29; i++ {
+		d1 := float64(i) / 29
+
+		d2 := 1.0
+		var s1 float64
+		for j := 1; j < len(x); j++ {
+			s1 += float64(j) * d2 * x[j]
+			d2 *= d1
+		}
+
+		d2 = 1.0
+		var s2 float64
+		for _, v := range x {
+			s2 += d2 * v
+			d2 *= d1
+		}
+
+		t := s1 - s2*s2 - 1
+		s3 := 2 * d1 * s2
+		d2 = 2 / d1
+		for j := range x {
+			grad[j] += d2 * (float64(j) - s3) * t
+			d2 *= d1
+		}
+	}
+	t := x[1] - x[0]*x[0] - 1
+	grad[0] += x[0] * (2 - 4*t)
+	grad[1] += 2 * t
+}
+
+func (Watson) Hess(hess mat.MutableSymmetric, x []float64) {
+	dim := len(x)
+	if dim != hess.Symmetric() {
+		panic("incorrect size of the Hessian")
+	}
+
+	for j := 0; j < dim; j++ {
+		for k := j; k < dim; k++ {
+			hess.SetSym(j, k, 0)
+		}
+	}
+	for i := 1; i <= 29; i++ {
+		d1 := float64(i) / 29
+		d2 := 1.0
+		var s1 float64
+		for j := 1; j < dim; j++ {
+			s1 += float64(j) * d2 * x[j]
+			d2 *= d1
+		}
+
+		d2 = 1.0
+		var s2 float64
+		for _, v := range x {
+			s2 += d2 * v
+			d2 *= d1
+		}
+
+		t := s1 - s2*s2 - 1
+		s3 := 2 * d1 * s2
+		d2 = 2 / d1
+		th := 2 * d1 * d1 * t
+		for j := 0; j < dim; j++ {
+			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))
+				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)
+}