Recursive Projection

Description of the algorithm

Recursive projection is an algorithm designed to compute the projection of a convex set. In general, it computes the approximation of a smooth set \(P\). To do so, it generates converging polyhedral approximations of \(P\), \(P_{inner}\) and \(P_{outer}\). The algorithm only needs an “oracle”, i.e. an algorithm or optimization problem that yields the extremal point of \(P\) in a given direction \(d\). As a generic optimization problem:

\[\begin{split}max. \qquad &d^T x\\ s.t. \qquad &x \in P\end{split}\]

The solution \(x^*\) of this problem is an extremal point in the direction \(d\). By solving repeatedly the above problem, we obtain:

  • The convex hull of all \(x^*\) is \(P_{inner}\).
  • The intersection of all halfspaces \(\{x \in \mathbb{R}^n | d^T x \leq d^T x^*\}\) forms \(P_{outer}\)

Now, the important point is how to choose the approriate sequence of directions \(d\). To do so, we compute the uncertainty volumes, i.e. the cuts of \(P_{outer}\) by the supporting hyperplanes of \(P_{inner}\). The direction \(d\) is chosen to be perpendicular to the supporting hyperplane that forms the largest uncertainty volume.

This is very useful to compute projections. Consider a convex body \(P\) in \(\mathbb{R}^{n+m}\). Computing the projection of this body onto \(\mathbb{R}^n\) can be done by specifying the following optimization problem:

\[\begin{split}max \qquad &d^T x\\ s.t. \qquad &(x, y) \in P\end{split}\]

This is particularly interesting when \(m >> n\). In this case, computing the direct projection (see for example this page) is prohibitively expensive as the complexity is exponential in \(m+n\). In our case, the complexity depends on the class of optimization problem being solved, but is typically polynomial in the dimension.

For more information, please refer to this paper.

How to use this class

This class is intended for developpers and researchers who wish to implement a new class of problems. If you are looking to compute stability or robust stability polygons and polyhedrons, please use stabilipy.StabilityPolygon. If you wish to compute the projection of a set of linear inequalities, please use stabilipy.LinearProjection. In general, one only needs to override the stabilipy.RecursiveProjectionProblem.solve() method.

Let us have a look at an example (available in sphere.py):

import stabilipy as stab


class SphereProjection(stab.RecursiveProjectionProblem):
  """Try to approximate a sphere of radius r"""

  def __init__(self, radius):
    """:param radius: Radius of the sphere
       :type radius: double"""
    stab.RecursiveProjectionProblem.__init__(self, dimension=3)
    self.radius = radius

  def solve(self, d):
    """We are computing a sphere so the extremal point in direction d is just r*d"""
    return self.radius*d

if __name__ == '__main__':
  sphere = SphereProjection(1.0)
  sphere.compute(solver='cdd', mode=stab.Mode.iteration, maxIter=50)

In this example we:

  • Implement a class that extends RecursiveProjectionProblem
  • Override solve to return the extremal point in the provided direction
  • Instanciate that object and compute the approximation using cdd as our double-description package
  • By default, this call will print the precision at each iteration and plot the result

See below for details of the API.

Class API

class stabilipy.RecursiveProjectionProblem(dimension, verbosity=<Verbosity.info: 2>)

Base class that encapsulates the recursive projection algorithm. To use it, you need to specify your problem by implementing the solve method. Then, this class will actually perform the projection.

Construct a projection problem.

Parameters:
  • dimension (int) – Dimension of the space on which you project.
  • verbosity (Verbosity) – Verbosity of the output. Default to info.
clearAlgo()

Resets internal state

compute(mode, maxIter=100, epsilon=0.0001, solver='cdd', plot_error=False, plot_init=False, plot_step=False, plot_direction=False, record_anim=False, plot_final=True, fname_polys=None)

Compute the polygon/polyhedron at a given precision or for a number of iterations or at best.

Parameters:
  • mode (stabilipy.Mode) – Stopping criterion.
  • maxIter (int) – Maximum number of iterations.
  • epsilon – Precision target.
  • solver (stabilipy.backends) – Backend that will be used.
  • plot_error – Make a running plot of the error during computation.
  • plot_init – Plot the initial state of the algorithm.
  • plot_step – Plot the state of the algorithm at each iteration.
  • plot_direction – Plot the direction found for the next iteration.
  • record_anim – Record all steps as images in files.
  • plot_final – Plot the final polygon/polyhedron.
  • fname_polys – Record successive iterations as text files.
make_problem()

This method is called upon launching the computation. Use it to build the complete problem from user-defined quantities. Does nothing by default.

outer_polyhedron()

Return the outer polyhedron as a set of vertices

plot()

Plot the current solution and polyhedrons

polyhedron()

Return the inner polyhedron as a set of vertices

save_outer(fname)

Save the outer polyhedron as a set of vertices :param fname: Filename to which the polyhedron is saved :type fname: string

save_polyhedron(fname)

Save the inner polyhedron as a set of vertices :param fname: Filename to which the polyhedron is saved :type fname: string

solve(d)

This method should return an extremal point in the given direction d, or None in case of error. You must reimplement this function to perform a computation or use one of the pre-implemented instances

Parameters:d (np.array((dim, 1))) – Search direction
class stabilipy.Mode

All polygon computations should select a mode of operation.

  • precision: will try to reach desired precision, however many iterations are required.
  • iteration: will iterate for a number of iterations and stop, regardless of accuracy.
  • best: will try to reach desired precision under the given number of iterations.