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

  1. 1. The problem statement, all variables and given/known data
    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).

    2. Relevant 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);

    3. 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
     
    Last edited: Mar 16, 2012
  2. jcsd
  3. NascentOxygen

    Staff: Mentor

    Hi NickF! [​IMG]

    Does the function have only one dip and that's the minimum?
     
    Last edited: Mar 16, 2012
  4. NascentOxygen

    Staff: Mentor

    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...&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.
     
    Last edited: Mar 16, 2012
  5. 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.

    Best regards,
    Nicolae
     
  6. 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
     
  7. NascentOxygen

    Staff: Mentor

    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?
     
  8. This sounds better, I will try to implement it. Thanks.

    Nicolae
     
  9. Ray Vickson

    Ray Vickson 6,000
    Science Advisor
    Homework Helper

    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 Search.pdf (Fibonacci) and
    http://en.wikipedia.org/wiki/Golden_section_search (Golden Section).

    RGV
     
  10. 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.

    Best regards,
    Nicolae
     
  11. Ray Vickson

    Ray Vickson 6,000
    Science Advisor
    Homework Helper

    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
     
  12. 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.

    Best regards,
    Nicolae
     
  13. Ray Vickson

    Ray Vickson 6,000
    Science Advisor
    Homework Helper

    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
     
  14. 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
     
    Last edited: Mar 17, 2012
  15. 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;
    ya = 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;
    }
    //---------------------------------------------------------------------------
     
  16. Ray Vickson

    Ray Vickson 6,000
    Science Advisor
    Homework Helper

    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
     
Know someone interested in this topic? Share a link to this question via email, Google+, Twitter, or Facebook

Have something to add?