package/gonum
1
0
Fork 0
mirror of https://github.com/gonum/gonum.git synced 2025-10-06 15:47:01 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

optimize: make Problem.Hess take a *mat.SymDense

Browse Source
This commit is contained in:
Dan Kortschak
2019-04-14 08:39:17 +09:30
parent 3e852a4e0b
commit ae324d6d48
3 changed files with 70 additions and 85 deletions
Download Patch File Download Diff File Expand all files Collapse all files

139
optimize/functions/functions.go
View File

@@ -60,14 +60,13 @@ func (Beale) Grad(grad, x []float64) []float64 {
return grad return grad
} }
func (Beale) Hess(hess mat.Symmetric, x []float64) mat.Symmetric { func (Beale) Hess(dst *mat.SymDense, x []float64) {
if len(x) != 2 { if len(x) != 2 {
panic("dimension of the problem must be 2") panic("dimension of the problem must be 2")
} }
if hess == nil { if dst.IsZero() {
hess = mat.NewSymDense(len(x), nil) *dst = *(dst.GrowSym(len(x)).(*mat.SymDense))
} } else if len(x) != dst.Symmetric() {
if len(x) != hess.Symmetric() {
panic("incorrect size of the Hessian") panic("incorrect size of the Hessian")
} }
@@ -81,11 +80,9 @@ func (Beale) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
h00 := 2 * (t1*t1 + t2*t2 + t3*t3) 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))) 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]))) h11 := 2 * x[0] * (x[0] + 2*f2 + x[1]*(6*f3+x[0]*x[1]*(4+9*x[1]*x[1])))
h := hess.(*mat.SymDense) dst.SetSym(0, 0, h00)
h.SetSym(0, 0, h00) dst.SetSym(0, 1, h01)
h.SetSym(0, 1, h01) dst.SetSym(1, 1, h11)
h.SetSym(1, 1, h11)
return h
} }
func (Beale) Minima() []Minimum { func (Beale) Minima() []Minimum {
@@ -595,25 +592,22 @@ func (BrownBadlyScaled) Grad(grad, x []float64) []float64 {
return grad return grad
} }
func (BrownBadlyScaled) Hess(hess mat.Symmetric, x []float64) mat.Symmetric { func (BrownBadlyScaled) Hess(dst *mat.SymDense, x []float64) {
if len(x) != 2 { if len(x) != 2 {
panic("dimension of the problem must be 2") panic("dimension of the problem must be 2")
} }
if hess == nil { if dst.IsZero() {
hess = mat.NewSymDense(len(x), nil) *dst = *(dst.GrowSym(len(x)).(*mat.SymDense))
} } else if len(x) != dst.Symmetric() {
if len(x) != hess.Symmetric() {
panic("incorrect size of the Hessian") panic("incorrect size of the Hessian")
} }
h00 := 2 + 2*x[1]*x[1] h00 := 2 + 2*x[1]*x[1]
h01 := 4*x[0]*x[1] - 4 h01 := 4*x[0]*x[1] - 4
h11 := 2 + 2*x[0]*x[0] h11 := 2 + 2*x[0]*x[0]
h := hess.(*mat.SymDense) dst.SetSym(0, 0, h00)
h.SetSym(0, 0, h00) dst.SetSym(0, 1, h01)
h.SetSym(0, 1, h01) dst.SetSym(1, 1, h11)
h.SetSym(1, 1, h11)
return h
} }
func (BrownBadlyScaled) Minima() []Minimum { func (BrownBadlyScaled) Minima() []Minimum {
@@ -681,21 +675,19 @@ func (BrownAndDennis) Grad(grad, x []float64) []float64 {
return grad return grad
} }
func (BrownAndDennis) Hess(hess mat.Symmetric, x []float64) mat.Symmetric { func (BrownAndDennis) Hess(dst *mat.SymDense, x []float64) {
if len(x) != 4 { if len(x) != 4 {
panic("dimension of the problem must be 4") panic("dimension of the problem must be 4")
} }
if hess == nil { if dst.IsZero() {
hess = mat.NewSymDense(len(x), nil) *dst = *(dst.GrowSym(len(x)).(*mat.SymDense))
} } else if len(x) != dst.Symmetric() {
if len(x) != hess.Symmetric() {
panic("incorrect size of the Hessian") panic("incorrect size of the Hessian")
} }
h := hess.(*mat.SymDense)
for i := 0; i < 4; i++ { for i := 0; i < 4; i++ {
for j := i; j < 4; j++ { for j := i; j < 4; j++ {
h.SetSym(i, j, 0) dst.SetSym(i, j, 0)
} }
} }
for i := 1; i <= 20; i++ { for i := 1; i <= 20; i++ {
@@ -707,23 +699,22 @@ func (BrownAndDennis) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
s3 := 2 * t1 * t2 s3 := 2 * t1 * t2
r1 := t + 2*t1*t1 r1 := t + 2*t1*t1
r2 := t + 2*t2*t2 r2 := t + 2*t2*t2
h.SetSym(0, 0, h.At(0, 0)+r1) dst.SetSym(0, 0, dst.At(0, 0)+r1)
h.SetSym(0, 1, h.At(0, 1)+d1*r1) dst.SetSym(0, 1, dst.At(0, 1)+d1*r1)
h.SetSym(1, 1, h.At(1, 1)+d1*d1*r1) dst.SetSym(1, 1, dst.At(1, 1)+d1*d1*r1)
h.SetSym(0, 2, h.At(0, 2)+s3) dst.SetSym(0, 2, dst.At(0, 2)+s3)
h.SetSym(1, 2, h.At(1, 2)+d1*s3) dst.SetSym(1, 2, dst.At(1, 2)+d1*s3)
h.SetSym(2, 2, h.At(2, 2)+r2) dst.SetSym(2, 2, dst.At(2, 2)+r2)
h.SetSym(0, 3, h.At(0, 3)+d2*s3) dst.SetSym(0, 3, dst.At(0, 3)+d2*s3)
h.SetSym(1, 3, h.At(1, 3)+d1*d2*s3) dst.SetSym(1, 3, dst.At(1, 3)+d1*d2*s3)
h.SetSym(2, 3, h.At(2, 3)+d2*r2) dst.SetSym(2, 3, dst.At(2, 3)+d2*r2)
h.SetSym(3, 3, h.At(3, 3)+d2*d2*r2) dst.SetSym(3, 3, dst.At(3, 3)+d2*d2*r2)
} }
for i := 0; i < 4; i++ { for i := 0; i < 4; i++ {
for j := i; j < 4; j++ { for j := i; j < 4; j++ {
h.SetSym(i, j, 4*h.At(i, j)) dst.SetSym(i, j, 4*dst.At(i, j))
} }
} }
return h
} }
func (BrownAndDennis) Minima() []Minimum { func (BrownAndDennis) Minima() []Minimum {
@@ -1352,14 +1343,13 @@ func (PowellBadlyScaled) Grad(grad, x []float64) []float64 {
return grad return grad
} }
func (PowellBadlyScaled) Hess(hess mat.Symmetric, x []float64) mat.Symmetric { func (PowellBadlyScaled) Hess(dst *mat.SymDense, x []float64) {
if len(x) != 2 { if len(x) != 2 {
panic("dimension of the problem must be 2") panic("dimension of the problem must be 2")
} }
if hess == nil { if dst.IsZero() {
hess = mat.NewSymDense(len(x), nil) *dst = *(dst.GrowSym(len(x)).(*mat.SymDense))
} } else if len(x) != dst.Symmetric() {
if len(x) != hess.Symmetric() {
panic("incorrect size of the Hessian") panic("incorrect size of the Hessian")
} }
@@ -1368,14 +1358,12 @@ func (PowellBadlyScaled) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
s2 := math.Exp(-x[1]) s2 := math.Exp(-x[1])
t2 := s1 + s2 - 1.0001 t2 := s1 + s2 - 1.0001
h := hess.(*mat.SymDense)
h00 := 2 * (1e8*x[1]*x[1] + s1*(s1+t2)) h00 := 2 * (1e8*x[1]*x[1] + s1*(s1+t2))
h01 := 2 * (1e4*(1+2*t1) + s1*s2) h01 := 2 * (1e4*(1+2*t1) + s1*s2)
h11 := 2 * (1e8*x[0]*x[0] + s2*(s2+t2)) h11 := 2 * (1e8*x[0]*x[0] + s2*(s2+t2))
h.SetSym(0, 0, h00) dst.SetSym(0, 0, h00)
h.SetSym(0, 1, h01) dst.SetSym(0, 1, h01)
h.SetSym(1, 1, h11) dst.SetSym(1, 1, h11)
return h
} }
func (PowellBadlyScaled) Minima() []Minimum { func (PowellBadlyScaled) Minima() []Minimum {
@@ -1619,18 +1607,17 @@ func (Watson) Grad(grad, x []float64) []float64 {
return grad return grad
} }
func (Watson) Hess(hess mat.Symmetric, x []float64) mat.Symmetric { func (Watson) Hess(dst *mat.SymDense, x []float64) {
dim := len(x) dim := len(x)
if hess == nil { if dst.IsZero() {
hess = mat.NewSymDense(len(x), nil) *dst = *(dst.GrowSym(len(x)).(*mat.SymDense))
} } else if len(x) != dst.Symmetric() {
if dim != hess.Symmetric() {
panic("incorrect size of the Hessian") panic("incorrect size of the Hessian")
} }
h := hess.(*mat.SymDense)
for j := 0; j < dim; j++ { for j := 0; j < dim; j++ {
for k := j; k < dim; k++ { for k := j; k < dim; k++ {
h.SetSym(j, k, 0) dst.SetSym(j, k, 0)
} }
} }
for i := 1; i <= 29; i++ { for i := 1; i <= 29; i++ {
@@ -1657,17 +1644,16 @@ func (Watson) Hess(hess mat.Symmetric, x []float64) mat.Symmetric {
v := float64(j) - s3 v := float64(j) - s3
d3 := 1 / d1 d3 := 1 / d1
for k := 0; k <= j; k++ { for k := 0; k <= j; k++ {
h.SetSym(k, j, h.At(k, j)+d2*d3*(v*(float64(k)-s3)-th)) dst.SetSym(k, j, dst.At(k, j)+d2*d3*(v*(float64(k)-s3)-th))
d3 *= d1 d3 *= d1
} }
d2 *= d1 d2 *= d1
} }
} }
t1 := x[1] - x[0]*x[0] - 1 t1 := x[1] - x[0]*x[0] - 1
h.SetSym(0, 0, h.At(0, 0)+8*x[0]*x[0]+2-4*t1) dst.SetSym(0, 0, dst.At(0, 0)+8*x[0]*x[0]+2-4*t1)
h.SetSym(0, 1, h.At(0, 1)-4*x[0]) dst.SetSym(0, 1, dst.At(0, 1)-4*x[0])
h.SetSym(1, 1, h.At(1, 1)+2) dst.SetSym(1, 1, dst.At(1, 1)+2)
return h
} }
func (Watson) Minima() []Minimum { func (Watson) Minima() []Minimum {
@@ -1747,29 +1733,26 @@ func (Wood) Grad(grad, x []float64) []float64 {
return grad return grad
} }
func (Wood) Hess(hess mat.Symmetric, x []float64) mat.Symmetric { func (Wood) Hess(dst *mat.SymDense, x []float64) {
if len(x) != 4 { if len(x) != 4 {
panic("dimension of the problem must be 4") panic("dimension of the problem must be 4")
} }
if hess == nil { if dst.IsZero() {
hess = mat.NewSymDense(len(x), nil) *dst = *(dst.GrowSym(len(x)).(*mat.SymDense))
} } else if len(x) != dst.Symmetric() {
if len(x) != hess.Symmetric() {
panic("incorrect size of the Hessian") panic("incorrect size of the Hessian")
} }
h := hess.(*mat.SymDense)
h.SetSym(0, 0, 400*(3*x[0]*x[0]-x[1])+2) dst.SetSym(0, 0, 400*(3*x[0]*x[0]-x[1])+2)
h.SetSym(0, 1, -400*x[0]) dst.SetSym(0, 1, -400*x[0])
h.SetSym(1, 1, 220.2) dst.SetSym(1, 1, 220.2)
h.SetSym(0, 2, 0) dst.SetSym(0, 2, 0)
h.SetSym(1, 2, 0) dst.SetSym(1, 2, 0)
h.SetSym(2, 2, 360*(3*x[2]*x[2]-x[3])+2) dst.SetSym(2, 2, 360*(3*x[2]*x[2]-x[3])+2)
h.SetSym(0, 3, 0) dst.SetSym(0, 3, 0)
h.SetSym(1, 3, 19.8) dst.SetSym(1, 3, 19.8)
h.SetSym(2, 3, -360*x[2]) dst.SetSym(2, 3, -360*x[2])
h.SetSym(3, 3, 200.2) dst.SetSym(3, 3, 200.2)
return h
} }
func (Wood) Minima() []Minimum { func (Wood) Minima() []Minimum {

8
optimize/minimize.go
View File

@@ -38,7 +38,7 @@ type Location struct {
X []float64 X []float64
F float64 F float64
Gradient []float64 Gradient []float64
Hessian mat.Symmetric Hessian *mat.SymDense
} }
// Method is a type which can search for an optimum of an objective function. // Method is a type which can search for an optimum of an objective function.
@@ -480,7 +480,11 @@ func evaluate(p *Problem, loc *Location, op Operation, x []float64) {
loc.Gradient = p.Grad(loc.Gradient, x) loc.Gradient = p.Grad(loc.Gradient, x)
} }
if op&HessEvaluation != 0 { if op&HessEvaluation != 0 {
loc.Hessian = p.Hess(loc.Hessian, x) // Make sure we have a destination in which to place the Hessian.
if loc.Hessian == nil {
loc.Hessian = &mat.SymDense{}
}
p.Hess(loc.Hessian, x)
} }
} }

8
optimize/types.go
View File

@@ -129,11 +129,9 @@ type Problem struct {
Grad func(grad []float64, x []float64) []float64 Grad func(grad []float64, x []float64) []float64
// Hess evaluates the Hessian at x and stores the result in-place in hess. // Hess evaluates the Hessian at x and stores the result in-place in hess.
// Hess must not modify x. Hess may use (and return) the provided Symmetric // Hess must not modify x. Hess must use the provided mat.SymDense, and
// if it is non-nil, or must allocate a new Symmetric otherwise. Minimize // must resize it if it is zero-sized.
// will 'give back' the returned Symmetric where possible, allowing Hess Hess func(hess *mat.SymDense, x []float64)
// 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 // 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 // error. This can be used to terminate early, for example when the function is