# NMaximize seems to be returning the wrong result

I'm trying to use NMaximize to show that entropy is maximized when the distribution is uniform. It seems to work when the random variable is valued, but breaks at 4+ variables:

p = {p1, p2, p3, p4};
c = Map[0 <= # <= 1 &, p];
NMaximize[{Sum[-Log[p[[i]]]*p[[i]], {i, 1, Length[p]}], Total[p] == 1,
Splice[c]}, p]


If p is changed to {p1, p2, p3}, then the result is p1 = p2 = p3 = 1/3 as expected. But for 4 probabilities, I get the following:

{1.30778, {p1 -> 0.3323, p2 -> 0.269063, p3 -> 0.0975481,
p4 -> 0.301089}}


Which seems clearly wrong. Am I using NMaximize incorrectly?

• Weirdly this works using the $L_1$ norm i.e Abs[p1] + Abs[p2] + Abs[p3] + Abs[p4] by doing - NMaximize[{Total[-# Log[#] & /@ p], Splice[c], Norm[p, 1] == 1}, p] and gives the result: {1.38629, {p1 -> 0.250001, p2 -> 0.25, p3 -> 0.25, p4 -> 0.25}} but note the slight error of .000001 which violates the constraint. Sep 3 '20 at 23:56
• If you use the "RandomSearch" method then it gives the right result straight away: NMaximize[{Total[-# Log[#] & /@ p], Splice[c], Total[p] == 1}, p, Method -> "RandomSearch"] Sep 4 '20 at 0:02
• Just tell NMaximize about the symmetry NMaximize[{Sum[-Log[p[[i]]]*p[[i]], {i, 1, Length[p]}], Total[p] == 1, p1 == p2 == p3 == p4, c} // Flatten, p]  . This works also with Maximize Maximize[{Sum[-Log[p[[i]]]*p[[i]], {i, 1, Length[p]}], Total[p] == 1, p1 == p2 == p3 == p4, c} // Flatten,  yields {Log, {p1 -> 1/4, p2 -> 1/4, p3 -> 1/4, p4 -> 1/4}}  Sep 4 '20 at 4:25
• @flinty The L1 norm trick doesn't seem to generalize if I add more parameters, but RandomSearch seems to work well consistently. Sep 4 '20 at 4:31
• @Akku14 the symmetry is something I knew about the nature of the solution in this case. It's not going to help me in other situations where the probabilities are going to be different Sep 4 '20 at 4:31

Clear["Global*"]

p = {p1, p2, p3, 1 - p1 - p2 - p3};

sol = Solve[Thread[D[Total[-p*Log[p]], {Most@p}] == 0], Most@p]

(* {{p1 -> 1/4, p2 -> 1/4, p3 -> 1/4}} *)

p /. sol[]

(* {1/4, 1/4, 1/4, 1/4} *)

Total[-p*Log[p]] /. sol[]

(* Log *)


EDIT: For larger number of probabilities

n = 50;

p = pr /@ Range[n - 1];
p = Append[p, 1 - Total[p]];

sol = Solve[Thread[D[Total[-p*Log[p]], {Most@p}] == 0], Most@p];

(p /. sol[] // Union) == {1/n}

(* True *)

(Total[-p*Log[p]] /. sol[]) == Log[n]

(* True *)
`
• That works alright, but it's a little manual. I would prefer to use mathematica's native optimization functions. Sep 4 '20 at 4:35