NMinimize seems to be on crack

Question

Say I have a crazy function, f, defined like so:

util[x_, y_, c_] := 0.5*Log[c-x] + 0.5*Log[c-y]
cost[x_, y_, l_] := c /. First[NSolve[util[x, y, c+l] == Log[10+l], c]]
prof[x_, y_]   := 0.01*Norm[{x,y}, 2]
liquid[x_, y_] := 0.01*Norm[{x,y}, 2]
f[x_, y_, a_, b_] := cost[a, b, liquid[x,y] + liquid[a-x, b-y]] - Max[a,b] 
      - cost[0,0,0] + prof[x,y] + liquid[x,y] + prof[a-x, b-y] + liquid[a-x, b-y]

Now I call NMinimize like this:

NMinimize[{f[50, 50, k, j], k >= 49, k <= 51, j >= 49, j <= 51}, {j, k}]

Which tells me this:

{-21.0465, {j -> 51., k -> 49.}}

But then if I actually check what f[50,50,49,51] is, it's this:

0.489033

Which is pretty different from the -21.0465 that NMinimize said. Is this par for the course with NMinimize? Floating point errors compounding or whatnot? Any ideas for beating NMinimize (or some such function) into submission?

Answer 1

It certainly seems to be related to your function f not being restricted to numerical arguments, plus the symbolic preprocessing performed by NMinimize. Once you change the signature to

f[x_?NumericQ, y_?NumericQ, a_?NumericQ, b_?NumericQ]:=...

The result is as expected, although it takes considerably longer to get it.

EDIT

We can dig deeper to reveal the true reason. First, note that your f (the original one, args unrestricted) is quite a function:

In[1423]:= f[50,50,49.,51.]
Out[1423]= 0.489033

In[1392]:= f[50,50,k,j]/.{j->51.`,k->49.`}
Out[1392]= -21.0465

The real culprit is NSolve, which gives two ordered solutions:

In[1398]:= NSolve[util[x,y,c+l]==Log[10+l],c]
Out[1398]= {{c->0.5 (-2. l+1. x+1. y-2. Sqrt[100.+20. l+1. l^2+0.25 x^2-0.5 x y+0.25 y^2])},
 {c->0.5 (-2. l+1. x+1. y+2. Sqrt[100.+20. l+1. l^2+0.25 x^2-0.5 x y+0.25 y^2])}}

The problem is, what is the ordering. It turns out to be different for symbolic and numeric arguments to NSolve, because in the latter case we don't have any symbols around. This can be seen as:

In[1399]:= 
Block[{cost},
   cost[x_,y_,l_]:=c/.Last[NSolve[util[x,y,c+l]==Log[10+l],c]];
   f[50,50,k,j]/.{j->51.,k->49.}]

Out[1399]= 0.489033

So you really have to settle on what is the right ordering for you, and which solution you really want to pick.