Opened 6 years ago
Last modified 5 years 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 )
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 6 years ago by
comment:2 Changed 6 years ago by
- Dependencies set to #20406
comment:3 Changed 6 years ago by
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^{3}=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 6 years ago by
- 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=7^{3}, 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: ↓ 6 Changed 6 years ago by
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: ↓ 7 Changed 6 years ago by
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: ↓ 8 Changed 6 years ago by
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: ↓ 9 Changed 6 years ago by
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 askis_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 6 years ago by
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 askis_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 6 years ago by
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 6 years ago by
- 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 5 years ago by
So isn't this essentially a duplicate of #18765 then?
comment:13 Changed 5 years ago by
- 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 5 years ago by
- Status changed from needs_info to needs_work
(See #20406 where this was discovered.)