Linearity check with presolve returns an empty set

[schillic]

I was debugging a bigger example (here) that leads to a (probably) wrong result when using GLKP with presolve activated. (I cannot say for sure that this is wrong because the input set is 4D, but when using the vertex representation, the result is a proper set and all 2D projections are proper too.)

I wonder if the following check for !isapproxzero(β_primal) is correct because the error message indicates that β_primal should just be nonnegative.

Polyhedra.jl/src/linearity.jl

Lines 180 to 181 in 07567d4

	if rep isa HRepresentation && !isapproxzero((β_primal = MOI.get(model, MOI.ConstraintPrimal(), β_con);); kws...)
	verbose >= 1 && @info("The polyhedron is empty as $β_primal is negative.")

In my example I get:

julia> Polyhedra.detect_new_linearities(rep, solver; verbose=1)
[ Info: The polyhedron is empty as 3.0 is negative.

(I think 3.0 is not negative 😅)

So it would be useful to get a confirmation that the check is correct. If so, I guess I have to reduce the example to get a more helpful input.

[blegat]

Sorry for the delay on this. The output seems weird, as the constraint is LessThan:

Polyhedra.jl/src/linearity.jl

Line 170 in 07567d4

β_con = MOI.add_constraint(model, hull[end], MOI.LessThan(zero(T)))

The ConstraintPrimal should be nonpositive for a feasible solution so if it is not zero, it should be negative. A MWE would definitely help.

[schillic]

Below is a simplified version. It definitely has to do with presolve in GLPK, so this might be the wrong place to report. If you look at this plot,

you see that the set is a normal polygon with no numerical issues. The only thing to note is that there are two almost-orthogonal diagonal constraints, and the difference in orthogonality is probably below numerical precision. But I would say that this should not matter because they are not parallel but bound in different directions.

Some explanation: This comes from an elimination algorithm that first extends the set with a dummy dimension and then projects some dimensions away.

julia> using Polyhedra, CDDLib, MathOptInterface, GLPK
julia> Ab = [
 ([1.0, 0.0], 10.0)
 ([0.0, 1.0], 1.5)
 ([0.0, -1.0], 1.5)
 ([1.0, 1.5], 11.125)
 ([-1.0, -1.5000000000000002], -6.875)
 ([-4.0, -1.0], -24.5)
];
julia> M = [1. 0];
julia> m, n = size(M);
julia> ₋Id_m = hcat(-1.0);
julia> Ax_leq_b = [Polyhedra.HalfSpace(vcat(zeros(m), Abi[1]), Abi[2]) for Abi in Ab];
julia> y_eq_Mx = [Polyhedra.HyperPlane(vcat(₋Id_m[i, :], Vector(M[i, :])), 0.0) for i in 1:m];
julia> Phrep = Polyhedra.hrep(y_eq_Mx, Ax_leq_b);
julia> backend = CDDLib.Library(:float);
julia> Phrep = polyhedron(Phrep, backend);
julia> method = Polyhedra.BlockElimination();
julia> Peli_block = Polyhedra.eliminate(Phrep, (m+1):(m+n), method);

### With presolve, the result is not correct

julia> solver = MathOptInterface.OptimizerWithAttributes(GLPK.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.RawOptimizerAttribute("presolve") => 1]);
julia> Polyhedra.removeduplicates(Polyhedra.hrep(Peli_block), solver)
H-representation CDDInequalityMatrix{Float64, Float64}:
2-element iterator of HyperPlane{Float64, Vector{Float64}}:
 HyperPlane([1.0], 0.0)
 HyperPlane([0.0], 1.0)

### Without presolve, the result is correct (just contains redundant constraints)

julia> solver2 = MathOptInterface.OptimizerWithAttributes(GLPK.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[]);
julia> Polyhedra.removeduplicates(Polyhedra.hrep(Peli_block), solver2)
H-representation CDDInequalityMatrix{Float64, Float64}:
7-element iterator of HalfSpace{Float64, Vector{Float64}}:
 HalfSpace([-4.0], -23.0)
 HalfSpace([-1.0], -4.625)
 HalfSpace([-0.0], 3.0)
 HalfSpace([1.0], 13.375)
 HalfSpace([1.480297366166875e-16], 4.250000000000001)
 HalfSpace([-5.0], -25.625)
 HalfSpace([1.0], 10.0)

### The following part leads to the funny error message

julia> rep = Polyhedra.hrep(Peli_block);
julia> aff = Polyhedra._linearity_space(rep, true);
julia> Polyhedra._hasnonlinearity(rep)
true

julia> isnothing(Polyhedra.detect_new_linearities(rep, solver))
true

julia> Polyhedra.detect_new_linearities(rep, solver; verbose=1)
[ Info: The polyhedron is empty as 3.0 is negative.

[blegat]

Thanks for the detailed MWE, that should be easier to investigate

[schillic]

Here is another example. I first compute the convex hull of the two polygons below and then call detecthlinearity. Again this is problem does not occur without the presolve option.

julia> using Polyhedra, JuMP, GLPK
julia> solver = JuMP.optimizer_with_attributes(GLPK.Optimizer, "presolve" => GLPK.GLP_ON);
julia> backend = DefaultLibrary{Float64}(solver);
julia> A, b = ([0.5 1.5; 0.0 -1.0; 0.5 -0.5; -1.0 0.0], [4.0, -1.0, 0.0, 0.0]);
julia> P = polyhedron(hrep(A, b), backend);
julia> A, b = ([-1.0 0.0; 0.5 1.5; 0.0 -1.0; 1.0 0.0], [1.0, 4.0, -1.0, 0.0]);
julia> Q = polyhedron(hrep(A, b), backend);

julia> R = convexhull(P, Q)
Polyhedron DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}:
8-element iterator of Vector{Float64}:
 [0.0, 1.0]
 [0.0, 2.666666666666667]
 [1.0, 1.0]
 [2.0, 2.0]
 [0.0, 1.0]
 [0.0, 2.666666666666667]
 [-1.0, 1.0]
 [-1.0, 3.0000000000000004]

julia> detecthlinearity(hrep(R), solver)
H-representation MixedMatHRep{Float64, Matrix{Float64}}:
3-element iterator of HyperPlane{Float64, Vector{Float64}}:
 HyperPlane([1.0, 0.0], 0.0)
 HyperPlane([0.0, 1.0], 0.0)
 HyperPlane([0.0, 0.0], 1.0)