 |
| Golden Ratio seen in nature |
Enter the Golden Section Search, a one dimensional search tool with a number of benefits over gradient descent. First, it is a class of minimization techniques which work by reducing a bracket containing a minimum. This means we can preset the bracket based on prior information, which can be a major upside. Second, it is typically more efficient than a gradient descent search. The downside is that it is limited to a single input dimension, although it can be combined with a conjugate gradient search method, minimizing a single value in the current search direction.
 |
| Golden Section Search Diagram |
So how does this approach work? Lets start with the above graph, where we have some unknown function f(x), and three points (x1, f1), (x2, f2), (x3, f3) which we already know something about. x4 is our next candidate x value, and f4a and f4b are possible y values for this candidate. We remember from the extreme value theorem that there will be a local minimum within a bracket when there is a point inside the bracket that falls below the line connecting the bracketing points. In this case, f2 clearly falls below the line connecting f1 and f3. Now we pick our bracketing points in a way that there will always be a slightly bigger side, and our next x value will be in the side with the bigger boundary. In this case, b is bigger than a (along the x-axis in the diagram), so our next point will be to the right of x2. Now with your minds eye, you can think of a parabola passing through f1, f2, and f4a, and you can see that the minimum value is in between f4a and f1... where as there would be a maximum if you made a parabola between f2, f4a, and f3. Because of this, you would use x1, x2, and x4 as your next three bracket points. If instead we were to do the same thing with f4b, the bracket would instead work out to x2, x4, and x3. This process is recursively repeated until your bracket is small enough to achieve a desired precision.
Now some code.
public delegate double Function1D(double x);
public class GoldenSectionSearch
{
private readonly double _goldenRatio = 2.0 - (1.0 + Math.Sqrt(5.0))/2.0;
public double Threshold { get; set; }
public GoldenSectionSearch()
{
Threshold = 0.001;
}
public double Minimize(Function1D function, double lowerBound, double upperBound)
{
var f = CachedFunction(function);
var midpoint = StartingSearchPoint(lowerBound, upperBound);
return Minimize(f, lowerBound, midpoint, upperBound);
}
private double Minimize(Function1D f, double lower, double midpoint, double upper)
{
var x = NextSearchPoint(lower, midpoint, upper);
if (IsSearchDone(lower, midpoint, upper, x))
return FinalPosition(lower, upper);
if (IsMinimumInNewSearchBracket(f, midpoint, x))
{
if (IsUpperBracketLarger(lower, midpoint, upper))
return Minimize(f, midpoint, x, upper);
return Minimize(f, lower, x, upper);
}
else
{
if (IsUpperBracketLarger(lower, midpoint, upper))
return Minimize(f, lower, midpoint, x);
return Minimize(f, x, midpoint, upper);
}
}
private Function1D CachedFunction(Function1D function)
{
var cache = new Dictionary();
Function1D f = (x) =>
{
if (cache.ContainsKey(x))
{
return cache[x];
}
var result = function(x);
cache[x] = result;
return result;
};
return f;
}
private double StartingSearchPoint(double lower, double upper)
{
return (upper - lower)*_goldenRatio + lower;
}
private double NextSearchPoint(double lower, double midpoint, double upper)
{
if(IsUpperBracketLarger(lower, midpoint, upper))
return midpoint + _goldenRatio * (upper - midpoint);
return midpoint - _goldenRatio * (midpoint - lower);
}
private bool IsSearchDone(double lower, double midpoint, double upper, double x)
{
return Math.Abs(upper - lower) < Threshold*(Math.Abs(midpoint) + Math.Abs(x));
}
private double FinalPosition(double lower, double upper)
{
return (upper + lower)/2.0;
}
private bool IsUpperBracketLarger(double lower, double midpoint, double upper)
{
return (upper - midpoint > midpoint - lower);
}
private bool IsMinimumInNewSearchBracket(Function1D f, double midpoint, double x)
{
return f(x) < f(midpoint);
}
}
Some explanation here...
First we use a function delegate to force the user to provide a 1d function to minimize. It is less flexible than the gradient descent delegates, but a necessary feature of this algorithm.
Inside the public Minimize function, we wrap the input with a caching decorator, so that any value we have already evaluated is reused.
The search termination criteria in the IsSearchDone is somewhat confusing, but a recommended formula by the numerical recipes guys. The gist is we say a minimum is found when the width of the bracket is smaller than a precision that we care about.
We use the "golden ratio" to pick the size of our bracket steps. This gives us special properties when we are reducing the size of the bracket to make sure the ratios never change, and also minimizes the number of steps required to reach the minimum. It is a surprising result that it is actually more efficient than bisection.
The private Minimize helper function is the workhorse, where the recursive search happens. Each iteration, we check to see if we are done, and if not, reduce the bracket, and pick a lower, midpoint, and upper bounds for the next iteration based on the intuition described above the code.
So now a test to show it in action....
[Test]
public void GoldenSectionSearchWorks()
{
var gss = new GoldenSectionSearch { Threshold = 0.0001 };
const double lowerBounds = -10.0;
const double upperBounds = 10.0;
Function1D f = (x) => Math.Pow(x - 4.0, 2.0);
var solution = gss.Minimize(f, lowerBounds, upperBounds);
Assert.AreEqual(solution, 4.0, 0.001);
}
And we are done! A C# Golden Section Search algorithm.
