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

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

[Bruno Luong]

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]

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.

[Julian]

The new default 'dual-simplex-highs' algorithm used by linprog/intlinprog is certainly not flawed or buggy in general, as it drives the HiGHS mixed-integer programming solver that requires vast numbers of LPs to be solved and is justifiably very popular. This particular instance contains a very large range of values in the constraint matrix - even after the tiny values have been ignored. You can't have values of 1e-40 being meaningful relative to values of order one when solving LPs using double precision arithmetic. This isn't to say that HiGHS can't cope with a very large range of values in the constraint matrix. You're just increasing the probability of numerical failure.

[Bruno Luong]

" You can't have values of 1e-40 being meaningful relative to values of order one when solving LPs using double precision arithmetic"

Non sense the row of my constraints contain has l2 norm about unity. It happens to have small value of only certain columns (of that row). The matrix condition is about 10, there is nothing challenging to deal numerical with such matrix. The other solvers do not have any issue.

Background: this is L1 norm polynomial fitting. The constraint matrix contain polynomial value of a regression points, they are propertly centered and normalizeed that why I can have such a good condition number. However the absissa of some of my data points (after centering) can be small, but not exactly 0.

If one cannot work with arbitrary abscissa then numerical fitting would be useless. The way 'dual-simplex-highs' deal with resaling is flawed IMO.

[Julian]

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]

"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]

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

[Matt J]

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]

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

[Bruno Luong]

Hello Deyra,

Thank you very much for investigating and informmation.

Best regards

Answer 2

[Matt J]

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

[Bruno Luong]

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]

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]

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 =

     []

[Bruno Luong]

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.