Improve performance of Linear Optimization

Question

EDIT

I am reframing the question here by pulling out the relevant part about linear optimization, you can find the original question and the entire code down below.

I am trying to increase the performance of the following linear optimization problem.

(* Defining functions and paramtere values*)
precision = 20;
chernofFunc[Num_, p_, \[Epsilon]_] := 
  chernofFunc[Num, p, \[Epsilon]] = 
   N[-Log[\[Epsilon]] (1. + Sqrt[1. - (2. p Num/Log[\[Epsilon]])]), 
    precision];
plBFunc[l_, pmu_, mu_] := 
  plBFunc[l, pmu, mu] = 
   N[Sum[pmu[[j]] Exp[-mu[[j]]] mu[[j]]^l/l!, {j, 3}], precision];
nzvec = {371307946, 593752, 99756};
nz = Total[nzvec];
pmu = {1/2, 1/4, 1/4};
mu = {2/5, 1.2 10^-3, 2 10^-4};
qz = 1/4;
nb = 10^11;
m = 20;
epshoeff = 10^-15;
epschern = 10^-15;
(constraints)
zsingconstr = N[Join[
Table[
 sumsub = 
  Exp[-mu[[j]]] pmu[[j]] Sum[
    N[Subscript[nzvar, l]] mu[[j]]^l/(l! plBFunc[l, pmu, mu]), {l,
      0, m}]; 
 69.0776 + sumsub >= nzvec[[j]] + Subscript[\[Delta]\[Mu]nz, j] >= 
  sumsub, {j, 3}],

{Sum[N[Subscript[\[Delta]\[Mu]nz, j]], {j, 3}] == 0.},

Table[-Sqrt[-Log[epshoeff/2.] nz/2.] <= 
  N[Subscript[\[Delta]\[Mu]nz, j]] <= 
  Sqrt[-Log[epshoeff/2.] nz/2.], {j, 3}],

Table[
 0. <= N[Subscript[nzvar, l]] <= 
  Min[plBFunc[l, pmu, mu] qz nb + 
    chernofFunc[qz nb, plBFunc[l, pmu, mu], epschern], nz], {l, 0,
   m}]
]];




(variables)
zsingvars = 
  N[Join[Table[Subscript[nzvar, l], {l, 0, m}], 
    Table[Subscript[[Delta][Mu]nz, j], {j, 3}]]];
zsingres = 
 LinearOptimization[N[Subscript[nzvar, 1]], zsingconstr, zsingvars(,
   Tolerance->10^(-10)),(WorkingPrecision->MachinePrecision,)
   Method -> "MOSEK", PerformanceGoal -> "Speed"] // AbsoluteTiming

The variables are of the form Subscript[nzvar, l] for l from 0 to m=20 and Subscript[\[Delta]\[Mu]nz, j] for j from 1 to 3. The constraints zsingconstr is constructed by joining four separate lists of bounds using some linear combination of the variables.

Heres everything that I've tried so far

• Represent every number as a floating point
• Enclosed every variable and function output in N[...]
• Tweaked the options for LinearOptimization (In the code below I'm using Method->"MOSEK" but if MOSEK isnt available Method->"InteriorPoint" will also work)

I am not really sure what else could be done. The objective function of the linear optimization is one of the variable itself and hence I don't expect Compile to be of much use there.

I think that maybe an improvement could be made in the way I'm constructing the constraints; doing it in a way thats faster and more 'easier' for LinearOptimization?

Original Question

(*------------------1----------------------*)
SetOptions[$FrontEndSession, NotebookAutoSave -> True];
NotebookSave[];
$PreRead = (# /. 
     s_String /; 
       StringMatchQ[s, NumberString] && 
        Precision@ToExpression@s == MachinePrecision :> s <> "`20." &);
(Defining functions that will called in the main function)
precision = 20;
chernofFunc[Num_, p_, [Epsilon]] := 
  chernofFunc[Num, p, [Epsilon]] = 
   N[-Log[[Epsilon]] (1. + Sqrt[1. - (2. p Num/Log[[Epsilon]])]), 
    precision];
plBFunc[l, pmu_, mu_] := 
  plBFunc[l, pmu, mu] = 
   N[Sum[pmu[[j]] Exp[-mu[[j]]] mu[[j]]^l/l!, {j, 3}], precision];
(------------------2----------------------)
(The main function)
keyrateFuncpaper[dB_?NumericQ, qx_?NumericQ, pmu1_?NumericQ, 
   pmu2_?NumericQ, mu1_?NumericQ, mu2_?NumericQ] := 
  Block[{$MaxExtraPrecision = 50(,$MachinePrecision=
    60,$MinPrecision=$MachinePrecision,$MaxPrecision=60)}, 
   Module[{mu3 = 2. 10.^-4., Y0 = 6007379741251949. 10.^-22., 
     eps = 10^-15., ed = 10.^-2., fEC = 116. 10.^-2., nb = 10.^11., 
     step = 10.^-2., index = 0., m = 20.,
     epstrunc = 10.^-15.,
     epshoeff = 10.^-15.,
     epschern = 10.^-15., mu, pmu , qz, pmu3, etasys, nxmu1, nxmu2, 
     nxmu3, Emu1, Emu2, Emu3, mxmu1, mxmu2, mxmu3, nzmu1, nzmu2, 
     nzmu3, mzmu1, mzmu2, mzmu3, nx, mx, nz, mz, eobs, sumsubx, 
     sumsubmz, nzvec, mzvec, tao0, nxmu3m, nxmu2p, auxSx0, Sx0, 
     Sx0true, tao1, nxmu2m, nxmu3p, nxmu1p, auxSx1, Sx1, Sx1true, 
     nzmu3m, nzmu2p, auxSz0, Sz0, Sz0true, nzmu2m, nzmu3p, nzmu1p, 
     auxSz1, Sz1, Sz1true, mzmu2p, mzmu3m, auxVz1, Vz1, 
     nzmu1singleerror, nzmu2singleerror, nzmu3singleerror, 
     nzsingleerror, VzbySz, aux1, aux2, nxvec, auxphix, phixfin, 
     phixcorr, nx0, nx1, nz1, mz1, phix, nxcoeff, lambdaEC, keylength,
      keyrate, zsingconstr, zsingres, zsingvars, zerrconstr, zerrres, 
     zerrvars, xconstr, xres, xvars},
qz = 1. - qx;
mu = {mu1, mu2, mu3};

pmu3 = 1. - pmu1 - pmu2;

pmu = {pmu1, pmu2, pmu3};

etasys = 10.^(-dB/10.);

nxmu1 = Floor[nb*qx*qx*pmu1*(1. - (1. - Y0)*Exp[-etasys*mu1])]; 
nxmu2 = Floor[nb*qx*qx*pmu2*(1. - (1. - Y0)*Exp[-etasys*mu2])]; 
nxmu3 = Floor[nb*qx*qx*pmu3*(1. - (1. - Y0)*Exp[-etasys*mu3])];

Emu1 = ((Y0/2. - ed)*Exp[-etasys*mu1] + 
    ed)/(1. - (1. - Y0)*Exp[-etasys*mu1]); 
Emu2 = ((Y0/2. - ed)*Exp[-etasys*mu2] + 
    ed)/(1. - (1. - Y0)*Exp[-etasys*mu2]); 
Emu3 = ((Y0/2. - ed)*Exp[-etasys*mu3] + 
    ed)/(1. - (1. - Y0)*Exp[-etasys*mu3]); 
mxmu1 = Ceiling[nxmu1*Emu1]; mxmu2 = Ceiling[nxmu2*Emu2]; 
mxmu3 = Ceiling[nxmu3*Emu3];

nzmu1 = 
 Floor[nb*(1. - qx)*(1. - qx)*
   pmu1*(1. - (1. - Y0)*Exp[-etasys*mu1])]; 
nzmu2 = Floor[
  nb*(1. - qx)*(1. - qx)*pmu2*(1. - (1. - Y0)*Exp[-etasys*mu2])]; 
nzmu3 = Floor[
  nb*(1. - qx)*(1. - qx)*pmu3*(1. - (1. - Y0)*Exp[-etasys*mu3])];

Emu1 = ((Y0/2. - ed)*Exp[-etasys*mu1] + 
    ed)/(1. - (1. - Y0)*Exp[-etasys*mu1]);
Emu2 = ((Y0/2. - ed)*Exp[-etasys*mu2] + 
    ed)/(1. - (1. - Y0)*Exp[-etasys*mu2]);
Emu3 = ((Y0/2. - ed)*Exp[-etasys*mu3] + 
    ed)/(1. - (1. - Y0)*Exp[-etasys*mu3]);

mzmu1 = Ceiling[nzmu1*Emu1];
mzmu2 = Ceiling[nzmu2*Emu2];
mzmu3 = Ceiling[nzmu3*Emu3];

nx = nxmu1 + nxmu2 + nxmu3;
mx = mxmu1 + mxmu2 + mxmu3;
nz = nzmu1 + nzmu2 + nzmu3;
mz = mzmu1 + mzmu2 + mzmu3;

eobs = If[nx == 0., 0., mx/nx];




(*First Linear Optimization*)

nzvec = {nzmu1, nzmu2, nzmu3};
zsingconstr = 
 N[Join[Table[
    sumsub = 
     Exp[-mu[[j]]] pmu[[j]] Sum[
       N[Subscript[nzvar, 
         l]] mu[[j]]^l/(l! plBFunc[l, pmu, mu]), {l, 0, m}]; 
    69.0776 + sumsub &gt;= 
     nzvec[[j]] + Subscript[\[Delta]\[Mu]nz, j] &gt;= sumsub, {j, 3}],
   {Sum[N[Subscript[\[Delta]\[Mu]nz, j]], {j, 3}] == 0.},
   Table[-Sqrt[-Log[epshoeff/2.] nz/2.] &lt;= 
     N[Subscript[\[Delta]\[Mu]nz, j]] &lt;= 
     Sqrt[-Log[epshoeff/2.] nz/2.], {j, 3}],
   Table[
    0. &lt;= N[Subscript[nzvar, l]] &lt;= 
     Min[plBFunc[l, pmu, mu] qz nb + 
       chernofFunc[qz nb, plBFunc[l, pmu, mu], epschern], nz], {l,
      0, m}]]];
zsingvars = 
 N[Join[Table[Subscript[nzvar, l], {l, 0, m}], 
   Table[Subscript[\[Delta]\[Mu]nz, j], {j, 3}]]];
zsingres = 
 LinearOptimization[N[Subscript[nzvar, 1]], zsingconstr, 
  zsingvars(*,Tolerance-&gt;10^(-10)*),(*WorkingPrecision-&gt;
  MachinePrecision,*)Method -&gt; &quot;MOSEK&quot;, 
  PerformanceGoal -&gt; &quot;Speed&quot;];


(*Second Linear Optimization*)
mzvec = {mzmu1, mzmu2, mzmu3};
zerrconstr = 
 N[Join[Table[
    sumsubmz = 
     Exp[-mu[[j]]] pmu[[j]] Sum[
       N[Subscript[mzvar, 
         l]] mu[[j]]^l/(l! plBFunc[l, pmu, mu]), {l, 0, m}]; 
    69.0776 + sumsubmz &gt;= 
     mzvec[[j]] + N[Subscript[\[Delta]\[Mu]mz, j]] &gt;= 
     sumsubmz, {j, 3}],
   {Sum[N[Subscript[\[Delta]\[Mu]mz, j]], {j, 3}] == 0.},
   Table[
    -Sqrt[-Log[epshoeff/2.] mz/2.] &lt;= 
     N[Subscript[\[Delta]\[Mu]mz, j]] &lt;= 
     Sqrt[-Log[epshoeff/2.] mz/2.], {j, 3}],
   Table[
    0. &lt;= N[Subscript[mzvar, l]] &lt;= 
     Min[plBFunc[l, pmu, mu] qz nb + 
       chernofFunc[qz nb, plBFunc[l, pmu, mu], epschern], mz], {l,
      0, m}]]];
zerrvars = 
 N[Join[Table[Subscript[mzvar, l], {l, 0, m}], 
   Table[Subscript[\[Delta]\[Mu]mz, j], {j, 3}]]];
zerrres = 
 LinearOptimization[N[-Subscript[mzvar, 1]], zerrconstr, zerrvars(*,
  Tolerance-&gt;10^(-10)*),(*WorkingPrecision-&gt;MachinePrecision(*,
  Tolerance-&gt;10^(-10)*),*)Method -&gt; &quot;MOSEK&quot;, 
  PerformanceGoal -&gt; &quot;Speed&quot;];




mz1 = (Subscript[mzvar, 1] /. zerrres);
nz1 = (Subscript[nzvar, 1] /. zsingres);

phix = 
 If[nz1 != 0., mz1/nz1, 0.] + 
  If[Or[nx == 0., nz == 0.], 0., 
   Sqrt[(nx + nz) (nx + 1.) Log[1./eps]/(2. nx^2. nz)]];

nxcoeff =(*SetPrecision[*)(1. - If[0. &lt; phix &lt; 1.,
    -phix*Log[2., phix] - (1. - phix) Log[2., 1. - phix], 0.
    ])(*,60]*);
lambdaEC = fEC If[0. &lt; eobs &lt; 1.,
   -eobs*Log[2., eobs] - (1. - eobs) Log[2., 1. - eobs], 0.
   ];





(*Third Linear Optimization*)

nxvec = {nxmu1, nxmu2, nxmu3};
xconstr = 
 N[Join[Table[
    sumsubx = 
     Exp[-mu[[j]]] pmu[[j]] Sum[
       N[Subscript[nxvar, 
         l]] mu[[j]]^l/(l! plBFunc[l, pmu, mu]), {l, 0, m}]; 
    69.0776 + sumsubx &gt;= 
     nxvec[[j]] + N[Subscript[\[Delta]\[Mu]nx, j]] &gt;= sumsubx, {j,
      3}],
   {Sum[N[Subscript[\[Delta]\[Mu]nx, j]], {j, 3}] == 0.},
   Table[-Sqrt[-Log[epshoeff/2.] nx/2.] &lt;= 
     N[Subscript[\[Delta]\[Mu]nx, j]] &lt;= 
     Sqrt[-Log[epshoeff/2.] nx/2.], {j, 3}],
   Table[
    0 &lt;= N[Subscript[nxvar, l]] &lt;= 
     Min[plBFunc[l, pmu, mu] (qx) nb + 
       chernofFunc[(qx) nb, plBFunc[l, pmu, mu], epschern], 
      nx], {l, 0, m}]]];
xvars = 
 N[Join[Table[Subscript[nxvar, l], {l, 0, m}], 
   Table[Subscript[\[Delta]\[Mu]nx, j], {j, 3}]]];
xres = 
 LinearOptimization[
  N[(Subscript[nxvar, 0] + Subscript[nxvar, 1] nxcoeff )], 
  xconstr, xvars(*,Tolerance-&gt;10^(-10)*),(*WorkingPrecision-&gt;
  MachinePrecision,*)Method -&gt; &quot;MOSEK&quot;, 
  PerformanceGoal -&gt; &quot;Speed&quot;];


nx0 = (Subscript[nxvar, 0] /. xres);
nx1 = (Subscript[nxvar, 1] /. xres);




keylength = 
 Floor[Max[(nx0 + nx1 (nxcoeff) - nx lambdaEC - 
     Log[1/(2.*eps^4.)]), 0.]];



keyrate = N[keylength/nb]


]];

keyrateFuncpaper[10, 0.55, 1/2, 1/4, 4/10, 0.0012] // AbsoluteTiming
(------------------3----------------------)
stepmul = 
  50.;(The higher the value, the faster the code runs;Ideally should 

be equal to 1)
step = 10.^-2.;
bitsDecoy = 2.;
bitsBasis = 2.;
stepBasis = 1./2.^bitsBasis;
stepDecoy = 1./2.^bitsDecoy;
mu3 = 2. 10^-04;
NdBpoints = 200; 
dBmin = 0.; 
dBmax = 65.;
keytabpaper = 
   ParallelTable[(Quiet[)
    Max[Table[
      keyrateFuncpaper[dB, qx, pmu1, pmu2, mu1, mu2], {mu1, 4/10, 1, 
       stepmul step}, {mu2, mu3 + 10^-3, Min[3/10, mu1 - mu3 - 10^-3],
        stepmul step}, {pmu1, 1/4, 1 - 2stepDecoy, stepDecoy}, {pmu2,
        stepDecoy, Min[1 - pmu1 - stepDecoy, 1/2], stepDecoy}, {qx, 
       1/2, 1 - stepBasis, stepBasis}]](])(]*), {dB, 0, 
     65, (dBmax - dBmin)/(NdBpoints - 1)}, 
    Method -> "FinestGrained"]; // AbsoluteTiming
ListLogPlot[{keytabpaper}, PlotRange -> Full, 
 PlotLegends -> {"paper"}]

Once the main function keyrateFuncpaper is defined, I'll call it over a range of its inputs. In the last block of the code I've defined stepmul which is a factor for the step size taken for the range of (two of the) inputs. Ideally this value should be 1, but it takes far too long for my use and I cant even run this on my machine (16 GB Mac air) because the system runs out of memory.

For ease of explaining, I've divided the code into three blocks. In the first block I've defined some functions. These are called in keyrateFuncpaper, defined in the second block.

keyrateFuncpaper is a function which takes 6 numerical values as input. The main components of this function is the three linear optimizations sections (there are comments at the start of each one). Each of these sections has the same structure; first I define the constraints, then the variables and finally the LinearOptimization. These three optimizations are very similiar, in terms of the constraints used.

In the final block, the code creates a table using the function we just defined. Out of the six inputs, I vary five of them and take the maximum function value for every value of the input dB. Finally the result is plotted.

Let me know I'm missed anything relevant or if you need any more info or explanation of the code. Thanks in advance!

Answer 1

For machine precision numbers, the actual backend of LinearOptimization will be some compiled and heavily optimized library. Such libraries require the problem to be formulated in some standard way. There are various "standard forms" for linear problems in use, but all of them formulate the optimization problem in terms of vectors and matrices - without any symbols involved altogether. One quite general one is this:

$$\text{Minimize} \quad c^T x \quad \text{subject to} \quad A \, x \geq a \quad \text{and} \quad B \,x = b,$$ where $A$ and $B$ are matrices and where $a$, $b$, $c$ are vectors. Moreover the inequality $A \, x \geq a$ is to be read componentwise.

So these libraries require $A$, $B$, $a$, $b$, $c$ to be submitted as input. There are several reasons for this practice:

1. These libraries are written in some strongly typed, down-to-earth languages like FORTRAN, C or C++. They do not have the symbolic capabilities that Mathematica users enjoy.

2. The formulation in terms of matrices and vectors allows to use highly optimized linear algebra routines. That makes them so efficient. In constrast, symbolic manipulations are just slow.

Hence what LinearOptimization has to do is to (i) take your expressions for the objective function, the constraints, and variables, convert them to suitable double precision matrices $A$, $B$ and vectors $a$, $b$, $c$; (ii) call the backend library with these arguments; and finally (iii) process and return the result of that library.

Of course, steps (ii) and (iii) cannot be avoided, if you want to use LinearOptimization. (You could try to write a LibraryLink function, that contains the whole loop over half a million problems. That would get you rid of quite some overhead, I think.)

But the developers at Wolfram Research were clever enough to allow you to avoid step (i) by offering the following syntax:

LinearOptimization[c, {A, a}, {B, b}]

This puts the responsibility to create $A$, $B$, $a$, $b$, $c$ onto your shoulders. But as an expert for the problem, you might be able to generate them directly without going through that many symbolic calculations. And as you say, you have many similar problems to solve. So maybe you can just modify these arrays each time, instead of generating entirely new ones (the latter might involve costly allocation operations). Maybe some of the data can even be resued. There is a lot of potential for short cuts!

I would like to show you what differences it does make if step (i) can be avoided. The following is to convert the problem into $A$, $B$, $a$, $b$, $c$. Note however that this is certainly not a very efficient way:

c = Developer`ToPackedArray[N[D[Subscript[nzvar, 1], {zsingvars, 1}]]];
inequalities = Flatten[
   zsingconstr /. {
     GreaterEqual[x_, y_] :> x - y, 
     GreaterEqual[x_, y_, z_] :> {x - y, y - z},
     LessEqual[x_, y_] :> y - x, 
     LessEqual[x_, y_, z_] :> {y - x, z - y},
     Equal[x_, y_] :> Nothing
     }
   ];
A = DeveloperToPackedArray[N[D[inequalities, {zsingvars, 1}]]]; a = DeveloperToPackedArray[
   N[inequalities /. Thread[zsingvars -> 0.]]];
equalities = Flatten[
   zsingconstr /. {
     GreaterEqual[x_, y_] :> Nothing, 
     GreaterEqual[x_, y_, z_] :> Nothing,
     LessEqual[x_, y_] :> Nothing, LessEqual[x_, y_, z_] :> Nothing,
     Equal[x_, y_] :> x - y
     }
   ];
B = DeveloperToPackedArray[N[D[equalities, {zsingvars, 1}]]]; b = DeveloperToPackedArray[N[equalities /. Thread[zsingvars -> 0.]]];

Here, Developer`ToPackedArray is to guarantee that the arrays are in packed form (in the same storage mode that is used in the C programming languag), so that no further conversions are needed. Now let's see how long it takes to solve the problem:

LinearOptimization[c, {A, a}, {B, b}, Method -> "CLP", PerformanceGoal -> "Speed"] // RepeatedTiming

{0.0004165, {382141., 1.7814710^8, 1.9347310^8, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -80951.6, 80951.6}}

For comparison, the original code took 0.0025653 seconds on my machine.

Btw., I use RepeatedTiming here in order to run this several times and to get an averaged timing. For small runtimes, this is more accurate than a single AbsoluteTiming. (The results should still be taken with a grain of salt, though.)