Blocks are a powerful feature of Objective-C and are easy to use once you get your head around them. They are very useful in scenarios where you want to pass a function as a parameter to another function like a Golden Section Search (read the linked Wikipedia article for a detailed explanation of the algorithm). It is designed to find the minimum (or maximum) of a function, provided the region is it searching is unimodal.
Firstly we need to create a block definition that defines our one-dimensional function. Create a class called PWAlgorithm and open the header file.
typedef CGFloat(^PWBlockFunction1D)(CGFloat x);
The syntax for block definitions is... well... confusing. So confusing, in fact, that there is an entire webpage dedicated to it, appropriately called fuckingblocksyntax.com. We start with the return type, then we choose a custom block name (preceded by a caret ^) and finally list all the parameters in a C-style format. This definition says, "Blocks of type PWBlockFunction1D will take a single float x and return a float..."
Now for the Golden Section Search implementation, Objective-C style.
#define M_RESPHI 0.38196601125
+ (CGFloat)goldenSectionSearchWithFunction1D:(PWBlockFunction1D)function1D a:(CGFloat)a b:(CGFloat)b fb:(CGFloat)fb c:(CGFloat)c tau:(CGFloat)tau {
CGFloat x;
if (c - b > b - a) {
x = b + M_RESPHI * (c - b);
} else {
x = b - M_RESPHI * (b - a);
}
if(fabs(c - a) < tau * (fabs(b) + fabs(x))) {
return (c + a) / 2.0;
}
CGFloat fx = function1D(x);
if(fabs(fx-fb) < 0.0000001) {
return (c + a) / 2.0;
}
if(fx < fb) {
if (c - b > b - a) {
return [PWAlgorithm goldenSectionSearchWithFunction1D:function1D a:b b:x fb:fx c:c tau:tau];
} else {
return [PWAlgorithm goldenSectionSearchWithFunction1D:function1D a:a b:x fb:fx c:b tau:tau];
}
} else {
if (c - b > b - a) {
return [PWAlgorithm goldenSectionSearchWithFunction1D:function1D a:a b:b fb:fb c:x tau:tau];
} else {
return [PWAlgorithm goldenSectionSearchWithFunction1D:function1D a:x b:b fb:fb c:c tau:tau];
}
}
return (c + a) / 2.0;
}The first parameter, locally named function1D, is a block that can be executed simply by calling function1D(x). See how the block definition PWBlockFunction1D can be used as a parameter type, exactly like any primitive.
The other parameters are the normal Golden Section Search inputs - a, b, c and tau. The function of b, fb, is also a parameter for the sole purpose of efficiency. It caches the result of the block for the value of b, avoiding unnecessary computation.
We don't need all the parameters to kick-start the recursive algorithm, so let's write another static convenience method to act as the launch point.
+ (CGFloat)goldenSectionSearchWithFunction1D:(PWBlockFunction1D)function1D unimodalLowerX:(CGFloat)unimodalLowerX unimodalUpperX:(CGFloat)unimodalUpperX tau:(CGFloat)tau {
CGFloat midX = (unimodalLowerX + unimodalUpperX)/2.0;
return [PWAlgorithm goldenSectionSearchWithFunction1D:function1D a:unimodalLowerX b:midX fb:function1D(midX) c:unimodalUpperX tau:tau
];
}All we now need to provide is the block (function1D), the lower (a) and upper (c) bounds of the unimodal range and the tolerance (tau). Initially, the range midpoint is sufficient to act as the b parameter.
Let's test our code with a very simple quadratic equation.
PWBlockFunction1D function1D = ^CGFloat(CGFloat x) {
return (3.0*x*x - 3.0*x + 1.0);
};
CGFloat xAtMinimum = [PWAlgorithm goldenSectionSearchWithFunction1D:function1D unimodalLowerX:-55.0 unimodalUpperX:106.0 tau:0.01];
CGFloat yAtMinimum = function1D(xAtMinimum);After running the search within an arbitrary unimodal range of [-55, 106], xAtMinimum = 0.5007752153 and yAtMinimum = 0.2500018029. Looking at the graph, the minimum of this function is at (x = 0.5, y = 0.25). Success!
Remember how the Golden Section Search only works if it is initially provided with a unimodal range? Solving quadratic equations with the Golden Section Search is easy as they only have one single turning point, meaning any range containing that turning point will be unimodal. If we have a more complex function with multiple turning points, then we will need to work out our initial unimodal range before we can even start our search. There are many ways to do this, some much more reliable than others. Let's do the easiest way - uniformly sampling our function. We will look for the sample with the smallest value, and use the surrounding region as our unimodal range.
+ (CGFloat)goldenSectionSearchWithFunction1D:(PWBlockFunction1D)function1D unimodalSearchLowerX:(CGFloat)unimodalSearchLowerX unimodalSearchUpperX:(CGFloat)unimodalSearchUpperX unimodalSearchSamples:(NSUInteger)unimodalSearchSamples tau:(CGFloat)tau {
if(unimodalSearchSamples < 3) {
// 2 or less samples is not enough to do the search
// just use the provided search bounds
return [PWAlgorithm goldenSectionSearchWithFunction1D:function1D unimodalLowerX:unimodalSearchLowerX unimodalUpperX:unimodalSearchUpperX tau:tau];
}
// width of each sampled region
CGFloat sampledRegionWidth = (unimodalSearchUpperX-unimodalSearchLowerX)/(CGFloat)(unimodalSearchSamples-1);
CGFloat minimumSampledFunctionOutput = CGFLOAT_MAX;
CGFloat xAtSampledMinimum = 0.0;
CGFloat currentFunctionOutput, currentX;
for(NSInteger s = 0; s<unimodalSearchSamples; ++s) {
currentX = unimodalSearchLowerX + sampledRegionWidth*(CGFloat)s;
currentFunctionOutput = function1D(currentX);
if(currentFunctionOutput < minimumSampledFunctionOutput) {
minimumSampledFunctionOutput = currentFunctionOutput;
xAtSampledMinimum = currentX;
}
}
CGFloat unimodalLowerX = xAtSampledMinimum - sampledRegionWidth;
// ensure lower bound stays within search range
if(unimodalLowerX < unimodalSearchLowerX) {
unimodalLowerX = unimodalSearchLowerX;
}
CGFloat unimodalUpperX = xAtSampledMinimum + sampledRegionWidth;
// ensure upper bound stays within search range
if(unimodalUpperX > unimodalSearchUpperX) {
unimodalUpperX = unimodalSearchUpperX;
}
return [PWAlgorithm goldenSectionSearchWithFunction1D:function1D unimodalLowerX:unimodalLowerX unimodalUpperX:unimodalUpperX tau:tau];
}We calculate the width of the region between two adjacent samples by dividing the width of the entire search range by the number of regions. The x value of each sample is calculated by multiplying the current sample number (starting at 0 of course) with the region width, and adding that to the lower search x bounds.
We evaluate the function output at each iteration, sampling the function. We keep a record of the smallest function output so far in the variable minimumSampledFunctionOutput, initialised to the largest float possible. If the function output at the current iteration is smaller than the recorded minimum, we overwrite it and note the x value where it occurred. Our final unimodal range is the x value at the sampled minimum plus/minus the region width, bounded by the search limits.
Let's try it out on a higher order polynomial.
If we set the search range to [-8,8] and the number of samples to 17, the region width becomes 1. The polynomial will be sampled at x = -8, -7, -6 ... 0 ... 6, 7, 8.
The red crosses indicate the the function output at each x division, and the smallest is circled at -5. Our unimodal search range becomes -5 ± 1, or [-6, -4]. As you can see from the plot, the function is unimodal in this region, and it does contain the global minimum.
PWBlockFunction1D function1D = ^CGFloat(CGFloat x) {
CGFloat output = 1.0/1000.0;
output *= (x + 6.3);
output *= (x + 3.7);
output *= (x + 1.0);
output *= (x - 2.0);
output *= (x - 3.0);
output *= (x - 6.0);
return output;
};
CGFloat xAtMinimum = [PWAlgorithm goldenSectionSearchWithFunction1D:function1D unimodalSearchLowerX:-8.0 unimodalSearchUpperX:8.0 unimodalSearchSamples:17 tau:0.01];
CGFloat yAtMinimum = function1D(xAtMinimum);After running the code we have xAtMinimum = -5.4270510674 and yAtMinimum = -4.7734761238. The ever reliable wolframalpha.com tells us that the global minimum for our polynomial is at [x = -4.77348, y = -5.42716]. Great, that's good enough for now.
The method we used to determine the unimodal range is far from fool proof. You can see that if we didn't sample enough points along the curve we could easily get an incorrect result. For example, if we set the search range to [-7,8] and the number of samples to 6, the region width becomes 3. The polynomial will be sampled at x = -7, -4, -1, 2, 5 & 8.
The smallest sampled function output is at x = 5 meaning our unimodal search range becomes 5 ± 3, or [2, 8]. This region does not contain the global minimum and hence our Golden Section Search will not find it.
Our simple unimodal search algorithm is OK, but keep it's limitations in mind. It could easily fail if you do not have a fairly good idea of where you should be searching for the global minimum, or you don't know how bumpy the function is. You might even want to write a more complex unimodal search algorithm of your own.
And remember, you can maximise any function by minimising it's negative.