Chasing what is wrong with 'dual-simplex-highs' in linprog

Question

Bruno Luong 2024-9-12，22:27

2
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2152430-chasing-what-is-wrong-with-dual-simplex-highs-in-linprog

编辑： Bruno Luong 2024-9-16，17:49

linprog_test.mat

I try to see why 'dual-simplex-highs' algorithm fails and 'dual-simplex-legacy' works OK on this specific LP problem of size 467.

The linear programming involves only linear equality constraints, and lower bounds x >= 0 on some components of x (but not all of them).

The Aeq size is 211 x 467 and the condion number is not high at all IMO (about 10). So I consider it is not a difficult problem numerically (?).

'dual-simplex-legacy' able to find the solution, however not the default algorithm 'dual-simplex-highs', the output does not help much what is wrong.

Can someone tell me where I could investigate further to the cause?

load('linprog_test.mat')
size(Aeq)
ans = 1×2
   211   467
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
cond(full(Aeq))
ans = 10.3680
linprogopt = optimset('Algorithm', 'dual-simplex-legacy');
[lpsol, ~, exitflag, out] =  linprog(c, [], [], Aeq, beq, LB, UB, linprogopt)
Optimal solution found.
lpsol = 467×1
   -0.0599
    0.2955
   -0.4765
    4.5536
   -0.1640
   -3.2941
   -0.7875
   -0.8662
    0.2426
    1.5786
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
exitflag = 1
out = struct with fields:
         iterations: 371
          algorithm: 'dual-simplex-legacy'
    constrviolation: 2.2172e-08
            message: 'Optimal solution found.'
      firstorderopt: 2.4052e-07
linprogopt = optimset('Algorithm', 'dual-simplex-highs');
[lpsol, ~, exitflag, out] =  linprog(c, [], [], Aeq, beq, LB, UB, linprogopt)
Solver stopped prematurely.


lpsol =

     []
exitflag = 0
out = struct with fields:
         iterations: -1
          algorithm: 'dual-simplex-highs'
            message: 'Solver stopped prematurely.'
    constrviolation: []
      firstorderopt: []

2 个评论
显示无隐藏无

Bruno Luong 2024-9-13，17:43

编辑：Bruno Luong 2024-9-13，17:45

After some more investigation, I conclude that the new default 'dual-simplex-highs' algorithm used by linprog/intlinprog is flawed, if not buggy, at least in this specific test case in the above question.

Bruno Luong 2024-9-16，16:07

编辑：Bruno Luong 2024-9-16，17:49

Another BUG of this function is the wrong exitflag. In the given example the error message is "Solver stopped prematurely". and exitflag is 0. In the documentation page it states

"0 Number of iterations exceeded options.MaxIterations or solution time in seconds exceeded options.MaxTime."

Which is obviously incorrect since the number iteration perfomred retunred by linprog is -1.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2024-9-12，23:17

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2152430-chasing-what-is-wrong-with-dual-simplex-highs-in-linprog#answer_1515745

编辑：Matt J 2024-9-12，23:19

在 MATLAB Online 中打开

linprog_test.mat

It doesn't like your super-small Aeq values.

load('linprog_test.mat')
Aeq(abs(Aeq)<1e-8)=0;
linprogopt = optimset('Algorithm', 'dual-simplex-highs');
[lpsol, ~, exitflag, out] =  linprog(c, [], [], Aeq, beq, LB, UB, linprogopt)
Optimal solution found.
lpsol = 467×1
   -0.0599
    0.2955
   -0.4765
    4.5536
   -0.1640
   -3.2941
   -0.7875
   -0.8661
    0.2427
    1.5786
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
exitflag = 1
out = struct with fields:
         iterations: 305
          algorithm: 'dual-simplex-highs'
    constrviolation: 5.6899e-15
            message: 'Optimal solution found.'
      firstorderopt: 6.8412e-13

2 个评论
显示无隐藏无

Bruno Luong 2024-9-13，7:04

编辑：Bruno Luong 2024-9-13，8:26

Interesting, in principle the elements of my matrix is a 2D polynomial evaluation of a normalized coordinates (x,y), and it can get arbitrary small value when the coordinates is close to (0,0). Back ground I do 2D polynomiam L1 fitting. Not sure exactly why HIGS algorithm presolver fails due to that. Rather than changing Aeq, I might change (x,y) input coordinates to (0,0) and in turns make corresponding Aeq elemenrs 0s.

The legacy presolver seems to work OK under the wider dynamic range in Aeq. The 'interior-point' algorithm is also working.

Bruno Luong 2024-9-13，9:08

编辑：Bruno Luong 2024-9-13，9:29

在 MATLAB Online 中打开

linprog_test_2.mat

I test another case where there is element of Aeq that as small as same order of 1e-44 and HIGHS is still able to find solution. So it seems the dynamic range of Aeq is not a direct cause of the failure.

load('linprog_test_2.mat')
size(Aeq)
ans = 1×2
   279   603
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
cond(full(Aeq))
ans = 13.3418
[i,j,s] = find(Aeq); 
[smin,kmin] = min(abs(s));
imin = i(kmin), % row of the minimum element
imin = 15
jmin = j(kmin), % column of the minimum element
jmin = 45
smin % The smallest value of Aeq
smin = 6.2776e-44
max(abs(Aeq(imin,:))) % max value of the sale row of min value
ans = sparse double
   (1,1)        1
max(abs(Aeq(:,jmin)))  % max value of the same column of min value
ans = sparse double
   (1,1)       0.4686
% linprog on the wide dynamic range
linprogopt = optimset('Algorithm', 'dual-simplex-highs');
[lpsol, ~, exitflag, out] =  linprog(c, [], [], Aeq, beq, LB, UB, linprogopt)
Optimal solution found.
lpsol = 603×1
    0.0028
    0.1305
   -0.2128
    1.1160
   -0.0364
   -0.8465
   -0.0889
   -0.1347
    0.0176
    0.1980
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
exitflag = 1
out = struct with fields:
         iterations: 495
          algorithm: 'dual-simplex-highs'
    constrviolation: 1.4867e-09
            message: 'Optimal solution found.'
      firstorderopt: 1.1886e-09

请先登录，再进行评论。

Answer 2

Catalytic 2024-9-13，0:29

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2152430-chasing-what-is-wrong-with-dual-simplex-highs-in-linprog#answer_1515760

在 MATLAB Online 中打开

linprog_test.mat

intlinprog offers some additional diagnostic output -

load('linprog_test.mat')
intlinprog(c,[], [], [], Aeq, beq, LB, UB)
Running HiGHS 1.7.0: Copyright (c) 2024 HiGHS under MIT licence terms
WARNING: LP matrix packed vector contains 1176 |values| in [1.20981e-70, 9.7101e-10] less than or equal to 1e-09: ignored
Coefficient ranges:
  Matrix [1e-09, 1e+00]
  Cost   [1e+00, 1e+00]
  Bound  [0e+00, 0e+00]
  RHS    [1e-03, 1e+00]
Presolving model
211 rows, 467 cols, 8741 nonzeros  0s
211 rows, 467 cols, 8741 nonzeros  0s
Presolve : Reductions: rows 211(-0); columns 467(-0); elements 8741(-0) - Not reduced
Problem not reduced by presolve: solving the LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 211(29.8113); Du: 0(3.82173e-08) 0s
Model   status      : Not Set
HiGHS run time      :          0.01

Solver stopped prematurely.

No integer variables specified.


ans =

     []

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Bruno Luong 2024-9-13，6:34

编辑：Bruno Luong 2024-9-13，7:21

Thanks it helps. INTLINPROG uses the same so callled "highs" algorithm (see also wiki) and possibly uses the same presolve. I still uncertain if it is due to my data or algorithm.

I'll keep this algo un observation, and I might switch to legacy if the issue occurs again with different data.

请先登录，再进行评论。