Fast find of all local maximums in C++

Question

Problem

I have a formula for calculation of 1D polynomial, joint function. I want to find all local maximums of that function within a given range.

My approach

My current solution is that i evaluate my function in a certain number of points from the range and then I go through these points and remember points where function changed from rising to decline. Of cause I can change number of samples within the interval, but I want to find all maximums with as lowest number of samples as possible.

Question

Can you suggest any effetive algorithm to me?

Answer 1

Finding all the maxima of an unknown function is hard. You can never be sure that a maximum you found is really just one maximum or that you have not overlooked a maximum somewhere.

However, if something is known about the function, you can try to exploit that. The simplest one is, of course, is if the function is known to be rational and bounded in grade. Up to a rational function of grade five it is possible to derive all four extrema from a closed formula, see http://en.wikipedia.org/wiki/Quartic_equation#General_formula_for_roots for details. Most likely, you don't want to implement that, but for linear, square, and cubic roots, the closed formula is feasible and can be used to find maxima of a quartic function.

That is only the most simple information that might be known, other interesting information is whether you can give a bound to the second derivative. This would allow you to reduce the sampling density when you find a strong slope.

You may also be able to exploit information from how you intend to use the maxima you found. It can give you clues about how much precision you need. Is it sufficient to know that a point is near a maximum? Or that a point is flat? Is it really a problem if a saddle point is classified as a maximum? Or if a maximum right next to a turning point is overlooked? And how much is the allowable error margin?

If you cannot exploit information like this, you are thrown back to sampling your function in small steps and hoping you don't make too much of an error.

---

Edit:
You mention in the comments that your function is in fact a kernel density estimation. This gives you at least the following information:

• Unless the kernel is not limited in extend, your estimated function will be a piecewise function: Any point on it will only be influenced by a precisely calculable number of measurement points.

• If the kernel is based on a rational function, the resulting estimated function will be piecewise rational. And it will be of the same grade as the kernel!

  • If the kernel is the uniform kernel, your estimated function will be a step function.
    This case needs special handling because there won't be any maxima in the mathematical sense. However, it also makes your job really easy.

  • If the kernel is the triangular kernel, your estimated function will be a piecewise linear function.

  • If the kernel is the Epanechnikov kernel, your estimated function will be a piecewise quadratic function.

  In all these cases it is next to trivial to produce the piecewise functions and to find their maxima.

• If the kernel is of too high grade or transcendental, you still know the measurements that your estimation is based on, and you know the kernel properties. This allows you to derive a heuristic on how dense your maxima can get.

• At the very least, you know the first and second derivative of the kernel.

  • In principle, this allows you to calculate the first and second derivative of the estimated function at any point.

  • In the case of a local kernel, it might be more prudent to calculate the first derivative and an upper bound to the second derivative of the estimated function at any point.

  With this information, it should be possible to constrain the search to the regions where there are maxima and avoid oversampling of the slopes.

As you see, there is a lot of useful information that you can derive from the knowledge of your function, and which you can use to your advantage.

Answer 2

The local maxima are among the roots of the first derivative. To isolate those roots in your working interval you can use the Sturm theorem, and proceed by dichotomy. In theory (using exact arithmetic) it gives you all real roots.

An equivalent approach is to express your polynomial in the Bezier/Bernstein basis and look for changes of signs of the coefficients (hull property). Dichotomic search can be efficiently implemented by recursive subdivision of the Bezier.

There are several classical algorithms available for polynomials, such as Laguerre, that usually look for the complex roots as well.