[Octave-bug-tracker] [bug #66010] Function glpk produces incorrect output

[Hendrik K]

Follow-up Comment #6, bug #66010 (group octave):

I looked at this and key issue seems to be that (some of the) data elements in
A and in b are smaller than eps, so essentially it is a dodgy input data
issue.
 

find(and(abs(A)<= eps(), A!= 0))'
ans =
 Columns 1 through 26:
   156   187   218   249   280   311   342   373   404   435   456   461   487
  492   518   523   549   554   580   585   606   611   616   637   642   647
 Columns 27 through 35:
   668   673   678   699   704   709   730   735   740

find(and(abs(b) < eps(), b!=0))'
ans =
    8    9   10   14   15   20



Very funny things can happen in numerical mathematics in such cases.

For example the glpk prescaling (set automatically by octave calling the api
and available via log info printing (set param.msglev = 3)) shows ultra large
prescaling ratios due to the values in A smaller than eps which are prone to
cause subsequent issues (e.g. a very small number < eps multiplied with such a
very large number will have a material error compared to the analytical
correct result...)
 

 A: min|aij| =  7.704e-34  max|aij| =  1.000e+00  ratio =  1.298e+33
GM: min|aij| =  5.308e-09  max|aij| =  1.884e+08  ratio =  3.549e+16
EQ: min|aij| =  2.817e-17  max|aij| =  1.000e+00  ratio =  3.549e+16


Possible solutions
==================
(1) Data cleansing
If one sets all A values smaller than eps() to zero, the issue disappears

A(abs(A) <= eps()) = 0;
leading to 
(A*xmin)(24) =  0.1412


(2) Appropriate glpk configuration
For example and depending on the use case, it may be possible to deactivate
the presolver (param.presol = 0) as in general this is risky according to the
glpk documentation: 
"Note that the presolver is able to recover only optimal solutions. If a
computed solution is infeasible or non-optimal, the corresponding solution of
the original problem cannot be recovered and therefore remains undefined. If
you need to know a basic solution even if it is infeasible or non-optimal, the
presolver should be disabled." 

This also leads to a solution without this issue.

param.presol = 0;
[xmin,fmin,status,extra]=glpk(c,A,b,lb,ub,ctype,vartype,sense,param);
(A*xmin)(24) = 0.1412


Note that I am not very skilled in linear programming and neither very
knowledgeable in glpk, but there are quite a few control parameter in glpk, so
I suggest having a deeper look which are the most appropriate one(s) for the
data and use case in question.   


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?66010>

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/

signature.asc
Description: PGP signature