Convex polytopes defined by equality constraints return wrong results

[LongPham7]

trafficstars

In my application domain, convex polytopes are defined by not only inequalities but also equality constraints (e.g. x = y). When I run sample_points on such convex polytopes, I get unsatisfactory or even wrong results.

For example, consider the following R code where I sample from a convex polytope defined by 0 <= x <= 5 and x = y.

library("volesti")
A <- matrix(rep(0, 4*2), nrow=4, ncol=2)
A[1,1] <- 1 # x <= 5
A[2,1] <- -1 # -x <= 0
A[3,1] <- 1; A[3,2] <- -1 # x - y <= 0
A[4,1] <- -1; A[4,2] <- 1 # -x + y <= 0
b <- c(5,0,0,0)
P <- Hpolytope$new(A, b)
sample_points(P, n=1, random_walk = list("walk" = "RDHR", walk_length=10))

The code returns NaN for both x and y. This is unsatisfactory because I expect to obtain at least some valid point from the convex polytope.

Furthermore, if I replace the last line of the above code with sample_points(P, n=1, random_walk = list("walk" = "BiW", walk_length=10, starting_point=c(1,1), nburns=0, L=0.5)), I get a wrong sample where x and y are unequal.

I use volesti v1.1.2. Does anyone know how to work around this issue? Thanks!

[vissarion]

@LongPham7 thanks for opening this issue and your example. We use to support low dimensional polytopes but this is dropped in https://github.com/GeomScale/volesti/pull/194 (see the files with names "get_full_dimensional_polytope") because the method was very unstable. See also https://github.com/GeomScale/volesti/discussions/213

However, that method might work in your example.

We are currently working on a new method that can scale in high dimensions with stability. This is already implemented in the python package https://github.com/GeomScale/dingo that uses volesti as a backend for sampling.

Thus, I am marking this issue as a feature request.

[vissarion]

Similar issue (feature request) https://github.com/GeomScale/volesti/issues/246