# NMinimize go out of bounds on some linear constraints

I'm trying to minimize a non-linear function of four variables with some linear constraints. NMinimize on Mathematica 8 seems to go out of bounds on some constraints giving complex values of the function at some point in the iteration. The function to minimize is

ff[lxw_, lwz_, c_, d_] := - J1 (lxw + lwz) - 2 J2 c +
T (-Log[2] - 1/2 (1 - lxw) Log[(1 - lxw)/4] -
1/2 (1 + lxw) Log[(1 + lxw)/4] -
1/2 (1 - lwz) Log[(1 - lwz)/4] -
1/2 (1 + lwz) Log[(1 + lwz)/4] + 1/2 (1 - d) Log[(1 - d)/16] +
1/8 (1 + 2 c + d - 2 lwz - 2 lxw) Log[
1/16 (1 + 2 c + d - 2 lwz - 2 lxw)])


where

T = 10;
J1 = 1;
J2 = -0.2;


are constant parameters. The constraints are

var = {lxw, lwz, c, d};
cons = And @@ Cases[ff @@ var, Log[x_] -> x > 0, Infinity] // Simplify


which amounts to

d < 1 && lwz < 1 && 1 + lwz > 0 && 1 + 2 c + d > 2 (lwz + lxw) &&
lxw < 1 && 1 + lxw > 0


Then

NMinimize[{ff @@ var, cons}, var, Method -> "DifferentialEvolution"]


gives

NMinimize::nrnum: The function value 7.60939\[VeryThinSpace]-3.87314 I
is not a real number at {c,d,lwz,lxw} =  {-0.267966,0.319033,0.899803,-0.0151082}. >>
{-12.5741, {c -> -0.236255, d -> -0.978425, lwz -> -0.681637,
lxw -> -0.939595}}

-
I suspect you want to get rid of the complex value and solve the problem ? If not - add more details. –  Sektor Feb 26 at 22:53
@belisarius No, you probably made the same mistake I did of using Thread[var->theBadSet]. The message states that it is using a reordering of var and for that the value is indeed complex. –  Daniel Lichtblau Feb 27 at 0:07
@DanielLichtblau Epic failure. I never check that! Thanks a lot. –  belisarius Feb 27 at 0:21
C.B. DeMille would never have made that mistake ("epic failure" indeed). –  Daniel Lichtblau Feb 27 at 0:26
@Sektor the problem is that the constraints, which are given in line below which amounts to , are not satisfied, i.e. at some point in the iteration some logarithm is evaluated with a negative argument, and this gives the complex result. If the constraints would be satisfied the computation should proceed with strictly positive arguments in all logarithms in the function to be minimized. Btw, I don't think there is any reordering problem, as stated by Daniel Lichblau . –  danielstariolo Feb 27 at 12:05

There are two things that can be done to improve matters. One is to use constraints with weak inequalities that keep the log arguments above a minmium threshold. The other is to use an altered logarithm that allows argumens less-equal to zero and simply returns a suitable large negative.

T = 10;
J1 = 1;
J2 = -0.2;
bigValue = 10^6;
myLog[x_?NumberQ] := If[x <= 0, -bigValue, Log[x]]
vars = {lxw, lwz, c, d};
ff[lxw_, lwz_, c_,
d_] := -J1 (lxw + lwz) - 2 J2 c +
T (-Log[2] - 1/2 (1 - lxw) Log[(1 - lxw)/4] -
1/2 (1 + lxw) Log[(1 + lxw)/4] -
1/2 (1 - lwz) Log[(1 - lwz)/4] -
1/2 (1 + lwz) Log[(1 + lwz)/4] + 1/2 (1 - d) Log[(1 - d)/16] +
1/8 (1 + 2 c + d - 2 lwz - 2 lxw) Log[
1/16 (1 + 2 c + d - 2 lwz - 2 lxw)]) /. Log -> myLog
eps = 1/1000;


So the constraints are now as follows.

cons = And @@ Cases[ff @@ var, myLog[x_] :> x >= eps, Infinity];

(* ut[57]= (1 - d)/16 >= 1/1000 && (1 - lwz)/4 >= 1/1000 && (1 + lwz)/
4 >= 1/1000 &&
1/16 (1 + 2 c + d - 2 lwz - 2 lxw) >= 1/1000 && (1 - lxw)/4 >= 1/
1000 && (1 + lxw)/4 >= 1/1000 *)


Here is the optimization.

{min, vals} =
NMinimize[Evaluate[{ff @@ vars, cons}], vars,
Method -> "DifferentialEvolution"]

(* Out[55]= {-27.7261, {lxw -> -0.996, lwz -> 0.996, c -> 4.33532,
d -> -4.65512}} *)


Notice that the variables lxw and lwz really want to live at the boundary of the allowed region. This makes me suspect there is an issue with the forumlation, that maybe it's not solving the problem you want it to handle? Or maybe you want to let them be -1 and 1 respectively, and remove terms with logs that would then contain zero? This seems plausible since they have factors that are also zero and, in a limiting sense the values "should" be zero. When I make those adjustments, (skipping the code changes and going right to the results) I get

{min, vals} =
NMinimize[Evaluate[{ff @@ vars, cons}], vars,
Method -> "DifferentialEvolution"]

(* Out[99]= {-28.0146, {c -> 4.33553, d -> -4.65527}} *)

-
thank you very much for your help, indeed now the algorithm is returning the results I expected. The proposed workaround is very good, simply changing a strong constraint with a penalty in the function which is being minimized, nice trick. –  danielstariolo Feb 28 at 13:57