3
$\begingroup$

I have been facing the problem of finding the true minimum of a rapidly varying functions. The commands like FindMinimum and NMinimize often lead to the minimum value which is the global minimum for that function. I try to present an example here:

y[x_] = Sin[10 x^2] + 5 Cos[20 x];  

When I plot this function say in region {x,0,5}, I can see the minima which are lower in value than returned by the FindMinimum and NMinimize commands. Can anyone help to figure out what is going wrong here? Thanks.

$\endgroup$
  • 2
    $\begingroup$ Have you tried the options to either NMinimize or FindMinimum? Take a look at the tutorial here. For example, try NMinimize[Sin[10 x^2] + 5 Cos[20 x], {x}, Method -> {"RandomSearch", "SearchPoints" -> 500}]. $\endgroup$ – dr.blochwave Oct 13 '16 at 7:25
  • 2
    $\begingroup$ I notice you haven't accepted any answers yet to any of your questions. It's not much, but it's the way to show gratitude on the site to those who have given you help voluntarily, without charge or recompense, and to help other visitors find answers to their questions. If there are answers that solve your problems, please select the one that in your opinion is the best. $\endgroup$ – Michael E2 Mar 9 '17 at 13:01
  • $\begingroup$ Oh I am sorry Sir. The answers have really helped me a lot. $\endgroup$ – Jee Mar 9 '17 at 17:39
  • $\begingroup$ But could you please tell me what do you mean by accepting the answer? I would definitely like to thank all the people who have helped me on this site in one or the other ways. $\endgroup$ – Jee Mar 9 '17 at 17:41
  • $\begingroup$ @Jee See my first comment: click the checkmark sign. There's a link that explains it, if you have trouble finding the sign to click. $\endgroup$ – Michael E2 Mar 9 '17 at 17:43
6
$\begingroup$

MinimalBy is an easy way.

y[x_] = Sin[10 x^2] + 5 Cos[20 x];
d = Table[{x, y[x]}, {x, 1, 5, 0.0001}];
{X, Y} = First@MinimalBy[d, Last]

{1.7278, -5.99996}

Edit

For a better overview one can do the same with a plot.

plot = Plot[y[x], {x, 0, 5}, PlotPoints -> 500];

Determine the plot points.

p1 = Join @@ Cases[Normal@plot, Line[x1__] :> x1, Infinity];
{X, Y} = First@MinimalBy[p1, Last]
{1.72783, -5.99996}

Length@p1
**8635**

Edit 2

plt = Plot[y[x], {x, 0, 5}, Mesh -> {{0}}, MeshFunctions -> {y'[#] &}, PlotPoints -> 500];
p2 = Cases[Normal[plt], Point[pt_] :> pt, Infinity];
{X, Y} = First@MinimalBy[p2, Last]
{1.72776, -5.99996}

Length@p2
**61**

To compare the number of points p1 and p2, I chose the same number of plot points. @J.M. 's hint to add Mesh->{{0}},MeshFunctions->{y'[#]&} in the plot greatly reduces the number of points in p2.

Plot[y[x], {x, 0, 5}, Mesh -> {{0}}, MeshFunctions -> {y'[#] &}, 
 MeshStyle -> PointSize[Medium], Epilog -> {Red, PointSize@Large, Point@{X, Y}}]

enter image description here

Thanks to @J.M. for pointing this out!

$\endgroup$
  • 1
    $\begingroup$ Since you're plotting anyway, use MeshFunctions: f[x_] := Sin[10 x^2] + 5 Cos[20 x]; plt = Plot[f[x], {x, 0, 5}, Mesh -> {{0}}, MeshFunctions -> {f'[#] &}, PlotPoints -> 75]; MinimalBy[Cases[Normal[plt], Point[pt_] :> pt, Infinity], Last]. $\endgroup$ – J. M. will be back soon Oct 13 '16 at 12:43
  • $\begingroup$ @ J. M. Hmm.., I want to get the red point for a visually test having taken the minimal point. $\endgroup$ – user36273 Oct 13 '16 at 15:40
  • $\begingroup$ @ J. M. Why to add Mesh->{{0}},MeshFunctions->{f'[#]&} in the plot only to replace Line[x1__] :> x1 with Point[pt_] :> pt ? Of course I have to add Join@@ at the beginning, but thats all. $\endgroup$ – user36273 Oct 13 '16 at 16:09
  • $\begingroup$ Because that extracts only the extrema of the function as opposed to all the points. $\endgroup$ – J. M. will be back soon Oct 13 '16 at 16:19
  • $\begingroup$ @ J. M. I know MeshFunctions->{f'[#]& . Thanks for the hint. $\endgroup$ – user36273 Oct 13 '16 at 16:24
5
$\begingroup$

Too long for a comment.


y[x_] = Sin[10 x^2] + 5 Cos[20 x];
d = Table[{x, y[x]}, {x, 1, 5, 0.0001}];

min = Min@d[[All, 2]]

-5.99996

pos = Position[d, min]

{{7279, 2}}

d[[pos[[1, 1]]]]

{1.7278, -5.99996}

(In an instant.)


On the other hand:

NMinimize[{y[x], 5 > x > 0}, x]

{-5.99954, {x -> 4.87009}}

NMinimize[{Sin[10 x^2] + 5 Cos[20 x], 5 > x > 0}, {x}, 
 Method -> {"RandomSearch", "SearchPoints" -> 500}]

{-5.99954, {x -> 4.87009}}


Also

fmin = Quiet@
    Table[FindMinimum[{y[x], 5 > x > 0}, {x, i}], {i, 0, 5, 
      0.01}]; // AbsoluteTiming

{4.29685, Null}

fmin[[Position[fmin, Min@fmin[[All, 1]]][[1, 1]]]]

{-5.99997, {x -> 1.72776}}

so the first method is fastest.

$\endgroup$
4
$\begingroup$

On such a univariate function, a direct calculus approach is feasible:

With[{f = Sin[10 x^2] + 5 Cos[20 x]},
 Min[f /. NSolve[{D[f, x], 0 <= x <= 5}, x]]
 ]
(*  -5.99997  *)

The numerical minimizers tend to involve some sort of search method, and generally only return the least local minimum found. From the "Details" of NMinimize:

  • If f and cons are linear, NMinimize can always find global minima, over both real and integer values.
  • Otherwise, NMinimize may sometimes find only a local minimum.
$\endgroup$
  • 1
    $\begingroup$ It is also possible to use Solve[D[f,x]==0 && 0<=x<=5, x] to get "exact" results. $\endgroup$ – Carl Woll Mar 9 '17 at 15:30
4
$\begingroup$

Using FindRoot[] to find the zeros of the first derivative of the function in a given interval {a,b} with a certain granularity of inspection d we arrive at:

y[x_] = Sin[10 x^2] + 5 Cos[20 x];

a = -1; (* lower endpoint of search interval *)
b = 3; (* upper endpoint of search interval *)
d = 0.01; (* sampling interval *)
psize = 0.015; (* point size of extrema in plot *)
data = {#, y[#]} & /@ 
  Union[Flatten[
    Table[x /. FindRoot[y'[x] == 0, {x, {z}}], {z, a, b, d}] // Chop],
    SameTest -> (Abs[#1 - #2] < d/2 &)];
datamax = Select[data, #[[2]] > 0 &]; 
datamin = Select[data, #[[2]] < 0 &];

Plot[y[x], {x, a, b}, 
 Epilog -> {PointSize -> psize, {Black, Point /@ datamax}, {Red, 
    Point /@ datamin}}]

enter image description here

$\endgroup$
4
$\begingroup$

One view of the problem focuses on the issue of getting a sample point in the one valley where the absolute minimum lies. If you know something about the minimum width of the valleys of your objective function(s), then one can sample at regular intervals of less than this width. To write a generic solver, the minimum width would be unknown. In this answer, I will present an approach that should work on a fairly smooth objective function.

The method is based on making a Chebyshev approximation to the function. When the approximation is fairly good, there will usually be several sample points per oscillation, assuming the amplitude of the oscillations are fairly large. The function values at the sample points are computed Chebyshev approximation. These can be used to locate valleys using PeakDetect[-values]. Finally FindMinimum is used to polish the minima, and the ones with the least function value are returned. (It's set up to return multiple solutions, if the absolute minimum occurs at several locations. It's possible that some locations are missed in some cases.)

Example uses

OP's example:

min1 = chebMinimize[Sin[10 x^2] + 5 Cos[20 x], {x, 0, 5}]  // AbsoluteTiming
(*  {0.020218,  {-5.99997, {{x -> 1.72776}}} }  *)

Minima at end points:

chebMinimize[-x^2, {x, -1, 1}]
(*  {-1., {{x -> 1.}, {x -> -1.}}}  *)

Example with almost 10000 local minima:

obj = Sin[50000 x + 1]/1000 + Sin[200 Cos[Pi x/2]^8 x];
min1 = chebMinimize[obj, {x, -1, 1}, MaxRecursion -> 20]  // AbsoluteTiming
(*  {4.79622,  {-1.001, {{x -> -0.605373}}} }  *)

The underlying Chebyshev approximation is much faster than 10000 FindMinima[] calls. One can limit the number of local minima and speed things up:

min1 = chebMinimize[obj, {x, -1, 1}, MaxRecursion -> 20, "MaxMinima" -> 100]  //
  AbsoluteTiming
(*  {0.08998,  {-1.001, {{x -> -0.605373}}} }  *)

The least local minima (in the discrete list of values of obj at the sample points) are chosen, but it's possible to undersample. The success depends on where the sample points fall with respect to the absolute minimum. Note that the result below, while close, is not quite right.

min2 = chebMinimize[obj, {x, -1, 1}, MaxRecursion -> 20, "MaxMinima" -> 10]  //
  AbsoluteTiming
(*  {0.046868,  {-1.001, {{x -> -0.340726}}} }  *)

min1[[2, 1]] - min2[[2, 1]]
(*  -1.07694*10^-7  *)

Code dump

There are two helper functions. The function cheb`approx[] computes a Chebyshev series approximation to a function f over an interval {x, a, b}; it returns the sample points/abscissae and function values in the form {xvals, fvals}. The function cheb`minpos[] returns the positions of the local minima in a list of values.

There is an irrelevant data structure cheb`data that I've been using to inspect the workings of the algorithm. For example, Length@cheb`data["Y"] will give the number of function values, and Length@cheb`data["M"] gives the number of local minima. (It gets overwritten with every call to chebMinimize[].)

ClearAll[chebMinimize];
Options[chebMinimize] = {MinRecursion -> 4, MaxRecursion -> 10, 
   Tolerance -> 1.*^-8, 
   "SuppressMessages" :> {FindMinimum::lstol, FindMinimum::reged}, 
   "MaxMinima" -> Infinity};
chebMinimize::ncv = 
  "chebMinimize failed to converge to prescribed accuracy after `1` \
recursive bisections of {`2`, `3`}. chebMinimize obtained `4` for the \
relative error estimate.";

Clear[cheb`minpos, cheb`approx];
(* Returns the indices of the search starting point of the bracketing sample points.
* Note that the Chebyshev abscissae are reversed *)
cheb`minpos[v_] := Flatten[SparseArray[PeakDetect[-v, 0, 0, Infinity]]["NonzeroPositions"]];

(* Returns the abscissae and function values of a
* Chebyshev approximation of f of absolute error less than tolerance *)
cheb`$keys = {
   "X", (* x-values/abscissae *)
   "Y", (* y-values/function values *)
   "F", (* the function f *)
   "I", (* the interval {x, a, b} *)
   "E", (* error *)
   "C", (* Chebyshev coefficients *)
   "M", (* positions of the local minima *)
   "FM" (* FindMinimum results *)
   };
cheb`data = AssociationMap[None &, cheb`$keys];  (* not needed; used for diagnostics *)
cheb`approx[f_, {x_, a_, b_}, tolerance_, minR_: 4, maxR_: 10] := 
  Module[{vf, xvals, fvals, error, r = minR,
    coeffs},  (* coeffs are not needed; saved in cheb`data *)
   If[Quiet@ListQ[f /. x -> {N@a}],    (* vf = Listable-ized f *)
    vf = Evaluate[f /. x -> #] &,
    vf = Map[Evaluate[f /. x -> #] &, #] &;
    ];
   {xvals, fvals} = NestWhile[
     Function[{xx, ff}, 
        (* recursive bisection of abscissae *)
        With[{new = Rescale[Cos[Pi Range[1., 2^r, 2]/2^r], {-1, 1}, {a, b}]},
         r++;
         {Riffle[xx, new], Riffle[ff, vf[new]]}   (* riffle new values *)
         ]
        ] @@ # &,
     (r++; {#, vf[#]}) &@     (* initial subdivision of {a,b} *)
      Rescale[Cos[Pi Range[0., 2^r]/2^r], {-1, 1}, {a, b}],
     (error =                 (* calculate error from tail of Chebyshev series *)
        Total@Abs@#[[-Max[2, r] ;;]]/Max@Abs[#] &[
          coeffs = FourierDCT[Last@#, 1]]) > tolerance &,
     1,
     1 + maxR - minR
     ];
   If[error > tolerance,      (* check error *)
    Message[chebMinimize::ncv, maxR, a, b, error]
    ];
   cheb`data = 
    AssociationThread[cheb`$keys -> {xvals, fvals, f, {x, a, b}, error, coeffs, {}, {}}];
   {xvals, fvals}
   ];

chebMinimize[f_, {x_, a_, b_}, opts : OptionsPattern[]] :=
  Module[{xvals, fvals, minima, res},
   (* get sample points & values *)
   {xvals, fvals} = 
    cheb`approx[f, {x, a, b}, 
     OptionValue[Tolerance],
     OptionValue[MinRecursion],
     OptionValue[MaxRecursion]];

   (* get positions of minima in fvals *)
   minima = cheb`minpos[fvals];
   If[NumericQ[OptionValue["MaxMinima"]] && 
       Length@minima > OptionValue["MaxMinima"],
    minima = minima[[
       Ordering[fvals[[minima]], OptionValue["MaxMinima"]]
       ]]
    ];
   cheb`data["M"] = minima;

   (* polish minima with FindMinimum *)
   res = With[{msgs = OptionValue[chebMinimize, {opts}, "SuppressMessages", Hold] /.
       {All | Automatic :> Sequence[], None -> {}}},
     msgs /. Hold[m___] :> MinimalBy[First]@
        Quiet[
         cheb`data["FM"] = FindMinimum[f,
             {x,
              xvals[[#]],
              xvals[[Clip[# + 1, {1, Length@xvals}]]], 
              xvals[[Clip[# - 1, {1, Length@xvals}]]]}
             ] & /@ minima,
         m]];
   MapAt[First, Transpose@res, 1]
   ];

Some analysis of the algorithm

The accuracy of a Chebyshev series approximation can be estimated from the tail of the series: $$\left|f(x) - \sum_{k=0}^N c_k T_k(x)\right| = \left|\sum_{k=N+1}^\infty c_k T_k(x)\right| \le \sum_{k=N+1}^\infty \left|c_k \right|$$ since $|T_k(x)| = |\cos(k \arccos x)| \le 1$. It takes $N+1$ Chebyshev sample points to get a series of order $N$. The accuracy of a Chebyshev series approximation to an oscillatory function generally bounces around the amplitude of the oscillations until the number is sufficiently large. The approximation cannot really get a grip on the function until there is at least one point per extrema, but often, especially for irregular oscillations, it has will be considerable more. When convergence starts, it will be rapid for an analytic function, and the error can estimated from the last few coefficients of the Chebyshev approximation.

The term Sin[200 Cos[Pi x/2]^8 x] used in the extremely oscillatory example above has a "maximum frequency" of 200 / (2 Pi) near x == 0. That would translate to about 127 extrema. We can see that the Chebyshev coefficients $c_k$ start to settle down around $k = 127$, but they do not start their convergence to zero until $k$ is more than $200$. At $k = 460$, the approximation reaches the limit of machine precision. (Note also that the function is odd, which means the odd-order Chebyshev coefficients are theoretically zero. The function also has some very flat minima, which causes FindMinimum to complain.)

obj = Sin[200 Cos[Pi x/2]^8 x];
chebMinimize[obj, {x, -1, 1}, "SuppressMessages" -> All];  (* load cheb`data *)
cdata = RealExponent@cheb`data["C"];
ListPlot[cdata, PlotRange -> {-18, 2}, GridLines -> {{127, 200, 460}, None}]

Mathematica graphics

A similar thing happens with the other term Sin[50000 x + 1]/1000, except that convergence does not start until we have about 50000 points. Note that the level of the pre-convergent phase is about the amplitude of the function, 1/1000.

obj = Sin[50000 x + 1]/1000;
chebMinimize[obj, {x, -1, 1}, MaxRecursion -> 20, "MaxMinima" -> 100]; (* load cheb`data *)
cdata = RealExponent@cheb`data["C"];
ListPlot[cdata, PlotRange -> {-18, 2}]

Mathematica graphics

The full example is a combination of these two processes. Note that oddness of the second term is reflected in the plateau at 10^-3, which is magnitude of the first term (lower left). In general, one should expect that for the algorithm to be guaranteed to be successful, the approximation should have an error much less than this minimum amplitude of the oscillations. When this is not known, one can see that a machine-precision approximation is possible. If the oscillations of this order, then it will hopeless anyway, unless the working precision is increased.

obj = Sin[50000 x + 1]/1000 + Sin[200 Cos[Pi x/2]^8 x];
chebMinimize[obj, {x, -1, 1}, MaxRecursion -> 20, "MaxMinima" -> 100]; (* load cheb`data *)
cdata = RealExponent@cheb`data["C"];
ListPlot[cdata, PlotRange -> {-18, 2}]
GraphicsRow[
 Module[{i, j, plot},
  plot := ListPlot[cdata[[i ;; j]], DataRange -> {i, j}, PlotRange -> {-18, 2}];
  {i = 1; j = 600; plot,
   i = 49500; j = Length@cdata - 15000; plot}
  ], ImageSize -> 500]

Mathematica graphics

Here is a closer look at the sampling and local minima. The sample points are shown in green and red, with the red one being the seeds for FindMinimum. The local minima returned by FindMinimum are in cyan.

With[{samples = Transpose@{cheb`data["X"], cheb`data["Y"]}},
 Plot[obj, {x, -1, 1},
  PlotPoints -> 60, PlotStyle -> {Thickness[0.003], Lighter@Blue}, 
  Background -> Black,
  Epilog -> {PointSize[0.004],
    Green, Point[samples],
    PointSize[0.012], Red, Point@samples[[cheb`data["M"]]],
    PointSize[0.005], Cyan, 
    Point[{x, obj} /. cheb`data["FM"][[All, 2]]]},
  Frame -> True, PlotRangePadding -> {0, Scaled[.05]}]
 ]

Mathematica graphics

Here is a close-up of the region where the absolute minimum was found. We can see that there are at least three sample points per oscillation, just enough to identify good seeds for FindMinimum, and that every local minimum has been found.

With[{samples = Transpose@{cheb`data["X"], cheb`data["Y"]}},
 Plot[obj, {x, -0.607, -0.6024},
  PlotPoints -> 60, PlotStyle -> {Thickness[0.003], Lighter@Blue}, 
  Background -> Black,
  Epilog -> {PointSize[0.008],
    Green, Point[samples],
    PointSize[0.015], Red, Point@samples[[cheb`data["M"]]],
    PointSize[0.008], Cyan, 
    Point[{x, obj} /. cheb`data["FM"][[All, 2]]]},
  Frame -> True, PlotRangePadding -> {0, Scaled[.05]}]
 ]

Mathematica graphics

$\endgroup$
  • $\begingroup$ I wonder: once you have the approximate locations of the extrema, why not differentiate the Chebyshev series and then use Newton-Raphson? $\endgroup$ – J. M. will be back soon Oct 16 '16 at 5:33
  • $\begingroup$ @J.M. Several reasons, I suppose. (1) Keep the code shorter for the post. (2) The original function & its derivative are probably faster. (3) FindMinimum does all sorts of checks, in case the initial points are not as good as imagined. It does in fact give warnings (hence Quiet[] above), and the reason for bracketing the roots is because in one test it gave a minimum well outside interval, past several local minima. (4) I don't know what to do when the initial point lies on a concave down segment. (I'm sure I could look it up. Someday.) $\endgroup$ – Michael E2 Oct 16 '16 at 11:54
  • 1
    $\begingroup$ @J.M. The better way, I think, it to use a piecewise Chebyshev approximation, differentiate, use the colleague matrices of the pieces, and pick out the minimum. $\endgroup$ – Michael E2 Oct 16 '16 at 17:02
2
$\begingroup$

Yet another way: A direct search for local minima with NDSolve, polished with NMinimize, which essentially uses FindMinimum on the seeds found (this approach was alluded to at the end this answer).

{seeds} = Last@Reap@NDSolve[{t'[x] == 1, t[0] == 0, f[x] == y[x],
      WhenEvent[f'[x] > 0, Sow[{x}]]},
     {}, {x, 0, 5}];

{min, sol} = NMinimize[{y[x], 0 <= x <= 5}, x, 
  Method -> {"RandomSearch", "InitialPoints" -> pts}]

NMinimize::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations.

(*  {-5.99997, {x -> 1.72776}}  *)

The problem with convergence is insufficient working precision. A fix is to use WorkingPrecision -> 16 or higher. At machine precision, the Newton update, which is appproximately $-y'(x)/y''(x) \approx 7\times10^{-10}$, results in no change in value of $y$:

Differences@y[{x, x - y'[x]/y''[x]} /. sol]
(*  {0.}  *)

Compare at a higher precision:

Differences@y[{x, x - y'[x]/y''[x]} /. SetPrecision[sol, 32]]
(*  {-7.896393357752822*10^-16}  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.