//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** @author John Miller * @version 1.3 * @date Wed Aug 24 19:53:22 EDT 2011 * @see LICENSE (MIT style license file). */ package scalation.maxima import scala.math.{abs, max, pow, sqrt} import scala.util.control.Breaks.{breakable, break} import scalation.calculus.Differential.{gradient, gradientD} import scalation.linalgebra.VectorD import scalation.math.FunctionS2S //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `GoldenSectionLS` class performs a line search on f(x) to find a maximal * value for 'f'. It requires no derivatives and only one functional evaluation * per iteration. A search is conducted from 'x1' (often 0) to 'xmax'. A guess * for 'xmax' must be given, but can be made larger during the expansion phase, * that occurs before the recursive golden section search is called. It works on * scalar functions (see `GoldenSectionLSTest`). If starting with a vector function * 'f(x)', simply define a new function 'g(y) = x0 + direction * y' (see `GoldenSectionLSTest2`). * @param f the scalar objective function to maximize */ class GoldenSectionLS (f: FunctionS2S) { private val DEBUG = false // debug flag private val EPSILON = 1E-7 // number close to zero private val MAX_ITER = 10 // maximum number of expansion iterations private val G_RATIO = (1.0 + sqrt (5.0)) / 2.0 // the golden ratio (1.618033988749895) private val G_SECTION = G_RATIO / (1.0 + G_RATIO) // the golden section number (0.6180339887498949) //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** A recursive golden section search requiring only one functional evaluation * per call. It works by comparing two center points x2 (given) and x4 computed. * @param left whether to search left (true) or right (false) side of last interval * @param x1 the left-most point * @param x2 the center point (.618 across for left and .382 across for right) * @param x3 the right-most point * @param f2 the functional value for the x2 center point */ def gsection (left: Boolean, x1: Double, x2: Double, x3: Double, f2: Double): Double = { var x4 = x2 var f4 = 0.0 if (DEBUG) println ("gsection: left = " + left + ", x1 = " + x1 + ", x2 = " + x2 + ", x3 = " + x3 + ", f2 = " + f2) val dist = x3 - x1 if (dist < EPSILON) return (x1 + x3) / 2.0 // mid point if (left) { x4 = x3 - G_SECTION * dist // search left: x1 < x4 < x2 < x3 f4 = f(x4) if (f4 > f2) gsection (true, x1, x4, x2, f4) else gsection (false, x4, x2, x3, f2) } else { x4 = x1 + G_SECTION * dist // search right: x1 < x2 < x4 < x3 f4 = f(x4) if (f2 > f4) gsection (true, x1, x2, x4, f2) else gsection (false, x2, x4, x3, f4) } // if } // gsection //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Perform a Line Search (LS) using the Golden Search Algorithm. Two * phases are used: an expansion phase (moving the end-point) to find a * down-up pattern, followed by a traditional golden section search. * @param xmax a rough guess for the right end-point of the line search * @param x1 the left (smallest) anchor point for the search (usually 0) */ def search (xmax: Double = 2.0, x1: Double = 0.0): Double = { val f1 = f(x1) var x2 = x1 var f2 = 0.0 var x3 = xmax var f3 = 0.0 var matched = false breakable { for (k <- 1 to MAX_ITER) { // expand right to try to find a down-up pattern val dist = x3 - x1 x2 = x1 + G_SECTION * dist f2 = f(x2) f3 = f(x3) if (DEBUG) println ("search: x1 = " + x1 + ", f1 = " + f1 + ", x2 = " + x2 + ", f2 = " + f2 + ", x3 = " + x3 + ", f3 = " + f3) if (f1 < f2 && f2 > f3) { matched = true; break } // found up-down pattern, e.g., 10, 30, 20 x2 = x3 x3 = x1 + G_RATIO * dist // increase upper bound }} // for if (matched) { gsection (true, x1, x2, x3, f2) // apply golden section search on expanded interval } else { x2 = x1 + G_SECTION * (xmax - x1) f2 = f(x2) gsection (true, x1, x2, xmax, f2) // apply golden section search on original interval } // if } // search //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Print the golden ratio and the golden section. */ def printGolden () { println ("GOLDEN_RATIO = " + G_RATIO) println ("GOLDEN_SECTION = " + G_SECTION) } // printGolden } // GoldenSectionLS class //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `GoldenSectionLSTest` object is used to test the `GoldenSectionLS` class * on scalar functions. */ object GoldenSectionLSTest extends App { // def f (x: Double): Double = 10.0 - (x - 4.0) * (x - 4.0) // no expansion phase def f (x: Double): Double = 10.0 - (x - 40.0) * (x - 40.0) // requires expansion phase val solver = new GoldenSectionLS (f) println ("\nProblem 1: 10 - (x - 4)^2") println ("optimal solution = " + solver.search (10.0)) } // GoldenSectionLSTest object //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `GoldenSectionLSTest2` object is used to test the `GoldenSectionLS` class * on vector functions. */ object GoldenSectionLSTest2 extends App { val zo = VectorD (0.0, 0.0) // zero vector, the origin val dir = VectorD (1.0, 1.0) // direction to search in val ymax = 5.0 var y = 0.0 var x = zo def f (x: VectorD): Double = 10.0 - (x(0) - 2.0) * (x(0) - 2.0) - (x(1) - 3.0) * (x(1) - 3.0) def g (y: Double): Double = f(zo + dir * y) def f2 (x: VectorD): Double = 10.0 - x(0)/4.0 + 5.0*x(0)*x(0) + pow(x(0),4) - 9.0*x(0)*x(0)*x(1) + 3.0*x(1)*x(1) - 2.0*pow(x(1),4) def g2 (y: Double): Double = f2(zo + dir * y) val solver = new GoldenSectionLS (g) val solver2 = new GoldenSectionLS (g2) println ("\nProblem 1: 10 - (x_0 - 2)^2 + (x_1 - 3)^2") y = solver.search (ymax) println ("optimal y solution = " + y) x = zo + dir * y println ("optimal x solution = " + x) println ("optimal f solution = " + f(x)) println ("\nProblem 4: 10 - x_0/4 + 5x_0^2 + x_0^4 - 9x_0^2 x_1 + 3x_1^2 - 2x_1^4") // @see http://math.fullerton.edu/mathews/n2003/gradientsearch/GradientSearchMod/Links/GradientSearchMod_lnk_5.html y = solver2.search (ymax) println ("optimal y solution = " + y) x = zo + dir * y println ("optimal x solution = " + x) println ("optimal f solution = " + f(x)) } // GoldenSectionLSTest2 object