Speed up linear programming computation?

Question

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.

Kagaratsch · Answer

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.

Speed up linear programming computation?

One Answer

Add your own answers!

Ask a Question