fibsearch.mws
Fibonacci Search in Optimization of Unimodal Functions
by Dr. William P. Fox, Department of Mathematics
Francis Marion University, Florence, SC 29501, wfox@fmarion.edu
(c) 2002 William P. Fox
This program performs the Fibonacci Line Search algorithm to find the maximum of a unimodal function,
f(x)
, over an interval,
a
<
x
<
b
. The program calculates the number of iterations required to insure the final interval is within the
user-specified
tolerance. This is found by solving for the smallest value of n that makes this inequality true: Fn >(b-a)/tolerance, where n is the Fibonacci number from the sequence {F0,F1, F2, ...}. The Fibonacci search concept involves placing two experiments between [a,b] using the ratios of Fibonacci numbers. (The limit of the ratio of Fibonacci numbers is the golden section 0.618 but the Fibonacci method converges quicker.) One experiment is placed at position,
, and the other at position,
. The function, to be maximized, is evaluated at these two points and the functional values are compared. We want to keep the larger functional value (in our maximization problem) and its corresponding opposite end -point. At the end of the required iterations, the final interval is the answer. At times when the final answer must be a single point and not an interval, the convention of selecting the midpoint is provided. This program works when the function is not differentiable and we are looking for a solution. If you need to minimize a function, then multiple the function by (-1) and find the maximum.
To find a maximum solution to given a function, f(x), on the interval [a, b] where the function, f(x), is unimodal.
Overview of the algorithm
INPUT
: endpoints a, b; tolerance, t, Fibonacci sequence
OUTPUT
: final interval [ai, bi], f(midpoint)
Step 1
. Initialize the tolerance, t >0.
Step 2
. Set Fn > (b-a)/t as the smallest Fn and define the test points
Step 3
. Calculate f(x1) and f(x2)
Step 4
. Compare f(x1) and f(x2)
(a.) If f(x1) < f(x2), then the new interval is [x1, b]:
a becomes the previous x1
b does not change
x1 becomes the previous x2
n = n-1
Find the new x2 using the formula in Step 2.
(b.) If f(x1) > f(x2), then the new interval is [a, x2]:
a does not change
b becomes the previous x2
x2 becomes the previous x1
n = n-1
Find the new x1 using the formula in Step 2.
Step 5
. If the length of the new interval from Step 4 is less than the tolerance specified, the stop. Otherwise go back to Step 3.
Step 6
. Estimate x* as the midpoint of the final interval and compute, f(x*), the estimated maximum of the function.
STOP
To utilize the Fibonacci routine you need the following:
User inputs are:
The user enters the function
f
using
f:=x-><enter the expression in x>
Type
FIBSearch(f, a, b, tolerance)
for specific values of
a, b,
and the
tolerance
.
The output is the iterative process (each step).
The last output provided is the midpoint of the final interval and the value of
f(x)
at that point.
Program Code
>
FIBSearch:=proc(f::procedure,a::numeric,b::numeric,T::numeric)
>
local L, M, j, C, x1, x2, N, val;
>
L:=(b-a)/T;
>
M:=round(L);
>
j:=0;
>
label_2;
>
C:=fib(j);
>
if C<L then
>
j:=j+1;
>
goto(label_2);
>
else
>
N:=j;
>
end if;
>
>
x1:=evalf(a+(fib(N-2)/fib(N))*(b-a));
>
x2:=evalf(a+(fib(N-1)/fib(N))*(b-a));
>
printf("The interval [a,b] is [% 4.2f,% 4.2f]and user specified tolerance level is% 6.5f.\n",a,b,T);
>
### WARNING: %x or %X format should be %y or %Y if used with floating point arguments
printf("The first 2 experimental endpoints are x1= % 6.3f and x2 = % 6.3f. \n",x1,x2);
>
printf(" \n");
>
printf(" \n");
>
>
printf(" Iteration x(1) x(2) f(x1) f(x2) Interval \n");
>
iterate(f,a,b,N,x1,x2);
>
>
val:=f(mdpt);
>
printf(" \n");
>
printf(" \n");
>
printf("The midpoint of the final interval is% 9.6f and f(midpoint) = % 7.3f. \n",mdpt, val);
>
printf(" \n");
>
printf(" \n");
>
### WARNING: %x or %X format should be %y or %Y if used with floating point arguments
printf("The maximum of the function is % 7.3f and the x value = % 9.6f \n",fkeep,xkeep);
>
printf(" \n");
>
printf(" \n");
>
end:
>
iterate:=proc(f::procedure,a::numeric,b::numeric, N::posint,x1::numeric,x2::numeric)
>
local x1n,x2n,an,bn,i,fx1,fx2,j,f1,f2,fmid;
>
global mdpt,fkeep,xkeep,fib;
>
i:=1;
>
x1n(1):=x1;
>
x2n(1):=x2;
>
an(1):=a;
>
bn(1):=b;
>
i:=1;
>
for j from 1 to N do
>
fx1(i):=evalf(f(x1n(i)));
>
fx2(i):=evalf(f(x2n(i)));
>
if fx1(i)<=fx2(i) then
>
an(i+1):=x1n(i);
>
bn(i+1):=bn(i);
>
x1n(i+1):=x2n(i);
>
x2n(i+1):=an(i+1)+(fib(N-i-1)/fib(N-i))*(bn(i+1)-an(i+1));
>
else
>
an(i+1):=an(i);
>
bn(i+1):=x2n(i);
>
x2n(i+1):=x1n(i);
>
x1n(i+1):=an(i+1)+(fib(N-i-2)/fib(N-i))*(bn(i+1)-an(i+1));
>
fi;
>
i:=i+1;
>
>
printf(" % 3.0f % 11.4f % 10.4f % 10.4f %10.4f [% 6.4f, % 6.4f]\n",i,x1n(i),x2n(i),fx1(i-1),fx2(i-1),an(i),bn(i));
>
mdpt := (an(i) + bn(i))/2;
>
>
>
if ((i+2)=N) then
>
if evalf(f(an(i))) > evalf(f(bn(i))) or evalf(f(an(i))) > evalf(f(mdpt)) then
>
>
fkeep := f(an(i)); xkeep := an(i);
>
else
>
if evalf(f(bn(i))) > evalf(f(mdpt)) then
>
fkeep := f(bn(i)); xkeep := bn(i);
>
else
>
fkeep := f(mdpt); xkeep := mdpt;
>
end if;
>
end if;
>
end if;
>
if(N-i-2)<0 then: goto(label_3):end if;
>
end do;
>
label_3;
>
end:
>
>
fib:=proc(n::numeric)
>
option rememer;
>
fib(0):=1:fib(1):=1;
>
if n<2 then
>
n;
>
else
>
fib(n-1)+fib(n-2);
>
end if;
>
end:
Examples
Example 1
Maximize the function
f(x)=-x^2+21.6*x+3
over the interval [0,20]
.
. Via calculus the maximum is at 10.8. Let's see how the Fibonacci Search does.
>
f:= x->-x^2+21.6*x+3;
>
FIBSearch(f,0,20,.001);
The interval [a,b] is [ 0.00, 20.00]and user specified tolerance level is .00100.
The first 2 experimental endpoints are x1= 7.639 and x2 = 12.361.
Iteration x(1) x(2) f(x1) f(x2) Interval
2 12.3607 15.2786 109.6501 117.2043 [ 7.6393, 20.0000]
3 10.5573 12.3607 117.2043 99.5818 [ 7.6393, 15.2786]
4 9.4427 10.5573 119.5811 117.2043 [ 7.6393, 12.3607]
5 10.5573 11.2461 117.7978 119.5811 [ 9.4427, 12.3607]
6 10.1316 10.5573 119.5811 119.4410 [ 9.4427, 11.2461]
7 10.5573 10.8204 119.1932 119.5811 [ 10.1316, 11.2461]
8 10.8204 10.9830 119.5811 119.6396 [ 10.5573, 11.2461]
9 10.7199 10.8204 119.6396 119.6065 [ 10.5573, 10.9830]
10 10.8204 10.8825 119.6336 119.6396 [ 10.7199, 10.9830]
11 10.7820 10.8204 119.6396 119.6332 [ 10.7199, 10.8825]
12 10.7583 10.7820 119.6397 119.6396 [ 10.7199, 10.8204]
13 10.7820 10.7967 119.6383 119.6397 [ 10.7583, 10.8204]
14 10.7967 10.8057 119.6397 119.6400 [ 10.7820, 10.8204]
15 10.7911 10.7967 119.6400 119.6400 [ 10.7820, 10.8057]
16 10.7967 10.8002 119.6399 119.6400 [ 10.7911, 10.8057]
17 10.8002 10.8022 119.6400 119.6400 [ 10.7967, 10.8057]
18 10.7988 10.8002 119.6400 119.6400 [ 10.7967, 10.8022]
19 10.8002 10.8009 119.6400 119.6400 [ 10.7988, 10.8022]
20 10.7995 10.8002 119.6400 119.6400 [ 10.7988, 10.8009]
21 10.8002 10.8002 119.6400 119.6400 [ 10.7995, 10.8009]
The midpoint of the final interval is 10.800154 and f(midpoint) = 119.640.
The maximum of the function is 119.640 and the x value = 10.799805
Example 2
Maximize a function whose derivative can not be solved explictly for 0.
Maximize the function,
f(x) = x-exp(x)
over the interval [-1,3].
>
f:= x->x^2-exp(x);
>
FIBSearch(f,0,3,.1);
The interval [a,b] is [ 0.00, 3.00]and user specified tolerance level is .10000.
The first 2 experimental endpoints are x1= 1.147 and x2 = 1.853.
Iteration x(1) x(2) f(x1) f(x2) Interval
2 1.8529 2.2941 1.1482 1.1937 [ 1.1471, 3.0000]
3 2.2941 2.5588 1.1937 1.2027 [ 1.8529, 3.0000]
4 2.5588 2.7353 1.2027 1.2036 [ 2.2941, 3.0000]
5 2.4706 2.5588 1.2036 1.2028 [ 2.2941, 2.7353]
6 2.3824 2.4706 1.2036 1.2036 [ 2.2941, 2.5588]
7 2.4706 2.4706 1.2033 1.2036 [ 2.3824, 2.5588]
The midpoint of the final interval is 2.470588 and f(midpoint) = 1.204.
The maximum of the function is 1.204 and the x value = 2.558824
Example 3.
Maximizing a function that does not have a derivative.
>
f:= x->-(abs(2-x)+abs(5-4*x)+abs(8-9*x));
>
FIBSearch(f,0,3,.1);
The interval [a,b] is [ 0.00, 3.00]and user specified tolerance level is .10000.
The first 2 experimental endpoints are x1= 1.147 and x2 = 1.853.
Iteration x(1) x(2) f(x1) f(x2) Interval
2 .7059 1.1471 -3.5882 -11.2353 [ 0.0000, 1.8529]
3 1.1471 1.4118 -5.1176 -3.5882 [ .7059, 1.8529]
4 .9706 1.1471 -3.5882 -5.9412 [ .7059, 1.4118]
5 .8824 .9706 -2.8824 -3.5882 [ .7059, 1.1471]
6 .7941 .8824 -2.6471 -2.8824 [ .7059, .9706]
7 .8824 .8824 -3.8824 -2.6471 [ .7941, .9706]
The midpoint of the final interval is .882353 and f(midpoint) = -2.647.
The maximum of the function is -2.882 and the x value = .970588
Example 5.
Minimize the function,
.
To minimize
,
we maximize
.
We keep the interval the same.
>
f:= x->-2*x^2+4*x;
>
FIBSearch(f,-1,2,.1);
The interval [a,b] is [-1.00, 2.00]and user specified tolerance level is .10000.
The first 2 experimental endpoints are x1= .147 and x2 = .853.
Iteration x(1) x(2) f(x1) f(x2) Interval
2 .8529 1.2941 .5450 1.9567 [ .1471, 2.0000]
3 .5882 .8529 1.9567 1.8270 [ .1471, 1.2941]
4 .8529 1.0294 1.6609 1.9567 [ .5882, 1.2941]
5 1.0294 1.1176 1.9567 1.9983 [ .8529, 1.2941]
6 .9412 1.0294 1.9983 1.9723 [ .8529, 1.1176]
7 1.0294 1.0294 1.9931 1.9983 [ .9412, 1.1176]
The midpoint of the final interval is 1.029412 and f(midpoint) = 1.998.
The maximum of the function is 2.000 and the x value = .985294
>
f:=(x)->1-exp(-x)+1/(1+x);
>
FIBSearch(f,0,20,.1);
The interval [a,b] is [ 0.00, 20.00]and user specified tolerance level is .10000.
The first 2 experimental endpoints are x1= 7.639 and x2 = 12.361.
Iteration x(1) x(2) f(x1) f(x2) Interval
2 4.7210 7.6395 1.1153 1.0748 [ 0.0000, 12.3605]
3 2.9185 4.7210 1.1659 1.1153 [ 0.0000, 7.6395]
4 1.8026 2.9185 1.2012 1.1659 [ 0.0000, 4.7210]
5 2.9185 3.6052 1.1919 1.2012 [ 1.8026, 4.7210]
6 2.4893 2.9185 1.2012 1.1900 [ 1.8026, 3.6052]
7 2.2318 2.4893 1.2036 1.2012 [ 1.8026, 2.9185]
8 2.4893 2.6609 1.2021 1.2036 [ 2.2318, 2.9185]
9 2.4034 2.4893 1.2036 1.2033 [ 2.2318, 2.6609]
10 2.4893 2.5751 1.2034 1.2036 [ 2.4034, 2.6609]
11 2.4893 2.4893 1.2036 1.2036 [ 2.4034, 2.5751]
The midpoint of the final interval is 2.489270 and f(midpoint) = 1.204.
The maximum of the function is 1.203 and the x value = 2.403433
