Convex Polyhedra and Constrained Convex Hulls

Suppose we have a set of nn points

V={viRdi=1,,n}.\begin{equation*} \mathcal{V} = \{ \bm{v}_i\in\mathbb{R}^d \mid i=1,\dots,n\}. \end{equation*}

A convex combination of these points is a weighted sum x=iwivi\bm{x}=\sum_iw_i\bm{v}_i with non-negative weights wiw_i satisfying iwi=1\sum_iw_i=1. The convex hull of V\mathcal{V} is the set of all convex combinations of points in V\mathcal{V}. Gathering the points in V\mathcal{V} into the columns of the matrix V=[v1,,vn]\bm{V}=[\bm{v}_1,\dots,\bm{v}_n], we can write the convex hull as the set

C=conv(V)={Vww0,eTw=1},\begin{equation}\label{1} \mathcal{C} = \mathrm{conv}(\mathcal{V}) = \{ \bm{V}\bm{w} \mid \bm{w}\geq\bm{0}, \bm{e}^T\bm{w}=1 \}, \end{equation}

where e\bm{e} is a vector of ones, w=[w1,,wn]T\bm{w}=[w_1,\dots,w_n]^T is the vector of weights, and the inequality is taken elementwise.

In this post, we explore how to put additional constraints on the weights w\bm{w} in (1)\eqref{1} and what the resulting set looks like.

Polyhedral Double Description

Let’s begin with some background on polyhedral geometry. A polyhedron is a geometric shape with flat sides (called faces). In this post we will assume we are talking about convex polyhedra; a polyhedron is convex when the line between any two points in the polyhedron is also completely contained in it. A bounded polyhedron is called a polytope and a two-dimensional polyhedron is called a polygon, some examples of which are shown in Figure 1 below.

Examples of two polygons.

Figure 1: Examples of two (convex) polygons, each of which can be defined as the convex hull of its vertices.

By the Minkowski-Weyl theorem, any polyhedron P\mathcal{P} can be equivalently represented in two different ways, together known as the double description of P\mathcal{P}. These two representations are:

Vertex Representation

The vertex representation (V-rep) can be written as

P={Vw+Uzw,z0,eTw=1},\begin{equation*} \mathcal{P} = \{ \bm{V}\bm{w} + \bm{U}\bm{z} \mid \bm{w},\bm{z}\geq\bm{0}, \bm{e}^T\bm{w}=1 \}, \end{equation*}

where each column of V\bm{V} is a vertex and each column of U\bm{U} is a ray (together these are referred to as generators). When U=0\bm{U}=\bm{0}, P\mathcal{P} is just a convex hull like (1)\eqref{1}. When V=0\bm{V}=\bm{0}, P\mathcal{P} instead reduces to a conical hull; however, we will only consider cases with U=0\bm{U}=\bm{0} in this post.

Halfspace Representation

The halfspace representation (H-rep) can be written as

P={xRdAxb},\begin{equation} \mathcal{P} = \{ \bm{x}\in\mathbb{R}^d \mid \bm{A}\bm{x}\leq\bm{b} \},\label{2} \end{equation}

which is a finite intersection of halfspaces. Equality constraints of the form Ax=b\bm{A}\bm{x}=\bm{b}, which serve to restrict to the polyhedron to a lower-dimensional subspace of Rd,\mathbb{R}^d, can also be included by equivalently using two inequalities:

[AA]x[bb].\begin{equation*} \begin{bmatrix} \bm{A} \\ -\bm{A} \end{bmatrix}\bm{x} \leq \begin{bmatrix} \bm{b} \\ -\bm{b} \end{bmatrix}. \end{equation*}

Converting from H-rep to V-rep is called vertex enumeration, while converting from V-rep to H-rep is called facet enumeration. No provably polynomial-time algorithms are known for these conversions, but nevertheless they can be computed fairly quickly when the dimension and number of generators and faces is small, using a library like cddlib or its Python bindings pycddlib.

You can read more about polyhedral double description here and see some examples of its application in robotics in this paper as well as this more recent paper of mine.

Computing Constrained Convex Hulls

Now that we’re equipped with some background on polyhedra, you may notice that the set of weights we used to define the convex hull (1)\eqref{1} is itself the H-rep of the polyhedron

W={wRnw0,eTw=1}.\begin{equation*} \mathcal{W} = \{\bm{w}\in\mathbb{R}^n\mid \bm{w}\geq\bm{0},\bm{e}^T\bm{w}=1\}. \end{equation*}

This particular polyhedron is known as a simplex and is a special case of the general form (2)\eqref{2} with

A=[1neTeT],b=[011],\begin{align*} \bm{A} &= \begin{bmatrix} -\bm{1}_n \\ \bm{e}^T \\ -\bm{e}^T \end{bmatrix}, & \bm{b} &= \begin{bmatrix} \bm{0} \\ 1 \\ -1 \end{bmatrix}, \end{align*}

where 1n\bm{1}_n is the n×nn\times n identity matrix. To incorporate additional constraints wu\bm{\ell}\leq\bm{w}\leq\bm{u}, where 0ue\bm{0}\leq\bm{\ell}\leq\bm{u}\leq\bm{e}, eT1\bm{e}^T\bm{\ell}\leq1, and eTu1\bm{e}^T\bm{u}\geq1 (without these latter two conditions, the constraints could be inconsistent because eTw=1\bm{e}^T\bm{w}=1 is impossible to realize if eT>1\bm{e}^T\bm{\ell}>1 or eTu<1\bm{e}^T\bm{u}<1), we add them to A\bm{A} and b\bm{b} to obtain

A=[1neTeT1n1n],b=[011u].\begin{align*} \bm{A} &= \begin{bmatrix} -\bm{1}_n \\ \bm{e}^T \\ -\bm{e}^T \\ -\bm{1}_n \\ \bm{1}_n \end{bmatrix}, & \bm{b} &= \begin{bmatrix} \bm{0} \\ 1 \\ -1 \\ \bm{\ell} \\ \bm{u} \end{bmatrix}. \end{align*}

Converting this H-rep to V-rep, we get the set

W={Qyy0,eTy=1}\begin{equation*} \mathcal{W}' = \{ \bm{Q}\bm{y} \mid \bm{y}\geq\bm{0},\bm{e}^T\bm{y}=1 \} \end{equation*}

with some matrix Q\bm{Q} (note W\mathcal{W}' is no longer a simplex). Substituting w=Qy\bm{w}=\bm{Q}\bm{y} into (1)\eqref{1}, we get the constrained convex hull

C={VQyy0,eTy=1},\begin{equation*} \mathcal{C}' = \{ \bm{V}\bm{Q}\bm{y} \mid \bm{y}\geq\bm{0},\bm{e}^T\bm{y}=1 \}, \end{equation*}

which is another V-rep polyhedron with weights y\bm{y} and vertices given by the columns of VQ\bm{V}\bm{Q}.

Visualizing Constrained Convex Hulls

Now that we can compute the vertices of a constrained convex hull, we can plot it. Figure 2 below shows the constrained convex hulls for different values of =e\bm{\ell}=\ell\bm{e} and u=ue\bm{u}=u\bm{e}. You can play with this yourself using the Python code here.

Examples of constrained convex hulls.

Figure 2: Examples of constrained convex hulls. The original convex hull is shown in blue, while each of the red polygons is the convex hull with different bounds applied. From top to bottom we have bounds: (,u)=(0.05,1)(\ell,u)=(0.05, 1), (,u)=(0,0.6)(\ell,u)=(0, 0.6), and (,u)=(0.05,0.6)(\ell,u)=(0.05, 0.6).

As we can see, the lower bound \bm{\ell} has the effect of “pulling” the hull in from the faces of the original polytope, while the upper bound u\bm{u} pulls the hull away from the vertices. We only consider bounds in this post, but in general we can add any desired affine constraints to w\bm{w}.

Analytical Solutions

It turns out that when have only lower or upper bounds (not both) then we can calculate the vertices of the constrained convex hulls analytically (that is, without converting between V-rep and H-rep).

Lower Bounds

Let us first consider lower bounds \bm{\ell}. The iith vertex vi\bm{v}_i, which corresponds to weight wiiw_i\geq\ell_i, pulls the convex hull toward it from each other vertex by a factor of i\ell_i. An example with the shape from Figure 2 with a lower bound on only one of the vertex weights and no (additional) upper bounds is shown in Figure 3 below.

A constrained convex hull with a single lower bound.

Figure 3: The constrained convex hull (red) of the polygon from Figure 2 (blue) with bounds =[,0,0,0,0]T\bm{\ell}=[\ell,0,0,0,0]^T, =0.2\ell=0.2, and u=e\bm{u}=\bm{e}. The vertex corresponding to the non-zero lower bound is on the far left side. The lower bound has the effect of “pulling” each of the other vertices toward the bounded vertex, proportionally to the value of the bound. To see this, a line has been drawn between the bounded vertex and each other vertex, with a fraction proportional to \ell coloured yellow and the remainder, proportional to (1)(1-\ell), in green (entire figure shown on gray background for better contrast).

Each vertex vi\bm{v}_i' of the constrained convex hull depends on the bounds at all vertices of the original shape. To compute it, we have the formula

vi=vi+j=1nj(vjvi).\begin{equation*} \bm{v}_i' = \bm{v}_i + \sum_{j=1}^n \ell_j(\bm{v}_j - \bm{v}_i). \end{equation*}

Implemented in Python code, it looks like

# V is the array of original vertices, one per row, shape (n, d) (yes, this is
# the transpose of the the definition of V in eq. (1) - sorry)
# l is the array of lower bounds, shape (n,)

V_new = []  # vertices of the lower-bounded convex hull
for i in range(n):
    v_new = V[i, :] + sum([l[j] * (V[j, :] - V[i, :]) for j in range(n)])
    V_new.append(v_new)

Upper Bounds

The upper bounds are a bit different. Each lower bound affected all other vertices, while each upper bound only affects the current vertex (i.e., the one corresponding to the bound), and this effect is a function of only the adjacent (i.e., neighbouring) vertices. Furthermore, unlike the lower-bounded case, the upper bound often results in a constrained convex hull with a different number of vertices than the original polytope. An example with the shape from Figure 2 with only a single upper bound on the vertex weights and no (additional) lower bounds is shown in Figure 4 below.

A constrained convex hull with a single upper bound.

Figure 4: The constrained convex hull (red) of the polygon from Figure 2 (blue) with bounds =0\bm{\ell}=\bm{0} and u=[u,1,1,1,1]T\bm{u}=[u,1,1,1,1]^T with u=0.6u=0.6. The vertex corresponding to the non-unity upper bound is on the far left side. The upper bound has the effect of pulling the hull toward the neighbouring vertices, proportionally to (1u)(1-u), creating a new face (dashed line). To see this, a line has been drawn between the bounded vertex and its two neighbours, with a fraction proportional to uu coloured yellow and the remainder, proportional to (1u)(1-u), in green.

Let V\mathcal{V}' be the set of vertices of the constrained convex hull, initially empty. Let A(vi)\mathcal{A}(\bm{v}_i) be the set of indices of the vertices adjacent to vertex vi\bm{v}_i. For each vertex vi\bm{v}_i and each of its adjacent vertices vj,jA(vi)\bm{v}_j,j\in\mathcal{A}(\bm{v}_i), we compute

v=vj+ui(vivj)\begin{equation*} \bm{v}' = \bm{v}_j + u_i(\bm{v}_i - \bm{v}_j) \end{equation*}

and add it to the set of new vertices:

VV{v}.\begin{equation*} \mathcal{V}'\leftarrow\mathcal{V}'\cap\{\bm{v}'\}. \end{equation*}

Once we’ve gone through all vertices and their neighbours, we are left with the complete set of vertices V\mathcal{V}' of the upper-bounded convex hull. In Python, we have

# V is the array of original vertices, one per row, shape (n, d)
# u is the array of upper bounds, shape (n,)

# note I assume we have a function adj(V, i) that returns the indices of
# vertices adjacent to V[i, :], but this can be computed fairly efficiently (in
# my actual code I use pycddlib, which uses linear programming to compute
# adjacency)

V_new = []  # vertices of the upper-bounded convex hull
for i in range(n):
    for j in adj(V, i):
        v_new = V[j, :] + u[i] * (V[i, :] - V[j, :])
        V_new.append(v_new)

Note that this approach for computing the upper-bounded convex hull only works when ui+uj1u_i+u_j\geq1 for all i=1,,ni=1,\dots,n and jA(vi)j\in\mathcal{A}(\bm{v}_i). Otherwise, the new faces added for each upper bound can intersect with each other, which complicates the solution considerably.

Both Lower and Upper Bounds

If we have both lower and upper bounds, we can independently compute the resulting constrained convex hulls using the analytical formulas above. To combine them into a single constrained convex hull, we need to take their intersection—unfortunately, I am not aware of an analytical solution for this in the vertex representation. Instead, we must again turn to polyhedral double description.

Let V\mathcal{V}_\ell and Vu\mathcal{V}_u be the sets of vertices of the lower- and upper-bounded convex hulls, respectively. Then:

  1. Convert both of these V-rep polytopes to the H-reps {xRdAxb}\{\bm{x}\in\mathbb{R}^d\mid\bm{A}_\ell\bm{x}\leq\bm{b}_\ell\} and {xRdAuxbu}\{\bm{x}\in\mathbb{R}^d\mid\bm{A}_u\bm{x}\leq\bm{b}_u\}.
  2. Stack the H-reps to combine them into the single polytope {xRdAxb},\{\bm{x}\in\mathbb{R}^d\mid\bm{A}\bm{x}\leq\bm{b}\}, whereA=[AAu],b=[bbu].\begin{align*} \bm{A} &= \begin{bmatrix} \bm{A}_\ell \\ \bm{A}_u \end{bmatrix}, & \bm{b} &= \begin{bmatrix} \bm{b}_\ell \\ \bm{b}_u \end{bmatrix}. \end{align*}
  3. Convert the combined H-rep polytope back into V-rep if desired (e.g., for plotting).

Uniform Padding

Finally, notice that the width of the padding between the boundary of the constrained convex hull and that of the original shape produced by the lower bound in Figure 2 is not uniform. If you do want uniform padding, do the following:

  1. Convert the polytope to the H-rep {xRdAxb}\{\bm{x}\in\mathbb{R}^d\mid\bm{A}\bm{x}\leq\bm{b}\}.
  2. For each row aiT\bm{a}_i^T of A\bm{A} and corresponding element bib_i of b\bm{b}, normalize and subtract the desired padding width rr:bibi/ai2r,aiai/ai2.\begin{align*} b_i &\leftarrow b_i/\|\bm{a}_i\|_2 - r, \\ \bm{a}_i &\leftarrow\bm{a}_i/\|\bm{a}_i\|_2. \end{align*}
  3. Convert back to the V-rep for plotting.

The normalization in the second step makes each row ai\bm{a}_i a unit vector without changing the polytope defined by the inequalities. This ensures each inequality has the same scale, so that subtracting rr from each results in the same width of padding for each face. An example of uniform padding is shown in Figure 5 below.

A polygon inset with uniform padding.

Figure 5: The same polygon from Figure 2 inset with uniform padding of width r=50r=50 pixels.

Thanks to Abhishek Goudar for reviewing a draft of this post.