Finding the minimum of an unknown function by the least amount of tries

[NickF]

Homework Statement

I am working on a program that queries an electronic instrument in order to find a value x for which the instrument will return a minimum value. I can provide any value for x that I want and the instrument will return the result f(x). The function has only one minimum. In order to test my algorithm, I created a function and I would like to know how can find the minimum of this function using as few tries as possible. Once I am pretty happy with the algorithm, I will implement it in the real program. I thought of a simple algorithm (which I admit is pretty dumb and it may not converge for a slightly different function or for any start value of x). With my algorithm, if I start with x = 0, it will take me 35 tries to identify x for which f(x) is minimum (when x = -4.8671, f(x) = 0.9426).

Homework Equations

f(x) = 0.5*(0.29*x*(x+2)+0.22*x+3-3*sqrt(x*x-3*x+9)+12.88-x/1.33);

The Attempt at a Solution

a. choose a random start value
b. choose a random step
c. query the function what is the result yleft = f(xstart - step)
d. query the function what is the result yright = f(xstart + step)
e. query the function what is the result of y = f(xstart)
f. if yleft < yright, then x must be somwhere to the left, so xleft = x-2*step, xright = x
g. if yright < yleft, xleft = x, xright = x+2*step
h. More details in the sample program provided below.

x = 0;
xStep = 2.0;
y = Function(x);
xl = x - xStep;
xr = x + xStep;

do {
n++;
yl = Function(xl);
yr = Function(xr);
x = (xl + xr) / 2;
y = Function(x);

if(yl <= yr && yl <= y) {
Dir = Left; //minimum is probably to the left of X, but not necessarily!
x = (xl + xr) / 2;
xr = x;
xl = xl - 2 * xStep;
}
else if(yr <= yl && yr <= y) {
Dir = Right; //minimum is probably to the right of X, but not necessarily!
x = (xl + xr) / 2;
xl = x;
xr = xr + 2 * xStep;
}
else if(y <= yl && y <= yr) {
Dir = Center;
xStep = xStep * 0.70; // I found by trial and error that 0.70 is a value that converges decently fast
xl = xl + xStep;
xr = xr - xStep;
}
}
while(fabs(yr - yl) > 0.0001 && n < 100);

Certainly there are much better algorithms than this one, so any suggestions are welcome.

Kind regards,
Nicolae

 

[NascentOxygen]

Hi NickF! http://img96.imageshack.us/img96/5725/red5e5etimes5e5e45e5e25.gif

Does the function have only one dip and that's the minimum?

 

[NascentOxygen]

Iff it has only one dip, wouldn't bisection be simplest? Get the y values for xmin, xmax, and the mid x value. Then move your attention to one of those halves, repeatedly.

The number of evaluations needed will depend on the accuracy you require.

Just noticed your polynomial test function, so graphed it here http://fooplot.com/index.php?&type0=0&type1=0&type2=0&type3=0&type4=0&y0=&y1=0.5*%280.29*x*%28x%2B2%29%2B0.22*x%2B3-3*sqrt%28x*x-3*x%2B9%29%2B12.88-x/1.33%29&y2=&y3=&y4=&r0=&r1=&r2=&r3=&r4=&px0=&px1=&px2=&px3=&px4=&py0=&py1=&py2=&py3=&py4=&smin0=0&smin1=0&smin2=0&smin3=0&smin4=0&smax0=2pi&smax1=2pi&smax2=2pi&smax3=2pi&smax4=2pi&thetamin0=0&thetamin1=0&thetamin2=0&thetamin3=0&thetamin4=0&thetamax0=2pi&thetamax1=2pi&thetamax2=2pi&thetamax3=2pi&thetamax4=2pi&ipw=1&ixmin=-40&ixmax=40&iymin=0&iymax=100&igx=10&igy=10&igl=1&igs=0&iax=1&ila=1&xmin=-40&xmax=40&ymin=0&ymax=100

It has only one dip, so if this is representative, maybe try bisection.

 

[NickF]

Hi NickF! http://img96.imageshack.us/img96/5725/red5e5etimes5e5e45e5e25.gif

Does the function have only one dip and that's the minimum?

Hi, the function has only one dip. Even the function I defined in there has only one dip (despite its complexity).
I will try what you said and report the results. Thanks for help.


Nicolae

 

[NickF]

Iff it has only one dip, wouldn't bisection be simplest? Get the y values for xmin, xmax, and the mid x value. Then move your attention to one of those halves, repeatedly.

The number of evaluations needed will depend on the accuracy you require.

Just noticed your polynomial test function, I guess it has more than one dip, so forget about bisection.

Could you explain please how to select which half contains my minimum? Let's say f(xmin) = 4, f(xavg) = 3, f(xmax)= 8. I don't see how can I guess from here if the minimum is in the left half or in the right half. There can be cases when the minimum is in either half (if the right side of the function is very steep,the minimum can be in the right half).

Regards,
Nicolae

 

[NascentOxygen]

Could you explain please how to select which half contains my minimum?

How does this sound: you have two halves, if you evaluate y for the midpoints of those two regions, you now have 5 function values. The minimum of these indicates the midpoint of the new two adjacent regions to examine. (Might need to make special provision for those occasions when there are equal function values.)

Can you see any problems with this?

 

[NickF]

How does this sound: you have two halves, if you evaluate y for the midpoints of those two regions, you now have 5 function values. The minimum of these indicates the midpoint of the new two adjacent regions to examine. (Might need to make special provision for those occasions when there are equal function values.)

Can you see any problems with this?

This sounds better, I will try to implement it. Thanks.

Nicolae

 

[Ray Vickson]

Homework Statement

I am working on a program that queries an electronic instrument in order to find a value x for which the instrument will return a minimum value. I can provide any value for x that I want and the instrument will return the result f(x). The function has only one minimum. In order to test my algorithm, I created a function and I would like to know how can find the minimum of this function using as few tries as possible. Once I am pretty happy with the algorithm, I will implement it in the real program. I thought of a simple algorithm (which I admit is pretty dumb and it may not converge for a slightly different function or for any start value of x). With my algorithm, if I start with x = 0, it will take me 35 tries to identify x for which f(x) is minimum (when x = -4.8671, f(x) = 0.9426).

Homework Equations

f(x) = 0.5*(0.29*x*(x+2)+0.22*x+3-3*sqrt(x*x-3*x+9)+12.88-x/1.33);

The Attempt at a Solution

a. choose a random start value
b. choose a random step
c. query the function what is the result yleft = f(xstart - step)
d. query the function what is the result yright = f(xstart + step)
e. query the function what is the result of y = f(xstart)
f. if yleft < yright, then x must be somwhere to the left, so xleft = x-2*step, xright = x
g. if yright < yleft, xleft = x, xright = x+2*step
h. More details in the sample program provided below.

x = 0;
xStep = 2.0;
y = Function(x);
xl = x - xStep;
xr = x + xStep;

do {
n++;
yl = Function(xl);
yr = Function(xr);
x = (xl + xr) / 2;
y = Function(x);

if(yl <= yr && yl <= y) {
Dir = Left; //minimum is probably to the left of X, but not necessarily!
x = (xl + xr) / 2;
xr = x;
xl = xl - 2 * xStep;
}
else if(yr <= yl && yr <= y) {
Dir = Right; //minimum is probably to the right of X, but not necessarily!
x = (xl + xr) / 2;
xl = x;
xr = xr + 2 * xStep;
}
else if(y <= yl && y <= yr) {
Dir = Center;
xStep = xStep * 0.70; // I found by trial and error that 0.70 is a value that converges decently fast
xl = xl + xStep;
xr = xr - xStep;
}
}
while(fabs(yr - yl) > 0.0001 && n < 100);

Certainly there are much better algorithms than this one, so any suggestions are welcome.

Kind regards,
Nicolae

Google "Fibonacci Search", which is optimal in the sense of minimizing the maximum size of the "uncertainty interval" (where the min is known to be) for a given number of steps. Somewhat easier (and at most one function evaluation worse) is "Golden Section Search". See, eg. http://math.fullerton.edu/mathews/n2003/fibonaccisearchmod.html or
http://www.math.ucf.edu/~xli/Fibonacci%20Search.pdf (Fibonacci) and
http://en.wikipedia.org/wiki/Golden_section_search (Golden Section).

RGV

 

[NickF]

How does this sound: you have two halves, if you evaluate y for the midpoints of those two regions, you now have 5 function values. The minimum of these indicates the midpoint of the new two adjacent regions to examine. (Might need to make special provision for those occasions when there are equal function values.)
Can you see any problems with this?

Hi again,
I implemented the suggested algorithm and it works quite well. Probably it is enough for my needs. Thank you very much for help. I got some more replies, I might try a different method, just for fun.


Nicolae

 

[Ray Vickson]

Iff it has only one dip, wouldn't bisection be simplest? Get the y values for xmin, xmax, and the mid x value. Then move your attention to one of those halves, repeatedly.

The number of evaluations needed will depend on the accuracy you require.

Just noticed your polynomial test function, so graphed it here http://fooplot.com/index.php?&type0=0&type1=0&type2=0&type3=0&type4=0&y0=&y1=0.5*%280.29*x*%28x%2B2%29%2B0.22*x%2B3-3*sqrt%28x*x-3*x%2B9%29%2B12.88-x/1.33%29&y2=&y3=&y4=&r0=&r1=&r2=&r3=&r4=&px0=&px1=&px2=&px3=&px4=&py0=&py1=&py2=&py3=&py4=&smin0=0&smin1=0&smin2=0&smin3=0&smin4=0&smax0=2pi&smax1=2pi&smax2=2pi&smax3=2pi&smax4=2pi&thetamin0=0&thetamin1=0&thetamin2=0&thetamin3=0&thetamin4=0&thetamax0=2pi&thetamax1=2pi&thetamax2=2pi&thetamax3=2pi&thetamax4=2pi&ipw=1&ixmin=-40&ixmax=40&iymin=0&iymax=100&igx=10&igy=10&igl=1&igs=0&iax=1&ila=1&xmin=-40&xmax=40&ymin=0&ymax=100

It has only one dip, so if this is representative, maybe try bisection.

Generally, if all one has access to are function evaluations, then the bisection scheme you propose will not necessarily work. For example, if we have x1 < x2 < x3 and f(x1) > f(x2) < f(x3), we cannot figure out where the min. will be. It could be in [x1,x2] or in [x2,x3]. You _could_ do a quadratic fit (and that is one of the standard minimization algorithms), by fitting a quadratic to (x1,f(x)), (x2,f(x)) and (x3,f(x3)), then finding the min x4 of the quadratic. You could pick the best three of the four points x1, x2, x3, x4 and go on from there, but without some form of "safeguarding" the method can (and in some examples, does) fail to converge, or if converges, does so very slowly. More precisely: sometimes the quadratic fit method works very well and is really fast, but sometimes it does not work well at all.

There is a vast literature devoted to univariate function minimization using only function values; there is no good reason to re-invent the wheel in practical problem-solving.

RGV

 

[NickF]

Google "Fibonacci Search", which is optimal in the sense of minimizing the maximum size of the "uncertainty interval" (where the min is known to be) for a given number of steps. Somewhat easier (and at most one function evaluation worse) is "Golden Section Search". See, eg. http://math.fullerton.edu/mathews/n2003/fibonaccisearchmod.html or
http://www.math.ucf.edu/~xli/Fibonacci%20Search.pdf (Fibonacci) and
http://en.wikipedia.org/wiki/Golden_section_search (Golden Section).

RGV

Hi Ray,

The Fibonacci algorithm seems a very good choice, I will try to implement it as well. I am curious how many steps I can save between the algorithm suggested by NascentOxygen and Fibonacci's algorithm on my test function. I will post the comparison results.


Nicolae

 

[Ray Vickson]

Hi Ray,

The Fibonacci algorithm seems a very good choice, I will try to implement it as well. I am curious how many steps I can save between the algorithm suggested by NascentOxygen and Fibonacci's algorithm on my test function. I will post the comparison results.


Nicolae

Fibonacci has some optimal properties, provided that you already have a starting interval [a,b] that contains your minimum. As I said, Golden Section search is at most one f(x) evaluation worse for the same performance, and is a bit easier.

However, suppose you do not have a starting interval; then what? For the Golden section method there are starting routines that are quite efficient. If r = (sqrt(5)-1)/2 = 0.618034 is the Golden Section, then: start with x0 = 0 (for example), and let d be some step size. Set x1 = x0+d, x2 = x0 + (1+r)*d, x3 = x0 + (1+r)^2*d, etc. Stop as soon as you find f(x_i) < f(x_{i+1}) : take a = x_{i_1} and b = x_{i+1} as the initial uncertainty interval. (In the special case f(x0) < f(x1), take a=x0 and b = x1.) Note that if i > 1 you already have one correctly-placed point (x_i) in this interval, so you need only one more f(x) evaluation to reduce length of the interval [a,b] by factor r. Each additional f evaluation reduces the interval by factor r, so after n function evaluations in [a,b] the new uncertainty interval's length will be r^n*(b-a). For example, for n = 10 we have r^10 = 0.00813.

RGV

 

[NickF]

Fibonacci has some optimal properties, provided that you already have a starting interval [a,b] that contains your minimum. As I said, Golden Section search is at most one f(x) evaluation worse for the same performance, and is a bit easier.

However, suppose you do not have a starting interval; then what? For the Golden section method there are starting routines that are quite efficient. If r = (sqrt(5)-1)/2 = 0.618034 is the Golden Section, then: start with x0 = 0 (for example), and let d be some step size. Set x1 = x0+d, x2 = x0 + (1+r)*d, x3 = x0 + (1+r)^2*d, etc. Stop as soon as you find f(x_i) < f(x_{i+1}) : take a = x_{i_1} and b = x_{i+1} as the initial uncertainty interval. (In the special case f(x0) < f(x1), take a=x0 and b = x1.) Note that if i > 1 you already have one correctly-placed point (x_i) in this interval, so you need only one more f(x) evaluation to reduce length of the interval [a,b] by factor r. Each additional f evaluation reduces the interval by factor r, so after n function evaluations in [a,b] the new uncertainty interval's length will be r^n*(b-a). For example, for n = 10 we have r^10 = 0.00813.

RGV

In my real life case, I can easily provide a domain that contains the minimum.
If I understand correctly, the Golden Section method is possibly one step slower than Fibonacci Search method, but it has the advantage of not requiring a known starting interval that contains the minimum? If this is the case, I would prefer this method (especially if it is easier to implement - although I had a brief look to Fibonacci Search and it does not seem too difficult).

Regards,
Nicolae

 

[NickF]

Hi All,

I completed the Fibonacci Search program. Because I could not find a C/C++ sample program, I will add my work in here, maybe it will help somebody else. I used Borland C++ Builder, but it is simple enough to be ported to other platforms.

Regards,
Nicolae

void __fastcall TForm1::btnFibonacciSearchClick(TObject *Sender)
{
FibonacciSearch();
}
//---------------------------------------------------------------------------

double __fastcall TForm1::Function(double x)
{
double y;
// y = 0.5*(0.29*x*(x+2)+0.22*x+3-3*sqrt(x*x-3*x+9)+12.88-x/1.33);
// y = x*x*x*x*x -(4.0/5.0) * x + 0.5;
y = x * x - sin(x);
return y;
}
//---------------------------------------------------------------------------

void __fastcall TForm1::FibonacciSearch()
{
double y, x;
double a, b, c, d;

double eps = 1e-4;
double e = 0.01;
double Fn;
double r;
double ya, yb, yc, yd;
int n, k;

a = 0;
b = 1;
Fn = (b - a) / eps;
n = ReverseFibonacci(Fn);
for(k = 1; k <= n; k++) {
r = (double) Fibonacci(n - k - 1) / (double) Fibonacci(n - k);
c = a + (1 - r) * (b - a);
d = a + r * (b - a);
yc = Function(c);
yd = Function(d);
if(yc < yd) { //remove the point outside the higher of the two center points - d is higher, remove b -> d bcomes b
b = d;
yb = yd;
d = c;
yd = yc;
c = a + (1.0 - (double) Fibonacci(n - 1 - k) / (double) Fibonacci(n - k)) * (b - a);
yc = Function(c);
}
else if(yc > yd) {
a = c;
you = yc;
c = d;
yc = yd;
d = a + ((double) Fibonacci(n - 1 - k) / (double) Fibonacci(n - k)) * (b - a);
yd = Function(d);
}
DispFibMsg(k, a, b, c, d, yc, yd);
}
}
//---------------------------------------------------------------------------

void __fastcall TForm1::DispFibMsg(int k, double a, double b, double c, double d, double yc, double yd)
{
mHistory->Lines->Add("k = " + IntToStr(k));
mHistory->Lines->Add("");
mHistory->Lines->Add("a = " + FloatToStrF(a, ffFixed, 8, 6) + ";\t b = " + FloatToStrF(b, ffFixed, 8, 6));
mHistory->Lines->Add("c = " + FloatToStrF(c, ffFixed, 8, 6) + ";\t d = " + FloatToStrF(d, ffFixed, 8, 6));
mHistory->Lines->Add("yc = " + FloatToStrF(yc, ffFixed, 10, 9) + ";\t yd = " + FloatToStrF(yd, ffFixed, 10, 9));
mHistory->Lines->Add("");
Application->ProcessMessages();
}
//---------------------------------------------------------------------------

long __fastcall TForm1::Fibonacci(int n) // n >= 3
{
long fibn, fibpn, temp; //fibonacci[n] and fibonacci[n-1]
int i;
fibn = 2;
fibpn = 1;
for(i = 1; i < n - 2; i++) {
temp = fibn;
fibn = fibn + fibpn;
fibpn = temp;
}
return fibn;
}
//---------------------------------------------------------------------------

int __fastcall TForm1::ReverseFibonacci(double Fn)
{ //returns which term of the Fibonacci series has a value greater than Fn
long fibn, fibpn, temp; //fibonacci[n] and fibonacci[n-1]
int i;
fibn = 2;
fibpn = 1;
for(i = 1; i < 1000; i++) {
temp = fibn;
fibn = fibn + fibpn;
fibpn = temp;
if(fibn >= Fn)
return i + 3;
}
return -1;
}
//---------------------------------------------------------------------------

 

[Ray Vickson]

Hi All,

I completed the Fibonacci Search program. Because I could not find a C/C++ sample program, I will add my work in here, maybe it will help somebody else. I used Borland C++ Builder, but it is simple enough to be ported to other platforms.

Regards,
Nicolae

void __fastcall TForm1::btnFibonacciSearchClick(TObject *Sender)
{
FibonacciSearch();
}
//---------------------------------------------------------------------------

double __fastcall TForm1::Function(double x)
{
double y;
// y = 0.5*(0.29*x*(x+2)+0.22*x+3-3*sqrt(x*x-3*x+9)+12.88-x/1.33);
// y = x*x*x*x*x -(4.0/5.0) * x + 0.5;
y = x * x - sin(x);
return y;
}
//---------------------------------------------------------------------------

void __fastcall TForm1::FibonacciSearch()
{
double y, x;
double a, b, c, d;

double eps = 1e-4;
double e = 0.01;
double Fn;
double r;
double ya, yb, yc, yd;
int n, k;

a = 0;
b = 1;
Fn = (b - a) / eps;
n = ReverseFibonacci(Fn);
for(k = 1; k <= n; k++) {
r = (double) Fibonacci(n - k - 1) / (double) Fibonacci(n - k);
c = a + (1 - r) * (b - a);
d = a + r * (b - a);
yc = Function(c);
yd = Function(d);
if(yc < yd) { //remove the point outside the higher of the two center points - d is higher, remove b -> d bcomes b
b = d;
yb = yd;
d = c;
yd = yc;
c = a + (1.0 - (double) Fibonacci(n - 1 - k) / (double) Fibonacci(n - k)) * (b - a);
yc = Function(c);
}
else if(yc > yd) {
a = c;
you = yc;
c = d;
yc = yd;
d = a + ((double) Fibonacci(n - 1 - k) / (double) Fibonacci(n - k)) * (b - a);
yd = Function(d);
}
DispFibMsg(k, a, b, c, d, yc, yd);
}
}
//---------------------------------------------------------------------------

void __fastcall TForm1::DispFibMsg(int k, double a, double b, double c, double d, double yc, double yd)
{
mHistory->Lines->Add("k = " + IntToStr(k));
mHistory->Lines->Add("");
mHistory->Lines->Add("a = " + FloatToStrF(a, ffFixed, 8, 6) + ";\t b = " + FloatToStrF(b, ffFixed, 8, 6));
mHistory->Lines->Add("c = " + FloatToStrF(c, ffFixed, 8, 6) + ";\t d = " + FloatToStrF(d, ffFixed, 8, 6));
mHistory->Lines->Add("yc = " + FloatToStrF(yc, ffFixed, 10, 9) + ";\t yd = " + FloatToStrF(yd, ffFixed, 10, 9));
mHistory->Lines->Add("");
Application->ProcessMessages();
}
//---------------------------------------------------------------------------

long __fastcall TForm1::Fibonacci(int n) // n >= 3
{
long fibn, fibpn, temp; //fibonacci[n] and fibonacci[n-1]
int i;
fibn = 2;
fibpn = 1;
for(i = 1; i < n - 2; i++) {
temp = fibn;
fibn = fibn + fibpn;
fibpn = temp;
}
return fibn;
}
//---------------------------------------------------------------------------

int __fastcall TForm1::ReverseFibonacci(double Fn)
{ //returns which term of the Fibonacci series has a value greater than Fn
long fibn, fibpn, temp; //fibonacci[n] and fibonacci[n-1]
int i;
fibn = 2;
fibpn = 1;
for(i = 1; i < 1000; i++) {
temp = fibn;
fibn = fibn + fibpn;
fibpn = temp;
if(fibn >= Fn)
return i + 3;
}
return -1;
}
//---------------------------------------------------------------------------

Have you looked at "Numerical Recipes in C" or "Numerical Recipes in C++"? They have routines for Golden Section search (but not Fibonacci?). The NRC book has a slightly obsolete version that is freely downloadable; just Google Numerical Recipes". If you actually purchase one of the books you get a CD or DVD (maybe nowadays a downloadable file) containing all the programs. Professional numerical analysts have complained that the NR routines are not always the best, or most efficient, or most robust, but they are almost always better than what an average user could put together in a few hours.

RGV