An issue concerning the NMaximize command in Mathematica

Question

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.

Answer 1

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}}

Answer 2

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}}

Addition. 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.

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

Answer 3

To long for a comment:

mutated into an answer

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]}]]