github.com/gopherd/gonum@v0.0.4/optimize/functions/functions.go (about) 1 // Copyright ©2015 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 functions 6 7 import ( 8 "math" 9 10 "github.com/gopherd/gonum/floats" 11 "github.com/gopherd/gonum/mat" 12 ) 13 14 // Beale implements the Beale's function. 15 // 16 // Standard starting points: 17 // Easy: [1, 1] 18 // Hard: [1, 4] 19 // 20 // References: 21 // - Beale, E.: On an Iterative Method for Finding a Local Minimum of a 22 // Function of More than One Variable. Technical Report 25, Statistical 23 // Techniques Research Group, Princeton University (1958) 24 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 25 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 26 type Beale struct{} 27 28 func (Beale) Func(x []float64) float64 { 29 if len(x) != 2 { 30 panic("dimension of the problem must be 2") 31 } 32 33 f1 := 1.5 - x[0]*(1-x[1]) 34 f2 := 2.25 - x[0]*(1-x[1]*x[1]) 35 f3 := 2.625 - x[0]*(1-x[1]*x[1]*x[1]) 36 return f1*f1 + f2*f2 + f3*f3 37 } 38 39 func (Beale) Grad(grad, x []float64) { 40 if len(x) != 2 { 41 panic("dimension of the problem must be 2") 42 } 43 if len(x) != len(grad) { 44 panic("incorrect size of the gradient") 45 } 46 47 t1 := 1 - x[1] 48 t2 := 1 - x[1]*x[1] 49 t3 := 1 - x[1]*x[1]*x[1] 50 51 f1 := 1.5 - x[0]*t1 52 f2 := 2.25 - x[0]*t2 53 f3 := 2.625 - x[0]*t3 54 55 grad[0] = -2 * (f1*t1 + f2*t2 + f3*t3) 56 grad[1] = 2 * x[0] * (f1 + 2*f2*x[1] + 3*f3*x[1]*x[1]) 57 } 58 59 func (Beale) Hess(dst *mat.SymDense, x []float64) { 60 if len(x) != 2 { 61 panic("dimension of the problem must be 2") 62 } 63 if len(x) != dst.SymmetricDim() { 64 panic("incorrect size of the Hessian") 65 } 66 67 t1 := 1 - x[1] 68 t2 := 1 - x[1]*x[1] 69 t3 := 1 - x[1]*x[1]*x[1] 70 f1 := 1.5 - x[1]*t1 71 f2 := 2.25 - x[1]*t2 72 f3 := 2.625 - x[1]*t3 73 74 h00 := 2 * (t1*t1 + t2*t2 + t3*t3) 75 h01 := 2 * (f1 + x[1]*(2*f2+3*x[1]*f3) - x[0]*(t1+x[1]*(2*t2+3*x[1]*t3))) 76 h11 := 2 * x[0] * (x[0] + 2*f2 + x[1]*(6*f3+x[0]*x[1]*(4+9*x[1]*x[1]))) 77 dst.SetSym(0, 0, h00) 78 dst.SetSym(0, 1, h01) 79 dst.SetSym(1, 1, h11) 80 } 81 82 func (Beale) Minima() []Minimum { 83 return []Minimum{ 84 { 85 X: []float64{3, 0.5}, 86 F: 0, 87 Global: true, 88 }, 89 } 90 } 91 92 // BiggsEXP2 implements the Biggs' EXP2 function. 93 // 94 // Standard starting point: 95 // [1, 2] 96 // 97 // Reference: 98 // Biggs, M.C.: Minimization algorithms making use of non-quadratic properties 99 // of the objective function. IMA J Appl Math 8 (1971), 315-327; doi:10.1093/imamat/8.3.315 100 type BiggsEXP2 struct{} 101 102 func (BiggsEXP2) Func(x []float64) (sum float64) { 103 if len(x) != 2 { 104 panic("dimension of the problem must be 2") 105 } 106 107 for i := 1; i <= 10; i++ { 108 z := float64(i) / 10 109 y := math.Exp(-z) - 5*math.Exp(-10*z) 110 f := math.Exp(-x[0]*z) - 5*math.Exp(-x[1]*z) - y 111 sum += f * f 112 } 113 return sum 114 } 115 116 func (BiggsEXP2) Grad(grad, x []float64) { 117 if len(x) != 2 { 118 panic("dimension of the problem must be 2") 119 } 120 if len(x) != len(grad) { 121 panic("incorrect size of the gradient") 122 } 123 124 for i := range grad { 125 grad[i] = 0 126 } 127 for i := 1; i <= 10; i++ { 128 z := float64(i) / 10 129 y := math.Exp(-z) - 5*math.Exp(-10*z) 130 f := math.Exp(-x[0]*z) - 5*math.Exp(-x[1]*z) - y 131 132 dfdx0 := -z * math.Exp(-x[0]*z) 133 dfdx1 := 5 * z * math.Exp(-x[1]*z) 134 135 grad[0] += 2 * f * dfdx0 136 grad[1] += 2 * f * dfdx1 137 } 138 } 139 140 func (BiggsEXP2) Minima() []Minimum { 141 return []Minimum{ 142 { 143 X: []float64{1, 10}, 144 F: 0, 145 Global: true, 146 }, 147 } 148 } 149 150 // BiggsEXP3 implements the Biggs' EXP3 function. 151 // 152 // Standard starting point: 153 // [1, 2, 1] 154 // 155 // Reference: 156 // Biggs, M.C.: Minimization algorithms making use of non-quadratic properties 157 // of the objective function. IMA J Appl Math 8 (1971), 315-327; doi:10.1093/imamat/8.3.315 158 type BiggsEXP3 struct{} 159 160 func (BiggsEXP3) Func(x []float64) (sum float64) { 161 if len(x) != 3 { 162 panic("dimension of the problem must be 3") 163 } 164 165 for i := 1; i <= 10; i++ { 166 z := float64(i) / 10 167 y := math.Exp(-z) - 5*math.Exp(-10*z) 168 f := math.Exp(-x[0]*z) - x[2]*math.Exp(-x[1]*z) - y 169 sum += f * f 170 } 171 return sum 172 } 173 174 func (BiggsEXP3) Grad(grad, x []float64) { 175 if len(x) != 3 { 176 panic("dimension of the problem must be 3") 177 } 178 if len(x) != len(grad) { 179 panic("incorrect size of the gradient") 180 } 181 182 for i := range grad { 183 grad[i] = 0 184 } 185 for i := 1; i <= 10; i++ { 186 z := float64(i) / 10 187 y := math.Exp(-z) - 5*math.Exp(-10*z) 188 f := math.Exp(-x[0]*z) - x[2]*math.Exp(-x[1]*z) - y 189 190 dfdx0 := -z * math.Exp(-x[0]*z) 191 dfdx1 := x[2] * z * math.Exp(-x[1]*z) 192 dfdx2 := -math.Exp(-x[1] * z) 193 194 grad[0] += 2 * f * dfdx0 195 grad[1] += 2 * f * dfdx1 196 grad[2] += 2 * f * dfdx2 197 } 198 } 199 200 func (BiggsEXP3) Minima() []Minimum { 201 return []Minimum{ 202 { 203 X: []float64{1, 10, 5}, 204 F: 0, 205 Global: true, 206 }, 207 } 208 } 209 210 // BiggsEXP4 implements the Biggs' EXP4 function. 211 // 212 // Standard starting point: 213 // [1, 2, 1, 1] 214 // 215 // Reference: 216 // Biggs, M.C.: Minimization algorithms making use of non-quadratic properties 217 // of the objective function. IMA J Appl Math 8 (1971), 315-327; doi:10.1093/imamat/8.3.315 218 type BiggsEXP4 struct{} 219 220 func (BiggsEXP4) Func(x []float64) (sum float64) { 221 if len(x) != 4 { 222 panic("dimension of the problem must be 4") 223 } 224 225 for i := 1; i <= 10; i++ { 226 z := float64(i) / 10 227 y := math.Exp(-z) - 5*math.Exp(-10*z) 228 f := x[2]*math.Exp(-x[0]*z) - x[3]*math.Exp(-x[1]*z) - y 229 sum += f * f 230 } 231 return sum 232 } 233 234 func (BiggsEXP4) Grad(grad, x []float64) { 235 if len(x) != 4 { 236 panic("dimension of the problem must be 4") 237 } 238 if len(x) != len(grad) { 239 panic("incorrect size of the gradient") 240 } 241 242 for i := range grad { 243 grad[i] = 0 244 } 245 for i := 1; i <= 10; i++ { 246 z := float64(i) / 10 247 y := math.Exp(-z) - 5*math.Exp(-10*z) 248 f := x[2]*math.Exp(-x[0]*z) - x[3]*math.Exp(-x[1]*z) - y 249 250 dfdx0 := -z * x[2] * math.Exp(-x[0]*z) 251 dfdx1 := z * x[3] * math.Exp(-x[1]*z) 252 dfdx2 := math.Exp(-x[0] * z) 253 dfdx3 := -math.Exp(-x[1] * z) 254 255 grad[0] += 2 * f * dfdx0 256 grad[1] += 2 * f * dfdx1 257 grad[2] += 2 * f * dfdx2 258 grad[3] += 2 * f * dfdx3 259 } 260 } 261 262 func (BiggsEXP4) Minima() []Minimum { 263 return []Minimum{ 264 { 265 X: []float64{1, 10, 1, 5}, 266 F: 0, 267 Global: true, 268 }, 269 } 270 } 271 272 // BiggsEXP5 implements the Biggs' EXP5 function. 273 // 274 // Standard starting point: 275 // [1, 2, 1, 1, 1] 276 // 277 // Reference: 278 // Biggs, M.C.: Minimization algorithms making use of non-quadratic properties 279 // of the objective function. IMA J Appl Math 8 (1971), 315-327; doi:10.1093/imamat/8.3.315 280 type BiggsEXP5 struct{} 281 282 func (BiggsEXP5) Func(x []float64) (sum float64) { 283 if len(x) != 5 { 284 panic("dimension of the problem must be 5") 285 } 286 287 for i := 1; i <= 11; i++ { 288 z := float64(i) / 10 289 y := math.Exp(-z) - 5*math.Exp(-10*z) + 3*math.Exp(-4*z) 290 f := x[2]*math.Exp(-x[0]*z) - x[3]*math.Exp(-x[1]*z) + 3*math.Exp(-x[4]*z) - y 291 sum += f * f 292 } 293 return sum 294 } 295 296 func (BiggsEXP5) Grad(grad, x []float64) { 297 if len(x) != 5 { 298 panic("dimension of the problem must be 5") 299 } 300 if len(x) != len(grad) { 301 panic("incorrect size of the gradient") 302 } 303 304 for i := range grad { 305 grad[i] = 0 306 } 307 for i := 1; i <= 11; i++ { 308 z := float64(i) / 10 309 y := math.Exp(-z) - 5*math.Exp(-10*z) + 3*math.Exp(-4*z) 310 f := x[2]*math.Exp(-x[0]*z) - x[3]*math.Exp(-x[1]*z) + 3*math.Exp(-x[4]*z) - y 311 312 dfdx0 := -z * x[2] * math.Exp(-x[0]*z) 313 dfdx1 := z * x[3] * math.Exp(-x[1]*z) 314 dfdx2 := math.Exp(-x[0] * z) 315 dfdx3 := -math.Exp(-x[1] * z) 316 dfdx4 := -3 * z * math.Exp(-x[4]*z) 317 318 grad[0] += 2 * f * dfdx0 319 grad[1] += 2 * f * dfdx1 320 grad[2] += 2 * f * dfdx2 321 grad[3] += 2 * f * dfdx3 322 grad[4] += 2 * f * dfdx4 323 } 324 } 325 326 func (BiggsEXP5) Minima() []Minimum { 327 return []Minimum{ 328 { 329 X: []float64{1, 10, 1, 5, 4}, 330 F: 0, 331 Global: true, 332 }, 333 } 334 } 335 336 // BiggsEXP6 implements the Biggs' EXP6 function. 337 // 338 // Standard starting point: 339 // [1, 2, 1, 1, 1, 1] 340 // 341 // References: 342 // - Biggs, M.C.: Minimization algorithms making use of non-quadratic 343 // properties of the objective function. IMA J Appl Math 8 (1971), 315-327; 344 // doi:10.1093/imamat/8.3.315 345 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 346 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 347 type BiggsEXP6 struct{} 348 349 func (BiggsEXP6) Func(x []float64) (sum float64) { 350 if len(x) != 6 { 351 panic("dimension of the problem must be 6") 352 } 353 354 for i := 1; i <= 13; i++ { 355 z := float64(i) / 10 356 y := math.Exp(-z) - 5*math.Exp(-10*z) + 3*math.Exp(-4*z) 357 f := x[2]*math.Exp(-x[0]*z) - x[3]*math.Exp(-x[1]*z) + x[5]*math.Exp(-x[4]*z) - y 358 sum += f * f 359 } 360 return sum 361 } 362 363 func (BiggsEXP6) Grad(grad, x []float64) { 364 if len(x) != 6 { 365 panic("dimension of the problem must be 6") 366 } 367 if len(x) != len(grad) { 368 panic("incorrect size of the gradient") 369 } 370 371 for i := range grad { 372 grad[i] = 0 373 } 374 for i := 1; i <= 13; i++ { 375 z := float64(i) / 10 376 y := math.Exp(-z) - 5*math.Exp(-10*z) + 3*math.Exp(-4*z) 377 f := x[2]*math.Exp(-x[0]*z) - x[3]*math.Exp(-x[1]*z) + x[5]*math.Exp(-x[4]*z) - y 378 379 dfdx0 := -z * x[2] * math.Exp(-x[0]*z) 380 dfdx1 := z * x[3] * math.Exp(-x[1]*z) 381 dfdx2 := math.Exp(-x[0] * z) 382 dfdx3 := -math.Exp(-x[1] * z) 383 dfdx4 := -z * x[5] * math.Exp(-x[4]*z) 384 dfdx5 := math.Exp(-x[4] * z) 385 386 grad[0] += 2 * f * dfdx0 387 grad[1] += 2 * f * dfdx1 388 grad[2] += 2 * f * dfdx2 389 grad[3] += 2 * f * dfdx3 390 grad[4] += 2 * f * dfdx4 391 grad[5] += 2 * f * dfdx5 392 } 393 } 394 395 func (BiggsEXP6) Minima() []Minimum { 396 return []Minimum{ 397 { 398 X: []float64{1, 10, 1, 5, 4, 3}, 399 F: 0, 400 Global: true, 401 }, 402 { 403 X: []float64{1.7114159947956764, 17.68319817846745, 1.1631436609697268, 404 5.1865615510738605, 1.7114159947949301, 1.1631436609697998}, 405 F: 0.005655649925499929, 406 Global: false, 407 }, 408 { 409 // X: []float64{1.22755594752403, X[1] >> 0, 0.83270306333466, X[3] << 0, X[4] = X[0], X[5] = X[2]}, 410 X: []float64{1.22755594752403, 1000, 0.83270306333466, -1000, 1.22755594752403, 0.83270306333466}, 411 F: 0.306366772624790, 412 Global: false, 413 }, 414 } 415 } 416 417 // Box3D implements the Box' three-dimensional function. 418 // 419 // Standard starting point: 420 // [0, 10, 20] 421 // 422 // References: 423 // - Box, M.J.: A comparison of several current optimization methods, and the 424 // use of transformations in constrained problems. Comput J 9 (1966), 67-77 425 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 426 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 427 type Box3D struct{} 428 429 func (Box3D) Func(x []float64) (sum float64) { 430 if len(x) != 3 { 431 panic("dimension of the problem must be 3") 432 } 433 434 for i := 1; i <= 10; i++ { 435 c := -float64(i) / 10 436 y := math.Exp(c) - math.Exp(10*c) 437 f := math.Exp(c*x[0]) - math.Exp(c*x[1]) - x[2]*y 438 sum += f * f 439 } 440 return sum 441 } 442 443 func (Box3D) Grad(grad, x []float64) { 444 if len(x) != 3 { 445 panic("dimension of the problem must be 3") 446 } 447 if len(x) != len(grad) { 448 panic("incorrect size of the gradient") 449 } 450 451 grad[0] = 0 452 grad[1] = 0 453 grad[2] = 0 454 for i := 1; i <= 10; i++ { 455 c := -float64(i) / 10 456 y := math.Exp(c) - math.Exp(10*c) 457 f := math.Exp(c*x[0]) - math.Exp(c*x[1]) - x[2]*y 458 grad[0] += 2 * f * c * math.Exp(c*x[0]) 459 grad[1] += -2 * f * c * math.Exp(c*x[1]) 460 grad[2] += -2 * f * y 461 } 462 } 463 464 func (Box3D) Minima() []Minimum { 465 return []Minimum{ 466 { 467 X: []float64{1, 10, 1}, 468 F: 0, 469 Global: true, 470 }, 471 { 472 X: []float64{10, 1, -1}, 473 F: 0, 474 Global: true, 475 }, 476 { 477 // Any point at the line {a, a, 0}. 478 X: []float64{1, 1, 0}, 479 F: 0, 480 Global: true, 481 }, 482 } 483 } 484 485 // BraninHoo implements the Branin-Hoo function. BraninHoo is a 2-dimensional 486 // test function with three global minima. It is typically evaluated in the domain 487 // x_0 ∈ [-5, 10], x_1 ∈ [0, 15]. 488 // f(x) = (x_1 - (5.1/(4π^2))*x_0^2 + (5/π)*x_0 - 6)^2 + 10*(1-1/(8π))cos(x_0) + 10 489 // It has a minimum value of 0.397887 at x^* = {(-π, 12.275), (π, 2.275), (9.424778, 2.475)} 490 // 491 // Reference: 492 // https://www.sfu.ca/~ssurjano/branin.html (obtained June 2017) 493 type BraninHoo struct{} 494 495 func (BraninHoo) Func(x []float64) float64 { 496 if len(x) != 2 { 497 panic("functions: dimension of the problem must be 2") 498 } 499 a, b, c, r, s, t := 1.0, 5.1/(4*math.Pi*math.Pi), 5/math.Pi, 6.0, 10.0, 1/(8*math.Pi) 500 501 term := x[1] - b*x[0]*x[0] + c*x[0] - r 502 return a*term*term + s*(1-t)*math.Cos(x[0]) + s 503 } 504 505 func (BraninHoo) Minima() []Minimum { 506 return []Minimum{ 507 { 508 X: []float64{-math.Pi, 12.275}, 509 F: 0.397887, 510 Global: true, 511 }, 512 { 513 X: []float64{math.Pi, 2.275}, 514 F: 0.397887, 515 Global: true, 516 }, 517 { 518 X: []float64{9.424778, 2.475}, 519 F: 0.397887, 520 Global: true, 521 }, 522 } 523 } 524 525 // BrownBadlyScaled implements the Brown's badly scaled function. 526 // 527 // Standard starting point: 528 // [1, 1] 529 // 530 // References: 531 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 532 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 533 type BrownBadlyScaled struct{} 534 535 func (BrownBadlyScaled) Func(x []float64) float64 { 536 if len(x) != 2 { 537 panic("dimension of the problem must be 2") 538 } 539 540 f1 := x[0] - 1e6 541 f2 := x[1] - 2e-6 542 f3 := x[0]*x[1] - 2 543 return f1*f1 + f2*f2 + f3*f3 544 } 545 546 func (BrownBadlyScaled) Grad(grad, x []float64) { 547 if len(x) != 2 { 548 panic("dimension of the problem must be 2") 549 } 550 if len(x) != len(grad) { 551 panic("incorrect size of the gradient") 552 } 553 554 f1 := x[0] - 1e6 555 f2 := x[1] - 2e-6 556 f3 := float64(x[0]*x[1]) - 2 // Prevent fused multiply subtract. 557 grad[0] = 2*f1 + 2*f3*x[1] 558 grad[1] = 2*f2 + 2*f3*x[0] 559 } 560 561 func (BrownBadlyScaled) Hess(dst *mat.SymDense, x []float64) { 562 if len(x) != 2 { 563 panic("dimension of the problem must be 2") 564 } 565 if len(x) != dst.SymmetricDim() { 566 panic("incorrect size of the Hessian") 567 } 568 569 h00 := 2 + 2*x[1]*x[1] 570 h01 := 4*x[0]*x[1] - 4 571 h11 := 2 + 2*x[0]*x[0] 572 dst.SetSym(0, 0, h00) 573 dst.SetSym(0, 1, h01) 574 dst.SetSym(1, 1, h11) 575 } 576 577 func (BrownBadlyScaled) Minima() []Minimum { 578 return []Minimum{ 579 { 580 X: []float64{1e6, 2e-6}, 581 F: 0, 582 Global: true, 583 }, 584 } 585 } 586 587 // BrownAndDennis implements the Brown and Dennis function. 588 // 589 // Standard starting point: 590 // [25, 5, -5, -1] 591 // 592 // References: 593 // - Brown, K.M., Dennis, J.E.: New computational algorithms for minimizing a 594 // sum of squares of nonlinear functions. Research Report Number 71-6, Yale 595 // University (1971) 596 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 597 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 598 type BrownAndDennis struct{} 599 600 func (BrownAndDennis) Func(x []float64) (sum float64) { 601 if len(x) != 4 { 602 panic("dimension of the problem must be 4") 603 } 604 605 for i := 1; i <= 20; i++ { 606 c := float64(i) / 5 607 f1 := x[0] + c*x[1] - math.Exp(c) 608 f2 := x[2] + x[3]*math.Sin(c) - math.Cos(c) 609 f := f1*f1 + f2*f2 610 sum += f * f 611 } 612 return sum 613 } 614 615 func (BrownAndDennis) Grad(grad, x []float64) { 616 if len(x) != 4 { 617 panic("dimension of the problem must be 4") 618 } 619 if len(x) != len(grad) { 620 panic("incorrect size of the gradient") 621 } 622 623 for i := range grad { 624 grad[i] = 0 625 } 626 for i := 1; i <= 20; i++ { 627 c := float64(i) / 5 628 f1 := x[0] + c*x[1] - math.Exp(c) 629 f2 := x[2] + x[3]*math.Sin(c) - math.Cos(c) 630 f := f1*f1 + f2*f2 631 grad[0] += 4 * f * f1 632 grad[1] += 4 * f * f1 * c 633 grad[2] += 4 * f * f2 634 grad[3] += 4 * f * f2 * math.Sin(c) 635 } 636 } 637 638 func (BrownAndDennis) Hess(dst *mat.SymDense, x []float64) { 639 if len(x) != 4 { 640 panic("dimension of the problem must be 4") 641 } 642 if len(x) != dst.SymmetricDim() { 643 panic("incorrect size of the Hessian") 644 } 645 646 for i := 0; i < 4; i++ { 647 for j := i; j < 4; j++ { 648 dst.SetSym(i, j, 0) 649 } 650 } 651 for i := 1; i <= 20; i++ { 652 d1 := float64(i) / 5 653 d2 := math.Sin(d1) 654 t1 := x[0] + d1*x[1] - math.Exp(d1) 655 t2 := x[2] + d2*x[3] - math.Cos(d1) 656 t := t1*t1 + t2*t2 657 s3 := 2 * t1 * t2 658 r1 := t + 2*t1*t1 659 r2 := t + 2*t2*t2 660 dst.SetSym(0, 0, dst.At(0, 0)+r1) 661 dst.SetSym(0, 1, dst.At(0, 1)+d1*r1) 662 dst.SetSym(1, 1, dst.At(1, 1)+d1*d1*r1) 663 dst.SetSym(0, 2, dst.At(0, 2)+s3) 664 dst.SetSym(1, 2, dst.At(1, 2)+d1*s3) 665 dst.SetSym(2, 2, dst.At(2, 2)+r2) 666 dst.SetSym(0, 3, dst.At(0, 3)+d2*s3) 667 dst.SetSym(1, 3, dst.At(1, 3)+d1*d2*s3) 668 dst.SetSym(2, 3, dst.At(2, 3)+d2*r2) 669 dst.SetSym(3, 3, dst.At(3, 3)+d2*d2*r2) 670 } 671 for i := 0; i < 4; i++ { 672 for j := i; j < 4; j++ { 673 dst.SetSym(i, j, 4*dst.At(i, j)) 674 } 675 } 676 } 677 678 func (BrownAndDennis) Minima() []Minimum { 679 return []Minimum{ 680 { 681 X: []float64{-11.594439904762162, 13.203630051207202, -0.4034394881768612, 0.2367787744557347}, 682 F: 85822.20162635634, 683 Global: true, 684 }, 685 } 686 } 687 688 // ExtendedPowellSingular implements the extended Powell's function. 689 // Its Hessian matrix is singular at the minimizer. 690 // 691 // Standard starting point: 692 // [3, -1, 0, 3, 3, -1, 0, 3, ..., 3, -1, 0, 3] 693 // 694 // References: 695 // - Spedicato E.: Computational experience with quasi-Newton algorithms for 696 // minimization problems of moderatly large size. Towards Global 697 // Optimization 2 (1978), 209-219 698 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 699 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 700 type ExtendedPowellSingular struct{} 701 702 func (ExtendedPowellSingular) Func(x []float64) (sum float64) { 703 if len(x)%4 != 0 { 704 panic("dimension of the problem must be a multiple of 4") 705 } 706 707 for i := 0; i < len(x); i += 4 { 708 f1 := x[i] + 10*x[i+1] 709 f2 := x[i+2] - x[i+3] 710 t := x[i+1] - 2*x[i+2] 711 f3 := t * t 712 t = x[i] - x[i+3] 713 f4 := t * t 714 sum += f1*f1 + 5*f2*f2 + f3*f3 + 10*f4*f4 715 } 716 return sum 717 } 718 719 func (ExtendedPowellSingular) Grad(grad, x []float64) { 720 if len(x)%4 != 0 { 721 panic("dimension of the problem must be a multiple of 4") 722 } 723 if len(x) != len(grad) { 724 panic("incorrect size of the gradient") 725 } 726 727 for i := 0; i < len(x); i += 4 { 728 f1 := x[i] + 10*x[i+1] 729 f2 := x[i+2] - x[i+3] 730 t1 := x[i+1] - 2*x[i+2] 731 f3 := t1 * t1 732 t2 := x[i] - x[i+3] 733 f4 := t2 * t2 734 grad[i] = 2*f1 + 40*f4*t2 735 grad[i+1] = 20*f1 + 4*f3*t1 736 grad[i+2] = 10*f2 - 8*f3*t1 737 grad[i+3] = -10*f2 - 40*f4*t2 738 } 739 } 740 741 func (ExtendedPowellSingular) Minima() []Minimum { 742 return []Minimum{ 743 { 744 X: []float64{0, 0, 0, 0}, 745 F: 0, 746 Global: true, 747 }, 748 { 749 X: []float64{0, 0, 0, 0, 0, 0, 0, 0}, 750 F: 0, 751 Global: true, 752 }, 753 { 754 X: []float64{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 755 F: 0, 756 Global: true, 757 }, 758 } 759 } 760 761 // ExtendedRosenbrock implements the extended, multidimensional Rosenbrock 762 // function. 763 // 764 // Standard starting point: 765 // Easy: [-1.2, 1, -1.2, 1, ...] 766 // Hard: any point far from the minimum 767 // 768 // References: 769 // - Rosenbrock, H.H.: An Automatic Method for Finding the Greatest or Least 770 // Value of a Function. Computer J 3 (1960), 175-184 771 // - http://en.wikipedia.org/wiki/Rosenbrock_function 772 type ExtendedRosenbrock struct{} 773 774 func (ExtendedRosenbrock) Func(x []float64) (sum float64) { 775 for i := 0; i < len(x)-1; i++ { 776 a := 1 - x[i] 777 b := x[i+1] - x[i]*x[i] 778 sum += a*a + 100*b*b 779 } 780 return sum 781 } 782 783 func (ExtendedRosenbrock) Grad(grad, x []float64) { 784 if len(x) != len(grad) { 785 panic("incorrect size of the gradient") 786 } 787 788 dim := len(x) 789 for i := range grad { 790 grad[i] = 0 791 } 792 // Prevent fused multiply add and fused multiply subtract. 793 for i := 0; i < dim-1; i++ { 794 grad[i] -= float64(2 * (1 - x[i])) 795 grad[i] -= float64(400 * (x[i+1] - float64(x[i]*x[i])) * x[i]) 796 } 797 for i := 1; i < dim; i++ { 798 grad[i] += float64(200 * (x[i] - float64(x[i-1]*x[i-1]))) 799 } 800 } 801 802 func (ExtendedRosenbrock) Minima() []Minimum { 803 return []Minimum{ 804 { 805 X: []float64{1, 1}, 806 F: 0, 807 Global: true, 808 }, 809 { 810 X: []float64{1, 1, 1}, 811 F: 0, 812 Global: true, 813 }, 814 { 815 X: []float64{1, 1, 1, 1}, 816 F: 0, 817 Global: true, 818 }, 819 { 820 X: []float64{-0.7756592265653526, 0.6130933654850433, 821 0.38206284633839305, 0.14597201855219452}, 822 F: 3.701428610430017, 823 Global: false, 824 }, 825 { 826 X: []float64{1, 1, 1, 1, 1}, 827 F: 0, 828 Global: true, 829 }, 830 { 831 X: []float64{-0.9620510206947502, 0.9357393959767103, 832 0.8807136041943204, 0.7778776758544063, 0.6050936785926526}, 833 F: 3.930839434133027, 834 Global: false, 835 }, 836 { 837 X: []float64{-0.9865749795709938, 0.9833982288361819, 0.972106670053092, 838 0.9474374368264362, 0.8986511848517299, 0.8075739520354182}, 839 F: 3.973940500930295, 840 Global: false, 841 }, 842 { 843 X: []float64{-0.9917225725614055, 0.9935553935033712, 0.992173321594692, 844 0.9868987626903134, 0.975164756608872, 0.9514319827049906, 0.9052228177139495}, 845 F: 3.9836005364248543, 846 Global: false, 847 }, 848 { 849 X: []float64{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, 850 F: 0, 851 Global: true, 852 }, 853 } 854 } 855 856 // Gaussian implements the Gaussian function. 857 // The function has one global minimum and a number of false local minima 858 // caused by the finite floating point precision. 859 // 860 // Standard starting point: 861 // [0.4, 1, 0] 862 // 863 // Reference: 864 // More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained optimization 865 // software. ACM Trans Math Softw 7 (1981), 17-41 866 type Gaussian struct{} 867 868 func (Gaussian) y(i int) (yi float64) { 869 switch i { 870 case 1, 15: 871 yi = 0.0009 872 case 2, 14: 873 yi = 0.0044 874 case 3, 13: 875 yi = 0.0175 876 case 4, 12: 877 yi = 0.0540 878 case 5, 11: 879 yi = 0.1295 880 case 6, 10: 881 yi = 0.2420 882 case 7, 9: 883 yi = 0.3521 884 case 8: 885 yi = 0.3989 886 } 887 return yi 888 } 889 890 func (g Gaussian) Func(x []float64) (sum float64) { 891 if len(x) != 3 { 892 panic("dimension of the problem must be 3") 893 } 894 895 for i := 1; i <= 15; i++ { 896 c := 0.5 * float64(8-i) 897 b := c - x[2] 898 d := b * b 899 e := math.Exp(-0.5 * x[1] * d) 900 f := x[0]*e - g.y(i) 901 sum += f * f 902 } 903 return sum 904 } 905 906 func (g Gaussian) Grad(grad, x []float64) { 907 if len(x) != 3 { 908 panic("dimension of the problem must be 3") 909 } 910 if len(x) != len(grad) { 911 panic("incorrect size of the gradient") 912 } 913 914 grad[0] = 0 915 grad[1] = 0 916 grad[2] = 0 917 for i := 1; i <= 15; i++ { 918 c := 0.5 * float64(8-i) 919 b := c - x[2] 920 d := b * b 921 e := math.Exp(-0.5 * x[1] * d) 922 f := x[0]*e - g.y(i) 923 grad[0] += 2 * f * e 924 grad[1] -= f * e * d * x[0] 925 grad[2] += 2 * f * e * x[0] * x[1] * b 926 } 927 } 928 929 func (Gaussian) Minima() []Minimum { 930 return []Minimum{ 931 { 932 X: []float64{0.398956137837997, 1.0000190844805048, 0}, 933 F: 1.12793276961912e-08, 934 Global: true, 935 }, 936 } 937 } 938 939 // GulfResearchAndDevelopment implements the Gulf Research and Development function. 940 // 941 // Standard starting point: 942 // [5, 2.5, 0.15] 943 // 944 // References: 945 // - Cox, R.A.: Comparison of the performance of seven optimization algorithms 946 // on twelve unconstrained minimization problems. Ref. 1335CNO4, Gulf 947 // Research and Development Company, Pittsburg (1969) 948 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 949 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 950 type GulfResearchAndDevelopment struct{} 951 952 func (GulfResearchAndDevelopment) Func(x []float64) (sum float64) { 953 if len(x) != 3 { 954 panic("dimension of the problem must be 3") 955 } 956 957 for i := 1; i <= 99; i++ { 958 arg := float64(i) / 100 959 r := math.Pow(-50*math.Log(arg), 2.0/3.0) + 25 - x[1] 960 t1 := math.Pow(math.Abs(r), x[2]) / x[0] 961 t2 := math.Exp(-t1) 962 t := t2 - arg 963 sum += t * t 964 } 965 return sum 966 } 967 968 func (GulfResearchAndDevelopment) Grad(grad, x []float64) { 969 if len(x) != 3 { 970 panic("dimension of the problem must be 3") 971 } 972 if len(x) != len(grad) { 973 panic("incorrect size of the gradient") 974 } 975 976 for i := range grad { 977 grad[i] = 0 978 } 979 for i := 1; i <= 99; i++ { 980 arg := float64(i) / 100 981 r := math.Pow(-50*math.Log(arg), 2.0/3.0) + 25 - x[1] 982 t1 := math.Pow(math.Abs(r), x[2]) / x[0] 983 t2 := math.Exp(-t1) 984 t := t2 - arg 985 s1 := t1 * t2 * t 986 grad[0] += s1 987 grad[1] += s1 / r 988 grad[2] -= s1 * math.Log(math.Abs(r)) 989 } 990 grad[0] *= 2 / x[0] 991 grad[1] *= 2 * x[2] 992 grad[2] *= 2 993 } 994 995 func (GulfResearchAndDevelopment) Minima() []Minimum { 996 return []Minimum{ 997 { 998 X: []float64{50, 25, 1.5}, 999 F: 0, 1000 Global: true, 1001 }, 1002 { 1003 X: []float64{99.89529935174151, 60.61453902799833, 9.161242695144592}, 1004 F: 32.8345, 1005 Global: false, 1006 }, 1007 { 1008 X: []float64{201.662589489426, 60.61633150468155, 10.224891158488965}, 1009 F: 32.8345, 1010 Global: false, 1011 }, 1012 } 1013 } 1014 1015 // HelicalValley implements the helical valley function of Fletcher and Powell. 1016 // Function is not defined at x[0] = 0. 1017 // 1018 // Standard starting point: 1019 // [-1, 0, 0] 1020 // 1021 // References: 1022 // - Fletcher, R., Powell, M.J.D.: A rapidly convergent descent method for 1023 // minimization. Comput J 6 (1963), 163-168 1024 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 1025 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 1026 type HelicalValley struct{} 1027 1028 func (HelicalValley) Func(x []float64) float64 { 1029 if len(x) != 3 { 1030 panic("dimension of the problem must be 3") 1031 } 1032 if x[0] == 0 { 1033 panic("function not defined at x[0] = 0") 1034 } 1035 1036 theta := 0.5 * math.Atan(x[1]/x[0]) / math.Pi 1037 if x[0] < 0 { 1038 theta += 0.5 1039 } 1040 f1 := 10 * (x[2] - 10*theta) 1041 f2 := 10 * (math.Hypot(x[0], x[1]) - 1) 1042 f3 := x[2] 1043 return f1*f1 + f2*f2 + f3*f3 1044 } 1045 1046 func (HelicalValley) Grad(grad, x []float64) { 1047 if len(x) != 3 { 1048 panic("dimension of the problem must be 3") 1049 } 1050 if len(x) != len(grad) { 1051 panic("incorrect size of the gradient") 1052 } 1053 if x[0] == 0 { 1054 panic("function not defined at x[0] = 0") 1055 } 1056 1057 theta := 0.5 * math.Atan(x[1]/x[0]) / math.Pi 1058 if x[0] < 0 { 1059 theta += 0.5 1060 } 1061 h := math.Hypot(x[0], x[1]) 1062 r := 1 / h 1063 q := r * r / math.Pi 1064 s := x[2] - 10*theta 1065 grad[0] = 200 * (5*s*q*x[1] + (h-1)*r*x[0]) 1066 grad[1] = 200 * (-5*s*q*x[0] + (h-1)*r*x[1]) 1067 grad[2] = 2 * (100*s + x[2]) 1068 } 1069 1070 func (HelicalValley) Minima() []Minimum { 1071 return []Minimum{ 1072 { 1073 X: []float64{1, 0, 0}, 1074 F: 0, 1075 Global: true, 1076 }, 1077 } 1078 } 1079 1080 // Linear implements a linear function. 1081 type Linear struct{} 1082 1083 func (Linear) Func(x []float64) float64 { 1084 return floats.Sum(x) 1085 } 1086 1087 func (Linear) Grad(grad, x []float64) []float64 { 1088 if len(x) != len(grad) { 1089 panic("incorrect size of the gradient") 1090 } 1091 1092 for i := range grad { 1093 grad[i] = 1 1094 } 1095 return grad 1096 } 1097 1098 // PenaltyI implements the first penalty function by Gill, Murray and Pitfield. 1099 // 1100 // Standard starting point: 1101 // [1, ..., n] 1102 // 1103 // References: 1104 // - Gill, P.E., Murray, W., Pitfield, R.A.: The implementation of two revised 1105 // quasi-Newton algorithms for unconstrained optimization. Report NAC 11, 1106 // National Phys Lab (1972), 82-83 1107 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 1108 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 1109 type PenaltyI struct{} 1110 1111 func (PenaltyI) Func(x []float64) (sum float64) { 1112 for _, v := range x { 1113 sum += (v - 1) * (v - 1) 1114 } 1115 sum *= 1e-5 1116 1117 var s float64 1118 for _, v := range x { 1119 s += v * v 1120 } 1121 sum += (s - 0.25) * (s - 0.25) 1122 return sum 1123 } 1124 1125 func (PenaltyI) Grad(grad, x []float64) { 1126 if len(x) != len(grad) { 1127 panic("incorrect size of the gradient") 1128 } 1129 1130 s := -0.25 1131 for _, v := range x { 1132 s += v * v 1133 } 1134 for i, v := range x { 1135 grad[i] = 2 * (2*s*v + 1e-5*(v-1)) 1136 } 1137 } 1138 1139 func (PenaltyI) Minima() []Minimum { 1140 return []Minimum{ 1141 { 1142 X: []float64{0.2500074995875379, 0.2500074995875379, 0.2500074995875379, 0.2500074995875379}, 1143 F: 2.2499775008999372e-05, 1144 Global: true, 1145 }, 1146 { 1147 X: []float64{0.15812230111311634, 0.15812230111311634, 0.15812230111311634, 1148 0.15812230111311634, 0.15812230111311634, 0.15812230111311634, 1149 0.15812230111311634, 0.15812230111311634, 0.15812230111311634, 0.15812230111311634}, 1150 F: 7.087651467090369e-05, 1151 Global: true, 1152 }, 1153 } 1154 } 1155 1156 // PenaltyII implements the second penalty function by Gill, Murray and Pitfield. 1157 // 1158 // Standard starting point: 1159 // [0.5, ..., 0.5] 1160 // 1161 // References: 1162 // - Gill, P.E., Murray, W., Pitfield, R.A.: The implementation of two revised 1163 // quasi-Newton algorithms for unconstrained optimization. Report NAC 11, 1164 // National Phys Lab (1972), 82-83 1165 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 1166 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 1167 type PenaltyII struct{} 1168 1169 func (PenaltyII) Func(x []float64) (sum float64) { 1170 dim := len(x) 1171 s := -1.0 1172 for i, v := range x { 1173 s += float64(dim-i) * v * v 1174 } 1175 for i := 1; i < dim; i++ { 1176 yi := math.Exp(float64(i+1)/10) + math.Exp(float64(i)/10) 1177 f := math.Exp(x[i]/10) + math.Exp(x[i-1]/10) - yi 1178 sum += f * f 1179 } 1180 for i := 1; i < dim; i++ { 1181 f := math.Exp(x[i]/10) - math.Exp(-1.0/10) 1182 sum += f * f 1183 } 1184 sum *= 1e-5 1185 sum += (x[0] - 0.2) * (x[0] - 0.2) 1186 sum += s * s 1187 return sum 1188 } 1189 1190 func (PenaltyII) Grad(grad, x []float64) { 1191 if len(x) != len(grad) { 1192 panic("incorrect size of the gradient") 1193 } 1194 1195 dim := len(x) 1196 s := -1.0 1197 for i, v := range x { 1198 s += float64(dim-i) * v * v 1199 } 1200 for i, v := range x { 1201 grad[i] = 4 * s * float64(dim-i) * v 1202 } 1203 for i := 1; i < dim; i++ { 1204 yi := math.Exp(float64(i+1)/10) + math.Exp(float64(i)/10) 1205 f := math.Exp(x[i]/10) + math.Exp(x[i-1]/10) - yi 1206 grad[i] += 1e-5 * f * math.Exp(x[i]/10) / 5 1207 grad[i-1] += 1e-5 * f * math.Exp(x[i-1]/10) / 5 1208 } 1209 for i := 1; i < dim; i++ { 1210 f := math.Exp(x[i]/10) - math.Exp(-1.0/10) 1211 grad[i] += 1e-5 * f * math.Exp(x[i]/10) / 5 1212 } 1213 grad[0] += 2 * (x[0] - 0.2) 1214 } 1215 1216 func (PenaltyII) Minima() []Minimum { 1217 return []Minimum{ 1218 { 1219 X: []float64{0.19999933335, 0.19131670128566283, 0.4801014860897, 0.5188454026659}, 1220 F: 9.376293007355449e-06, 1221 Global: true, 1222 }, 1223 { 1224 X: []float64{0.19998360520892217, 0.010350644318663525, 1225 0.01960493546891094, 0.03208906550305253, 0.04993267593895693, 1226 0.07651399534454084, 0.11862407118600789, 0.1921448731780023, 1227 0.3473205862372022, 0.36916437893066273}, 1228 F: 0.00029366053745674594, 1229 Global: true, 1230 }, 1231 } 1232 } 1233 1234 // PowellBadlyScaled implements the Powell's badly scaled function. 1235 // The function is very flat near the minimum. A satisfactory solution is one 1236 // that gives f(x) ≅ 1e-13. 1237 // 1238 // Standard starting point: 1239 // [0, 1] 1240 // 1241 // References: 1242 // - Powell, M.J.D.: A Hybrid Method for Nonlinear Equations. Numerical 1243 // Methods for Nonlinear Algebraic Equations, P. Rabinowitz (ed.), Gordon 1244 // and Breach (1970) 1245 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 1246 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 1247 type PowellBadlyScaled struct{} 1248 1249 func (PowellBadlyScaled) Func(x []float64) float64 { 1250 if len(x) != 2 { 1251 panic("dimension of the problem must be 2") 1252 } 1253 1254 f1 := 1e4*x[0]*x[1] - 1 1255 f2 := math.Exp(-x[0]) + math.Exp(-x[1]) - 1.0001 1256 return f1*f1 + f2*f2 1257 } 1258 1259 func (PowellBadlyScaled) Grad(grad, x []float64) { 1260 if len(x) != 2 { 1261 panic("dimension of the problem must be 2") 1262 } 1263 if len(x) != len(grad) { 1264 panic("incorrect size of the gradient") 1265 } 1266 1267 f1 := 1e4*x[0]*x[1] - 1 1268 f2 := math.Exp(-x[0]) + math.Exp(-x[1]) - 1.0001 1269 grad[0] = 2 * (1e4*f1*x[1] - f2*math.Exp(-x[0])) 1270 grad[1] = 2 * (1e4*f1*x[0] - f2*math.Exp(-x[1])) 1271 } 1272 1273 func (PowellBadlyScaled) Hess(dst *mat.SymDense, x []float64) { 1274 if len(x) != 2 { 1275 panic("dimension of the problem must be 2") 1276 } 1277 if len(x) != dst.SymmetricDim() { 1278 panic("incorrect size of the Hessian") 1279 } 1280 1281 t1 := 1e4*x[0]*x[1] - 1 1282 s1 := math.Exp(-x[0]) 1283 s2 := math.Exp(-x[1]) 1284 t2 := s1 + s2 - 1.0001 1285 1286 h00 := 2 * (1e8*x[1]*x[1] + s1*(s1+t2)) 1287 h01 := 2 * (1e4*(1+2*t1) + s1*s2) 1288 h11 := 2 * (1e8*x[0]*x[0] + s2*(s2+t2)) 1289 dst.SetSym(0, 0, h00) 1290 dst.SetSym(0, 1, h01) 1291 dst.SetSym(1, 1, h11) 1292 } 1293 1294 func (PowellBadlyScaled) Minima() []Minimum { 1295 return []Minimum{ 1296 { 1297 X: []float64{1.0981593296997149e-05, 9.106146739867375}, 1298 F: 0, 1299 Global: true, 1300 }, 1301 } 1302 } 1303 1304 // Trigonometric implements the trigonometric function. 1305 // 1306 // Standard starting point: 1307 // [1/dim, ..., 1/dim] 1308 // 1309 // References: 1310 // - Spedicato E.: Computational experience with quasi-Newton algorithms for 1311 // minimization problems of moderatly large size. Towards Global 1312 // Optimization 2 (1978), 209-219 1313 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 1314 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 1315 type Trigonometric struct{} 1316 1317 func (Trigonometric) Func(x []float64) (sum float64) { 1318 var s1 float64 1319 for _, v := range x { 1320 s1 += math.Cos(v) 1321 } 1322 for i, v := range x { 1323 f := float64(len(x)+i+1) - float64(i+1)*math.Cos(v) - math.Sin(v) - s1 1324 sum += f * f 1325 } 1326 return sum 1327 } 1328 1329 func (Trigonometric) Grad(grad, x []float64) { 1330 if len(x) != len(grad) { 1331 panic("incorrect size of the gradient") 1332 } 1333 1334 var s1 float64 1335 for _, v := range x { 1336 s1 += math.Cos(v) 1337 } 1338 1339 var s2 float64 1340 for i, v := range x { 1341 f := float64(len(x)+i+1) - float64(i+1)*math.Cos(v) - math.Sin(v) - s1 1342 s2 += f 1343 grad[i] = 2 * f * (float64(i+1)*math.Sin(v) - math.Cos(v)) 1344 } 1345 1346 for i, v := range x { 1347 grad[i] += 2 * s2 * math.Sin(v) 1348 } 1349 } 1350 1351 func (Trigonometric) Minima() []Minimum { 1352 return []Minimum{ 1353 { 1354 X: []float64{0.04296456438227447, 0.043976287478192246, 1355 0.045093397949095684, 0.04633891624617569, 0.047744381782831, 1356 0.04935473251330618, 0.05123734850076505, 0.19520946391410446, 1357 0.1649776652761741, 0.06014857783799575}, 1358 F: 0, 1359 Global: true, 1360 }, 1361 { 1362 // TODO(vladimir-ch): If we knew the location of this minimum more 1363 // accurately, we could decrease defaultGradTol. 1364 X: []float64{0.05515090434047145, 0.05684061730812344, 1365 0.05876400231100774, 0.060990608903034337, 0.06362621381044778, 1366 0.06684318087364617, 0.2081615177172172, 0.16436309604419047, 1367 0.08500689695564931, 0.09143145386293675}, 1368 F: 2.795056121876575e-05, 1369 Global: false, 1370 }, 1371 } 1372 } 1373 1374 // VariablyDimensioned implements a variably dimensioned function. 1375 // 1376 // Standard starting point: 1377 // [..., (dim-i)/dim, ...], i=1,...,dim 1378 // 1379 // References: 1380 // More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained optimization 1381 // software. ACM Trans Math Softw 7 (1981), 17-41 1382 type VariablyDimensioned struct{} 1383 1384 func (VariablyDimensioned) Func(x []float64) (sum float64) { 1385 for _, v := range x { 1386 t := v - 1 1387 sum += t * t 1388 } 1389 1390 var s float64 1391 for i, v := range x { 1392 s += float64(i+1) * (v - 1) 1393 } 1394 s *= s 1395 sum += s 1396 s *= s 1397 sum += s 1398 return sum 1399 } 1400 1401 func (VariablyDimensioned) Grad(grad, x []float64) { 1402 if len(x) != len(grad) { 1403 panic("incorrect size of the gradient") 1404 } 1405 1406 var s float64 1407 for i, v := range x { 1408 s += float64(i+1) * (v - 1) 1409 } 1410 for i, v := range x { 1411 grad[i] = 2 * (v - 1 + s*float64(i+1)*(1+2*s*s)) 1412 } 1413 } 1414 1415 func (VariablyDimensioned) Minima() []Minimum { 1416 return []Minimum{ 1417 { 1418 X: []float64{1, 1}, 1419 F: 0, 1420 Global: true, 1421 }, 1422 { 1423 X: []float64{1, 1, 1}, 1424 F: 0, 1425 Global: true, 1426 }, 1427 { 1428 X: []float64{1, 1, 1, 1}, 1429 F: 0, 1430 Global: true, 1431 }, 1432 { 1433 X: []float64{1, 1, 1, 1, 1}, 1434 F: 0, 1435 Global: true, 1436 }, 1437 { 1438 X: []float64{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, 1439 F: 0, 1440 Global: true, 1441 }, 1442 } 1443 } 1444 1445 // Watson implements the Watson's function. 1446 // Dimension of the problem should be 2 <= dim <= 31. For dim == 9, the problem 1447 // of minimizing the function is very ill conditioned. 1448 // 1449 // Standard starting point: 1450 // [0, ..., 0] 1451 // 1452 // References: 1453 // - Kowalik, J.S., Osborne, M.R.: Methods for Unconstrained Optimization 1454 // Problems. Elsevier North-Holland, New York, 1968 1455 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 1456 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 1457 type Watson struct{} 1458 1459 func (Watson) Func(x []float64) (sum float64) { 1460 for i := 1; i <= 29; i++ { 1461 d1 := float64(i) / 29 1462 1463 d2 := 1.0 1464 var s1 float64 1465 for j := 1; j < len(x); j++ { 1466 s1 += float64(j) * d2 * x[j] 1467 d2 *= d1 1468 } 1469 1470 d2 = 1.0 1471 var s2 float64 1472 for _, v := range x { 1473 s2 += d2 * v 1474 d2 *= d1 1475 } 1476 1477 t := s1 - s2*s2 - 1 1478 sum += t * t 1479 } 1480 t := x[1] - x[0]*x[0] - 1 1481 sum += x[0]*x[0] + t*t 1482 return sum 1483 } 1484 1485 func (Watson) Grad(grad, x []float64) { 1486 if len(x) != len(grad) { 1487 panic("incorrect size of the gradient") 1488 } 1489 1490 for i := range grad { 1491 grad[i] = 0 1492 } 1493 for i := 1; i <= 29; i++ { 1494 d1 := float64(i) / 29 1495 1496 d2 := 1.0 1497 var s1 float64 1498 for j := 1; j < len(x); j++ { 1499 s1 += float64(j) * d2 * x[j] 1500 d2 *= d1 1501 } 1502 1503 d2 = 1.0 1504 var s2 float64 1505 for _, v := range x { 1506 s2 += d2 * v 1507 d2 *= d1 1508 } 1509 1510 t := s1 - s2*s2 - 1 1511 s3 := 2 * d1 * s2 1512 d2 = 2 / d1 1513 for j := range x { 1514 grad[j] += d2 * (float64(j) - s3) * t 1515 d2 *= d1 1516 } 1517 } 1518 t := x[1] - x[0]*x[0] - 1 1519 grad[0] += x[0] * (2 - 4*t) 1520 grad[1] += 2 * t 1521 } 1522 1523 func (Watson) Hess(dst *mat.SymDense, x []float64) { 1524 dim := len(x) 1525 if len(x) != dst.SymmetricDim() { 1526 panic("incorrect size of the Hessian") 1527 } 1528 1529 for j := 0; j < dim; j++ { 1530 for k := j; k < dim; k++ { 1531 dst.SetSym(j, k, 0) 1532 } 1533 } 1534 for i := 1; i <= 29; i++ { 1535 d1 := float64(i) / 29 1536 d2 := 1.0 1537 var s1 float64 1538 for j := 1; j < dim; j++ { 1539 s1 += float64(j) * d2 * x[j] 1540 d2 *= d1 1541 } 1542 1543 d2 = 1.0 1544 var s2 float64 1545 for _, v := range x { 1546 s2 += d2 * v 1547 d2 *= d1 1548 } 1549 1550 t := s1 - s2*s2 - 1 1551 s3 := 2 * d1 * s2 1552 d2 = 2 / d1 1553 th := 2 * d1 * d1 * t 1554 for j := 0; j < dim; j++ { 1555 v := float64(j) - s3 1556 d3 := 1 / d1 1557 for k := 0; k <= j; k++ { 1558 dst.SetSym(k, j, dst.At(k, j)+d2*d3*(v*(float64(k)-s3)-th)) 1559 d3 *= d1 1560 } 1561 d2 *= d1 1562 } 1563 } 1564 t1 := x[1] - x[0]*x[0] - 1 1565 dst.SetSym(0, 0, dst.At(0, 0)+8*x[0]*x[0]+2-4*t1) 1566 dst.SetSym(0, 1, dst.At(0, 1)-4*x[0]) 1567 dst.SetSym(1, 1, dst.At(1, 1)+2) 1568 } 1569 1570 func (Watson) Minima() []Minimum { 1571 return []Minimum{ 1572 { 1573 X: []float64{-0.01572508644590686, 1.012434869244884, -0.23299162372002916, 1574 1.2604300800978554, -1.51372891341701, 0.9929964286340117}, 1575 F: 0.0022876700535523838, 1576 Global: true, 1577 }, 1578 { 1579 X: []float64{-1.5307036521992127e-05, 0.9997897039319495, 0.01476396369355022, 1580 0.14634232829939883, 1.0008211030046426, -2.617731140519101, 4.104403164479245, 1581 -3.1436122785568514, 1.0526264080103074}, 1582 F: 1.399760138096796e-06, 1583 Global: true, 1584 }, 1585 // TODO(vladimir-ch): More, Garbow, Hillstrom list just the value, but 1586 // not the location. Our minimizers find a minimum, but the value is 1587 // different. 1588 // { 1589 // // For dim == 12 1590 // F: 4.72238e-10, 1591 // Global: true, 1592 // }, 1593 // TODO(vladimir-ch): netlib/uncon report a value of 2.48631d-20 for dim == 20. 1594 } 1595 } 1596 1597 // Wood implements the Wood's function. 1598 // 1599 // Standard starting point: 1600 // [-3, -1, -3, -1] 1601 // 1602 // References: 1603 // - Colville, A.R.: A comparative study of nonlinear programming codes. 1604 // Report 320-2949, IBM New York Scientific Center (1968) 1605 // - More, J., Garbow, B.S., Hillstrom, K.E.: Testing unconstrained 1606 // optimization software. ACM Trans Math Softw 7 (1981), 17-41 1607 type Wood struct{} 1608 1609 func (Wood) Func(x []float64) (sum float64) { 1610 if len(x) != 4 { 1611 panic("dimension of the problem must be 4") 1612 } 1613 1614 f1 := x[1] - x[0]*x[0] 1615 f2 := 1 - x[0] 1616 f3 := x[3] - x[2]*x[2] 1617 f4 := 1 - x[2] 1618 f5 := x[1] + x[3] - 2 1619 f6 := x[1] - x[3] 1620 return 100*f1*f1 + f2*f2 + 90*f3*f3 + f4*f4 + 10*f5*f5 + 0.1*f6*f6 1621 } 1622 1623 func (Wood) Grad(grad, x []float64) { 1624 if len(x) != 4 { 1625 panic("dimension of the problem must be 4") 1626 } 1627 if len(x) != len(grad) { 1628 panic("incorrect size of the gradient") 1629 } 1630 1631 f1 := x[1] - x[0]*x[0] 1632 f2 := 1 - x[0] 1633 f3 := x[3] - x[2]*x[2] 1634 f4 := 1 - x[2] 1635 f5 := x[1] + x[3] - 2 1636 f6 := x[1] - x[3] 1637 grad[0] = -2 * (200*f1*x[0] + f2) 1638 grad[1] = 2 * (100*f1 + 10*f5 + 0.1*f6) 1639 grad[2] = -2 * (180*f3*x[2] + f4) 1640 grad[3] = 2 * (90*f3 + 10*f5 - 0.1*f6) 1641 } 1642 1643 func (Wood) Hess(dst *mat.SymDense, x []float64) { 1644 if len(x) != 4 { 1645 panic("dimension of the problem must be 4") 1646 } 1647 if len(x) != dst.SymmetricDim() { 1648 panic("incorrect size of the Hessian") 1649 } 1650 1651 dst.SetSym(0, 0, 400*(3*x[0]*x[0]-x[1])+2) 1652 dst.SetSym(0, 1, -400*x[0]) 1653 dst.SetSym(1, 1, 220.2) 1654 dst.SetSym(0, 2, 0) 1655 dst.SetSym(1, 2, 0) 1656 dst.SetSym(2, 2, 360*(3*x[2]*x[2]-x[3])+2) 1657 dst.SetSym(0, 3, 0) 1658 dst.SetSym(1, 3, 19.8) 1659 dst.SetSym(2, 3, -360*x[2]) 1660 dst.SetSym(3, 3, 200.2) 1661 } 1662 1663 func (Wood) Minima() []Minimum { 1664 return []Minimum{ 1665 { 1666 X: []float64{1, 1, 1, 1}, 1667 F: 0, 1668 Global: true, 1669 }, 1670 } 1671 } 1672 1673 // ConcaveRight implements an univariate function that is concave to the right 1674 // of the minimizer which is located at x=sqrt(2). 1675 // 1676 // References: 1677 // More, J.J., and Thuente, D.J.: Line Search Algorithms with Guaranteed Sufficient Decrease. 1678 // ACM Transactions on Mathematical Software 20(3) (1994), 286–307, eq. (5.1) 1679 type ConcaveRight struct{} 1680 1681 func (ConcaveRight) Func(x []float64) float64 { 1682 if len(x) != 1 { 1683 panic("dimension of the problem must be 1") 1684 } 1685 return -x[0] / (x[0]*x[0] + 2) 1686 } 1687 1688 func (ConcaveRight) Grad(grad, x []float64) { 1689 if len(x) != 1 { 1690 panic("dimension of the problem must be 1") 1691 } 1692 if len(x) != len(grad) { 1693 panic("incorrect size of the gradient") 1694 } 1695 xSqr := x[0] * x[0] 1696 grad[0] = (xSqr - 2) / (xSqr + 2) / (xSqr + 2) 1697 } 1698 1699 // ConcaveLeft implements an univariate function that is concave to the left of 1700 // the minimizer which is located at x=399/250=1.596. 1701 // 1702 // References: 1703 // More, J.J., and Thuente, D.J.: Line Search Algorithms with Guaranteed Sufficient Decrease. 1704 // ACM Transactions on Mathematical Software 20(3) (1994), 286–307, eq. (5.2) 1705 type ConcaveLeft struct{} 1706 1707 func (ConcaveLeft) Func(x []float64) float64 { 1708 if len(x) != 1 { 1709 panic("dimension of the problem must be 1") 1710 } 1711 return math.Pow(x[0]+0.004, 4) * (x[0] - 1.996) 1712 } 1713 1714 func (ConcaveLeft) Grad(grad, x []float64) { 1715 if len(x) != 1 { 1716 panic("dimension of the problem must be 1") 1717 } 1718 if len(x) != len(grad) { 1719 panic("incorrect size of the gradient") 1720 } 1721 grad[0] = math.Pow(x[0]+0.004, 3) * (5*x[0] - 7.98) 1722 } 1723 1724 // Plassmann implements an univariate oscillatory function where the value of L 1725 // controls the number of oscillations. The value of Beta controls the size of 1726 // the derivative at zero and the size of the interval where the strong Wolfe 1727 // conditions can hold. For small values of Beta this function represents a 1728 // difficult test problem for linesearchers also because the information based 1729 // on the derivative is unreliable due to the oscillations. 1730 // 1731 // References: 1732 // More, J.J., and Thuente, D.J.: Line Search Algorithms with Guaranteed Sufficient Decrease. 1733 // ACM Transactions on Mathematical Software 20(3) (1994), 286–307, eq. (5.3) 1734 type Plassmann struct { 1735 L float64 // Number of oscillations for |x-1| ≥ Beta. 1736 Beta float64 // Size of the derivative at zero, f'(0) = -Beta. 1737 } 1738 1739 func (f Plassmann) Func(x []float64) float64 { 1740 if len(x) != 1 { 1741 panic("dimension of the problem must be 1") 1742 } 1743 a := x[0] 1744 b := f.Beta 1745 l := f.L 1746 r := 2 * (1 - b) / l / math.Pi * math.Sin(l*math.Pi/2*a) 1747 switch { 1748 case a <= 1-b: 1749 r += 1 - a 1750 case 1-b < a && a <= 1+b: 1751 r += 0.5 * ((a-1)*(a-1)/b + b) 1752 default: // a > 1+b 1753 r += a - 1 1754 } 1755 return r 1756 } 1757 1758 func (f Plassmann) Grad(grad, x []float64) { 1759 if len(x) != 1 { 1760 panic("dimension of the problem must be 1") 1761 } 1762 if len(x) != len(grad) { 1763 panic("incorrect size of the gradient") 1764 } 1765 a := x[0] 1766 b := f.Beta 1767 l := f.L 1768 grad[0] = (1 - b) * math.Cos(l*math.Pi/2*a) 1769 switch { 1770 case a <= 1-b: 1771 grad[0]-- 1772 case 1-b < a && a <= 1+b: 1773 grad[0] += (a - 1) / b 1774 default: // a > 1+b 1775 grad[0]++ 1776 } 1777 } 1778 1779 // YanaiOzawaKaneko is an univariate convex function where the values of Beta1 1780 // and Beta2 control the curvature around the minimum. Far away from the 1781 // minimum the function approximates an absolute value function. Near the 1782 // minimum, the function can either be sharply curved or flat, controlled by 1783 // the parameter values. 1784 // 1785 // References: 1786 // - More, J.J., and Thuente, D.J.: Line Search Algorithms with Guaranteed Sufficient Decrease. 1787 // ACM Transactions on Mathematical Software 20(3) (1994), 286–307, eq. (5.4) 1788 // - Yanai, H., Ozawa, M., and Kaneko, S.: Interpolation methods in one dimensional 1789 // optimization. Computing 27 (1981), 155–163 1790 type YanaiOzawaKaneko struct { 1791 Beta1 float64 1792 Beta2 float64 1793 } 1794 1795 func (f YanaiOzawaKaneko) Func(x []float64) float64 { 1796 if len(x) != 1 { 1797 panic("dimension of the problem must be 1") 1798 } 1799 a := x[0] 1800 b1 := f.Beta1 1801 b2 := f.Beta2 1802 g1 := math.Sqrt(1+b1*b1) - b1 1803 g2 := math.Sqrt(1+b2*b2) - b2 1804 return g1*math.Sqrt((a-1)*(a-1)+b2*b2) + g2*math.Sqrt(a*a+b1*b1) 1805 } 1806 1807 func (f YanaiOzawaKaneko) Grad(grad, x []float64) { 1808 if len(x) != 1 { 1809 panic("dimension of the problem must be 1") 1810 } 1811 if len(x) != len(grad) { 1812 panic("incorrect size of the gradient") 1813 } 1814 a := x[0] 1815 b1 := f.Beta1 1816 b2 := f.Beta2 1817 g1 := math.Sqrt(1+b1*b1) - b1 1818 g2 := math.Sqrt(1+b2*b2) - b2 1819 grad[0] = g1*(a-1)/math.Sqrt(b2*b2+(a-1)*(a-1)) + g2*a/math.Sqrt(b1*b1+a*a) 1820 }