12
$\begingroup$

I'm trying to solve a maximization problem that apparently is too complicated (it's a convex function) and NMaximize just runs endlessly.

I'd like to have an approximate result, though. How can I tell NMaximize to just give up after $n$ seconds and give me the best it has found so far?

$\endgroup$
  • 2
    $\begingroup$ There's MaxIterations, but it's not a constraint based on time. $\endgroup$ – Michael E2 Nov 28 '16 at 3:06
  • 2
    $\begingroup$ I have often wanted such an option $\endgroup$ – M.R. Nov 28 '16 at 5:22
11
$\begingroup$

This is how I usually deal with this kind of problem. Keywords to the solution are TimeConstrained, AbortProtect, Throw and Catch.

Consider the two target functions:

fun[x_] := (Pause[1]; x^4 - 3 x^2 - x);
fun2[x_, y_] := (Pause[1]; x + y);

Now we define our own optimization routine, that is TimeConstrained by some timeLimit.

timedOptimization[timeLimit_, vars_, target_, constraints___] :=
  Module[{objective, bestVal = 100, symbolicVars = HoldForm /@ vars, tempVal},
   objective[v_ /; VectorQ[v, NumericQ]] := Module[{val}, 
     val = target@@v; If[val < bestVal, bestVal = val; tempVal = ReleaseHold@v;];
     val
     ];
   AbortProtect[Catch@TimeConstrained[
      NMinimize[
       {objective[vars],constraints},
       vars,
       MaxIterations -> 400
      ],
      timeLimit,
      Print["Too slow! Aborting with value ", bestVal, " for parameter ", Thread[symbolicVars -> tempVal]]; 
      Throw[{bestVal, Flatten[Thread[symbolicVars -> tempVal]]}]
      ]
    ]
   ];

This should be quite clear. Mainly, AbortProtect is needed to protect against the Abort generated by TimeConstrained as soon as timeLimit is reached. The current best value and the parameter is Thrown within the fail expression for TimeConstrained and the outer Catch is needed to, well, catch these values.

Due to the BlankNullSequence, constraints is an optional argument which does not need to be specified. varsis a list with all parameters, even if it is a single one. See below examples of how to use with one/multiple arguments and with/without constraints. Please note that the MaxIterations option is actually not needed and just included for demonstration purposes (it exists).

AbsoluteTiming@timedOptimization[10, {x}, fun]

Too slow! Aborting with value -2.8916 for parameter x->0.965034

{10.001186, {-2.8916, x -> 0.965034}}

AbsoluteTiming@timedOptimization[10, {x}, fun, x > 1]

Too slow! Aborting with value -3.24616 for parameter {x->1.09146}

{10.001160, {-3.24616, {x -> 1.09146}}}

AbsoluteTiming@timedOptimization[10, {x,y}, fun2]

Too slow! Aborting with value -0.672799 for parameter {x->-0.535769,y->-0.13703}

{10.000125, {-0.672799, {x -> -0.535769, y -> -0.13703}}}

AbsoluteTiming@timedOptimization[10, {x,y}, fun2, {x, y} \[Element] Disk[]]

Too slow! Aborting with value -0.70015 for parameter {x->-0.351705,y->-0.348445} {10.000125, {-0.70015, {x -> -0.351705, y ->-0.348445}}}

$\endgroup$
  • $\begingroup$ Okay, I like where this is going. How would I add multiple variables and constraints? Also, in the If statement, wouldn't it make more sense to say If[val<bestVal,...]? $\endgroup$ – Ziofil Nov 28 '16 at 15:19
  • $\begingroup$ @Ziofil Right, indeed there was a typo in the If statement. With respect to the rest, please see the updated answer. $\endgroup$ – Lukas Nov 28 '16 at 15:46
  • $\begingroup$ Great! Do you also mean val = target[v] and fun2[{x_,y_}]:=...? $\endgroup$ – Ziofil Nov 28 '16 at 15:52
  • $\begingroup$ @Ziofil Glad I could help. Thanks again for pointing out the mistake again. I should copy-paste the whole notebook after editing it instead of fixing the lines manually here (and forgetting crucial stuff...). The syntax for fun2 was correct as it is, however there was a mistake for val=... in the objective function. I Apply the target function to the list of parameters. To see why, evaluate f@@{a,b,c}. $\endgroup$ – Lukas Nov 28 '16 at 15:58
  • $\begingroup$ Okay. One last weird behaviour: if fun2 is just the sum of the inputs, why in your third and fourth example is the value never the sum of the inputs? $\endgroup$ – Ziofil Nov 28 '16 at 16:04
8
$\begingroup$

Each of the four methods in NMinimize has built-in hooks you can use to get the current values through StepMonitor. There are some great advantages to this approach over ones that hijack the user's objective function and only minor drawbacks:

  • The regular NMinimize interface can be used.
  • The user's objective function will be used as is, so that any analysis of the function normally done by NMinimize is not prevented by wrapping the objective function in shell that turns the function into a numeric black box.
  • The hooks give direct access to the state of the method algorithm being used. One could hardly ask for more.
  • The hooks are easy to access. Some post-processing is often needed, depending on the method. One might want further to apply FindMinimum[] to the results (not shown).
  • The hooks are undocumented AFAIK, and perhaps they are subject to change. OTOH, the code is open to inspection and this approach can be adapted should new methods be added or old methods improved. I believe the current code has been stable for a fairly long time.

Note that NMaximize[f[x], x] basically calls NMinimize[-f[x], x], so I will speak primarily in terms of minimization. The raw values you get with the following approaches will also be of -f[x], when using NMaximize[]. In most examples below, which call NMaximize[], this is accounted for that the maximum is returned.

In NMinimize the objective function consists of the user's function plus a penalty function. There are two values that it keeps track of, val and fval. The value initially optimized is val, which equals fval plus a penalty (often 0.), where fval is the value of the function. At the end of the method, there is post-processing of the results, sometimes using the equivalent of FindMinimum to polish the results.

Each algorithm is different and there is not a uniform user-interface to them. Here are examples of each:

"DifferentialEvolution"

In this example foo contains the most recent pools of points (vecs) and values. One can use linked lists to keep track of each step (see "SimulatedAnnealing" at the end).

(*  "DifferentialEvolution"  *)
TimeConstrained[
 NMaximize[
  {7 x - 4 x^2 + y - y^2, {x, y} ∈ Disk[]}, {x, y},
  Method -> "DifferentialEvolution", 
  StepMonitor :> 
   If[ValueQ@Optimization`NMinimizeDump`fvals, 
    foo = {Optimization`NMinimizeDump`vals, Optimization`NMinimizeDump`vecs, 
      Optimization`NMinimizeDump`fvals}]],
 0.1]

tolerance = 10^-7;
Last@ Sort@ Pick[Transpose[{-foo[[1]], Thread[{x, y} -> #] & /@ foo[[2]]}], 
   UnitStep[foo[[1]] - foo[[3]] - tolerance], 0]
(*
$Aborted
{3.31235, {x -> 0.870649, y -> 0.491204}}
*)

"NelderMead"

Like with "DifferentialEvolution", foo contains the most recent pools of points (vecs, the "simplex") and values.

(*  "NelderMead"  *)
TimeConstrained[
 NMaximize[
  {7 x - 4 x^2 + y - y^2, {x, y} ∈ Disk[]}, {x, y},
  Method -> {"NelderMead",
    "RandomSeed" -> 1 (* for reproducibility *)}, 
  StepMonitor :> 
   If[ValueQ@Optimization`NMinimizeDump`fvals, 
    foo = {Optimization`NMinimizeDump`vals, Optimization`NMinimizeDump`vecs, 
      Optimization`NMinimizeDump`fvals}]],
 0.05]

tolerance = 10^-7;
Last@ Sort@ Pick[Transpose[{-foo[[1]], Thread[{x, y} -> #] & /@ foo[[2]]}], 
   UnitStep[tolerance + foo[[3]] - foo[[1]]], 1]
(*
$Aborted
{3.31236, {x -> 0.871075, y -> 0.491156}}
*)

Note that each value val has a penalty. Hence the need for a positive tolerance, or all points would be rejected. (This would be cleaned up in post-processing, which was aborted by the time constraint.) As one can see below, the selected point above does not satisfy {x, y} ∈ Disk[]. The user will have to decide how to treat the results in their particular case. (This applies to all methods, in fact.)

Norm[{x, y}] /. Last[%] // InputForm
(*  1.0000027734542223  *)

foo[[1]] - foo[[3]]
(*  {5.21808*10^-8, 3.6511*10^-8, 5.48056*10^-8}  *)

"RandomSearch"

In "RandomSearch" results is initialized to the pool of initial points. Each point is replaced by the result of a local minimizer (which may be specified with the "Method" suboption to Method). The post-processing code chooses the best result of the minimizer.

(*  "RandomSearch"  *)
TimeConstrained[
 NMaximize[{7 x - 4 x^2 + y - y^2, {x, y} ∈ Disk[]}, {x, y}, 
  Method -> "RandomSearch", 
  StepMonitor :> 
   If[ValueQ@Optimization`NMinimizeDump`results, 
    foo = Optimization`NMinimizeDump`results]],
 0.2]
Select[foo, ! FreeQ[#, "Converged"] &]
(*
$Aborted
{{{-3.31236, {0.871072 -> 0.871072, 0.491155 -> 0.491155}}, {True, "Converged"}},
 {{-3.31236, {0.871072 -> 0.871072, 0.491155 -> 0.491155}}, {True, "Converged"}},
 {{-3.31236, {0.871072 -> 0.871072, 0.491155 -> 0.491155}}, {True, "Converged"}},
 {{-3.31236, {0.871072 -> 0.871072, 0.491155 -> 0.491155}}, {True, "Converged"}},
 {{-3.31236, {0.871072 -> 0.871072, 0.491155 -> 0.491155}}, {True, "Converged"}},
 {{-3.31236, {0.871072 -> 0.871072, 0.491155 -> 0.491155}}, {True, "Converged"}}}
*)

"SimulatedAnnealing"

Simulated annealing works somewhat like the random search method in that it starts with a pool of initial points and processes each individually. It then does some post-processing and returns the best result found. As it processes each point, it keeps track of the best result of processing that point in the form

 Optimization`NMinimizeDump`best = {val, vec, fval}

To accumulate all the results efficiently, I used link lists. I wrapped the list Optimization`NMinimizeDump`best in an undefined head called hold to make flattening the linked list in post-processing easier. Note it does not Hold[] the results.

(*  "SimulatedAnnealing"  *)
TimeConstrained[
 ClearAll[hold];
 foo = {};
 NMaximize[{7 x - 4 x^2 + y - y^2, {x, y} ∈ Disk[]}, {x, y}, 
  Method -> "SimulatedAnnealing", 
  StepMonitor :> 
   If[ValueQ@Optimization`NMinimizeDump`best, 
    foo = {hold[Optimization`NMinimizeDump`best], foo}]],
 0.02]
foo = Flatten@foo /. hold -> Identity;

tolerance = 0;
First@ Sort@ Pick[
   Transpose[{-foo[[All, 1]], Thread[{x, y} -> #] & /@ foo[[All, 2]]}], 
   UnitStep[tolerance + foo[[All, 3]] - foo[[All, 1]]], 1]
(*
$Aborted
{3.25174, {x -> 0.753606, y -> 0.457427}}
*)
$\endgroup$
  • $\begingroup$ I ended using the other solution, but I'm also playing around with your code. At the moment, for some reason, foo is not assigned a value when I use my objective function. I'll keep playing. $\endgroup$ – Ziofil Dec 1 '16 at 14:16
  • $\begingroup$ @Ziofil Can you tell me the function? Or you might try StepMonitor without the If[ValueQ@.., part. NMinimize does one symbolic evaluation for some reason, which messed up the post-processing. If foo is not assigned a value, it seems to imply StepMonitor is not evaluated. $\endgroup$ – Michael E2 Dec 1 '16 at 16:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.