Time constrained optimization?

Question

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?

Answer 1

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}}}

Answer 2

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}}
*)