Opened 3 years ago

Last modified 23 months ago

#20447 needs_work defect

Wrong result from delsarte_bound_additive_hamming_space with GLPK exact simplex

Reported by: mkoeppe Owned by:
Priority: major Milestone: sage-wishlist
Component: coding theory Keywords: lp
Cc: dimpase Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #20406 Stopgaps:

Description (last modified by mkoeppe)

delsarte_bound_additive_hamming_space should be guarded, using epsilons, against floating point fuzz that will appear with numerical solvers and even the GLPK exact solver because of the nature of its interface.

sage: delsarte_bound_additive_hamming_space(19,15,7,isinteger=False)
3
sage: from sage.numerical.backends.generic_backend import get_solver
sage: def glpk_exact_solver():                                               
        b = get_solver(solver="GLPK")
        b.solver_parameter("simplex_or_intopt", "exact_simplex_only")
        return b
sage: delsarte_bound_additive_hamming_space(19,15,7,solver=glpk_exact_solver,isinteger=False)
glp_exact: 54 rows, 20 columns, 795 non-zeros
...
2

Change History (14)

comment:1 Changed 3 years ago by mkoeppe

(See #20406 where this was discovered.)

comment:2 Changed 3 years ago by dimpase

  • Dependencies set to #20406

comment:3 Changed 3 years ago by dimpase

the function performs the following, in a loop: maximises the sum of variables, gets the value of this maximum, say M, then rounds M down to the closest m:=7d, adds a constraint that the sum of the variables is at most m, and repeats. In the case of rational variables, this process obviously must stop after one iteration (it need not stop if all the variables are integers). When it stops, it returns the latest d.

glp_exact: 54 rows, 20 columns, 795 non-zeros
GNU MP bignum library is being used
      0:   infsum =                      1   (14)
      1:   infsum =                      0   (14)
*     1:   objval =                      1   (14)
*    19:   objval =       1671.30573248408   (0)
OPTIMAL SOLUTION FOUND
glp_exact: 54 rows, 20 columns, 795 non-zeros
GNU MP bignum library is being used
      0:   infsum =                      1   (14)
      1:   infsum =                      0   (14)
*     1:   objval =                      1   (14)
*    17:   objval =                    343   (0)
OPTIMAL SOLUTION FOUND

here it should have been done, and return 3 (as 73=343), but it is not!

glp_exact: 54 rows, 20 columns, 795 non-zeros
GNU MP bignum library is being used
      0:   infsum =                      1   (14)
      1:   infsum =                      0   (14)
*     1:   objval =                      1   (14)
*    16:   objval =                     49   (0)
OPTIMAL SOLUTION FOUND
2

comment:4 Changed 3 years ago by dimpase

  • Report Upstream changed from N/A to Not yet reported upstream; Will do shortly.
  • Status changed from new to needs_info

in its infinite wisdom, GLPK returns the result of an exact computation as a float. And then the line

      if q_base**(m+1) == bd:

compares, for equality, the exact value, 343=73, with the float 342.999999999999943, and founds them unequal. Whereas for the correct functioning of the algorithm, they must be equal.

Not sure whether this is an upstream bug, or an upstream feature...

comment:5 follow-up: Changed 3 years ago by mkoeppe

How come you never run into trouble with the floating-point based solvers with this kind of code?

It is certainly ironic that the exact solver has more floating-point fuzz than the floating-point solvers. But your code cannot be correct if solver is a floating-point solver.

And there is an apparent upstream issue with GLPK. I wouldn't call it a bug, but given the design decision to return the exact values as double-floats, it should at least be improved to return all integers that have an exact representation in double-float as such.

comment:6 in reply to: ↑ 5 ; follow-up: Changed 3 years ago by dimpase

Replying to mkoeppe:

How come you never run into trouble with the floating-point based solvers with this kind of code?

of course I did, and docs explicitly say that as soon as you specify another solver, you are on your own.

- solver – the LP/ILP solver to be used. Defaults to PPL. It is arbitrary precision, thus there will be no rounding errors. With other solvers (see MixedIntegerLinearProgram for the list), you are on your own!

It is certainly ironic that the exact solver has more floating-point fuzz than the floating-point solvers. But your code cannot be correct if solver is a floating-point solver.

I don't think there is an API to check whether a solver is exact, and so I never bothered to check this in the code.

And there is an apparent upstream issue with GLPK. I wouldn't call it a bug, but given the design decision to return the exact values as double-floats, it should at least be improved to return all integers that have an exact representation in double-float as such.

this design is called lazyness, in CS, and not only :-)

comment:7 in reply to: ↑ 6 ; follow-up: Changed 3 years ago by mkoeppe

Replying to dimpase:

I don't think there is an API to check whether a solver is exact, and so I never bothered to check this in the code.

You can query the base_ring of the MIP and then ask is_exact.

comment:8 in reply to: ↑ 7 ; follow-up: Changed 3 years ago by dimpase

Replying to mkoeppe:

Replying to dimpase:

I don't think there is an API to check whether a solver is exact, and so I never bothered to check this in the code.

You can query the base_ring of the MIP and then ask is_exact.

well, is_exact is a bit too much to ask, one merely needs extra precision. (one can figure out how much, in fact). So you can use a backend that allows you to set the base ring to e.g. RealField(2000).

Besides, I don't think there are places in Sage that forbid doing things cause they are "risky".

comment:9 in reply to: ↑ 8 Changed 3 years ago by mkoeppe

Replying to dimpase:

Replying to mkoeppe:

Replying to dimpase:

I don't think there is an API to check whether a solver is exact, and so I never bothered to check this in the code.

You can query the base_ring of the MIP and then ask is_exact.

well, is_exact is a bit too much to ask, one merely needs extra precision. (one can figure out how much, in fact). So you can use a backend that allows you to set the base ring to e.g. RealField(2000).

Perhaps you mean a RealIntervalField here, because certainly RealField(2000) does not guarantee that the result of some unspecified numerical algorithm such as the implementation of the simplex method in the solver leads to results with 2000 correct bits, just like double-floats don't guarantee 53 correct bits.

Besides, I don't think there are places in Sage that forbid doing things cause they are "risky".

I think there more places should forbid things like that, or at least display a warning. For example, the polyhedral code in Sage is written in a way that it assumes exact arithmetic -- and if fed with floating point numbers, leads to mysterious errors. One would usually assume from algorithms that accept floating point numbers that they have some (however naive) accommodation for floating point fuzz, in the form of some epsilons.

comment:10 Changed 3 years ago by mkoeppe

The relevant code is in GLPK's glpapi07.c; it's using mpq_get_d to store the rational results in a double.

comment:11 Changed 3 years ago by mkoeppe

  • Milestone changed from sage-7.2 to sage-wishlist
  • Report Upstream changed from Not yet reported upstream; Will do shortly. to N/A

Actually, storing the rational results as doubles is fine, in fact I have a test in glpk_backend for that. The error is happening when the big coefficients are converted into doubles, and then reconstructed by GLPK's exact simplex.

So there's no upstream bug to be reported -- except that we still really want GLPK to make the interface to glpssx.h public (#18765).

Guarding your code against floating-point fuzz is a wishlist item -- I'm marking this ticket as such.

comment:12 Changed 23 months ago by jdemeyer

So isn't this essentially a duplicate of #18765 then?

comment:13 Changed 23 months ago by mkoeppe

  • Description modified (diff)

No, #18765 is about making a "proper" rational GLPK API available.

I've changed the description of the present ticket to say what should be done.

comment:14 Changed 23 months ago by mkoeppe

  • Status changed from needs_info to needs_work
Note: See TracTickets for help on using tickets.