very slow sampling from a 3-dimensional HPolytope (Julia 1.5.2)

[ko56]

I have the following 3-dimensional polytope:

julia> GXent.GPolyTope HPolytope{Float64,Array{Float64,1}}(HalfSpace{Float64,Array{Float64,1}}[HalfSpace{Float64,Array{Float64,1}}([1.0, 1.0, 0.0], 4.0), HalfSpace{Float64,Array{Float64,1}}([-1.0, -1.0, -0.0], -4.0), HalfSpace{Float64,Array{Float64,1}}([-1.0, 0.0, 0.0], -0.0), HalfSpace{Float64,Array{Float64,1}}([0.0, -1.0, 0.0], -0.0), HalfSpace{Float64,Array{Float64,1}}([0.0, 0.0, -1.0], -0.0), HalfSpace{Float64,Array{Float64,1}}([1.0, 0.0, 0.0], 10.0), HalfSpace{Float64,Array{Float64,1}}([0.0, 1.0, 0.0], 10.0), HalfSpace{Float64,Array{Float64,1}}([0.0, 0.0, 1.0], 10.0), HalfSpace{Float64,Array{Float64,1}}([0.0, 1.0, 1.0], 6.0)])

On a reasonably fast machine I see

julia> @time sample(GXent.GPolyTope) 5.632507 seconds (30.91 M allocations: 3.224 GiB, 7.10% gc time) 3-element Array{Float64,1}: 0.48320334527456144 3.5167966552486147 1.6576822969795617

julia>

For comparison, sampling from a random 3-dimensional polytope is almost instantaneous. Any help would be appreciated.

[mforets]

Hi @ko56,

I think the underlying issue here is that the set is flat in some dimensions. I don't know if you are familiar with the fact that the default sampling, sample(X, npoints), uses a method called rejection sampling. It is for that reason that random sampling from a flat set may take long time.

If you work only in three dimensions, it shouldn't be too hard to have a specialized method that checks if the set is flat on 2d projections and use a specialized method.

[ko56]

Hi Marcelo,

Thanks for replying.  By "flat" do you mean that it has the x-, y-, and z-planes as facets (first orthant)? Yes, I know that the default sampler does rejection sampling, using a bounding hyperrectangle. But if what you mean by "flat" is what I think, I still don't see why rejection sampling is slow.

I'm not only working in 3 dimensions, I chose 3 because I could produce a small example that  shows the problem.  But most of the polytopes I'm interested in have the constraints x >= 0, y >= 0, etc.

Kostas

On Mon, 2021-02-01 at 09:50 -0800, Marcelo Forets wrote:

Hi @ko56, I think the underlying issue here is that the set is flat in some dimensions. I don't know if you are familiar with the fact that the default sampling, sample(X, npoints), uses a method called rejection sampling. It is for that reason that random sampling from a flat set may take long time. If you work only in three dimensions, it shouldn't be too hard to have a specialized method that checks if the set is flat on 2d projections and use a specialized method. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[mforets]

Flat in the sense that the projection onto the 1-2 plane the set is a line segment. I think that's the reason why it takes time for the rejection sampling because it has to try many many times until it hits one point in the set.

We've opened a pull request that improves random sampling from line segments, https://github.com/JuliaReach/LazySets.jl/pull/2538 I think that we can combine such feature to sample more quickly from the set in your question, I'll upload a notebook shortly.

[mforets]

Flat in the sense that the projection onto the 1-2 plane the set is a line segment.

you can see this by doing

julia> plot(Projection(X, [1, 2]), 1e-2)

[mforets]

About my previous comment, please see https://nbviewer.jupyter.org/github/mforets/escritoire/blob/master/2021/Week5/SamplingFlat.ipynb which samples ~40 points from the set in ~1ms. the solution is indeed very particular to the problem and it's not useful in general -- but maybe you can get some inspiration! ;)

again i think choosing an appropriate algorithm depends on the number of dimensions (basically 2, 3, or more than 3). in LazySets there is only rejection sampling implemented, but there exist other methods in the literature such as those based on random walks on the vertices.

[mforets]

if you want to try that notebook you have to checkout the branch like so

] add LazySets#mforets/2537

[ko56]

Ok, I see what you mean.  The example polytope I provided is generated from a JuMP model, which involves the constraints

x1, x2, x3 >= 0,  x1 + x2 = 4, x2 + x3 <= 6

So maybe the general issue is sampling from a polytope that is not full- dimensional?

On Mon, 2021-02-01 at 11:09 -0800, Marcelo Forets wrote:

Flat in the sense that the projection onto the 1-2 plane the set is a line segment. you can see this by doing julia> plot(Projection(X, [1, 2]), 1e-2) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[mforets]

Exactly, i think that rejection sampling (at least in the naive/straightforward way which is currently implemented) is not a good choice for non full dimensional polytopes.

[ko56]

Ok, this situation is not uncommon when dealing with optimization models.   So let me look into reducing the polytope using some of the functions in Polyhedra. It seems to me this may be a better approach, rather than investing effort into improving rejection sampling.  What do you think?

(And also, thanks for the notebook. Very informative.)

On Mon, 2021-02-01 at 15:12 -0800, Marcelo Forets wrote:

Exactly, i think that rejection sampling (ar least in the naive/straightforward way which is currently implemented) is not a good choice for non full dimensional polytopes. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[mforets]

About sampling strategies -- I'm not an expert -- but there is a class of methods called "hit and run", see eg. the second answer here. In particular the method implemented in this Python package seems interesting! Unfortunately there is no license so we would have to explicitly ask the author(s) for use.

It seems to me this may be a better approach, rather than investing effort into improving rejection sampling. What do you think?

I don't know, it may depend in the "big picture" of your problem -- if the computation of vertices in X is an intermediate step there may be other workarounds; if it is the "final" step then finding a more robust way to sample seems relevant (or to modify the set such that it is full dimensional and use rejection sampling, that's what you propose in your last comment IIUC).

[ko56]

The "big picture" in what I'm doing is this: you solve a concave maximization problem with JuMP. The domain of the problem is a convex set C (polytope here). Now you want to test a conjecture on the set consisting of C minus a norm ball of a certain radius centered at the solution of the maximization problem. Let this set (non-convex) set be C'. The conjecture is that a certain inequality holds at all points of the set C'. You are working toward proving this conjecture, but you certainly want to see if there's any evidence that it's false.

I also am no expert on sampling. But at least in the setting I described, it seems to me that modifying the polytope is a more general approach.

[mforets]

Ok, I see. When I have some free time again I'll check whether the following idea works: (1) compute the relative interior of X, call it Xint; (2) compute an inner approximation of Xint; (3) sample from such set, eg. by rejection sampling.

[ko56]

This sounds good. The only problem I see, for my situation, is that there is some possibility that the inequality I'm trying to check may be true in the interior, but may fail on the boundary of the set.

[schillic]

Another option to sample from a polytope P is to observe that x ∈ P iff x is a convex combination of the vertices ∑ αᵢ vᵢ. So you would just sample random factors αᵢ. This is of course costly for sets with many vertices but will always work independent of the shape.

[mforets]

Unfortunately there is no license so we would have to explicitly ask the author(s) for use.

I have asked the author for permission and he replied that he would be fine that we use / translate to Julia his package.

[mforets]

sorry if it is a bit off-topic here (feel free to open new issues), about:

The conjecture is that a certain inequality holds at all points of the set C'. You are working toward proving this conjecture, but you certainly want to see if there's any evidence that it's false.

how do you compute the set difference? did you consider working only in the constraint representation to prove the conjecture, by intersecting with the inequality? below is a MWE of what i mean assuming that the inequality is linear (if it is not, you can consider converting it to one or more polytopes using ICP):

julia> using LazySets
[ Info: Precompiling LazySets [b4f0291d-fe17-52bc-9479-3d1a343d9043]

julia> C = rand(HPolygon);

julia> B(ε) = BallInf(zeros(2), ε);

julia> C′ = pontryagin_difference(C, B(0.1))
HPolytope{Float64,Array{Float64,1}}(HalfSpace{Float64,Array{Float64,1}}[HalfSpace{Float64,Array{Float64,1}}([0.9485726302597515, 0.10802315532423185], 0.3933841651889365), HalfSpace{Float64,Array{Float64,1}}([0.2414270181292819, 0.30071680177214843], 1.1867559370604959), HalfSpace{Float64,Array{Float64,1}}([0.6027322972581786, 0.9141764803480165], 3.6521947138374444), HalfSpace{Float64,Array{Float64,1}}([0.14182145934518875, 0.5654023544368563], 2.5482308687918622), HalfSpace{Float64,Array{Float64,1}}([-1.3255052032041248, 1.265353716335288], 8.425362910339212), HalfSpace{Float64,Array{Float64,1}}([-2.1099953382912515, -2.6733648546904627], -4.17431007148387), HalfSpace{Float64,Array{Float64,1}}([1.5009471365029754, -0.4803076535260784], -1.445878316520327)])

julia> μ = HalfSpace([1.0, 1.0], 0.0) # certain inequality
HalfSpace{Float64,Array{Float64,1}}([1.0, 1.0], 0.0)

julia> vcat(C′.constraints, μ) |> HPolytope |> isempty
true

[ko56]

On Wed, 2021-02-03 at 02:51 -0800, Marcelo Forets wrote:

sorry if it is a bit off-topic here (feel free to open new issues), about: The conjecture is that a certain inequality holds at all points of the set C'. You are working toward proving this conjecture, but you certainly want to see if there's any evidence that it's false. how do you compute the set difference?  Sorry, I was imprecise: when I said "minus a ball" I meant ordinary set difference, not difference in the Minkowski sense. So I just check if the point generated in the polytope is also in the ball. did you consider working only in the constraint representation to prove the conjecture, by intersecting with the inequality. below is a MWE of what i mean assuming that the inequality is linear (if it is not, you can consider converting it to one or more polytopes using ICP): julia> using LazySets [ Info: Precompiling LazySets [b4f0291d-fe17-52bc-9479-3d1a343d9043]

julia> C = rand(HPolygon);

julia> B(ε) = BallInf(zeros(2), ε);

julia> C′ = pontryagin_difference(C, B(0.1)) HPolytope{Float64,Array{Float64,1}}(HalfSpace{Float64,Array{Float64,1}}[HalfSp ace{Float64,Array{Float64,1}}([0.9485726302597515, 0.10802315532423185], 0.3933841651889365), HalfSpace{Float64,Array{Float64,1}}([0.2414270181292819, 0.30071680177214843], 1.1867559370604959), HalfSpace{Float64,Array{Float64,1}}([0.6027322972581786, 0.9141764803480165], 3.6521947138374444), HalfSpace{Float64,Array{Float64,1}}([0.14182145934518875, 0.5654023544368563], 2.5482308687918622), HalfSpace{Float64,Array{Float64,1}}([-1.3255052032041248, 1.265353716335288], 8.425362910339212), HalfSpace{Float64,Array{Float64,1}}([-2.1099953382912515, -2.6733648546904627], -4.17431007148387), HalfSpace{Float64,Array{Float64,1}}([1.5009471365029754, -0.4803076535260784], -1.445878316520327)])

julia> μ = HalfSpace([1.0, 1.0], 0.0) # certain inequality HalfSpace{Float64,Array{Float64,1}}([1.0, 1.0], 0.0)

julia> vcat(C′.constraints, μ) |> HPolytope |> isempty true — Nice example, thanks.  The inequality I'm working with involves the function maximized by JuMP, so it's non-linear. I didn't know about ICP though, seems interesting.   But in the end, improved random sampling, in one of the ways we discussed, would be a simpler approach. You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[mforets]

OK. I'll plan to implement the basic hit-and-run method soon-ish (it is now https://github.com/JuliaReach/LazySets.jl/issues/2543)

[schillic]

I implemented my idea from https://github.com/JuliaReach/LazySets.jl/issues/2536#issuecomment-771679414 in #2547. I found out that sampling like that is not uniform: you get points closer to the center. But at least you get them fast. If you want to try out the branch, you can sample with sample(X, m, sampler=LazySets.PolytopeSampler) where X is the set and m is the number of samples.

julia> @time V = sample(X, 1000, sampler=LazySets.PolytopeSampler);
  0.005130 seconds (37.68 k allocations: 1.727 MiB)

julia> plot(LazySets.project(X, 1:2))
julia> scatter!([v[1] for v in V], [v[2] for v in V])

[ko56]

It's definitely a nice idea.  But the sampling closer to the center is a problem in my case, since one of the  primary things I want to check what happens close to the boundary.

On Thu, 2021-02-04 at 12:41 -0800, Christian Schilling wrote:

I implemented my idea from #2536 (comment) in #2547. I found out that sampling like that is not uniform: you get points closer to the center. But at least you get them fast. If you want to try out the branch, you can sample with sample(X, m, sampler=LazySets.PolytopeSampler) where X is the set and m is the number of samples. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[mforets]

the hit and run algorithm also has an issue with sampling from a flat polytope X, because throwing a random direction d of the unit sphere does not necessarily hit any point in X (it always hits a hyperplane if the set is non flat). but i think we can go around this issue by solving an LP: take d at random but compute the support vector of X and build a new d' by taking the normalized difference d' - d; in this case it is guaranteed that it is contained in X.

[ko56]

I don't know enough about hit-and-run to comment intelligently, but I'm wondering if we are overcomplicating the issue. If we want to "sample uniformly" from a polytope, it seems to me that the only sensible definition is to say that we sample uniformly from its relative interior.   E.g. to sample "uniformly in R^2" from x1 >=0 , x2 >= 0, x1+x2=5, we sample from  {x1 | 0 <= x1 <= 5} in R^1.  If so, the problem boils down to reducing the polytope to its full-dimensional relative interior.  Thoughts?

On Mon, 2021-02-08 at 06:10 -0800, Marcelo Forets wrote:

the hit and run algorithm also has an issue with sampling from a flat polytope X, because throwing a random direction d of the unit sphere does not necessarily hit any point in X (it always hits a hyperplane if the set is non flat). but i think we can go around this issue by solving an LP: take d at random but compute the support vector of X and build a new d' by taking the normalized difference d' - d; in this case it is guaranteed that it is contained in X. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[mforets]

yes, i think that's true! but there is no functionality ready to use that will make such reduction for a general polytope.

if we know which variables are redundant we can apply variable elimination (see example use here), but blindly applying elimination to all possible subsets of variables seems costly.

[mforets]

E.g. to sample "uniformly in R^2" from x1 >=0 , x2 >= 0, x1+x2=5, we sample from {x1 | 0 <= x1 <= 5} in R^1

it shouldn't be too hard to apply Polyhedra.eliminate to obtain this transformation. that can be a good first step. but more generally we need a method to specify which variables to eliminate. linear algebra to the rescue :book:

[mforets]

seems like we're looking for https://faculty.math.illinois.edu/Macaulay2/doc/Macaulay2-1.16/share/doc/Macaulay2/Polyhedra/html/_affine__Hull.html

i'll ask to Polyhedra devs if this functionality is available (i searched but didn't find any)

[ko56]

Exactly, I was thinking of Polyhedra also. There is a function which computes the dimension of the polytope, (dimension of its affine hull).  Wouldn't that be helpful?

Also, the result would depend on which variables one eliminates, but I'm not clear on how much that matters.

On Tue, 2021-02-09 at 06:59 -0800, Marcelo Forets wrote:

E.g. to sample "uniformly in R^2" from x1 >=0 , x2 >= 0, x1+x2=5, we sample from {x1 | 0 <= x1 <= 5} in R^1 it shouldn't be too hard to apply Polyhedra.eliminate to obtain this transformation. that can be a good first step. but more generally we need a method to specify which variables to eliminate. linear algebra to the rescue 📖 — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[schillic]

I implemented a new idea in #2563 because you mentioned that you want a bias toward the edges. I am not sure if there is a formal argument but my intuition tells me that this approach satisfies that goal, and the plot seems to support that (green is the new sampling).

EDIT On second thought (and also from the plot), I think this approach is biased toward the side with more corners. So again not quite what you want :disappointed:

[schillic]

I updated the approach with a seemingly better distribution (although still higher frequency in the corner areas).

[ko56]

How do I try it?

On Tue, 2021-02-09 at 12:25 -0800, Christian Schilling wrote:

I updated the approach with a seemingly better distribution (although still higher frequency in the corner areas).

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[mforets]

you can try it by checking out the branch

julia> ] add LazySets#schillic/polytope_sampling

then V2 = sample(X, 1000, sampler=LazySets.PolytopeSampler2) as in this comment

[mforets]

seems like there are now a bunch of initiatives on the table. here is a detailed summary:

• (internal change) refactor samples.jl file such that it is easier to add algorithm options (cf. this comment)
• make PolytopeSampler accept different options (#2563)
• finish HitAndRun (#2548) using an_element as starting point << the problem is that it generally returns a vertex, not an inner point
• implement minover (minimum overlap) algorithm to find an inner point of a polytope in hrep
• make HitAndRun to use the point found with minover
• extend HitAndRun with an option that uses support vector projections (i expect this to work with flat polytopes)
• method to sample from the affine hull of a polytope << hjow to do this systematically?

[schillic]

@ko56 There is a new release of LazySets with a new default sampler that automatically falls back to the sampler for (potentially flat) polytopes if rejection sampling does not make progress. So I consider this issue closed.

julia> using LazySets
julia> import Polyhedra

P = HPolytope([HalfSpace([1.0, 1.0, 0.0], 4.0), HalfSpace([-1.0, -1.0, -0.0], -4.0), HalfSpace([-1.0, 0.0, 0.0], -0.0),
 HalfSpace([0.0, -1.0, 0.0], -0.0), HalfSpace([0.0, 0.0, -1.0], -0.0), HalfSpace([1.0, 0.0, 0.0], 10.0),
 HalfSpace([0.0, 1.0, 0.0], 10.0), HalfSpace([0.0, 0.0, 1.0], 10.0), HalfSpace([0.0, 1.0, 1.0], 6.0)])

julia> @time sample(P, 100)
glp_simplex: unable to recover undefined or non-optimal solution
  0.000884 seconds (3.67 k allocations: 266.484 KiB)
100-element Vector{Vector{Float64}}:
...

The only catch is that you need to load Polyhedra. I just sent a proposal to print a better error message to indicate that in #2846.

[ko56]

Thanks, this looks good. Issue closed.