# Maximize simple polynomial: Wrong answer (Gives the minimum!)

I want to maximize poly under the constraints $g1,g2,g3=0$

poly = x1^4 + 2 x1^2 x2^2 + 2 x1^2 x3^2 + 2 x1^2 x4^2 + 0.2 x2^4 +
0.4 x2^2 x3^2 + 0.4 x2^2 x4^2 + 0.2 x3^4 + 0.4 x3^2 x4^2 + 0.2 x4^4;

g1 = 6.4 x1 x2^3 + 6.4 x1 x2 x3^2 + 6.4 x1 x2 x4^2;
g2 = 6.4 x1 x2^2 x3 + 6.4 x1 x3^3 + 6.4 x1 x3 x4^2;
g3 = -1 + x1^2 + x2^2 + x3^2 + x4^2;


Therefore I used all these different methods, but none of them work.

Maximize[{poly, g1 == 0, g2 == 0, g3 == 0}, {x1, x2, x3, x4}]


{0.2, {x1 -> -1.11978*10^-9, x2 -> -0.20688, x3 -> 0.930234, x4 -> -0.303094}}

NMaximize[{poly, g1 == 0, g2 == 0, g3 == 0}, {x1, x2, x3, x4}]


{0.2, {x1 -> -1.11978*10^-9, x2 -> -0.20688, x3 -> 0.930234, x4 -> -0.303094}}

FindMaximum[{poly, g1 == 0, g2 == 0, g3 == 0}, {x1, x2, x3, x4}]


{0.2, {x1 -> -4.20087*10^-11, x2 -> -0.705858, x3 -> -0.705858, x4 -> -0.059409}}

FindMaxValue[{poly, g1 == 0, g2 == 0, g3 == 0}, {x1, x2, x3, x4}]


0.2

But there is the obvious solution $\vec{x}=({1,0,0,0})$

poly /. x1 -> 1 /. x2 -> 0 /. x3 -> 0 /. x4 -> 0


1.

which also fulfills the constraints:

g1 /. x1 -> 1 /. x2 -> 0 /. x3 -> 0 /. x4 -> 0
g2 /. x1 -> 1 /. x2 -> 0 /. x3 -> 0 /. x4 -> 0
g3 /. x1 -> 1 /. x2 -> 0 /. x3 -> 0 /. x4 -> 0


0.

0.

0.

Actually the value 0.2 is the global minimum of the polynomial under the constraints...

• You use (inexact) machine precision numbers, which makes Maximize switch over to numeric mode. Try Rationalize on your input (which delivers {1, {x1 -> -1, x2 -> 0, x3 -> 0, x4 -> 0}}). – Yves Klett Nov 24 '14 at 17:30
• Straight from the help: "If Maximize is given an expression containing approximate numbers, it automatically calls NMaximize." – Yves Klett Nov 24 '14 at 17:38
• @Silvia hehe we had parallel evolution - just wanted to augment on Method -> "DifferentialEvolution" ;-). NMaximize being stuck at the minimum is a surpring coincidence though. – Yves Klett Nov 24 '14 at 17:44
• @YvesKlett Still I feel surprised the default setting gives a global minimum. (Or does it?) – Silvia Nov 24 '14 at 17:46
• I would do it this way: Maximize[Rationalize[{poly, g1 == 0, g2 == 0, g3 == 0}], {x1, x2, x3, x4}]. NMaximize sometimes yields wrong answers when working with constraints, see e.g. Maximize violating constraints – Artes Nov 24 '14 at 17:51

NMaximize[{poly, g1 == 0, g2 == 0, g3 == 0}, {x1, x2, x3, x4}, Method -> "DifferentialEvolution"]

(* {1., {x1 -> 0.999999, x2 -> 0.0000693388, x3 -> -0.00106407, x4 -> 0.00102883}}*)


or

NMaximize[{poly, g1 == 0, g2 == 0, g3 == 0}, {x1, x2, x3, x4}, WorkingPrecision -> 15]
(* {1.00000006317100, {x1 -> -0.999999798232786,
x2 -> 0.0000936576222886767, x3 -> -0.000460488000721031,
x4 -> 0.000462924413626347}} *)


Curiously:

sol = FindInstance[ poly == 1 && And @@ {g1 == 0, g2 == 0, g3 == 0}, {x1, x2, x3, x4}, Reals, 10]

(* {{x1 -> -1., x2 -> 0, x3 -> 0, x4 -> 0}, {x1 -> 1., x2 -> 0, x3 -> 0,  x4 -> 0}}*)

{poly /. sol, {g1, g2, g3} /. sol}

(* {{1., 1.}, {{0., 0., 0.}, {0., 0., 0.}}} *)


Here's something I noticed.

The maximization works if you use exact rational numbers instead of real numbers as the coefficients:

poly = x1^4 + 2 x1^2 x2^2 + 2 x1^2 x3^2 + 2 x1^2 x4^2 + 2 /10 x2^4 +
4/10 x2^2 x3^2 + 4/10 x2^2 x4^2 + 2/10 x3^4 + 4/10 x3^2 x4^2 +
2/10 x4^4;
g1 = 64/10 x1 x2^3 + 64/10 x1 x2 x3^2 + 64/10 x1 x2 x4^2;
g2 = 64/10 x1 x2^2 x3 + 64/10 x1 x3^3 + 64/10 x1 x3 x4^2;
g3 = -1 + x1^2 + x2^2 + x3^2 + x4^2;


And the result you get:

In:= Maximize[{poly, {g1 == 0, g2 == 0, g3 == 0}}, {x1, x2, x3, x4}]
Out= {1, {x1 -> -1, x2 -> 0, x3 -> 0, x4 -> 0}}


I'm afraid I can't explain the wrong answer if you have real-valued coefficients, though...