1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
|
// 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
//
// The resulting H will be stored in dst. 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.
//
// If the dst matrix is empty it will be resized to the correct dimensions,
// otherwise the dimensions of dst must match the length of x or Hessian will panic.
// Hessian will panic if the derivative order of the formula is not 1.
func Hessian(dst *mat.SymDense, f func(x []float64) float64, x []float64, settings *Settings) {
n := len(x)
if dst.IsEmpty() {
*dst = *(dst.GrowSym(n).(*mat.SymDense))
} else if dst.SymmetricDim() != n {
panic("hessian: dst size mismatch")
}
dst.Zero()
// 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
}
hessianConcurrent(dst, nWorkers, evals, f, x, formula.Stencil, step, originKnown, originValue)
}
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 x in case it is modified during the call.
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 data 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.SymmetricDim()
type run struct {
i, j int
iIdx, jIdx int
result float64
}
send := make(chan run, evals)
ans := make(chan run, evals)
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)
}
}
|
