Wrong result from delsarte_bound_additive_hamming_space with GLPK exact simplex

[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

Depends on #20406

CC: @dimpase

Component: coding theory

Issue created by migration from https://trac.sagemath.org/ticket/20447

[mkoeppe]

comment:1

(See #20406 where this was discovered.)

[dimpase]

Dependencies: #20406

[dimpase]

comment:3

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:=7^d, 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 7³=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

[dimpase]

comment:4

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=7³, 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...

[dimpase]

Upstream: Not yet reported upstream; Will do shortly.

[mkoeppe]

comment:5

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.

[dimpase]

comment:6

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 :-)

[mkoeppe]

comment:7

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.

[dimpase]

comment:8

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

[mkoeppe]

comment:9

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.

[mkoeppe]

comment:10

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

[mkoeppe]

comment:11

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.

[mkoeppe]

Changed upstream from Not yet reported upstream; Will do shortly. to none

[jdemeyer]

comment:12

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

[mkoeppe]

comment:13

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.