Projecting polytopes

Consider a polytope given in halfspace representation by:

\begin{equation*} P = \{ x \in \mathbb{R}^n \ | \ A x \leq b, C x = d \} \end{equation*}

Define the affine projection of this polytope onto a lower-dimensional space by:

\begin{equation*} \pi(P) = \{ y \in \mathbb{R}^m \ | \ \exists x \in P, y = E x + f \} \end{equation*}

We are interested in finding the halfspace representation of \(\pi(P)\).

Double description method

Projecting is easy in vertex representation: for a polytope \(P = \mathrm{conv}(V)\), the projection is simply obtained by \(\pi(P) = \{ E v + f, v \in V\}\). One simple algorithm to our problem is then:

  • Enumerate vertices of the halfspace representation of \(P\) using the double description method
  • Project the vertices to \(\pi(V) = \{ E v + f, v \in V \}\)
  • Apply the double description method again to convert the vertex representation \(\pi(V)\) into a halfspace representation

The cdd library can be used for this purpose. One limitation of this approach is that it can be numerically unstable, and becomes significantly slower as the dimension of the polytope increases.

Open source software

An implementation of this method is available in the projection submodule of the pypoman Python library.

Ray shooting methods

Ray shooting methods provide a better alternative when the size of the \(\pi(P)\) is small, and are typically used when the dimension \(m\) of the projected polytope is 2 or 3. They are output-sensitive algorithms: their complexity depends on the number of facets of \(\pi(P)\), which is reasonably bounded in 2D or 3D from McMullen's Upper Bound Theorem.

Projecting to a 2D polygon

The 2D ray shooting algorithm was introduced in robotics by Bretl and Lall to compute the static-equilibrium polygon of a legged robots in multi-contact. It was used in subsequent humanoid works to compute other criteria such as multi-contact ZMP support areas and CoM acceleration cones.

The core idea of the algorithm goes as follows: let \(\theta\) denote an angle in the 2D plane, with \(u_\theta = [\cos(\theta) \ \sin(\theta)]\) the corresponding unit vector. We solve the following linear program (LP):

\begin{equation*} \begin{array}{rl} \underset{x \in \mathbb{R}^n, y \in \mathbb{R}^m}{\text{maximize}} & u_\theta^T y \\ \text{subject to} & A x \leq b \\ & C x = d \\ & y = E x + f \end{array} \end{equation*}

This LP computes simultaneously \(x \in P\) and its projection \(y = \pi(x)\) via its inequality and inequality constraints. With the choice of \(u_\theta\) in its objective function, the output \(y^*\) will be the point of \(\pi(P)\) furthest away in the direction of \(u_\theta\). This has two important consequences:

  • \(y^*\) is a vertex of \(\pi(P)\)
  • \(u_\theta\) is the normal vector of a supporting hyperplane of \(\pi(P)\) at \(y^*\)

A supporting hyperplane (or line, since we are in 2D) is a plane tangent to the polygon at a given point. Formally, it separates the plane in two: one side \(S_0\) that contains no point of \(\pi(P)\), and its complement that contains the polygon fully. In our present case,

\begin{equation*} S_0 = \left\{z \in \mathbb{R}^2 \ \right|\left. \ u_\theta^T (z - y^*) > 0 \right\} \end{equation*}

Supporting hyperplanes are members of the dual of a polytope. They are generated by conic combination of the halfspace representation: in \(A x \leq b\), lines of \(A\) correspond to normal vectors of supporting hyperplanes, while the corresponding entries in \(b\) correspond to their respective affine coordinates. Hence, using the output of the LP above, we have obtained one line of an H-rep of \(\pi(P)\): \(u_\theta^T y \leq u_\theta^T y^*\). By repeatedly solving LPs with different values of the angle \(\theta\), we will thus obtain a halfspace representation of an outer approximation of the projected polygon.

Meanwhile, \(y^*\) gives us one vertex of \(\pi(P)\). By repeatedly solving LPs with different angles \(\theta\), we obtain a subset \(\{ y_1^*, \ldots, y_k^*\} \subset \pi(V)\) which gives us a vertex representation for an inner approximation of the projected polygon. Bretl and Lall showed how, as the number LPs solved for different angles \(\theta\) increases, these two approximations converge toward \(\pi(P)\). Furthermore, the number of iterations \(k\) required to reach a desired precision can be computed exlicitly (by looking at the surface between the inner and outer approximations), and the algorithm is anytime, i.e., computations can be interrupted after solving any number of LPs and still return a valid solution.

Open source software

A Python implementation of this algorithm written by Quang-Cuong Pham is available in the project_polytope_bretl() function of the pypoman library.

Bretl and Lall's algorithm was also extended to projecting 3D polytopes by Audren and Kheddar, used this time to compute the 3D robust static-equilibrium polyhedron. A Python implementation is available in the StabiliPy library.

Fourier-Motzkin elimination

If the projection is done offline and once and for all, one may consider Fourier-Motzkin Elimination (FME). Contrary to the previous two approaches, this method works directly and exclusively on the halfspace representation of the input polytope. One of its benefit is that it can be applied analytically when the signs of the linear coefficients are known in advance. This feature was e.g. used to calculate the analytical formula of single-contact wrench cones.

Pivoting rule

FME relies on a pivoting rule similar to Gaussian elimination, it "eliminates" \(x\) dimensions from the combined polytope on \([x \ y]\):

\begin{equation*} P_{\times} = \left\{[x\ y] \ \right|\left. \ A x \leq b, C x = d, y = E x + f \right\} \end{equation*}

Its core idea goes as follows. First, use the equality constraints \(C x = d\) and \(y = E x + f\) to eliminate as many coordinates \(x_i\) from \(x\) as possible. This leaves us with a reduced polytope:

\begin{equation*} P' = \left\{[x_0 \ \cdots \ x_k \ y] \ \right|\left. \ A' [x_0 \ \cdots \ x_k \ y] \leq b'\right\} \end{equation*}

Now, consider the first coordinate \(x_0\). Each line of the inequality constraint can be written as:

\begin{equation*} A'_{i, 0} x_0 + \cdots + A'_{i, k} x_k + A'_{i,k+1} y_0 + A'_{i, k+2} y_1 \leq b'_i \end{equation*}

One of three outcomes can happen:

  • \(A'_{i, 0} < 0\): in this case, dividing by this number yields a lower bound \(x_0 \geq l_i\), with
\begin{equation*} l_i := \frac{b'_i}{A'_{i, 0}} - \sum_{j=1}^k \frac{A'_{i, j}}{A'_{i, 0}} x_j - \frac{A'_{i, k+1}}{A'_{i, 0}} y_0 - \frac{A'_{i, k+2}}{A'_{i, 0}} y_1 \end{equation*}
  • \(A'_{i, 0} > 0\): then, doing the same yields an upper bound \(x_0 \leq u_i\)
  • \(A'_{i, 0} = 0\): in this case, the inequality is independent from \(x_0\), i.e. it is already part of the projection where \(x_0\) has been eliminated.

By construction, \(x_0 \in [l^*, u^*]\) where \(l^*\) is the highest of all lower bounds obtained by the above process, and \(u^*\) is similarly the lowest of all upper bounds. Therefore, the system of linear inequalities has solutions if and only if this interval is nonempty, that is to say,

\begin{equation*} \forall i, j,\ l_i \leq u_j \end{equation*}

The key idea here is that the expressions of \(l_i\)'s and \(u_j\)'s only depend on \(x_1, \ldots, x_k, y_0, y_1\). To eliminate \(x_0\), we will therefore:

  • remove from \(A'\) and \(b'\) all lines \(i\) for which \(A'_{i, 0} \neq 0\)
  • add one new line for each pair (i, j) of lower and upper bounds on \(x_0\) corresponding to \(l_i - u_j \leq 0\).

The result from this operation is a new polytope over \([x_1 \ \cdots \ x_k \ y_0 \ y_1]\), from which the elimination can be repeated until the only remaining coordinates are \(y_0\) and \(y_1\).

As can be readily seen from the combinatorics of its update rule, Fourier-Motzkin elimination has a (doubly) exponential memory complexity. Imbert has proposed a set of syntactic rules to avoid considering pairs of lower-upper bounds that yield redundant inequalities, but without lowering the worst-case complexity of the overall algorithm.

Open source software

The pyfme toolbox provides a Python implementation of Fourier-Motzkin elimination, using Imbert's acceleration theorems and multiprocessing for parallel computations.

To go further

The paper Equality set projection: A new algorithm for the projection of polytopes in halfspace representation proposes an iterative projection algorithm from the ray shooting family that closely resembles Bretl & Lall's algorithm. Its Section 8 compares a number of alternative methods. A dual reformulation of this algorithm was also proposed by Del Prete et al. to compute static-equilibrium conditions for legged robots.

An extensive discussion of Fourier-Motzkin elimination is given in the reference book Lectures on polytopes by Ziegler. This algorithm was, for instance, used in robotic grasping by Ponce et al. to compute force-closure grasps of polyhedral objects.


There are no comments yet. Feel free to leave a reply using the form below.

Post a comment

You can use Markdown with $\LaTeX$ formulas in your comment.

You agree to the publication of your comment on this page under the CC BY 4.0 license.

Your email address will not be published.

© Stéphane Caron — All content on this website is under the CC BY 4.0 license.