# Maximize does not return the global maximum. Why?

Bug introduced in between 10.0.1 and 10.4.1, and fixed in 11.1

I have the following code

Maximize[{-(((1-3 p0^2-3 p1^2-3 p2^2) (18 p0 p1-27 p0^2 p1-27 p0 p1^2+18 p0    p2-27 p0^2 p2+18 p1 p2-27 p1^2 p2-27 p0 p2^2-27 p1 p2^2))/(9 Sqrt[3]))+((-2+9 p0^2-9 p0^3+9 p1^2-9 p1^3+9 p2^2-9 p2^3) (-p1 p2-p0 (p1+p2)))/(3 Sqrt[3]),p0>=0&&p0<=1&&p1>=0&&p1<=1&&p2>=0&&p2<=1&&p0+p1+p2==1&&1/27 (-2+9 p0^2-9 p0^3+9 p1^2-9 p1^3+9 p2^2-9 p2^3)<=0},{p0,p1,p2}]


The result I get is zero. But this is wrong. This triple

p0->0.284088318172241,p1->0.4318233636555181,p2->0.284088318172241

satisfies all the constraints and the objective expression is positive. It equals 0.0000921367. Curiously, if I multiply the objective function by $\sqrt{3}$ I get the presumed global maximum

$-{16\over(-50013-6049\sqrt{69})}>0.$

The documentation to Maximize says If f and cons are linear or polynomial, Maximize will always find a global maximum. But $\sqrt{3}$ is just a number in my expression and it should not matter that it is not rational.

• I get -((22228 Sqrt[3])/1953125) + (24196 Sqrt[23])/5859375, which equals 0.0000921367 approximately. V11.1 Commented Mar 27, 2017 at 0:20
• Hmm, my version is 11.0.1.0. Yes, 0.0000921367 is the answer one should get. Commented Mar 27, 2017 at 0:21
• Quite shocking to be honest that this high in the version numbering an error like this can happen. Commented Mar 27, 2017 at 0:29
• It works fine in V9.0.1. The bug was present in V10.4.1; I can't check the versions in between. Commented Mar 27, 2017 at 0:32
• It works fine in V10.0.1. Commented Mar 27, 2017 at 4:00

Amplifying on comment by Michael E2

\$Version

(*  "11.1.0 for Mac OS X x86 (64-bit) (March 16, 2017)"  *)

sys = {-(((1 - 3 p0^2 - 3 p1^2 - 3 p2^2) (18 p0 p1 - 27 p0^2 p1 -
27 p0 p1^2 + 18 p0 p2 - 27 p0^2 p2 + 18 p1 p2 - 27 p1^2 p2 -
27 p0 p2^2 - 27 p1 p2^2))/(9 Sqrt[3])) + ((-2 + 9 p0^2 - 9 p0^3 +
9 p1^2 - 9 p1^3 + 9 p2^2 - 9 p2^3) (-p1 p2 - p0 (p1 + p2)))/(3 Sqrt[
3]), p0 >= 0 && p0 <= 1 && p1 >= 0 && p1 <= 1 && p2 >= 0 && p2 <= 1 &&
p0 + p1 + p2 == 1 &&
1/27 (-2 + 9 p0^2 - 9 p0^3 + 9 p1^2 - 9 p1^3 + 9 p2^2 - 9 p2^3) <= 0};


The exact solution is

max1 = Maximize[sys, {p0, p1, p2}]

(*  {-((22228 Sqrt[3])/1953125) + (24196 Sqrt[23])/
5859375, {p0 ->
1 - Root[{-23 + #1^2 &, -3 + #2^2 &,
2613168 #1 + 24048128 #2 - 39197520 #1 #3 - 735721920 #2 #3 +
7382812500 #2 #3^2 - 30164062500 #2 #3^3 + 64072265625 #2 #3^4 -
71191406250 #2 #3^5 + 35595703125 #2 #3^6 &}, {2, 2, 2}] -
Root[{-23 + #1^2 &, -3 + #2^2 &,
2613168 #1 + 24048128 #2 - 39197520 #1 #3 - 735721920 #2 #3 +
7382812500 #2 #3^2 - 30164062500 #2 #3^3 + 64072265625 #2 #3^4 -
71191406250 #2 #3^5 + 35595703125 #2 #3^6 &,
72588 #1 - 200052 #2 + 7812500 #2 #3 - 25390625 #2 #3^2 +
35156250 #2 #3^3 - 17578125 #2 #3^4 + 7812500 #2 #4 -
148437500 #2 #3 #4 + 439453125 #2 #3^2 #4 - 562500000 #2 #3^3 #4 +
263671875 #2 #3^4 #4 - 25390625 #2 #4^2 + 439453125 #2 #3 #4^2 -
843750000 #2 #3^2 #4^2 + 527343750 #2 #3^3 #4^2 + 35156250 #2 #4^3 -
562500000 #2 #3 #4^3 + 527343750 #2 #3^2 #4^3 - 17578125 #2 #4^4 +
263671875 #2 #3 #4^4 &}, {2, 2, 2, 1}],
p1 -> Root[{-23 + #1^2 &, -3 + #2^2 &,
2613168 #1 + 24048128 #2 - 39197520 #1 #3 - 735721920 #2 #3 +
7382812500 #2 #3^2 - 30164062500 #2 #3^3 + 64072265625 #2 #3^4 -
71191406250 #2 #3^5 + 35595703125 #2 #3^6 &}, {2, 2, 2}],
p2 -> Root[{-23 + #1^2 &, -3 + #2^2 &,
2613168 #1 + 24048128 #2 - 39197520 #1 #3 - 735721920 #2 #3 +
7382812500 #2 #3^2 - 30164062500 #2 #3^3 + 64072265625 #2 #3^4 -
71191406250 #2 #3^5 + 35595703125 #2 #3^6 &,
72588 #1 - 200052 #2 + 7812500 #2 #3 - 25390625 #2 #3^2 +
35156250 #2 #3^3 - 17578125 #2 #3^4 + 7812500 #2 #4 -
148437500 #2 #3 #4 + 439453125 #2 #3^2 #4 - 562500000 #2 #3^3 #4 +
263671875 #2 #3^4 #4 - 25390625 #2 #4^2 + 439453125 #2 #3 #4^2 -
843750000 #2 #3^2 #4^2 + 527343750 #2 #3^3 #4^2 + 35156250 #2 #4^3 -
562500000 #2 #3 #4^3 + 527343750 #2 #3^2 #4^3 - 17578125 #2 #4^4 +
263671875 #2 #3 #4^4 &}, {2, 2, 2, 1}]}}  *)


Simplifying the forms of p0, p1, and p2

max2 = max1 // RootReduce // ToRadicals

(*  {8/3 Sqrt[2/(837672973 + 100842879 Sqrt[69])], {p0 -> 1/75 (49 - 2 Sqrt[69]),
p1 -> 1/75 (13 + Sqrt[69]), p2 -> 1/75 (13 + Sqrt[69])}}  *)


Verifying that all of the conditions are met

verify = sys /. max2[[2]] // Simplify

(*  {(4 (-50013 + 6049 Sqrt[69]))/(5859375 Sqrt[3]), True}  *)


The form of the maximum given in verify is simpler than the ones in max1 or max2

LeafCount /@ First /@ {max1, max2, verify}

(*  {19, 21, 18}  *)


The different forms of the maximum are equal

max1[[1]] == max2[[1]] == verify[[1]] // FullSimplify

(*  True  *)


Mathematica does not give you the whole truth.

Maximize does not recognize, that the system is totaly symmetric and there are three equal maxima,

corresponding to the six permutations of {p0,p1,p2} divided by two, because of two p's are equal (I corrected my first version).

Redefine sys, calculate the maximum and take the permutations (using MMA 8.0)

    sys[p0_, p1_, p2_] = sys

max1 = Maximize[sys[p0, p1, p2], {p0, p1, p2}]

(*   {-((22228 Sqrt[3])/1953125) + (24196 Sqrt[23])/
5859375, {p0 -> 1 + 2/75 (-13 - Sqrt[69]),
p1 -> 1/75 (13 + Sqrt[69]), p2 -> 1/75 (13 + Sqrt[69])}}    *)

sys[Sequence @@ #] & /@ Permutations[{p0, p1, p2}] /.
max1[[2]] // N

(*   {{0.0000921367, True}, {0.0000921367, True}, {0.0000921367,
True}, {0.0000921367, True}, {0.0000921367, True}, {0.0000921367,
True}}    *)