How to obtain better speed performance with NMinimize over a complicated objective function?

Question

Edit.1:

I've rewritten the question in order to make it more accessible for troubleshooting / debugging / maintenance purposes, as suggested in the comments. I apologize for any inconvenience.

I'm trying to find the best way to write the objective function so as to minimize the time it takes NMinimize to solve my constraint minimization problem.

The particulars of the problem are presented below.

I have tried defining the objective function using three different ways which I originally thought that would produce progressively better time gains. Unfortunately, I was mistaken.

The first two definitions use SetDelayed and Function while the last one uses a combination of Function and Compile.

Quoting from the original question:

How can I improve on either of the three versions or how to adopt a whole different approach in order to obtain sensible timing for the optimization problem at hand? Eg: using the first approach (simple definition-SetDelayed-for the objective function) with 100 observations and 48 variables was a 6.5 min long evaluation (using AbsoluteTiming) on a run. I'd hope to achieve half the time for double the variables if at all possible, but I'd settle for any changes that point towards an improving performance.

---

In what follows, I initialize the problem instance and then present the three different versions of the objective function used, along with the results obtained.

This initializes the problem instance:

vs = 4;
nObs = 100;
xs = makeVars["x", vs];
ys = makeVars["y", vs];
args = Pattern[#, Blank[]] & /@ Join[xs, ys];
vars = Join[xs, ys];

{data, parameters} = makeData[nObs, vs];
fAux = Variance[fAuxObjective[xs, ys, #] & /@ data];

constraints = makeConstraints[xs, ys, parameters];
init = initVariables[xs, ys, ConstantArray[1/vs, vs], parameters];

As an example, I'm using two sets of 4 variables each ie 8 variables in total and 100 observations. The number of observations and total variables used affect the number of terms in the objective function.

The formula for the objective function is contained in fAux.

Below are presented the three different cases I have already tried:

Case 1:

(* case 1. objective function defined with SetDelayed *)

Clear[fObjective]
Evaluate[Apply[fObjective][args]] := Evaluate[fAux]

fcons = Join[{Apply[fObjective][vars]}, constraints];
NMinimize[fcons, init] 

{0.0989053, {x1 -> 0.244909, x2 -> 0.311147, x3 -> 0.260978, x4 -> 0.182966, 
y1 -> 0.248345, y2 -> 0.662958, y3 -> -0.485397, y4 -> 0.130779}}

Case 2:

(* case 2. objective function defined with Function and Set *)

Clear[fObjective]
fObjective = Function[Evaluate[vars], Evaluate[fAux]];

fcons = Join[{Apply[fObjective][vars]}, constraints];
NMinimize[fcons, init]

{0.0989053, {x1 -> 0.244909, x2 -> 0.311147, x3 -> 0.260978, x4 -> 0.182966, 
y1 -> 0.248345, y2 -> 0.662958, y3 -> -0.485397, y4 -> 0.130779}}

Case 3:

(* case 3. objective function defined with Compile and Set *)

Clear[fObjective]
copts = {CompilationTarget -> "C"};
gAux = Function[Evaluate[vars], Evaluate[fAux]];
fComp = Compile[{{args, _Real, 1}}, Apply[gAux][args], {{_, _Real}}, Evaluate[copts]];
fObjective[args_?(VectorQ[#, NumericQ] &)] := fComp[args];

fcons = Join[{fObjective[vars]}, constraints];
NMinimize[fcons, init]

{0.0989053, {x1 -> 0.244909, x2 -> 0.311147, x3 -> 0.260978, x4 -> 0.182966, 
 y1 -> -1., y2 -> 0.0161671, y3 -> -0.0143411, y4 -> 0.187533}}

In what follows are auxiliary definitions used in the generation of the objective function and the constraints.

(* create symbols to use as 'variables' *)
(* takes as input the string of the symbol to use and the number of symbols needed 
   and outputs appropriately numbered symbols *)
(* eg. evaluating makeVars["x", 3] outputs {x1, x2, x3} *)
makeVars[symbolString_, nVariables_] := Table[
  ToExpression[StringJoin[symbolString, ToString[i]]], {i, 1, nVariables}];


(* generate appropriate data points used in the objective function *)
(* takes as input the number of observation points (rows) to generate per variable, 
   the number of variables needed (columns) and a seed for random number generation *)
(* it outputs a list containing a matrix 'nObservations' x 'nVariables' and 
   a list 'nVariables' long *)
makeData[nObservations_, nVariables_, seed_: 2145879] := BlockRandom[
  {RandomReal[{-1, 1}, {nObservations, nVariables}], 
    RandomReal[{0, 1}, nVariables]}, RandomSeeding -> seed]


(* auxiliary function used in the construction of the objective function *)
(* takes as arguments three equally long lists; the former two are lists of 
   symbols (variables) and the later is a list of values (observations) *)
(* outputs a list (vector) of values )*)
(* eg fAuxObjective[{x1, x2}, {y1, y2}, {-0.1, 0.325}] 
   outputs -0.1 x1 - x1 y1 + 0.325 x2 - x2 y2 *)
fAuxObjective = Total[#1 #3 - #1 #2 ] &;


(* generate constraints *)
makeConstraints[variables1_, variables2_, 
  hiBounds2_, lowBounds2_: - 1., hiBounds1_: 1.] := Join[
    {Plus @@ variables1 == 1.},
    Thread[lowBounds2 <= variables2 <= hiBounds2],
    Thread[0. <= variables1 <= hiBounds1]
 ]


(* init variables *)
initVariables[variables1_, variables2_, initValues1_, initValues2_, 
  factor_: 0.01] := MapAt[
    Apply[Sequence][{# (1 - factor), # (1 + factor)}] &,
      Join[Transpose[{variables1, initValues1}], 
        Transpose[{variables2, initValues2}]], {All, 2}]

---

Original Q:

I'm trying to solve an optimization problem of modest size and I'm having trouble saving time by modifying the objective function.

The first attempt uses an objective function defined with a simple SetDelayed:

f = Total[vars[[1]] # - Times @@ vars] &;
BlockRandom[
  With[{data = RandomReal[{-1, 1}, {n, vs}], params = RandomReal[{0, 1}, vs]},
    With[{fval = Variance[f /@ data]},
      (* objective function *)
      objf[args__] := Evaluate[fval];
     ];
   ], RandomSeeding -> seed]

Please note that the objective function is dependent on a set of input data and other numeric parameters (I'm using random numbers to simulate them).

The second attempt follows more or less the same path but uses a Function instead:

BlockRandom[
  With[{data = RandomReal[{-1, 1}, {n, vs}], params = RandomReal[{0, 1}, vs]},
    With[{fval = Variance[f /@ data]},
      cf = Function[Evaluate[Flatten@vars], fval]
     ];
    objf[args__] := cf[args];
   ], RandomSeeding -> seed]

The third attempt uses Compile to define the objective function:

BlockRandom[
  With[{data = RandomReal[{-1, 1}, {n, vs}], params = RandomReal[{0, 1}, vs]},
    With[{fval = Variance[f /@ data]},
      With[{g = Function[Evaluate[Flatten@vars], fval]},
        cf = Compile[{{args, _Real, 1}}, Apply[g][args], {{_, _Real}}, Evaluate[copts]]
       ]
     ];
    objf[args_?(VectorQ[#, NumericQ] &)] := cf[args];
   ], RandomSeeding -> seed]

Please note the third argument in Compile without which Mathematica complains with CompiledFunction::cfte and CompiledFunction::cfexe.

After (limited) fooling around with the code provided in the end of this question, it seems that the first version is roughly equivalent to the second one (the second is marginally faster than the first) while the third is way slower than the first one (by a factor of 8x approx.).

And here is the problem: I was expecting a progressive improvement in evaluation speed between the three versions, which fails to obtain; especially so, after the unexpected slowdown from the version using Compile.

How can I improve on either of the three versions or how to adopt a whole different approach in order to obtain sensible timing for the optimization problem at hand (please see code blocks at the end of this question)?

Eg: using the first approach (simple definition-SetDelayed-for the objective function) with 100 observations and 48 variables was a 6.5 min long evaluation (using AbsoluteTiming) on a run. I'd hope to achieve half the time for double the variables if at all possible, but I'd settle for any changes that point towards an improving performance.

---

Block 1

Block[{makeVars, f, objf},

  (* no. of variables, symbols to use, observations (used in 'data' - see below) *)
  With[{vs = 4, syms = {"x", "y"}, n = 100, seed = 3254780},
    (* creates variables' symbols *)
    makeVars = Table[ToExpression[StringJoin[#1, ToString[i]]], {i, 1, vs}] &;

    With[{vars = makeVars /@ syms},

      (* 'f' depends on input data *)
      f = Total[vars[[1]] # - Times @@ vars] &;

      (* using random numbers in place of data and parameters *)
      BlockRandom[
        With[{data = RandomReal[{-1, 1}, {n, vs}], params = RandomReal[{0, 1}, vs]},

          (* 'fval' is the objective function *)
          With[{fval = Variance[f /@ data]},
            (* use a definition *)  
            objf[args__] := Evaluate[fval];
           ];

          (* add constraints and minimize *)
          With[{csts = {Total[vars[[1]]] == 1., Thread[0. <= vars[[1]] <= 1.], Thread[-1. <= vars[[-1]] <= params]}, d = 0.01, init = Transpose[{Flatten[vars], Join[ConstantArray[1/vs, vs], params]}]},
            NMinimize[Join[{Apply[objf][Flatten@vars]}, Flatten@csts], MapAt[Apply[Sequence][{# (1 - d), # (1 + d)}] &, init, {All, 2}]] // AbsoluteTiming
           ]
         ], RandomSeeding -> seed]
     ]
   ]
 ]

Block 2

Block[{makeVars, f, cf, objf},

  (* no. of variables, symbols to use, observations (used in 'data' - see below) *)
  With[{vs = 4, syms = {"x", "y"}, n = 100, seed = 3254780},
    (* creates variables' symbols *)
    makeVars = Table[ToExpression[StringJoin[#1, ToString[i]]], {i, 1, vs}] &;

    With[{vars = makeVars /@ syms},

      (* 'f' depends on input data *)
      f = Total[vars[[1]] # - Times @@ vars] &;

      (* using random numbers in place of data and parameters *)
      BlockRandom[
        With[{data = RandomReal[{-1, 1}, {n, vs}], params = RandomReal[{0, 1}, vs]},

          (* 'fval' is the objective function *)
          With[{fval = Variance[f /@ data]},
            (* enclose it in a function *)  
            cf = Function[Evaluate[Flatten@vars], fval];
           ];
          objf[args__] := cf[args];

          (* add constraints and minimize *)
          With[{csts = {Total[vars[[1]]] == 1., Thread[0. <= vars[[1]] <= 1.], Thread[-1. <= vars[[-1]] <= params]}, d = 0.01, init = Transpose[{Flatten[vars], Join[ConstantArray[1/vs, vs], params]}]},
            NMinimize[Join[{Apply[objf][Flatten@vars]}, Flatten@csts], MapAt[Apply[Sequence][{# (1 - d), # (1 + d)}] &, init, {All, 2}]] // AbsoluteTiming
           ]
         ], RandomSeeding -> seed]
     ]
   ]
 ]

Block 3

Block[{makeVars, f, cf, objf},

  (* no. of variables, symbols to use, observations (used in 'data' - see below) *)
  With[{vs = 4, syms = {"x", "y"}, n = 100, seed = 3254780},
    (* creates variables' symbols *)
    makeVars = Table[ToExpression[StringJoin[#1, ToString[i]]], {i, 1, vs}] &;

    With[{vars = makeVars /@ syms},

      (* 'f' depends on input data *)
      f = Total[vars[[1]] # - Times @@ vars] &;

      (* using random numbers in place of data and parameters *)
      BlockRandom[
        With[{data = RandomReal[{-1, 1}, {n, vs}], params = RandomReal[{0, 1}, vs]},

          (* 'fval' is the objective function *)
          With[{fval = Variance[f /@ data]},
            (* enclose it in a function *)  
            cf = Compile[{{args, _Real, 1}}, Apply[g][args], {{_,_Real}}, Evaluate[copts]];
           ];
          objf[args_?(VectorQ[#, NumericQ] &)] := cf[args];

          (* add constraints and minimize *)
          With[{csts = {Total[vars[[1]]] == 1., Thread[0. <= vars[[1]] <= 1.], Thread[-1. <= vars[[-1]] <= params]}, d = 0.01, init = Transpose[{Flatten[vars], Join[ConstantArray[1/vs, vs], params]}]},
            NMinimize[Join[{objf[Flatten@vars]}, Flatten@csts], MapAt[Apply[Sequence][{# (1 - d), # (1 + d)}] &, init, {All, 2}]] // AbsoluteTiming
           ]
         ], RandomSeeding -> seed]
     ]
   ]
 ]

Answer 1

This should give you about a factor of 5:

Clear[fObjective]
copts = {
 CompilationTarget -> "WVM", 
 CompilationOptions -> {"InlineExternalDefinitions" -> True}
};
gAux = Function[Evaluate[vars], Evaluate[fAux]];
cargs = {#, _Real} & /@ vars;
cobj = gAux @@ vars;
fComp = Compile[
  Evaluate@cargs,
  Evaluate@cobj,
  Evaluate[copts]
 ]
];
fObjective[args_?(VectorQ[#, NumericQ] &)] := fComp @@ args;

fcons = Join[{fObjective[vars]}, constraints];
NMinimize[fcons, init,Method -> "NelderMead"]

There were a couple of issues:

1. There were some syntax errors in your compile command. As a result, the compiled function was just an external call back to the math kernel. (Just slowing things down.)
2. Similar to the above, you need to tell Compile to inline the external definitions (otherwise it is just going to end up calling back to the Kernel).
3. Without specifying the Method to be use by NMinimize, Mathematica was wasting a lot of time using its default methods. (Which probably inlcude determing the gradient of the object function, which it would have to do numerically for the compiled case.)

As a final note: I used the "Wolfram Virtual Machine" as target for the compile for two reasons. 1) It allows for greater debug easy using CompiledFuntionTools 2) The machine I am testing on did not have a "C" compiler installed. You should probably be able to gain more speed by switching the compilation target back to "C".