github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/diff/fd/watson_test.go (about)

     1  // Copyright ©2017 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package fd
     6  
     7  import "github.com/jingcheng-WU/gonum/mat"
     8  
     9  // Watson implements the Watson's function.
    10  // Dimension of the problem should be 2 <= dim <= 31. For dim == 9, the problem
    11  // of minimizing the function is very ill conditioned.
    12  //
    13  // This is copied from gonum.org/v1/optimize/functions for testing Hessian-like
    14  // derivative methods.
    15  //
    16  // References:
    17  //  - Kowalik, J.S., Osborne, M.R.: Methods for Unconstrained Optimization
    18  //    Problems. Elsevier North-Holland, New York, 1968
    19  //  - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained
    20  //    optimization software. ACM Trans Math Softw 7 (1981), 17-41
    21  type Watson struct{}
    22  
    23  func (Watson) Func(x []float64) (sum float64) {
    24  	for i := 1; i <= 29; i++ {
    25  		d1 := float64(i) / 29
    26  
    27  		d2 := 1.0
    28  		var s1 float64
    29  		for j := 1; j < len(x); j++ {
    30  			s1 += float64(j) * d2 * x[j]
    31  			d2 *= d1
    32  		}
    33  
    34  		d2 = 1.0
    35  		var s2 float64
    36  		for _, v := range x {
    37  			s2 += d2 * v
    38  			d2 *= d1
    39  		}
    40  
    41  		t := s1 - s2*s2 - 1
    42  		sum += t * t
    43  	}
    44  	t := x[1] - x[0]*x[0] - 1
    45  	sum += x[0]*x[0] + t*t
    46  	return sum
    47  }
    48  
    49  func (Watson) Grad(grad, x []float64) {
    50  	if len(x) != len(grad) {
    51  		panic("incorrect size of the gradient")
    52  	}
    53  
    54  	for i := range grad {
    55  		grad[i] = 0
    56  	}
    57  	for i := 1; i <= 29; i++ {
    58  		d1 := float64(i) / 29
    59  
    60  		d2 := 1.0
    61  		var s1 float64
    62  		for j := 1; j < len(x); j++ {
    63  			s1 += float64(j) * d2 * x[j]
    64  			d2 *= d1
    65  		}
    66  
    67  		d2 = 1.0
    68  		var s2 float64
    69  		for _, v := range x {
    70  			s2 += d2 * v
    71  			d2 *= d1
    72  		}
    73  
    74  		t := s1 - s2*s2 - 1
    75  		s3 := 2 * d1 * s2
    76  		d2 = 2 / d1
    77  		for j := range x {
    78  			grad[j] += d2 * (float64(j) - s3) * t
    79  			d2 *= d1
    80  		}
    81  	}
    82  	t := x[1] - x[0]*x[0] - 1
    83  	grad[0] += x[0] * (2 - 4*t)
    84  	grad[1] += 2 * t
    85  }
    86  
    87  func (Watson) Hess(hess mat.MutableSymmetric, x []float64) {
    88  	dim := len(x)
    89  	if dim != hess.Symmetric() {
    90  		panic("incorrect size of the Hessian")
    91  	}
    92  
    93  	for j := 0; j < dim; j++ {
    94  		for k := j; k < dim; k++ {
    95  			hess.SetSym(j, k, 0)
    96  		}
    97  	}
    98  	for i := 1; i <= 29; i++ {
    99  		d1 := float64(i) / 29
   100  		d2 := 1.0
   101  		var s1 float64
   102  		for j := 1; j < dim; j++ {
   103  			s1 += float64(j) * d2 * x[j]
   104  			d2 *= d1
   105  		}
   106  
   107  		d2 = 1.0
   108  		var s2 float64
   109  		for _, v := range x {
   110  			s2 += d2 * v
   111  			d2 *= d1
   112  		}
   113  
   114  		t := s1 - s2*s2 - 1
   115  		s3 := 2 * d1 * s2
   116  		d2 = 2 / d1
   117  		th := 2 * d1 * d1 * t
   118  		for j := 0; j < dim; j++ {
   119  			v := float64(j) - s3
   120  			d3 := 1 / d1
   121  			for k := 0; k <= j; k++ {
   122  				hess.SetSym(k, j, hess.At(k, j)+d2*d3*(v*(float64(k)-s3)-th))
   123  				d3 *= d1
   124  			}
   125  			d2 *= d1
   126  		}
   127  	}
   128  	t1 := x[1] - x[0]*x[0] - 1
   129  	hess.SetSym(0, 0, hess.At(0, 0)+8*x[0]*x[0]+2-4*t1)
   130  	hess.SetSym(0, 1, hess.At(0, 1)-4*x[0])
   131  	hess.SetSym(1, 1, hess.At(1, 1)+2)
   132  }