@JohnZ, to extend on my comment:
There is a classic method for approximating polyhedra [1], useful for applications that require operating with high-dimensional subsets of \mathbb{R}^n [1] (“high” meaning n > 10, i.e. when options 2 and 3 above are not a viable alternative). This is the kind of approach taken in our library LazySets.jl.
[1] Kamenev, G. K. (1996). An algorithm for approximating polyhedra. Computational Mathematics and Mathematical Physics, 36(4), 533–544,
Let me comment in passing that LazySets.jl does heavily rely on Polyhedra.jl: whenever your computation requires to work with H-rep and V-rep interchangeably, quantifier elimination and so on. The tools interface with each other nicely. On the other hand, there are many other things that can be done “on-the-fly” very efficiently! This is where the mathematical methods implemented in LazySets.jl shine.
The following code implements the example from the OP. Note that m had 12 times the same constraint, so I only kept one as the half-space C. I also changed the sign in front of Data2; otherwise plots were too boring (the result was always a square). This result shows that the feasible set is certainly contained in the triangle. Moreover, this triangle is tight by construction. But we can say more than that: the result is actually exact, because the underapproximation, coincides with the overapproximation for eps small enough (you can also specify the directions in the overapproximation too, but Kamenev’s method picks them automatically).
Overapproximation in 3D is also possible, though it seems like a good idea that you first get familiar with the 2D code. See the documentation of LazySets.jl for details or ask back if you need more help!
About runtime, X̂(0.01) takes < 1ms and the call to underapproximate is 3 times as fast. This approach scales well with the dimension, as it is internally just solving LPs where appropriate.
using LazySets, Plots
Data1 = rand(1.0:20,12)
Data2 = rand(1.0:15,12)
B = Hyperrectangle(low=zeros(24), high=ones(24))
C = HalfSpace(vcat(Data1, Data2), 0.3)
P = intersection(B, C)
X = Projection(P, [23, 24]) # lazy projection
ρ(vcat(Data1, zeros(12)), P) # compute optimum value max sum(Data1[i]*x[i])
X̂(ε) = overapproximate(X, HPolygon, ε); # epsilon-close approximation
Xunder = plot(underapproximate(X, BoxDirections(2)))
pover = plot()
plot!(pover, X̂(1.0), alpha=.5, title="Overapproximation", lab="ε = 1.0", lw=3.0)
plot!(pover, X̂(0.01), alpha=.5, lab="ε = 0.01")
punder = plot(Xunder, alpha=0.1, title="Underapproximation", lab="boxdirs")
plot(pover, punder, layout=(1, 2))
