Project to nearest point on convex polytope

I have a point $$y \in \mathbb{R}^d$$ and a convex polytope $$\mathcal{P}$$ given as the intersection of half-spaces:

$$\mathcal{P} = \{x \in \mathbb{R}^d \mid a_1 \cdot x \le b_1, \dots, a_n \cdot x \le b_n\}.$$

I would like to project $$y$$ onto the polytope, i.e., to find the nearest point $$z \in \mathcal{P}$$: in other words, to minimize $$\|y-z\|_2$$ subject to $$z \in \mathcal{P}$$. I know there are algorithms using quadratic programming, but I am hoping for a simple to implement method, even if it is not optimal.

Here is one possible incremental method: pick the halfspace that $$y$$ is furthest from, i.e., find the index $$i$$ that maximizes $$a_i \cdot y - b_i$$, then project $$y$$ onto that halfspace, i.e., replace $$y$$ with $$y' = y - (a_i \cdot y - b_i) a_i$$, and repeat. (I have assumed, without loss of generality, the inequalities have been normalized so $$\|a_i\|_2=1$$.) While this might not yield the optimal solution, I hope that after it a fixed number of iterations it will get close to the optimal solution.

Is this a good method? Is there a better method that is simple to implement and does reasonably well?

• Have you implemented and tested this? It sounds reasonable but sometimes playing around a little bit helps to get insight. Jul 22 '20 at 3:27
• Also potentially answered here: stackoverflow.com/questions/2433298/… . Gilbert–Johnson–Keerthi distance algorithm might be what you are looking for. Jul 22 '20 at 3:57
• @hardmath, good call, thank you! I've updated the question.
– D.W.
Dec 15 '21 at 3:17

Assume that $$y$$ is not in the polyhedron (it is easy to check whether it is, and we know that the distance is zero in that case). If $$y$$ is outside then the closest point will be on the surface of the polyhedron.

So I came up with the following (horrible) algorithm, which will give you an upper bound. Let $$y^0=y$$.

1. Find distance of point $$y^n$$ to all planes $$a_i\cdot x = b_i$$.
2. Pick the closest one, save the distance, project $$y$$ on the plane and get $$y^{n+1}$$.
3. Remove the plane from the list for the purposes of steps 1 and 2.
4. Repeat steps 1-3 with until $$y^{n+1}$$ is inside the polyhedron.

The sum of the distances you saved during the process will be the upper bound. Since it is a convex polyhedron, this algorithm should terminate in at most 5 iteration. I am not so sure about this last claim, so I am going to remove it.

You can also potentially compute the distance between $$y$$ and $$y^{n+1}$$ to get a better upper bound.

• Thank you for suggesting ideas. How well does this work? For instance, does it always terminate? Does it give a good approximation? If we stop it early, does it give a good approximation? The algorithm I listed in the question finds the farthest plane and projects to that, whereas your answer proposes to find the nearest plane, but is otherwise very similar. What's the rationale for picking the nearest plane rather than the farthest plane? If we stop early after a fixed number of iterations, I suspect what I proposed will be better, but I haven't tried to check that.
– D.W.
Jul 22 '20 at 5:56
• For most cases, this will terminate faster than your suggestion. Both our suggestions should terminate, at worst, in (number of faces of polyhedron) iterations. However, the thing I cannot prove, with either your suggestion or mine, is whether $y'$ will be on or in the polyhedron when the algorithm stops. I could not create a counter-example so I am going to hypothesize that it will be, but I don't know. Good approximation, maybe not but can't say, with either algorithm. Mine will (at least, should) give an upper bound always, yours might find a lower bound sometimes. That is my rationale. Jul 22 '20 at 16:58