Techniques to find all local minima of black box function with n continuous derivatives?

Question

I have a blackbox function f on [a,b] that is known to be continuous, and has n continuous derivatives. The function uses polynomial interpolations (but isn't a polynomial interpolation itself), so no guarantee is made about the (n+1)st derivative, though it's likely to be discontinuous.

Even though the function has n continuous derivatives, the derivatives cannot be evaluated directly. You can of course look at the limit of (f(x+h)-f(x))/h as h->0, but the accuracy depends on how you take the numerical limit.

Is there a good general technique for finding all local minima of f?

A technique that works surprisingly well is to Plot f, look at the points Plot has chosen (ie, the [[1,1,3,2,1]] part of what Plot returns), and find local minima in the list (which is easy). The local minima of the list are often the local minima of f.

If the function wiggles too much on [a,b], breaking [a,b] into equal size subintervals and Plot-ting works.

This sort of makes sense, because Plot "seeks out" minima (and maxima) to plot accurately, instead of plotting at uniform intervals.

However, this annoys me. Plot is supposed to be for plotting, and using it to find minima seems wrong.

What's a better method?

Answer 1

If you want to use Plot to find the extrema of a black-box function, here is an implementation. I decided to look for maxima here, but you can turn this into minima by flipping the sign of the function. It works by first applying MaxDetect to the plot points in a given interval, and then constructing brackets for FindMaximum from the neighboring points. I'll apply it to the same example function that I used in a related answer for root finding:

data = NDSolve[{1.09 x''[t] - 0.05 x'[t] + 1.1759 Sin[x[t]] == 0, 
    x[0] == Pi/3, x'[0] == 0}, x, {t, 0, 50}];

f[t_] = (x /. First[data])[t];

Clear[findAllMaxima]
SyntaxInformation[
   findAllMaxima] = {"LocalVariables" -> {"Plot", {2, 2}}, 
   "ArgumentsPattern" -> {_, _, OptionsPattern[]}};
SetAttributes[findAllMaxima, HoldAll];

Options[findAllMaxima] = 
  Join[{"ShowPlot" -> False, PlotRange -> All}, 
   FilterRules[Options[Plot], Except[PlotRange]]];

findAllMaxima[fn_, {l_, lmin_, lmax_}, opts : OptionsPattern[]] := 
 Module[{pl, p, x, px, localFunction, brackets}, 
  localFunction = ReleaseHold[Hold[fn] /. HoldPattern[l] :> x];
  If[lmin != lmax, 
   pl = Plot[localFunction, {x, lmin, lmax}, 
     Evaluate@
      FilterRules[Join[{opts}, Options[findAllMaxima]], Options[Plot]]];
   p = Cases[pl, Line[{x__}] :> x, Infinity];
   px = p[[All, 1]];
   If[OptionValue["ShowPlot"], 
    Print[Show[pl, PlotLabel -> "Finding maxima for this function", 
      ImageSize -> 200, BaseStyle -> {FontSize -> 8}]]], p = {}];
  brackets = 
   Transpose[
    Map[Pick[px, #, 1] &, {RotateLeft[#], RotateRight[#]} &[
      MapAt[0 &, MaxDetect[p[[All, 2]]], {{1}, {-1}}]]]];
  Flatten@Apply[FindArgMax[localFunction, {x, ##1}] &, brackets, {1}]
  ]

findAllMaxima[f[t], {t, 0, 50}, "ShowPlot" -> True]

{6.58301,13.4054,20.6403,28.7857}

The special case that a maximum is situated right at the interval boundary is handled by dropping the corresponding entry in MaxDetect. This is what the MapAt command is for.

Here are the minima of the above function:

findAllMaxima[-f[t], {t, 0, 50}]

{3.26812,9.95657,16.9538,24.533,34.2571}

Edit: comparison to NDSolve approach

For comparison with the approach based on NDSolve, I copied MichaelE2's test function and tried it with the above method:

ng[t_?NumericQ] = g[[3]];

findAllMaxima[-ng[t], {t, 0, 10}]

{1.35809,2.84282,3.79441,4.5406,5.18681,5.75923,6.27529,6.759,7.21652,7.63862,8.03896,8.41996,8.78808,9.1421,9.47799,9.80125}

The results agree within $10^{-7}$, and the Plot based approach is noticeably faster than the NDSolve approach.

Answer 2

Here's an adaptation of the NDSolve method used by Daniel Lichtblau that uses numeric (finite difference) derivatives. Just what to return and in what form is a choice. The following returns a list of FindMinimum style solutions {{m1, {x -> x1}}, {m2, {x -> x2}},...}. The solutions may be polished further with FindMinimum.

findAllMinima[f_, {x_, x1_, x2_}, opts : OptionsPattern[NDSolve]] := 
   Module[{nf, res, xx, df, y},
   nf = Experimental`CreateNumericalFunction[{x}, {f}, {1}, Jacobian -> FiniteDifference];
   df[x0_?NumericQ] := nf["Jacobian"[{x0}]][[1, 1]];
   res = Reap[NDSolve[
       {y'[xx] == df[xx], y[x1] == First@nf[{x1}],
        WhenEvent[y'[xx] > 0, Sow[{y[xx], {x -> xx}}]]}, 
       y, {xx, x1, x2}, opts]][[2, 1]]
   ];

Sample numeric functions (InterpolatingFunctions) obtained by integrating a slightly noisy, oscillatory function:

SeedRandom[1];
xdata = Sort@ DeleteDuplicates[Join[{0., 100.}, RandomReal[{0, 100}, 1000]], 
    Chop@Subtract[##] == 0 &];
ydata = xdata^2 (Sin[xdata] + Cos[3 xdata]) + 10 RandomReal[{-1, 1}, Length@xdata];
g0 = Interpolation[Transpose[{Sqrt[xdata], ydata}], InterpolationOrder -> 1];

g = NestList[Integrate[#, t] &, g0[t], 3]; (* four functions, increasing differentiable *)

We'll use the last two, which are relatively smooth (see plots below), and we'll wrap the functions in a NumericQ-protected black-box (just in case).

ng[t_?NumericQ] = g[[3]];
findAllMinima[ng[t], {t, 0, 10}]
(*
  {{-1.19007, {t -> 1.35809}}, {-5.11378, {t -> 2.84282}}, {-9.42955, {t -> 3.79441}},
   {-12.3119, {t -> 4.5406}}, {-13.9851, {t -> 5.18681}}, {-13.762, {t -> 5.75923}},
   {-10.335, {t -> 6.27529}}, {-1.29706, {t -> 6.759}}, {-1.39131, {t -> 7.21652}},
   {-4.57187, {t -> 7.63862}}, {-5.02481, {t -> 8.03896}}, {-4.8891, {t -> 8.41996}},
   {-2.75207, {t -> 8.78808}}, {-22.9376, {t -> 9.1421}}, {-39.9945, {t -> 9.47799}},
   {-54.2516, {t -> 9.80125}}}
*)

Probably, one should set the PrecisionGoal smaller by default, as it is quite a bit faster. And it is still accurate on this example.

findAllMinima[ng[t], {t, 0, 10}]; // AbsoluteTiming
findAllMinima[ng[t], {t, 0, 10}, PrecisionGoal -> 3]; // AbsoluteTiming
(*
  {1.31999, Null}
  {0.056883, Null}
*)

Examples:

ng[t_?NumericQ] = g[[4]];
mins = findAllMinima[ng[t], {t, 0, 10}, PrecisionGoal -> 3];
Plot[g[[3]], {t, 0, 10}, 
 Epilog -> {PointSize[Medium], Red, 
   Point[Transpose[{t /. mins[[All, 2]], mins[[All, 1]]}]]}, 
 PlotPoints -> 201, PlotRange -> All]

ng[t_?NumericQ] = g[[4]];
mins = findAllMinima[ng[t], {t, 0, 10}, PrecisionGoal -> 3];
Plot[g[[4]], {t, 0, 10}, 
 Epilog -> {PointSize[Medium], Red, 
   Point[Transpose[{t /. mins[[All, 2]], mins[[All, 1]]}]]}, 
 PlotPoints -> 201, PlotRange -> All]