NMinimize usage

Question

I would like to minimize functions including the following.

NMinimize[{1 - (1 - 1/n)^x - x/n, n > x, x > 1}, {n, x}]

However Mathematica complains about values not being real. How should I do this?

Answer 1

There are many techniques to optimize objectives with respect to constraints. Computing and enforcing constraints can be expensive. "Penalty methods" instead allow the constraints to be violated, but add largish values to the objective so that it increases quickly for arguments beyond the constraint boundaries. Clever methods of increasing the penalties can produce a sequence of solutions that converge to one that (just barely) satisfies the constraints.

As evidence that Mathematica is using such methods, read more deeply into the error message:

f[n_, x_] := 1 - (1 - 1/n)^x - x/n;
NMinimize[{f[n, x], n > x, x > 1}, {n, x}]

The function value $-0.00737136+0.00079843\ I$ is not a real number at $\ \{n,x\} = \{0.960381,1.00629\}$.

Notice that the argument values supplied do not satisfy all the constraints: although indeed $x\gt 1$, it is not the case that $n=0.960381 \gt 1.00629=x$.

In short, you cannot assume the constraints will hold during the search for an optimum.

A principle in successfully optimizing functions is to make them as smooth and convex as possible. This is not a place for discussing these issues--they would (and have) filled entire books. Suffice it to say that to the extent possible you want to avoid imposing constraints in the form of absolute values (which are not differentiable everywhere) but instead use squares (which are differentiable).

This question nicely illustrates the application of this principle (and a few others we will uncover as we go along). Motivated by the desire to re-express the variables in smooth convex ways, to enforce $x\ge 1$, we might choose to write $$x = (1+y)^2$$ (with no constraint on $y$) and to enforce $n\ge x$, write $$n = m^2 + x = m^2 + (1 + y)^2$$ (again with no constraint on $m$). (There are other ways to do this, but these two expressions are straightforward and turn out to work well.) We therefore attempt an unconstrained version of the original problem in this form:

replacements = {n -> m^2 + (1 + y)^2, x -> (1 + y)^2}
sol = NMinimize[f[n, x] /. replacements , {m, y}])

$\{-0.367879,\{m\to -0.194048,y\to 19133.8\}\}$

This is quick: it takes less than $0.04$ seconds here, compared to around a half second for the constrained version of the same problem shown below.

If you are interested in the arguments where the objective is optimized, solve for them:

Solve[Equal @@@ replacements /. Last@sol, {x, n}]

$\{\{x\to 3.66139\times 10^8,n\to 3.66139\times 10^8\}\}$

---

This problem happens to be a difficult one: the surface is very flat near the optimum. Accordingly, the exact location of the optimum is uncertain and different solution methods can produce wildly different answers. Compare with this one, for instance, which cleverly sidesteps the problem by allowing complex values:

alt = NMinimize[{Re[1 - (1 - 1/n)^x - x/n], Im[1 - (1 - 1/n)^x - x/n] == 0 
        && n > x && x > 1}, {x, n}]

$\{-0.367879,\{x\to 2.31153\times 10^7,n\to 2.31153\times 10^7\}\}$

Although nominally it produces the same value at the optimum, the location of the optimum is far from that previously produced: these values of $x$ and $n$ are less than one-tenth the previous ones! Which solution is better? Let's compare:

First@alt - First@sol

$1.48898\times 10^{-8}$

The difference is slight but real: the first solution is a better minimum. I do not claim it is the best! It is possible that increasing the precision of the calculation will further change the location of the minimum. What I am showing, though, is that often (not always) following these principles, where applicable, can lead to better solutions that are obtained in substantially less time:

1. Create an objective that is well-defined even beyond the constraints,

2. Use smooth functions to express the objective and the constraints, and

3. Replace constraints altogether through suitable smooth re-expressions of the variables.

---

Incidentally, one reason this is a hard problem is that the minimum is never actually achieved. The huge values of $x$ and $n$ returned by the various numerical minimization attempts attest to that. Noticing, though, that usually $x=n$ at a minimum, we might care to explore the values of $f$ along this line:

Limit[f[n, n], n -> Infinity]

$-\frac{1}{e}$

One way to handle situations like this is to re-express the variables so that they lie within a compact domain. For instance, we could let both $x$ and $n$ be the tangents of some angle between $0$ and $\pi/2$. (Obviously values near $\pi/2$ correspond to near-infinity tangents.) Let's try plotting the logarithm of the (negative) of $f$, because the values of $f$ can vary so much. By starting at angles of $\pi/4$, whose tangent is $1$, and constraining the plot to show only the region where the (angle for) $n$ exceeds (the angle for) $x$, we can survey the "landscape" of $f$ subject to its constraints in one glance:

ContourPlot[Log[-f[Tan[n], Tan[x]]], {n, \[Pi]/4, \[Pi]/2}, {x, \[Pi]/4, \[Pi]/2},
  RegionFunction -> Function[{n, x, z}, n > x]]

It's clear (from the fact that lighter colors represent more negative values of $f$ itself) that all the action (concerning a minimum) is happening along the diagonal. This time let's plot the values of $f$ itself, because there's not much variation evident along the diagonal:

Plot[f[Tan[x], Tan[x]], {x, \[Pi]/4, \[Pi]/2}, 
  PlotStyle -> Thick, AxesOrigin -> {\[Pi]/4, 0}, AxesLabel -> {x, f[z, z]}]

(For brevity, the vertical axis label uses "$z$" to refer to $\tan(x)$.)

Yep, the global minimum is approached in the limit as $n\to\infty$ and $x\to\infty$ in such a way that $n$ and $x$ stay close to one another, and must equal $\exp(-1)$. This isn't difficult to prove, because $f$ is relatively simple. In more complicated situations we often don't have the luxury of being able to analyze $f$ so well and have to be content with the clues offered by our numerical optimizations and our plots.

Answer 2

This

f[x_, n_] := 1 - (1 - 1/n)^x - x/n

NMinimize[{Re[1 - (1 - 1/n)^x - x/n], 
  Im[1 - (1 - 1/n)^x - x/n] == 0 && n > x && x > 1}, {x, n}]

does give an answer (I minimize the real part, constraining the imaginary part to vanish).

Plot3D[f[x, n], {x, 1, 100}, {n, x, 100}]

I'd say you need to re-think about your mathematical problem.

Answer 3

The error message you get is

NMinimize::nrnum: The function value -0.00737136+0.00079843\ I is not a real number at \ {n,x} = {0.960381,1.00629}.

Note that the constraint n > x && x > 1 is not strictly being observed. One can be more explicit and avoid that error:

f = 1 - (1 - 1/n)^x - x/n;
NMinimize[{f, n > x && x > 1 && n > 1}, {n, x}]
(* {-0.367879, {n -> 1.05422*10^7, x -> 1.05422*10^7}} *)

But you get more messages that indicate a problem:

NMinimize::cvdiv: Failed to converge to a solution. The function may be unbounded.

NMinimize::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations.

---

A plot of the function indicates the minimum might occur along the boundary of the region n > x && x > 1:

Plot3D[f, {x, 1, 10}, {n, 1, 10}, MeshFunctions -> {#1 - #2 &}, 
 Mesh -> {{0}}, MeshShading -> {Yellow, Opacity[0.5]}, 
 PlotRange -> {-1, 0}]

Along the boundary n == x, one gets a limit of -1/E == -0.367879....

Limit[f /. n -> x, x -> Infinity]
(* -1/E *)
N[-1/E]
(* -0.367879 *)

Also, the derivative along the boundary can be shown to be negative:

df = D[f /. n -> x, x];
FullSimplify[df < 0, n > x && x > 1 && n > 1]
(* 1 + (-1 + x) Log[(-1 + x)/x] > 0 *)

and the last inequality is easily shown to be true for x > 1. In other words, the function is decreasing along the boundary x -- n.

Finally, the Hessian, that is, the determinant of the second derivative (matrix of second partials), appears to be always negative. This would imply there are no local extrema away from the boundaries (the surface would have negative curvature; any critical point would be a saddle). Below we show a plot of Boole[hf < 0], under a coordinate transformation {x -> 1/(1 - s), n -> 1/(1 - t)} that remaps the intervals from 1 to Infinity to the interval from 0 to 1.

hf = Det@D[D[f, {{x, n}}], {{x, n}}];
Plot3D[Evaluate[Boole[hf < 0] /. {x -> 1/(1 - s), n -> 1/(1 - t)}],
 {s, 0, 1}, {t, 0, 1}, MeshFunctions -> {#1 - #2 &}, Mesh -> {{0}}, 
 MeshShading -> {Yellow, Opacity[0.5]}]

The inequality is bit too hard for me to solve quickly, or perhaps at all.

FullSimplify[hf < 0 /. n -> 1/(1 - y),  1/(1 - y) > x && x > 1]
(* x y^x (-2 y^2 + y^x (1 + y)) Log[y]^2 >
   (-1 + y) (-y + y^x) (-y + y^x + 2 x y^x Log[y]) *)

So it appears that the has no minima in the region, but approaches its greatest lower bound along the boundary x == n as x, n approach infinity.

Answer 4

There are a couple of problems: one sytax and one conceptual. First, here's one interpretation of the problem that does give an answer:

NMinimize[{1 - Abs[(1 - 1/n)]^x - x/n, n > x && x > 1}, {n, x}] 

with answer:

{-0.3678, {n -> 8.1993*^7, x -> 8.1993*^7}}

You can see the syntax error: comma instead of && signs. But the major problem is the raising to a fractional power can give complex results. Putting the Abs[] around this solves that.

Thinking about this a bit more, you can recognize the answer is nothing more than a numerical approximation to -1/Exp[1].