# 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. – J. M.'s technical difficulties May 27 at 2:30
• From the docs: "Otherwise, NMaximize may sometimes find only a local maximum." – Michael E2 May 27 at 2:31
• In this case, Maximize[{f, tmin <= t <= tmax}, t] works. – Michael E2 May 27 at 2:32
• @J. M., I see. Ill remove the tag "bug". – A. Kato May 27 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. – Michael E2 May 27 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,
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. – A. Kato May 27 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. – A. Kato May 27 at 10:30
• @A.Kato: Unfortunately,Method->"DifferentialEvolution" does not work here. – user64494 May 27 at 11:29

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]
] 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[]. – J. M.'s technical difficulties May 27 at 12:32

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... – A. Kato May 27 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. – J. M.'s technical difficulties May 27 at 12:08