Speed up linear programming computation?

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]][] // 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.

• 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? Nov 9 '20 at 11:36
• @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. Nov 9 '20 at 12:02
• 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. Nov 9 '20 at 18:19

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.