3
$\begingroup$

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 >= 
         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}]]];
    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"];
    
    
    (*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 >= 
         mzvec[[j]] + N[Subscript[\[Delta]\[Mu]mz, j]] >= 
         sumsubmz, {j, 3}],
       {Sum[N[Subscript[\[Delta]\[Mu]mz, j]], {j, 3}] == 0.},
       Table[
        -Sqrt[-Log[epshoeff/2.] mz/2.] <= 
         N[Subscript[\[Delta]\[Mu]mz, j]] <= 
         Sqrt[-Log[epshoeff/2.] mz/2.], {j, 3}],
       Table[
        0. <= N[Subscript[mzvar, l]] <= 
         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->10^(-10)*),(*WorkingPrecision->MachinePrecision(*,
      Tolerance->10^(-10)*),*)Method -> "MOSEK", 
      PerformanceGoal -> "Speed"];
    
    
    
    
    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. < phix < 1.,
        -phix*Log[2., phix] - (1. - phix) Log[2., 1. - phix], 0.
        ])(*,60]*);
    lambdaEC = fEC If[0. < eobs < 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 >= 
         nxvec[[j]] + N[Subscript[\[Delta]\[Mu]nx, j]] >= sumsubx, {j,
          3}],
       {Sum[N[Subscript[\[Delta]\[Mu]nx, j]], {j, 3}] == 0.},
       Table[-Sqrt[-Log[epshoeff/2.] nx/2.] <= 
         N[Subscript[\[Delta]\[Mu]nx, j]] <= 
         Sqrt[-Log[epshoeff/2.] nx/2.], {j, 3}],
       Table[
        0 <= N[Subscript[nxvar, l]] <= 
         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->10^(-10)*),(*WorkingPrecision->
      MachinePrecision,*)Method -> "MOSEK", 
      PerformanceGoal -> "Speed"];
    
    
    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 - 2*stepDecoy, 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!

$\endgroup$
8
  • 1
    $\begingroup$ That's a huge code dump with no explanation. You shouldn't expect us to wade through your code and understand it in order to even start helping. Please reduce this down to a simpler example if you can, or thoroughly explain what each portion of your code does, and what problem you are solving with it. $\endgroup$
    – MarcoB
    Aug 29, 2023 at 12:45
  • $\begingroup$ Alright,Ive updated the question $\endgroup$
    – Dotman
    Aug 29, 2023 at 13:22
  • 2
    $\begingroup$ Don't use subscript variable names! newConstr = zsingconstr /. {Subscript[\[Delta]\[Mu]nz, x_] :> \[Delta]\[Mu]nz[x], Subscript[nzvar, x_] :> nzvar[x]}; newVars = zsingvars /. {Subscript[\[Delta]\[Mu]nz, x_] :> \[Delta]\[Mu]nz[x], Subscript[nzvar, x_] :> nzvar[x]}; And now LinearOptimization[nzvar[1], newConstr, newVars] works :-D $\endgroup$
    – ydd
    Aug 29, 2023 at 15:20
  • 3
    $\begingroup$ 3. Avoid using subscripted symbols in your code $\endgroup$
    – Roman
    Aug 29, 2023 at 17:07
  • 1
    $\begingroup$ Your code works fine for me with subscripts. However, the LinearOptimization only takes 5ms for me. Is this really too slow? $\endgroup$
    – Carl Woll
    Aug 30, 2023 at 20:57

1 Answer 1

6
+50
$\begingroup$

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 = Developer`ToPackedArray[N[D[inequalities, {zsingvars, 1}]]];
a = Developer`ToPackedArray[
   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 = Developer`ToPackedArray[N[D[equalities, {zsingvars, 1}]]];
b = Developer`ToPackedArray[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.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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