Incorrect intersection (not sure if due to underlying optimizer)

[tomerarnon]

I was testing some code related to my recent suggestion in #1176 when I ran into this. In particular, I was plotting the intersection of an HPolytope with the positive orthant; see below.

using LazySets, Polyhedra, CDDLib, Optim, LinearAlgebra
import LazySets: dim, HalfSpace
using Combinatorics, Plots

plotproj!(p, x, v; kw...) = plot!(p, convert(HPolytope, LazySets.project(x, v)); kw...)
axesplot() = plot([0], seriestype = [:vline, :hline], l = (:black, :dash, 2))
positive_orthant(X) = positive_orthant = HPolytope(Diagonal(-ones(dim(X))), zeros(dim(X)))


# I created this particular hpolytope with rand(HPolytope, dim = 3), but since the 
# problem only happens some of the time, I'm pasting this specific one here. 
# Other polytopes also displayed this problem, but it isn't universal.
A = [ -0.325001    0.162046    0.0255077
 -0.0251797   0.083553   -0.0196931
  0.0509044  -0.10032    -0.0377418
 -0.0589181   0.0268046   0.0580025
  0.10905    -0.042451   -0.0435752
  0.0160815   0.0276424   0.0240291]
b = [ -0.009374874516494713
  0.014763719930503618
  0.13525877586303783
  0.0022307775191686076
  0.08895373078378208
  0.011333472249832692]

hp = HPolytope(A, b)

# both methods of intersection lead to this problem sometimes. This made me think the 
# problem was in the projection (which would be fine for me in this case). 
# But checking the support vector confirms there is a discrepancy.
hpp = intersection(hp, positive_orthant(hp))
# hpp = HPolytope([hp.constraints; [HalfSpace((a = zeros(3); a[i]=-1; a), 0.0) for i in 1:3]])

plots = []
for v in combinations(1:3, 2)
    p = axesplot()
    plotproj!(p, hp, v, title = string(v))
    plotproj!(p, hpp, v)
    push!(plots, p)
end

plot(plots..., legend = nothing)

# should be the same:
ρ(ones(3), hpp), ρ(ones(3), hp) # (0.7047530240332701, 0.9896255059528571)

Note that the projection operation in plotproj yields the warning

glp_simplex: unable to recover undefined or non-optimal solution

which is the same one mentioned in #2278 (comment). This made me think at first the problem is optimizer related, but I confirmed it occurs only during the projection, and therefore the support vector test demonstrates that the problem may be elsewhere (or overlapping?)

[schillic]

Quick comments without looking at the rest:

positive_orthant(X) = positive_orthant = HPolytope(Diagonal(-ones(dim(X))), zeros(dim(X)))

1. The positive_orthant in the middle is strange.
2. This should be an HPolyhedron because the set is unbounded.

[tomerarnon]

1. The positive_orthant in the middle is strange.

It was copy pasted from a line into a function, which is why the variable name is there. It doesn't make a difference though (confirmed)

2. This should be an HPolyhedron because the set is unbounded.

Fair enough (the constructor allows it, but I guess the finer point is that certain assumptions may fail later on). That change also does not make a difference (also confirmed)

[mforets]

Fair enough (the constructor allows it, but I guess the finer point is that certain assumptions may fail later on)

yes. we've recently added a more efficient method for checking boundedness, so we may reconsider to check that HPolytopes are bounded upon construction (with an optional flag, as usual).

[mforets]

intersection of HPolygon is buggy (see here). maybe related?

[tomerarnon]

Does that problem extend to hpolytopes in general? The set I'm working with was 3D (and later 4D) to avoid the pitfalls of simplifications that are only possible in 2D. After calling intersection and Intersection, I also tried manually concatenating the constraints (this should work too, right?) and all 3 methods led to this sort of issue.

Edit:
Although it seems truly impossible that the error is projection/plotting related, I'm going through the process now of plotting the 3D sets directly using MeshCat to confirm absolutely.

[mforets]

Does that problem extend to hpolytopes in general?

actually my comment is not relevant because the method used in the OP is intersection of HPolytope's which works fine.

and all 3 methods led to this sort of issue.

after thinking this further, i don't see a problem in the obtained result. when you plot plotproj!(p, hp, v, title = string(v)) that's projecting hp onto the plane, say [1, 2] in such a way that if the set hp has some parts that go below that plane, they are also projected back to that plane. on the other hand, if you plot after intersecting with the positive orthant, such pieces are chopped out. that explains that the green is always an overapproximation of the violet one (and it's worth mentioning that both computations, and the respective plots, are exact).

it would be nice to visualize this argument with a 3d plot but it's atm broken upstream :/

---

as a side note regarding the glp_simplex warnings, in this case at least it is a harmeless warning, i checked by using another backend. it's a warning issued upstream (from the default Polyhedra library); in general, i would recommend to use CDDLib.

currently the function project (defined here) doesn't let on specify the backend, though that should be easy to add.

calling directly the linear map,

z = linear_map(M, hpp, backend=CDDLib.Library())
z = linear_map(M, hpp, algorithm="vrep", backend=CDDLib.Library())

but these don't work as i'd expect, i think that's because we have to fix the backend argument, #2297

on the other hand, doing directly

z = linear_map(M, VPolytope(vertices_list(hpp, backend=CDDLib.Library())));

does return the same set as with the default backend.

[schillic]

on the other hand, if you plot after intersecting with the positive orthant, such pieces are chopped out

Just to clarify: in general project(intersection(X, Y)) = intersection(project(X), project(Y)) does not hold.

[mforets]

here is an example in 2D:

using LazySets: project, Interval

x = VPolygon([[-1.0, 1.0], [-0.5, -0.1], [2.0, -1.0], [0.5, 1.0]])
p = convert(Interval, project(x, [1]))
y = intersection(x, positive_orthant(x))
q = convert(Interval, project(y, [1]));

plot([0], seriestype = [:vline, :hline], l = (:black, :dash, 2), lab="")
plot!(x, lab="x", legend=:topleft, c=:lightblue)
ϵ = 0.03
lp = LineSegment([min(p), -ϵ], [max(p), -ϵ])
plot!(lp, lw=3.0, lab="proj x ", c=:green)
lq = LineSegment([min(q), ϵ], [max(q), ϵ])
plot!(lq, lw=3.0, lab="proj int x, O", c=:blue)

[tomerarnon]

in general project(intersection(X, Y)) = intersection(project(X), project(Y)) does not hold.

Indeed I know! I am confirming whether this is what is happening or if this issue is real.

I am actually trying to use precisely this idea to try to prune the number of sets in the concrete rectification. By first intersecting with the positive orthant, a lot of the set may be lost. Comparing the support vector before and after this intersection can reveal which dimensions are important for the rectification.

[mforets]

I am confirming whether this is what is happening or if this issue is real.

sorry, could you reformulate what is the issue?

should be the same:
ρ(ones(3), hpp), ρ(ones(3), hp) # (0.7047530240332701, 0.9896255059528571)

i disagree; but it should hold that ρ(ones(3), hpp) <= ρ(ones(3), hp)

[tomerarnon]

sorry, could you reformulate what is the issue?

If there was really a problem with intersection.

You are both right, however, that I am wrong! I must have expected to see the entire positive projection covered because I was viewing many plots of rectifications beforehand (and for a couple of days...), and that let me down a misleading rabbit hole. There is no issue as far as I can tell, including after checking with MeshCat in 3D.

[mforets]

alright! feel free to keep bringing issues even under doubt :) it's very useful for the library.