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

Question

Bruno Luong 2024-9-12

2
链接

此问题的直接链接

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

编辑： Matt J 2024-10-21，22:45

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: []

8 个评论
显示 6更早的评论隐藏 6更早的评论

Julian Hall 2024-10-20

I’m not saying that your LP is numerically challenging to solve: I know that it isn’t, so have no reason to doubt your statement that the matrix condition is about 10. Leaving aside your problem’s LP condition (which is not the same as its constraint matrix condition) omitting the very small values in your matrix (HiGHS drops 1176 |values| in [1.20981e-70, 9.7101e-10]) won’t affect the numerical linear algebra. Hence, in this sense, the small values aren’t meaningful.

If the constraint matrix of an LP has small values that are meaningful – because they lead to a correspondingly large matrix condition number and sensitivity in the optimal solution – then solvers may struggle to solve the problem. For values of “1e-40”, this could be hopelessly hard. Hence small numbers are dropped by LP solvers to reduce the risk of trying to solve a problem with excessive constraint matrix condition – and if the problem was well conditioned anyway, the small values are safe to drop.

For problems like yours, where the largest matrix entries are of order unity and there are many small values (even after dropping excessively small values) the default scaling scheme in the HiGHS LP solver can result in small values being scaled up significantly in the LP that is solved internally – before scaling is removed and the solution cleaned up if necessary. It was my analysis of this that informed Derya’s response.

LP conditioning goes beyond the conditioning of the constraint matrix. Since an LP with a perfectly conditioned constraint matrix can have a non-unique solution – a situation that’s obviously not possible with a well-conditioned linear system of equations – small changes in the objective, RHS and matrix coefficients can lead to significant changes in the optimal solution.

Incorporating the HiGHS solvers as default in the Mathworks optimization toolbox is a significant enhancement in general, and has been done with due consideration by Mathworks. They wouldn’t be advocating solvers that are “flawed, if not buggy”.

Bruno Luong 2024-10-21，9:16

编辑：Bruno Luong 2024-10-21，11:09

在 MATLAB Online 中打开

linprog_test.mat

"the small values aren’t meaningful"

They are my data, not sure what you mean by meaningful. No I don't care of the small value they have it just happens to take values what they are.

Look my problem is well conditioned. I can drop the small values of Aeq, the solution is little perturbed, according to dual-simplex-legacy (and it can be showed mathematically the welll poseness of the problem theorically):

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 = linprog(c, [], [], Aeq, beq, LB, UB, linprogopt);
Optimal solution found.
% Drop Aeq elements smaller than "tol"
tol = 1e-10;
[i,j,Aval] = find(Aeq);
issmall = abs(Aval) < tol;
issmall = sub2ind(size(Aeq),i(issmall),j(issmall));
Aeq_perturb = Aeq;
Aeq_perturb(issmall) = 0;
% Solve new problem with perturbed Aeq
lpsol_perturb =  linprog(c, [], [], Aeq_perturb, beq, LB, UB, linprogopt);
Optimal solution found.
% Compare the solutions
for p = [1 2 Inf]
   RekLpErr = norm(lpsol_perturb-lpsol, p) / norm(lpsol, p);
    fprintf("Relative l%d norm error of the solution = %e\n", p, RekLpErr)
end
Relative l1 norm error of the solution = 3.978788e-07
Relative l2 norm error of the solution = 5.040052e-07
Relative lInf norm error of the solution = 7.460135e-07

I could also shift smalll values aways from 0 when using HIGS, but why should I doing such dirty trick?

Clearly HiGHS mess with the scaling and do something odd with small values that should not matter numerrically. This problem has been implemented since 10 years and work reliable with legacy method.

The HIGH method get disrupted by small values, and there is no obvious reason it should be so. I stand by my claim. This method HiGHS is flawed.

Bruno Luong 2024-10-21，9:24

编辑：Bruno Luong 2024-10-21，11:20

My workaround is still using legacy method. Sorry but HIGHS as default is not a wise decision from TMW IMHO.

Matt J 2024-10-21，15:05

编辑：Matt J 2024-10-21，22:45

在 MATLAB Online 中打开

linprog_test.mat

This seems like a more applicable test of whether the problem is well-conditioned. There does appear to be some sensitivity to it, though of course it may depend on epsilon

load('linprog_test.mat')
epsilon=1e-6;
linprogopt = optimoptions('linprog','Display','none','Algorithm','dual-simplex-legacy');
Warning: The 'dual-simplex-legacy' algorithm will be removed in a future release. To avoid this warning or a future error, choose a different algorithm: 'dual-simplex-highs', 'interior-point' or 'interior-point-legacy'.
[lpsol0, ~, exitflag, out] =  linprog(c, [], [], Aeq, beq, LB, UB, linprogopt);
size(lpsol0)
ans = 1×2
   467     1
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
dev=inf(1,10);
for i=1:numel(dev)
 [lpsol,~,exitflag] =  linprog(c+randn(size(c))*epsilon, [], [], Aeq, beq, LB, UB, linprogopt);
 if exitflag~=1, continue; end
 dev(i)=norm(lpsol-lpsol0,inf);
end
dev
dev = 1×10
    2.4327    0.0000    3.7828    0.0000    1.2906    0.0000    2.2675    0.0000    0.0000    0.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

请先登录，再进行评论。

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

Answer 1

Derya 2024-10-8

1
链接

此回答的直接链接

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

Hello all,

Thank you for providing the data for the failing problem, Bruno. And I'm sorry for the inconvenience.

As you noted there are two issues here: exitflag/interations/message inconsistency and dual-simplex-highs not finding a solution. We'll resolve the inconsistency in the exit condition shortly. We'll also address the issue with the dual-simplex algorithm (of HiGHS) that gets stuck at the initial iteration probably due to some automatic scaling done during the algorithm. Note that presolve doesn't do reductions in this case and is not the culprit:

Presolve : Reductions: rows 211(-0); columns 467(-0); elements 8741(-0) - Not reduced

Problem not reduced by presolve: solving the LP

Kind Regards,

Derya

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

Bruno Luong 2024-10-8

Hello Deyra,

Thank you very much for investigating and informmation.

Best regards

请先登录，再进行评论。

Answer 2

Matt J 2024-9-12

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

在 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

编辑：Bruno Luong 2024-9-13

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

编辑：Bruno Luong 2024-9-13

在 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 3

Catalytic 2024-9-13

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

编辑：Bruno Luong 2024-9-13

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.

请先登录，再进行评论。

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

8 个评论
显示 6更早的评论隐藏 6更早的评论

采纳的回答

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

更多回答（2 个）

2 个评论
显示无隐藏无

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

8 个评论 显示 6更早的评论隐藏 6更早的评论

采纳的回答

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

更多回答（2 个）

2 个评论 显示 无隐藏 无

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

8 个评论
显示 6更早的评论隐藏 6更早的评论

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

2 个评论
显示无隐藏无

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