# An issue concerning the NMaximize command in Mathematica

I have two sets. One set consists of constants and the other one consists of expressions depending on two variables "a" and "b". I find absolute values of correlations between the aforementioned two sets of values. I get the maximum value of 0.96707 among these absolute values of correlations using the NMaximize command (utilizing the method "DifferentialEvolution") over -10<= a <= 10 and -10 <= b <= 10:

Input:  NMaximize[{RealAbs[Correlation[{12^(2 a + 3 b), 13^(1 + 2 a + 3 b), 2^(2 a + b), 7^(3 a + b), 11^(4 a + b), 13^(a + 3 b)}, {-49.82, -51.5, -50.82, -50.69, -50.4, -53.71}]], -10<= a <= 10 && -10 <= b <= 10}, {a, b}, Method -> "DifferentialEvolution"]

Output: {0.96707, {a -> -1.43894, b -> 10.}}


However, if the constraint "-10 <= a <= 10 && -10 <= b <= 10" is replaced with "-3 <= a <= 3 && -3 <= b <= 3", the correlation increases: that is, the code

NMaximize[{RealAbs[Correlation[{12^(2 a + 3 b), 13^(1 + 2 a + 3 b), 2^(2 a + b), 7^(3 a + b), 11^(4 a + b), 13^(a + 3 b)}, {-49.82, -51.5, -50.82, -50.69, -50.4, -53.71}]], -3 <= a <= 3 && -3 <= b <= 3}, {a, b}, Method -> "DifferentialEvolution"]


gives output

{0.972989, {a -> -1.41882, b -> 0.504752}}


Certainly, the maximum value in the second case should be less than or equal to 0.96707 (the one obtained in the first case) because of the constraints. I don't know where I am making a mistake. Any help will be much appreciated.

• These methods "RandomSearch","DifferentialEvolution", "SimulatedAnnealing" find the global maximum, which isn't unique in your example, too! Commented Feb 21 at 9:53

Your problem seems that there are 2 maxima close together, One being the global, the other only a local maximum.

If you want to use the method of "DifferentialEvolution" you have to look what is does, e.g. "tutorial/ConstrainedOptimizationGlobalNumerical". You see it uses a population of random points. Now if you have 2 maxima close together, the chance of getting the right one increases with more sample points. Therefore add the option: "Method -> {"DifferentialEvolution", "SearchPoints" -> 250}":

NMaximize[{RealAbs[
Correlation[{12^(2  a + 3  b), 13^(1 + 2  a + 3  b), 2^(2  a + b),
7^(3  a + b), 11^(4  a + b),
13^(a + 3  b)}, {-49.82, -51.5, -50.82, -50.69, -50.4, \
-53.71}]], -10 <= a <= 10 && -10 <= b <= 10}, {a, b},
Method -> {"DifferentialEvolution", "SearchPoints" -> 250}]

{0.972989, {a -> -1.41882, b -> 0.504752}}

• The result of Plot3D[RealAbs[ Correlation[{12^(2 a + 3 b), 13^(1 + 2 a + 3 b), 2^(2 a + b), 7^(3 a + b), 11^(4 a + b), 13^(a + 3 b)}, {-49.82, -51.5, -50.82, -50.69, -50.4, 53.71}]], {a, -10, 10}, {b, -10, 10}, PlotPoints -> 100] shows the real state of things, not your speculations. Commented Feb 21 at 11:24

Similar questions were asked and answered several times. Up to the documentation,

If f is linear or concave and cons are linear or convex, the result given by NMaximize will be the global maximum, over both real and integer values; otherwise, the result may sometimes only be a local maximum

Making use of another options, we obtain the same results.

NMaximize[{RealAbs[ Correlation[{12^(2  a + 3  b), 13^(1 + 2  a + 3  b), 2^(2  a + b),
7^(3  a + b), 11^(4  a + b),
13^(a + 3  b)}, {-49.82, -51.5, -50.82, -50.69, -50.4, -53.71}]],
-10 <= a <= 10 && -10 <= b <= 10}, {a, b}, Method -> {"RandomSearch", "SearchPoints" -> 250}]


{{0.972989, {a -> -1.41882, b -> 0.504752}}}

NMaximize[{RealAbs[ Correlation[{12^(2  a + 3  b), 13^(1 + 2  a + 3  b), 2^(2  a + b),
7^(3  a + b), 11^(4  a + b),
13^(a + 3  b)}, {-49.82, -51.5, -50.82, -50.69, -50.4, -53.71}]],
-3 <= a <= 3 && -3 <= b <= 3}, {a, b}, Method -> {"RandomSearch", "SearchPoints" -> 250}]


{0.972989, {a -> -1.41882, b -> 0.504752}}

Plot3D[RealAbs[Correlation[{12^(2   a + 3   b), 13^(1 + 2   a + 3   b),
2^(2   a + b), 7^(3   a + b), 11^(4   a + b),
13^(a + 3   b)}, {-49.82, -51.5, -50.82, -50.69, -50.4,

• 53.71}]], {a, -10, 10}, {b, -10, 10}, PlotPoints -> 100]

shows the real state of things.

Edit. A typo: 53.71 is relaced by -53.71. This produces the same results of NMaximize, but not in Plot3D.

• Response to user64494: Thanks for your kind help. Is there any method that gives the global maximum in every case; for example, if I change the numbers in the constraints given above, or if I add additional elements (of the previous type) in the considered two sets?
– Tona
Commented Feb 21 at 7:51
• @Tona: It is difficult to answer your request without a changed data. It's clear there is no method to find the global maximum "for all seasons". Commented Feb 21 at 8:03
• Thanks for your kind response. What about the above example? Does your modified code give the global maximum for this example?
– Tona
Commented Feb 21 at 9:40
• @Tona: See the plot in the addition to my answer concerning the real state of things. Commented Feb 21 at 11:27

To long for a comment:

The "real state of things" (thanks at @user64494) follows with rationalized function parameters

fun[a_, b_] :=
RealAbs[Correlation[{12^(2 a + 3 b), 13^(1 + 2 a + 3 b), 2^(2 a + b),
7^(3 a + b), 11^(4 a + b),
13^(a + 3 b)}, {-49.82, -51.5, -50.82, -50.69, -50.4, -53.71} //
Rationalize[#, 0] &]]
Plot3D[fun[a, b], {a, -10, 10}, {b, -10, 10}, PlotPoints -> 100,MaxRecursion -> 4]


plot shows a horizontal line (edge) at the maximum. Global maximum isn't unique as already correctly suspected by @DanielHuber!

Solution(three different methods) for parameter ranges -3<a<3,-3<b<3

pts = Map[(max =
NMaximize[{fun[a, b], -3 <= a <= 3 && -3 <= b <= 3}, {a, b},Method -> {# }];
Join[{a, b} /. max[[2]], {max[[1]]}]) &,
{"DifferentialEvolution","SimulatedAnnealing", "RandomSearch"}]
(*{{-1.38247, 3., 0.954796}, {-1.37254, 0.557541, 0.95574},{-1.37254,0.557543, 0.95574}}*)

Show[Plot3D[fun[a, b], {a, -10, 10}, {b, -10, 10},MeshFunctions -> (#3 &) , PlotPoints -> 100], Graphics3D[{Blue, PointSize[Large], Point[pts]}]]


Solution(three different methods) for parameter ranges -10<a<10,-10<b<10

pts = Map[(max =
NMaximize[{fun[a, b], -10 <= a <= 10 && -10 <= b <= 10}, {a, b},Method -> {# }];
Join[{a, b} /. max[[2]], {max[[1]]}]) &,
{"DifferentialEvolution","SimulatedAnnealing", "RandomSearch"}]
(*{{-1.37423, 10., 0.958818}, {-1.37423, 10., 0.958818}, {-1.37423,10.,0.958818}}*)

Show[Plot3D[fun[a, b], {a, -10, 10}, {b, -10, 10},MeshFunctions -> (#3 &) , PlotPoints -> 100], Graphics3D[{Blue, PointSize[Large], Point[pts]}]]


• Are you sure in the plot? The result of Plot3D[Correlation[{12^(2 a + 3 b), 13^(1 + 2 a + 3 b), 2^(2 a + b), 7^(3 a + b), 11^(4 a + b), 13^(a + 3 b)}, Rationalize[{-49.82, -51.5, -50.82, -50.69, -50.4, 53.71}, 0]], {a, -10, 10}, {b, -10, 10}, PlotPoints -> 100,MaxRecursion->4] is the same as in the addition to my answer. Commented Feb 21 at 12:23
• Yes I'm sure. Look at the edge of the max plateau which is different from your plot Commented Feb 21 at 12:25
• You are right. I typed 53.71 instead of -53.71. Commented Feb 21 at 12:30
• NMaximize[{RealAbs[ Correlation[{12^(2 a + 3 b), 13^(1 + 2 a + 3 b), 2^(2 a + b), 7^(3 a + b), 11^(4 a + b), 13^(a + 3 b)},Rationalize[ {-49.82, -51.5, -50.82, -50.69, -50.4, \ -53.71},0]]], -10 <= a <= 10 && -10 <= b <= 10}, {a, b}, Method -> {"RandomSearch", "SearchPoints" -> 250}] results in {0.972989,{a->-1.41882,b->0.504752}} which is better than your 0.958818. Commented Feb 21 at 17:48
• @user64494 I corrected my function definition fun[a_,b_] to Rationalize[#,0]& (instead of Rationalize[#,9]& which now gives 0.96707 . Value is a bit smaller than your maximum, perhaps because I didn't use the option "SearchPoints"->250 Commented Feb 21 at 18:31