# 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, 2020 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, 2020 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[4], {p1 -> 1/4, p2 -> 1/4, p3 -> 1/4, p4 -> 1/4}}  Sep 4, 2020 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, 2020 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, 2020 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]]

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

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

(* Log[4] *)


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[[1]] // Union) == {1/n}

(* True *)

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

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