# How to use a compiled function in numerical optimization

I'm trying to wrap my head around function Compile. So, this is a spin-off question from: a previous question on using compile[] correctly.

I have a function that has been expressed as a compiled function. All variables are known except $\beta_{g,n}$ which is an optimization variable and should remain in the compiled function. Here's what I've done to express it usingCompile.

compileCost =
Compile[{{P, _Real}, {Ns, _Integer}, {gh, _Real, 2},
{usrK, _Integer,1}, {G, _Integer}, {betha, _Integer,2}},
N[Sum[-Exp[usrK[[g]]/(P gh[[g, n]])] (usrK[[g]] betha[[g,n]])/
Log ExpIntegralEi[-(usrK[[g]]/(P gh[[g, n]]))]
, {g, 1, G}, {n, 1, Ns}]
]]


When evaluated with the input:

P = 5; Ns = 7; G = 3; usrK = RandomInteger[{2, 5}, G];
gh = RandomReal[{1, 3}, {G, Ns}];
betha = Table[1, {g, G}, {n, Ns}];

compileCost[P, Ns, gh, usrK, G, betha]

101.298


Then I tried using this an a numerical optimization problem (removing the constraint functions for simplicity).

With[{variable = compileCost},
fxnCompiledCost = Compile[{{P, _Real},{Ns, _Integer},{gh, _Real, 2},{usrK, _Integer, 1},
{G,_Integer},{betha,_Integer,2}}, variable[P,Ns,gh,usrK,G,betha],
"RuntimeOptions" -> {"EvaluateSymbolically" -> False},
CompilationOptions -> {"InlineCompiledFunctions" -> False},
Parallelization -> Automatic]
];

subcindexSet = Flatten[Table[\beta_{g, n}, {g, G}, {n, Ns}]];
NMaximize[fxnCompiledCost[P, Ns, gh, usrK, G, \beta_{g, n}], subcindexSet]


When I try executing this snippet, I get a not-a-number message.

NMaximize::nnum: The function value -CompiledFunction ... is not a number at {...}

How do I move forward from here please?

I've followed examples from here: Using compiled function inside NMinimize

===============Edit =========

Actually there are constraint but I thought I could remove them to simplify my question:

 intConstRelax = Thread[0 <= subcindexSet <= 1];
binaryConst = List[subcindexSet \[Element] Integers];
nonShareConstr = Table[Sum[\beta_{g,n}, {g, 1, G}] == 1, {n, Ns}];
constrSet = Join[nonShareConstr, intConstRelax, binaryConst];


===============Edit =========

Computation slows down considerably when I add one more constraint equation to the above:

 qosConstr = Table[Sum[\beta_{g,n} Log2[1 + (P/Ns) gh[[g, n]] ], {n,Ns}] >= constQ, {g, G}];


So constraint list now becomes:

 constrSet = Join[nonShareConstr, intConstRelax, binaryConst, qosConstr];


===============Edit =========

Thank you everyone for your help. Due to time difference I was not able to comment instantly on your questions. Some of the procedure are also a little bit over my head so I'm checking up on them. Thanks so much.

• Do your constraints mean that all $\beta_{g,n}$ are either 0 or 1, and nothing else? If so, I've updated my answer a bit. Sep 30, 2014 at 17:29
• Yes. the constraints means 0 or 1. Just binary integer. Oct 1, 2014 at 8:46
• Ok, well look through the answers below! Oct 1, 2014 at 8:48

## Introduction

The update of the question invalidates my original approach, which depended on the problem being a simple special case. The current problem consists of a linear objective function with linear constraints, to be maximized over the integers. Theoretically, one can solve the problem with LinearProgramming. In the OP's setup, the variables betha form a (rank-2) matrix instead of a rank-1 vector needed to use LinearProgramming. One can flatten the problem out and use integer linear programming. Alternatively, one can do a brute-force search, which turns out to be faster in this case.

## Brute force is a little faster

Define the cost in terms of a memoized coefficient function costCoeff, and it turns out an exhaustive search of the possible inputs is at least twice as fast as LinearProgramming. See further down for some discussion of the use of listability in the code for costCoeff etc. For instance, Times automatically multiplies corresponding parts of lists and matrices. Such properties allows us to rewrite the OP's code more efficiently in terms of the whole matrices and vectors instead of their parts.

constQ = 2;  (* not given -- made up *)
P = 5; Ns = 7; G = 3;
SeedRandom;
usrK = RandomInteger[{2, 5}, G];
gh = RandomReal[{1, 3}, {G, Ns}];

costCoeff[P_?NumericQ, Ns_?NumericQ, gh_?(MatrixQ[#, NumericQ] &),
usrK_?(VectorQ[#, IntegerQ] &), G_?NumericQ] :=
costCoeff[P, Ns, gh, usrK, G] = -Exp[usrK/(P gh)] (usrK) / Log ExpIntegralEi[-(usrK/(P gh))]

cost[P_?NumericQ, Ns_?NumericQ, gh_?(MatrixQ[#, NumericQ] &),
usrK_?(VectorQ[#, IntegerQ] &), G_?NumericQ, betha_] :=
Total[costCoeff[P, Ns, gh, usrK, G] betha, 2];


We use UnitStep[x - constQ] to implement the OP's qosConstr constraints. The values of UnitStep have be 1 on all rows for the constraints to be met; if one is not met, the minimum of the values will be 0. This criterion is used to Pick out of all the 0-1 matrices with exactly one 1 in each column, those that satisfy these constraints.

domain = Transpose[Tuples[Permutations[Append[ConstantArray[0, G - 1], 1]], 7], {1, 3, 2}];
domain2 =
Pick[
domain,
Min /@ Transpose @
UnitStep[Total[Transpose[domain, {3, 1, 2}] Log2[1 + (P/Ns) gh], {2}] - constQ],
1];
betamax = domain2[[
Last @ Ordering[cost[P, Ns, gh, usrK, G, Transpose[domain2, {3, 1, 2}]]]
]]
cost[P, Ns, gh, usrK, G, betamax]
(*
{{0, 0, 1, 0, 0, 1, 1}, {1, 0, 0, 0, 1, 0, 0}, {0, 1, 0, 1, 0, 0, 0}}
41.9077
*)


## Linear programming is conceptually more natural

LinearProgramming minimizes its objective function. To get a maximum we can take the negative of the OP's cost function. Below, costCoeff represents the matrix of coefficients. To maximize the cost function, I used -Flatten @ costCoeff[P, Ns, gh, usrK, G] for the objective function. Each linear constraints is specified in two parts, a list of coefficients of the variables and a pair {b, s} representing a boundary value b and a code s -- -1, 0, or 1 -- for the inequality as described in the documentation for LinearProgramming. The second-to-last argument specifies the bounds on each variable to be 0 and 1.

solFlat = LinearProgramming[
-Flatten @ costCoeff[P, Ns, gh, usrK, G],           (* objective function *)
Table[Flatten @ Transpose @                         (* one 1 in each column *)
ReplacePart[ConstantArray[0., {Ns, G}], i -> ConstantArray[1, G]], {i, Ns}] ~Join~
MapIndexed[                                        (* fourth constraint coeffs *)
Flatten @ ReplacePart[ ConstantArray[0., {G, Ns}], #2 -> #1] &,
Log2[1 + (P/Ns) gh]],
ConstantArray[{1, -1}, Ns] ~Join~                   (* constraint inequalities *)
ConstantArray[{constQ, 1}, G],
ConstantArray[{0, 1}, G*Ns],                        (* 0-1 constraint *)
Integers
];
sol = Partition[solFlat, Ns]                           (* reshape into G x Ns matrix *)
cost[P, Ns, gh, usrK, G, sol]                          (* maximum cost *)
(*
{{0, 0, 1, 0, 0, 1, 1}, {1, 0, 0, 0, 1, 0, 0}, {0, 1, 0, 1, 0, 0, 0}}
41.9077
*)


## Remarks and explanation

### Objective function - listability and efficiency

To use Mathematica efficiently with matrices, vectors, and other arrays, it pays to learn how theListable attribute (or more precisely, internal vectorization) and threading work with the numeric functions, including Plus and Times, as well as the importance of packed arrays (see What is a Mathematica packed array?). Here are a few ways to compute the OP's cost function, the OP's way costOP, the way from the original answer costME2, and the way above costFlat. They are equivalent.

costOP = Sum[-Exp[usrK[[g]]/(P gh[[g, n]])] (usrK[[g]] betha[g, n]) /
Log ExpIntegralEi[-(usrK[[g]]/(P gh[[g, n]]))], {g, 1, G}, {n, 1, Ns}];
costME2 = Total[costCoeff[P, Ns, gh, usrK, G] Array[betha, {G, Ns}], 2];
costFlat = Flatten@costCoeff[P, Ns, gh, usrK, G] . Flatten@Array[betha, {G, Ns}];
costOP == costME2 == costFlat
(* True *)


Remark: In computing the objective function, instead of using Sum and Part, one should use the listability of the built-in numeric functions, which Mathematica handles much more efficiently. This would be an important point if the objective function were to be evaluated many times, as it would if NMinimize were used. In the above costCoeff is evaluated only once. If we turn them into functions, we can compare the timings:

ClearAll[costOP, costME2, costFlat];
costOP[betha_] := Sum[-Exp[usrK[[g]]/(P gh[[g, n]])] (usrK[[g]] betha[[g, n]])/
Log ExpIntegralEi[-(usrK[[g]]/(P gh[[g, n]]))], {g, 1, G}, {n, 1, Ns}];
costFlat[betha_] := Flatten@costCoeff[P, Ns, gh, usrK, G] . Flatten@betha;
costME2[betha_] := Total[costCoeff[P, Ns, gh, usrK, G] betha, 2];

SeedRandom;
testOP = Table[costOP[RandomInteger[1, {G, Ns}]], {1000}]; // AbsoluteTiming
SeedRandom;
testFlat = Table[costFlat[RandomInteger[1, {G, Ns}]], {1000}]; // AbsoluteTiming
SeedRandom;
testME2 = costME2[Transpose[RandomInteger[1, {1000, G, Ns}], {3, 1, 2}]]; // AbsoluteTiming
(*
{0.372749, Null}
{0.130880, Null}
{0.000734, Null}  <-- Wow.  Compare with the brute-force method above.
*)

testOP == testFlat == testME2
(* True *)


### Constraints

The constraints are implemented in three pieces.

First, the lists

Table[Flatten @ Transpose @ ReplacePart[ConstantArray[0., {Ns, G}], i -> {1, 1, 1}], {i, Ns}]


together with

ConstantArray[{1, -1}, Ns]


specify the column sums should be at most 1.

Second, the lists

MapIndexed[
Flatten @ ReplacePart[ConstantArray[0., {G, Ns}], #2 -> #1] &,
Log2[1 + (P/Ns) gh]]


together with

ConstantArray[{constQ, 1}, G]


specify the row sums with the coefficients from the corresponding row of the matrix Log2[1 + (P/Ns) gh] should be at most constQ.

Third, the argument

ConstantArray[{0, 1}, G*Ns],                        (* 0-1 constraint *)


specifies that each variable should be between 0 and 1 (inclusive), and the final argument Integers specifies that the variables are to be integers.

### Timing

The AbsoluteTiming of the brute-force method is 0.001875 sec. The AbsoluteTiming of the LinearProgramming code is 0.004235, which is more than twice as long as the brute-force method.

First, your objective function is linear in your variables betha, which form a matrix in your setup.

Second, your constraints imply that an input is a 0-1 matrix with exactly one entry being 1 in each column. Edit note: the OP later added a constraint that makes the following not apply to the current question.

The fastest way to maximize the function is to pick among each column of coefficients, the position of the maximum coefficient and set the corresponding variable to 1 (and set the others in the column to zero, of course).

Transpose @
Map[ReplacePart[ConstantArray[0, G], # -> 1] &,
Last /@ Ordering /@ Transpose @ costCoeff[P, Ns, gh, usrK, G]
]
cost[P, Ns, gh, usrK, G, %]
(*
{{0, 1, 1, 0, 1, 1, 1}, {1, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}}
42.5
*)


A check

A simple check is to evaluate the function on all possible inputs and take the max.

domain = Transpose[
Tuples[Permutations[Append[ConstantArray[0, G - 1], 1]], 7], {1, 3, 2}];
domain[[
Last @ Ordering[cost[P, Ns, gh, usrK, G, Transpose[domain, {3, 1, 2}]]]
]]

(*
{{0, 1, 1, 0, 1, 1, 1}, {1, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}}
42.5
*)

• hadn't thought of this, nice! Sep 30, 2014 at 19:01
• Your heuristic approach is quite interesting and Yes, your understanding of the problem is spot-on. My current (extended modified version of the problem) is super slow and I'm guessing this approach should be faster. I notice your current suggestion does not include the constraints. Would you consider updating with the constraint please? I've added four constraints in my edited question. Oct 1, 2014 at 9:40
• wooow. Lots of heavy stuffs going on here. I'm taking 'em a bit at a time. Thanks for helping out. Oct 2, 2014 at 11:51
• Amazing stuffs. I don't know who you are but I'll buy someone a cup of coffee on your behalf! Thank you for putting so much work and explanation into this. Oct 2, 2014 at 12:58
• @Afloz, Thank you very much. Mathematica much more efficiently deals with numeric arrays as arrays. Learning to think in those terms takes time; at least it did for me. Here are some simple examples to play with. For threading, compare how Times works in {2, 3, 5} {{a, b}, {c, d}, {e, f}}, {2, 3} {{a, b, c}, {d, e, f}} and {{1, 2}, {3, 4}, {5, 6}} {{a, b}, {c, d}, {e, f}}. For transposing multiple levels, let mat = Table[a[i] b[j] c[k], {i, 2}, {j, 3}, {k, 4}]; matT = Transpose[mat, {3, 1, 2}], and figure out why these are the same: mat[[1, 3, 4]] and matT[[3, 4, 1]]. Oct 2, 2014 at 13:51

So we have a couple things going on here.

First, as mentioned in an answer to your previous question, ExpIntegralEi is not compilable (MainEvaluate called). So speedup from compiling this function will be negligible and it may be easier to not even compile in this case.

Second, to properly use NMaximize we need the function to accept numeric arguments only. Since some arguments are lists, we have to get a little creative on how we specify argument patterns. We can define the function as such:

compileCost[P_?NumericQ, Ns_?NumericQ, gh : {{_?NumericQ ..} ..}, usrK : {_?NumericQ..},
G_?NumericQ, betha : {{_?NumericQ ..} ..}] :=
N[Sum[-Exp[usrK[[g]]/(P gh[[g, n]])] (usrK[[g]] betha[[g, n]])/Log
ExpIntegralEi[-(usrK[[g]]/(P gh[[g, n]]))], {g, 1, G}, {n, 1, Ns}]]


This ensures that the function will only evaluate for lists of certain dimensions that contain only numeric values.

Third, NMaximize is not being set up properly. We can adjust like so:

subcindexSet = Table[Symbol["betha" <> ToString@g <> ToString@n], {g, G}, {n, Ns}];
NMaximize[compileCost[P, Ns, gh, usrK, G, subcindexSet], Flatten@subcindexSet]


Without constraints, the function seems to be unbounded. Also, a word of caution that global minimization/maximization becomes very unwieldy, especially when guessing so many variables.

If I reinterpret your constraints so they make sense,

constraints = Apply[And,
Flatten@{
Apply[LessEqual, {0, #, 1}] & /@ Flatten@subcindexSet,
Map[Total@# == 1 &, Transpose@subcindexSet, {1}]
}]

(*
0 <= betha11 <= 1 && 0 <= betha12 <= 1 && 0 <= betha13 <= 1 &&
0 <= betha14 <= 1 && 0 <= betha15 <= 1 && 0 <= betha16 <= 1 &&
0 <= betha17 <= 1 && 0 <= betha21 <= 1 && 0 <= betha22 <= 1 &&
0 <= betha23 <= 1 && 0 <= betha24 <= 1 && 0 <= betha25 <= 1 &&
0 <= betha26 <= 1 && 0 <= betha27 <= 1 && 0 <= betha31 <= 1 &&
0 <= betha32 <= 1 && 0 <= betha33 <= 1 && 0 <= betha34 <= 1 &&
0 <= betha35 <= 1 && 0 <= betha36 <= 1 && 0 <= betha37 <= 1 &&
betha11 + betha21 + betha31 == 1 &&
betha12 + betha22 + betha32 == 1 &&
betha13 + betha23 + betha33 == 1 &&
betha14 + betha24 + betha34 == 1 &&
betha15 + betha25 + betha35 == 1 &&
betha16 + betha26 + betha36 == 1 &&
betha17 + betha27 + betha37 == 1
*)


However, I'm struggling to understand what you mean by restricting subcindexSet to be between 0 and 1, and then setting them all to be integers. Do you mean they are all equal to 0 or 1, and can therefore take no other values?

I guess you could redefine the function with an _?IntegerQ for betha,

compileCost[P_?NumericQ, Ns_?NumericQ, gh : {{_?NumericQ ..} ..},
usrK : {_?NumericQ ..}, G_?NumericQ, betha : {{_?IntegerQ ..} ..}] :=
N[Sum[-Exp[usrK[[g]]/(P gh[[g, n]])] (usrK[[g]] betha[[g, n]])/
Log ExpIntegralEi[-(usrK[[g]]/(P gh[[g, n]]))],
{g, 1, G}, {n, 1, Ns}]]


Then run it with the same constraints as given above, and do a little bit of magic-trickery, because I can't seem to get it to work with $\in$ Integers properly.

Similarly, restricting to (betha11 == 0 || betha11 == 1) && ... in the constraints seems to take forever, compared to this fudge(?).

resultrules = FindMaximum[
{compileCost[P, Ns, gh, usrK, G, subcindexSet], constraints},
Flatten@subcindexSet, WorkingPrecision -> 18];

resultunrounded = subcindexSet /. Last@resultrules
(*
{{1.58489319246111726*10^-14, 1.58489319246111726*10^-14,
0.999999999999984151, 1.58489319246111726*10^-14,
0.999999999999984151, 0.999999999999984151,
0.999999999999984151}, {0.999999999999984151, 0.999999999999984151,
1.58489319246111726*10^-14, 0.999999999999984151,
1.58489319246111726*10^-14, 1.58489319246111726*10^-14,
1.58489319246111726*10^-14}, {1.58489319246111726*10^-14,
1.58489319246111726*10^-14, 1.58489319246111726*10^-14,
1.58489319246111726*10^-14, 1.58489319246111726*10^-14,
1.58489319246111726*10^-14, 1.58489319246111726*10^-14}}
*)

result = Round@subcindexSet /. Last@resultrules
(* {{0, 0, 1, 0, 1, 1, 1}, {1, 1, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}} *)


Does this seem like a correct solution?

As @kale mentions, you need constraints. Otherwise, look what happens when I use FindMaximum (rather than NMaximize, for this example).

P = 5; Ns = 7; G = 3; usrK = RandomInteger[{2, 5}, G];
gh = RandomReal[{1, 3}, {G, Ns}];
betha = Table[1, {g, G}, {n, Ns}];

compileCost[P_?NumericQ, Ns_?NumericQ, gh : {{_?NumericQ ..} ..},
usrK : {_?NumericQ ..}, G_?NumericQ, betha : {{_?NumericQ ..} ..}] :=
N[Sum[-Exp[usrK[[g]]/(P gh[[g, n]])] (usrK[[g]] betha[[g, n]])/
Log ExpIntegralEi[-(usrK[[g]]/(P gh[[g, n]]))], {g, 1, G}, {n, 1, Ns}]]

subcindexSet = Table[Symbol["betha" <> ToString@g <> ToString@n], {g, G}, {n, Ns}];

FindMaximum[compileCost[P, Ns, gh, usrK, G, subcindexSet], Flatten@subcindexSet]
(* {1.97295*10^59,
{betha11 -> 3.44624*10^57, betha12 -> 2.26929*10^57,
betha13 -> 3.28037*10^57, betha14 -> 2.17458*10^57,
betha15 -> -5.88154*10^56, betha16 -> 4.10653*10^57,
betha17 -> 2.87052*10^57, betha21 -> 1.21735*10^57,
betha22 -> 1.00252*10^57, betha23 -> 7.45909*10^56,
betha24 -> 2.4225*10^57, betha25 -> 2.67294*10^57,
betha26 -> 3.08175*10^57, betha27 -> 5.92319*10^56,
betha31 -> 3.30424*10^57, betha32 -> 2.08312*10^57,
betha33 -> 1.4828*10^57, betha34 -> 1.68492*10^57,
betha35 -> 1.46808*10^57, betha36 -> 1.67782*10^57,
betha37 -> 3.11295*10^56}} *)


And as a further example,

constraints = Apply[And, Apply[LessEqual, {#, 100.}] & /@ Flatten@subcindexSet]
(*
betha11 <= 100. && betha12 <= 100. && betha13 <= 100. &&
betha14 <= 100. && betha15 <= 100. && betha16 <= 100. &&
betha17 <= 100. && betha21 <= 100. && betha22 <= 100. &&
betha23 <= 100. && betha24 <= 100. && betha25 <= 100. &&
betha26 <= 100. && betha27 <= 100. && betha31 <= 100. &&
betha32 <= 100. && betha33 <= 100. && betha34 <= 100. &&
betha35 <= 100. && betha36 <= 100. && betha37 <= 100.
*)

FindMaximum[{compileCost[P, Ns, gh, usrK, G, subcindexSet], constraints},
Flatten@subcindexSet]
(* {9372.77,
{betha11 -> 100., betha12 -> 100., betha13 -> 100.,
betha14 -> 100., betha15 -> 100., betha16 -> 100., betha17 -> 100.,
betha21 -> 100., betha22 -> 100., betha23 -> 100., betha24 -> 100.,
betha25 -> 100., betha26 -> 100., betha27 -> 100., betha31 -> 100.,
betha32 -> 100., betha33 -> 100., betha34 -> 100., betha35 -> 100.,
betha36 -> 100., betha37 -> 100.}} *)


And as a comparison of timing,

FindMaximum[{compileCost[P, Ns, gh, usrK, G, subcindexSet],
constraints}, Flatten@subcindexSet]; // AbsoluteTiming
NMaximize[{compileCost[P, Ns, gh, usrK, G, subcindexSet],
constraints}, Flatten@subcindexSet]; // AbsoluteTiming

(* FindMaximum: 0.84 seconds *)
(* NMaximize: 2.86 seconds *)

• I guess the $\in$ operator, binary restrictions on $\beta_{g,n}$ and formatting of my constraints seems to be what makes my code run for forever. Appears FindMaximum is faster because is Local right? I'm wondering if you did a comparison on optimal value in addition to the time comparison in the last step. Maybe I can consider results of FindMaximum an approximation if they are close enough. Oct 1, 2014 at 9:07
• The optimal results came out the same with FindMaximum and NMaximize when I ran it, but of course this is not always the case because the former is a local optimizer. Oct 1, 2014 at 9:13