Convex Hulls

In this section we illustrate the convex hull operation. We give examples of the symbolic implementation, and the concrete convex hull in low dimensions. We show how to test if a point lies in the convex hull of a set of points in the plane using LazySets. Moreover, we give examples of creating the convex hull of sets whose vertices are represented as static vectors, which can dramatically improve performance in many use cases. Finally, we give an example of creating the convex hull of points in higher dimensions.

Convex Hulls

Symbolic convex hull

The lazy convex hull, ConvexHull, is the binary operator that implements the convex hull of the union between two convex sets.

using Plots, LazySets

A = 1/sqrt(2.) * [1 -1; 1 1]
Bn = n -> BallInf(ones(n), 0.2)

X = Bn(2)
Y = CH(X, exp(A) * X)

ConvexHull{Float64, BallInf{Float64, Vector{Float64}}, LinearMap{Float64, BallInf{Float64, Vector{Float64}}, Float64, Matrix{Float64}}}(BallInf{Float64, Vector{Float64}}([1.0, 1.0], 0.2), LinearMap{Float64, BallInf{Float64, Vector{Float64}}, Float64, Matrix{Float64}}([1.5418634570456318 -1.3175384087798812; 1.3175384087798812 1.5418634570456318], BallInf{Float64, Vector{Float64}}([1.0, 1.0], 0.2)))

The name CH is an alias for ConvexHull, so you can use both interchangeably. This type is parametric in the operands's types.

p = plot(X, 1e-2, color="blue")
plot!(p, exp(A) * X, color="green")
plot!(p, Y, color="red", alpha=0.2)

We can as well work with a 100-dimensional set:

using SparseArrays

X = Bn(100)
A = blockdiag([sparse(A) for i in 1:50]...)
Y = CH(X, exp(Matrix(A)) * X)

dim(Y)

To take the convex hull of a large number of sets, there is the n-ary type ConvexHullArray. For instance, below we create a collection of balls b via list comprehension, and pass them to create a new ConvexHullArray instance.

b = [Ball2([2*pi*i/100, sin(2*pi*i/100)], 0.05) for i in 1:100];
c = ConvexHullArray(b);

plot(c, alpha=0.1, color="blue")
plot!(b, alpha=0.5, color="red")

2D convex hull

In two dimensions the convex_hull function computes the concrete convex hull of a set of points.

points = N -> [randn(2) for i in 1:N]
v = points(30)
hull = convex_hull(v)
typeof(hull), length(v), length(hull)

(Vector{Vector{Float64}}, 30, 7)

Notice that the output is a vector of floating point numbers representing the coordinates of the points, and that the number of points in the convex hull has decreased.

We can plot both the random points and the polygon generated by the convex hull with the plot function:

p = plot([Singleton(vi) for vi in v])
plot!(p, VPolygon(hull), alpha=0.2)

Test point in convex hull

One can check whether a point lies inside or outside of a convex hull efficiently in two dimensions, using the fact that the output of convex_hull returns the points ordered in counter-clockwise fashion.

Note

To check if a point p::AbstractVector is in another set, e.g. a polygon in vertex representation V, use p ∈ V. However, if you are working with a Singleton, which is a set with one element, use set inclusion ⊆. The following example illustrates this difference.

julia> Singleton(v[1]) ∈ VPolygon(hull)

ERROR: cannot make a point-in-set check if the left-hand side is a set; either
check for set inclusion, as in `S ⊆ X`, or check for membership, as in
`element(S) ∈ X` (the results are equivalent but the implementations may differ)

As the error suggests, either use element to access the element of the singleton and check if it belongs to the right-hand side set:

element(Singleton(v[1])) ∈ VPolygon(hull)

true

Or use set inclusion between the singleton and the right-hand side set:

Singleton(v[1]) ⊆ VPolygon(hull)

true

Let us note that one can also make the point-in-convex-hull test by solving a feasibility problem; actually, this is the fallback implementation used for in any dimension. However, the specialized approach in 2D is more efficient.

Using static vectors

Usual vectors are such that you can push! and pop! without changing its type: the size is not a static property. Vectors of fixed size, among other types, are provided by the StaticArrays.jl package from the JuliaArrays ecosystem. Using static arrays for vectors of "small" dimension can dramatically improve performance.

Since the convex hull algorithm supports any AbstractVector, it can be applied with static vectors. The following example illustrates this fact.

v = points(1000)
convex_hull(points(3)) # warm-up

@time convex_hull(v)

11-element Vector{Vector{Float64}}:
 [-4.531752441957591, -0.12329522596832723]
 [-2.6406859092978037, -2.141028678495725]
 [-1.7642766649415607, -2.4218501941126793]
 [0.19674014665189077, -2.7951321956504755]
 [2.0432124300785195, -2.4869755401551554]
 [2.2636646804353004, -2.2811875779114734]
 [3.2444608111681617, -0.943070518702561]
 [2.7671219693996707, 2.0677326533026044]
 [0.30138255899627026, 2.574394501803201]
 [-0.43266703686873814, 2.7229028130426465]
 [-1.3776662919062785, 2.655904838537312]

Now working with static vectors:

using StaticArrays

convex_hull([@SVector(rand(2)) for i in 1:3]) # warm-up

v_static = [SVector{2, Float64}(vi) for vi in v]
@time convex_hull(v_static)

11-element Vector{StaticArraysCore.SVector{2, Float64}}:
 [-4.531752441957591, -0.12329522596832723]
 [-2.6406859092978037, -2.141028678495725]
 [-1.7642766649415607, -2.4218501941126793]
 [0.19674014665189077, -2.7951321956504755]
 [2.0432124300785195, -2.4869755401551554]
 [2.2636646804353004, -2.2811875779114734]
 [3.2444608111681617, -0.943070518702561]
 [2.7671219693996707, 2.0677326533026044]
 [0.30138255899627026, 2.574394501803201]
 [-0.43266703686873814, 2.7229028130426465]
 [-1.3776662919062785, 2.655904838537312]

Higher-dimensional convex hull

One can compute the convex hull of points in higher dimensions using convex_hull. The appropriate algorithm is decided based on the dimensionality of the given points.

using Polyhedra

v = [randn(3) for _ in 1:30]
hull = convex_hull(v)
typeof(hull), length(v), length(hull)

(Vector{Vector{Float64}}, 30, 15)

Here, convex_hull is now using the concrete polyhedra library Polyhedra, hence it needs to be loaded beforehand.

One can check whether a point belongs to the convex hull using ∈ as follows:

P = VPolytope(hull)
x = sum(hull)/length(hull)

x ∈ P

true

Here x ∈ P solves a feasibility problem; see the docs of ?∈ for details. Equivalently, using set inclusion:

Singleton(x) ⊆ P

true

If no additional arguments are passed, convex_hull uses the default polyhedra library from default_polyhedra_backend for the given input; different options can be passed through the backend keyword; see the Julia polyhedra website for all the available backends.