0
$\begingroup$

Consider setting up the following system for a linear programming optimization

sz = 50;
d = Table[Abs[i - j], {i, 1, sz}, {j, 1, sz}];
vd = (d // Flatten)[[2 ;;]] // N;
fs = Table[f[i, j], {i, 1, sz}, {j, 1, sz}];
vfs = (fs // Flatten)[[2 ;;]];
uv = Table[1, sz];
ClearAll[o]; ClearAll[t];
on = Table[o[i], {i, sz}];
tw = Table[t[i], {i, sz}];
eqs = {fs[[1, 1]], -on - fs.uv, -tw - uv.fs} /. Solve[-1. == fs /. List -> Plus, f[1, 1]][[1]] // Flatten;
{b, m} = CoefficientArrays[eqs, vfs];

so that, e.g., generating some random input

one = Table[RandomReal[], {i, 1, sz}];
one = one/Total[one];
two = Table[RandomReal[], {i, 1, sz}];
two = two/Total[two];

we can evaluate the optimization

AbsoluteTiming[
    Do[o[i] = one[[i]]; t[i] = two[[i]];, {i, sz}]; 
    out = LinearProgramming[vd, m, b, Method -> "InteriorPoint"]; 
    vd.out
]

{0.0143417, 2.13037}

As we can see, for systems of size sz = 50 one optimization takes about 0.015 seconds. My question is:

Is 0.015 seconds in this case a good performance, or may there be ways to speed up the optimization?

I need to evaluate 10^12 such optimization steps, which at that rate would take 475.647 years to complete. I wonder if I should be looking into performance improvements, or rather change my approach if no significant speedup can be expected to be achieved.

Improve this question
$\endgroup$
3
  • $\begingroup$ Why are you doing this: Do[o[i] = one[[i]]; t[i] = two[[i]];, {i, sz}] when you could just do this: o = one; t = two; and it seems like you don't appear to use o or t in the LinearProgramming? $\endgroup$
    – flinty
    Commented Nov 9, 2020 at 11:36
  • $\begingroup$ @flinty note that o and t are not lists. The function values o[i], t[i] enter the LinearProgramming through the predefined vector b and matrix m. Their values change dynamically when the components are assigned. $\endgroup$
    – Kagaratsch
    Commented Nov 9, 2020 at 12:02
  • 2
    $\begingroup$ It's an unusual way to write it. I can suggest some improvements here pastebin.com/QFN7JkPX where I've used Array, ConstantArray, and Normalize to put it in a more functional and less Table-like / loop style. Nevertheless, I don't think you can improve on 0.015 seconds. If you use a less strict Tolerance->.1 then it gets a 2x times faster at best. LP is $\mathcal{O}(n^{2})$ or $\mathcal{O}(n^{3})$ and I think Mathematica is already using a decent optimizer (COIN?) so there's not much you can do. $\endgroup$
    – flinty
    Commented Nov 9, 2020 at 18:19

1 Answer 1

Reset to default
1
$\begingroup$

Turns out, in this particular case one can find an equivalent solution via dynamic programming:

getMat = Compile[{{One, _Real, 1}, {Two, _Real, 1}},
   Module[{
     sz = Length[One],
     mat = DiagonalMatrix[Table[Min[One[[i]], Two[[i]]], {i, sz}]],
     dif = Two - One,
     ii = 1, jj = 2},
    While[ii + jj < 2 sz ,
     If[dif[[ii]] > 0 && dif[[jj]] < 0,
       If[dif[[ii]] <= -dif[[jj]], mat[[jj, ii]] = dif[[ii]]; 
         dif[[jj]] = dif[[jj]] + dif[[ii]]; dif[[ii]] = 0; 
         ii = ii + 1; jj = ii + 1;
         ,
         mat[[jj, ii]] = -dif[[jj]];
         dif[[ii]] = dif[[jj]] + dif[[ii]]; dif[[jj]] = 0; 
         If[jj < sz, jj = jj + 1;, ii = ii + 1; jj = ii + 1;];
         ];
       ,
       If[dif[[ii]] < 0 && dif[[jj]] > 0,
         If[-dif[[ii]] >= dif[[jj]], mat[[ii, jj]] = dif[[jj]]; 
           dif[[ii]] = dif[[jj]] + dif[[ii]]; dif[[jj]] = 0; 
           If[jj < sz, jj = jj + 1;, ii = ii + 1; jj = ii + 1;];
           ,
           mat[[ii, jj]] = -dif[[ii]]; 
           dif[[jj]] = dif[[jj]] + dif[[ii]];
           dif[[ii]] = 0; ii = ii + 1; jj = ii + 1;
           ];
         ,
         If[dif[[ii]] === 0, ii = ii + 1; jj = ii + 1;, 
           If[jj < sz, jj = jj + 1;, ii = ii + 1; jj = ii + 1;];];
         ];
       ];
     ]; mat], CompilationTarget -> "C"];

This evaluates about 100 times faster than the linear programming routine

AbsoluteTiming[
 mat = getMat[one, two];
 vd.(Flatten[mat][[2 ;;]])
 ]

{0.0002348, 1.31883}

And we can check that the output is the same up to numerical error

Table[one = Table[RandomReal[], {i, 1, sz}];
   one = one/Total[one];
   two = Table[RandomReal[], {i, 1, sz}];
   two = two/Total[two];
   Do[o[i] = one[[i]]; t[i] = two[[i]];, {i, sz}]; 
   out = LinearProgramming[vd, m, b, Method -> "InteriorPoint"];
   mat = getMat[one, two]; 
   vd.(out - Flatten[mat][[2 ;;]])
, {jjj, 1, 1000}]//Abs//Max

7.02253*10^-7

where 10^-6 is the default linear programming precision of Mathematica.

Improve this answer
$\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.