I have an A.X-B linear system which I have transformed into A.x-b to be solvable as a linear program using the code below:

Linearize[A_, B_] :=
 {m, n, id, Aid, xv, bv},
 {m, n} = Dimensions[A];
 id = SparseArray[{{i_, i_} -> 1}, {n, n}];
 Aid = SparseArray[
  Flatten[Map[Map[Function[t, Times[t, #]], id, {2}] &, 
  A], {{1, 2}, {3, 4}}]];
 bv = SparseArray[Flatten[B]];
 {Aid, bv}

I have implemented a 1- and an Infinity- norm minimization algorithm based on this article, with the difference that I pass the arguments as a list of two elements to make it compatible with the function above and (because of the crash described below) I have made all the arrays sparse:

L1Solve[{A_, b_}] :=
 {m, n, Aall, ball, c, x, bds, id},
 {m, n} = Dimensions[A];
 id = SparseArray[{{i_, i_} -> 1}, {m, m}];
 Aall = SparseArray[
  Join[Transpose[Join[id, Transpose[-A]]], 
  Transpose[Join[id, Transpose[A]]]]];
 ball = SparseArray[Join[-b, b]];
 c = SparseArray[Join[Table[1, {m}], Table[0, {n}]]];
 bds = SparseArray[Table[{-Infinity, Infinity}, {m + n}]];
 x = LinearProgramming[c, Aall, ball, bds];
 Take[x, -n]

The Infinity-norm algorithm works fine, and so does the 1-norm for a small test found in the article guide:

mat = {{1, 2}, {5, 6}, {4.5, 6}};
vec = {5, 6, 8};
L1Solve[{mat, vec}]

But when I run it for the full dataset, which is an 500x3 matrix, the kernel crashes with no error. Memory usage remains below 50% during this and I am not sure what else to try. What do you suggest trying? Thanks.

  • 2
    $\begingroup$ Might be a problem in the default method (which I am guessing is simplex). Could try setting Method->"InteriorPoint" or perhaps a different method. $\endgroup$ – Daniel Lichtblau Dec 30 '18 at 15:08
  • $\begingroup$ Thanks for the suggestion. There are three Methods available: Simplex, RevisedSimplex and InteriorPoint. InteriorPoint is the one that crashes on the large dataset, which is the default for problems of this size. I think Simplex is unable to work with Sparse Arrays and gives an error. I got a solution with RevisedSimplex, but it took a couple of minutes and I have an even bigger dataset after this one. All three methods work fine on the small example, even with Sparse Arrays. And LinearProgramming cannot be parallelized. Any other suggestions? PS Should I edit to include this info in the post? $\endgroup$ – sgkion Dec 30 '18 at 19:16
  • 2
    $\begingroup$ (1) Yes, please put the details in the post as to which methods behave in which ways. (2) Also there is an undocumented Method->"CLP" which uses the interior point code in the COIN-CBC library. (3) Possibly those sparse arrays can be densified, via Normal, so that simplex can be used. (4) If possible, add the actual arrays in the crashing example. They would be useful for purposes of diagnosing the issue. $\endgroup$ – Daniel Lichtblau Dec 30 '18 at 20:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.