Projection of a point onto a simple convex polytope

Question

I have a convex polytope described by the intersection of a hypercube and two parallel lower and upper hyperplanes cutting through the hypercube. I want to find the projection of a point onto the polytope.

I read the question Projection of a point onto a convex polyhedra and see that the general problem of finding a projection onto a convex polytope can be solved with quadratic programming.

However, taken alone, projection onto a hypercube or a hyperplane is trivial. For the hypercube, each dimension is defined by a lower+upper bounds and we just pull each dimension within those bounds, which take $O(n)$ time in $n$ dimensions. For the hyperplane just find the projection onto it.

So I'm wondering whether there is simpler method for finding such a projection, maybe somehow combining the two independent projections.

Answer 1

You may want to check out Dykstra's projection algorithm (not to be confused with Dijkstra's algorithm). It precisely does what you want: It is an iterative method to compute the projection onto an intersection of convex sets using the projections onto these sets.

Note, that simply projecting onto the sets alone does not work. Consider the hypercube $[-1,1]^2$ and the hyperplane $\{x: x_1+x_2=1\}$. The projection of $(2,1)$ onto the intersection is $(1,0)$. If one first projects onto the cube, then onto the plane yields $(1/2,1/2)$, which is not the wanted projection.

Answer 2

From your description, I believe the optimization problem describing your projection is:

$$ \begin{aligned} \min & \quad \frac{1}{2} ||x - y||_2^2 \\ \text{s.t.} & \quad \ell_i \leq x_i \leq u_i &\quad i = 1, \dots, n \\ & \quad \alpha \leq c^T x \leq \beta \end{aligned} $$

I believe that the best way would be a simple custom interior point method. The objective augmented with log-barriers for the constraints is:

$$ A(x; \mu) = \frac{1}{2} ||x - y||_2^2 - \mu \left[ \sum_{i=1}^n \ln(u_i - x_i) + \sum_{i=1}^n \ln(x_i - \ell_i) + \ln(c^T x - \alpha) + \ln(\beta - c^T x) \right] $$

Fortunately, the Hessian has a nice structure of a diagonal matrix + two rank-one matrices, and can be easily inverted in linear time by applying the well-known Sherman–Morrison formula twice. So you get very quick convergence, and each iteration takes only linear time.