Azzera filtri
Azzera filtri

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

86 visualizzazioni (ultimi 30 giorni)
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 Commenti
Bruno Luong
Bruno Luong il 21 Ott 2024 alle 9:24
Modificato: Bruno Luong il 21 Ott 2024 alle 11:20
My workaround is still using legacy method. Sorry but HIGHS as default is not a wise decision from TMW IMHO.
Matt J
Matt J il 21 Ott 2024 alle 15:05
Modificato: Matt J il 21 Ott 2024 alle 22:45
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>

Accedi per commentare.

Risposta accettata

Derya
Derya il 8 Ott 2024 alle 20:17
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 Commento
Bruno Luong
Bruno Luong il 8 Ott 2024 alle 20:34
Hello Deyra,
Thank you very much for investigating and informmation.
Best regards

Accedi per commentare.

Più risposte (2)

Matt J
Matt J il 12 Set 2024
Modificato: Matt J il 12 Set 2024
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 Commenti
Bruno Luong
Bruno Luong il 13 Set 2024
Modificato: Bruno Luong il 13 Set 2024
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
Bruno Luong il 13 Set 2024
Modificato: Bruno Luong il 13 Set 2024
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

Accedi per commentare.


Catalytic
Catalytic il 13 Set 2024
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 Commento
Bruno Luong
Bruno Luong il 13 Set 2024
Modificato: Bruno Luong il 13 Set 2024
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.

Accedi per commentare.

Categorie

Scopri di più su Systems of Nonlinear Equations in Help Center e File Exchange

Prodotti


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by