How do I efficiently find the maximum value in an array containing values of a smooth function?

Question

I have a function that takes a floating point number and returns a floating point number. It can be assumed that if you were to graph the output of this function it would be 'n' shaped, ie. there would be a single maximum point, and no other points on the function with a zero slope. We also know that input value that yields this maximum output will lie between two known points, perhaps 0.0 and 1.0.

I need to efficiently find the input value that yields the maximum output value to some degree of approximation, without doing an exhaustive search.

I'm looking for something similar to Newton's Method which finds the roots of a function, but since my function is opaque I can't get its derivative.

Answer 1

You could reduce it to a simple linear fit on the delta's, finding the place where it crosses the x axis. Linear fit can be done very quickly.

Or just take 3 points (left/top/right) and fix the parabola.

It depends mostly on the nature of the underlying relation between x and y, I think.

edit this is in case you have an array of values like the question's title states. When you have a function take Newton-Raphson.

Answer 2

Edit 2: The method described in Jive Dadson is a better way to go about this. I'm leaving my answer up since it's easier to implement, if speed isn't too much of an issue.

Use a form of binary search, combined with numeric derivative approximations.

Given the interval [a, b], let x = (a + b) /2 Let epsilon be something very small.

Is (f(x + epsilon) - f(x)) positive? If yes, the function is still growing at x, so you recursively search the interval [x, b] Otherwise, search the interval [a, x].

There might be a problem if the max lies between x and x + epsilon, but you might give this a try.

Edit: The advantage to this approach is that it exploits the known properties of the function in question. That is, I assumed by "n"-shaped, you meant, increasing-max-decreasing. Here's some Python code I wrote to test the algorithm:

def f(x):
    return -x * (x - 1.0)

def findMax(function, a, b, maxSlope):
    x = (a + b) / 2.0
    e = 0.0001
    slope = (function(x + e) - function(x)) / e
    if abs(slope) < maxSlope:
        return x
    if slope > 0:
        return findMax(function, x, b, maxSlope)
    else:
        return findMax(function, a, x, maxSlope)

Typing findMax(f, 0, 3, 0.01) should return 0.504, as desired.

Answer 3

The Levenberg-Marquardt algorithm is a Newton's method like optimizer. It has a C/C++ implementation levmar that doesn't require you to define the derivative function. Instead it will evaluate the objective function in the current neighborhood to move to the maximum.

BTW: this website appears to be updated since I last visited it, hope it's even the same one I remembered. Apparently it now also support other languages.

Answer 4

Given that it's only a function of a single variable and has one extremum in the interval, you don't really need Newton's method. Some sort of line search algorithm should suffice. This wikipedia article is actually not a bad starting point, if short on details. Note in particular that you could just use the method described under "direct search", starting with the end points of your interval as your two points.

I'm not sure if you'd consider that an "exhaustive search", but it should actually be pretty fast I think for this sort of function (that is, a continuous, smooth function with only one local extremum in the given interval).

Answer 5

For optimizing a concave function, which is the type of function you are talking about, without evaluating the derivative I would use the secant method.

Given the two initial values x[0]=0.0 and x[1]=1.0 I would proceed to compute the next approximations as:

def next_x(x, xprev):
    return x - f(x) * (x - xprev) / (f(x) - f(xprev))

and thus compute x[2], x[3], ... until the change in x becomes small enough.

Edit: As Jive explains, this solution is for root finding which is not the question posed. For optimization the proper solution is the Brent minimizer as explained in his answer.

Answer 6

I would like to down-thumb all the other answers so far, for various reasons, but I won't.

An excellent and efficient method for minimizing (or maximizing) smooth functions when derivatives are not available is parabolic interpolation. It is common to write the algorithm so it temporarily switches to the golden-section search (Brent's minimizer) when parabolic interpolation does not progress as fast as golden-section would.

I wrote such an algorithm in C++. Any offers?

UPDATE: There is a C version of the Brent minimizer in GSL. The archives are here: ftp://ftp.club.cc.cmu.edu/gnu/gsl/ Note that it will be covered by some flavor of GNU "copyleft."

As I write this, the latest-and-greatest appears to be gsl-1.14.tar.gz. The minimizer is located in the file gsl-1.14/min/brent.c. It appears to have termination criteria similar to what I implemented. I have not studied how it decides to switch to golden section, but for the OP, that is probably moot.

UPDATE 2: I googled up a public domain java version, translated from FORTRAN. I cannot vouch for its quality. http://www1.fpl.fs.fed.us/Fmin.java I notice that the hard-coded machine efficiency ("machine precision" in the comments) is 1/2 the value for a typical PC today. Change the value of eps to 2.22045e-16.