# Projecting polytopes

Consider a polytope given in halfspace representation by:

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

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):

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,

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]\):

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:

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

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

- \(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,

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.