# Why does NMaximize miss this global maximum?

I am having trouble maximizing a function which appears as a curvature of a planar curve.

{tmin, tmax} = {0, 2 Pi}

f = -((6-3 Cos[t] - Cos[3 t])/((-11+6 Cos[t] + 8 Cos[2 t] - 6 Cos[3 t] + Cos[4 t])
Sqrt[Cos[t]^2 + 9 Sin[t]^2 - 12 Cos[t] Sin[t]^2 + 4 Cos[t]^2 Sin[t]^2]));

NMaximize[{f, tmin <= t <= tmax}, t]


says that the maximum of $$f$$ is attained at

{1.37888, {t -> 5.78352}}


But,

Plot[f, {t, tmin, tmax}, PlotRange -> Full]


indicates that the true maximum is attained at $$t=\pi$$.

Why is this happening? I'm using Mathematica version 12.0.0 for Microsoft Windows (64-bit).

• Do not use bugs as a tag until other people have confirmed what you see is a bug. In this case, it definitely isn't; NMaximize[] isn't always guaranteed to give a global optimum. May 27, 2020 at 2:30
• From the docs: "Otherwise, NMaximize may sometimes find only a local maximum." May 27, 2020 at 2:31
• In this case, Maximize[{f, tmin <= t <= tmax}, t] works. May 27, 2020 at 2:32
• @J. M., I see. Ill remove the tag "bug". May 27, 2020 at 2:35
• You can also use NMaximize[{f, tmin <= t <= tmax}, t, Method -> {"DifferentialEvolution", "SearchPoints" -> 15}], NMaximize[{f, tmin <= t <= tmax}, t, Method -> {"NelderMead", "InitialPoints" -> List /@ Subdivide[tmin, tmax, 5]}] and so forth in addition to simulated annealing. My question is how do you know it's more robust? How do you know it's even giving a correct answer? (In this example, you know the answer ahead of time, which is an unrealistic use-case — if you know the answer already, you wouldn't use NMaximize to find it.) In general global optimization is difficult. May 27, 2020 at 12:36

This kind of problem — smooth, univariate function over a finite and relatively small domain — can be handled numerically by using NDSolve to locate the relative maxima, polishing them with FindMaximum, and then selecting the greatest one:

MaximalBy[First]@
With[{df2 = D[f, {t, 2}]},
FindMaximum[{f, tmin <= t <= tmax}, {t, #}] & /@
First@Last@Reap@NDSolve[
{y'[t] == D[f, t], y[0] == 0,
WhenEvent[y'[t] == 0 && df2 < 0, Sow[t]]},
y, {t, tmin, tmax}]
]

(*  {{5., {t -> 3.14159}}}  *)


[I'm sure this has been shown elsewhere on site, probably by me and several others. This problem can in fact be done exactly by Maximize, but the OP suggests there are other cases that might need a numerical approach.]

• Thank you very much for your neat code. I've been using MMA for a long time, but I have never used WhenEvent. May 27, 2020 at 3:02

Making use of Method, one obtains

NMaximize[{f, tmin <= t <= tmax}, t, Method -> "RandomSearch"]
(*{5., {t -> 3.14159}}*)


So does Method -> "SimulatedAnnealing".

• Thank you very much! That seems to be a very simple but robust way to get the right answer. May 27, 2020 at 10:30
• @A.Kato: Unfortunately,Method->"DifferentialEvolution" does not work here. May 27, 2020 at 11:29

Since your problem is single variable, we can also use Grid Search.

grid = Subdivide[2 π, 1000] // N;
val = f /@ grid;

Extract[#, Ordering[val, -1]] & /@ {val, grid}


{5., 3.14159}

Alternatively, as suggested by @J.M., we can use PeakDetect

plot = Plot[f[t], {t, tmin, tmax}, PlotPoints -> 1000, PlotRange -> All];
points = Join @@ Cases[Normal@plot, Line[x_] :> x, ∞];
peaks = Pick[points, PeakDetect[points[[All, 2]]], 1];
MaximalBy[peaks, Last]


{{3.14162, 5.}}

ListPlot[points, Epilog -> {Red, Point[peaks]}, PlotRange -> All]


• You could systematize this by using Plot[]'s adaptive sampling along with PeakDetect[]. May 27, 2020 at 12:32

Another option is to do as we do in calculus class. Find derivative, set to zero, find roots, find hessian, check sign. (not checking for saddle point :)

ClearAll["Global*"];
{tmin, tmax} = {0, 2 Pi};
f = -((6 - 3 Cos[t] -
Cos[3 t])/((-11 + 6 Cos[t] + 8 Cos[2 t] - 6 Cos[3 t] +
Cos[4 t]) Sqrt[
Cos[t]^2 + 9 Sin[t]^2 - 12 Cos[t] Sin[t]^2 +
4 Cos[t]^2 Sin[t]^2]));

diff    = D[f, t];
roots   = NSolve[diff == 0 && tmin <= t <= tmax, t]
hessian = D[f, {t, 2}] /. roots;
pts     = MapThread[{If[#2 > 0, Red, Blue], PointSize[0.02],
Point[{#1, f /. t -> #1}]} &, {t /. roots, hessian}];

Plot[f, {t, tmin, tmax}, PlotRange -> All, Epilog -> pts,
GridLines -> Automatic, GridLinesStyle -> LightGray,
PlotLabel->Row[{"Blue is local max, red is local min"}],BaseStyle->12]
]


When f is not linear then NMaximize may return a local maximum.

{tmin, tmax} = {0, 2 Pi};

f = -((6 - 3 Cos[t] -
Cos[3 t])/((-11 + 6 Cos[t] + 8 Cos[2 t] - 6 Cos[3 t] + Cos[4 t]) Sqrt[
Cos[t]^2 + 9 Sin[t]^2 - 12 Cos[t] Sin[t]^2 + 4 Cos[t]^2 Sin[t]^2]));


Find all of the maximum in the interval and select the largest.

max = SortBy[{f /. #, #} & /@
NSolve[{D[f, t] == 0, D[f, {t, 2}] < 0, tmin <= t <= tmax}, t],
First] // Last

(* {5., {t -> 3.14159}} *)

• Thank you for your answer. I now understand that NMaximize is very close to FindMaximum. I wonder why they have different function names... May 27, 2020 at 10:36
• @A.Kato FindMaximum[] does not work as hard as NMaximize[], since it's only intended to find a local optimum near the starting value you gave. By contrast, NMaximize[]` works harder, because it tries (note the emphasis) to find a global optimum. Sometimes it works out, sometimes it doesn't. Here, you were just unlucky. May 27, 2020 at 12:08