Given a set of vectors $V = \{ \mathbf{v}_1, \ldots, \mathbf{v}_n \} \subset \mathbb{R}^d$, I want to project a point $\mathbf{x}_0 \in \mathbb{R}^d$ onto the convex hull $\text{conv}(V)$ of the vectors in $V$.

I know this is a quadratic program, to find $\mathbf{z}^*$ that minimizes $\frac{1}{2}\|\mathbf{x}_0 - \mathbf{z}\|^2$ subject to $\mathbf{z} \in \text{conv}(V)$.

I also know that $\text{conv}(V)$ is a polytope and expressible as a set $S = \{ \mathbf{x} : A\mathbf{x} \le \mathbf{b} \}$.

However, I don't know how to derive the constraint matrix and constraint vector $(A,\mathbf{b})$ from the vectors in $V$.

Secondarily, I'm wondering if there are simple and fast algorithms to solve this problem. The number of vectors in $V$ will be less than $250$ and the dimensionality will be less than $50$.

And finally, I am hoping to express the solution $\mathbf{z}^* = \text{Proj}(\mathbf{x}_0)$ in something like barycentric coordinates with respect to the vectors in $V$. In other words, I'd like to express $\mathbf{z}^*$ as the vector $(\alpha_1, \ldots, \alpha_n)$ such that $\mathbf{z}^* = \sum_i \alpha_i \mathbf{v}_i$ with $\alpha_i \ge 0$ and $\sum_i \alpha_i = 1$. Given that the set $V$ won't be linearly independent (because $n > d$), I know that such barycentric coordinates are not well defined. I'm hoping to use something akin to a least norm solution $\|\alpha\|^2$ here.


You do not need to express the convex hull explicitly as a system of inequalities $\mathbf{A x} \leq \mathbf{b}$. The convex hull, by definition, is $$ \operatorname{conv}(\mathbf{v}_1, \dots, \mathbf{v}_m) = \left \{ \sum_{i=1}^m w_i \mathbf{v}_i : \sum_{i=1}^m w_i=1, w_i \geq 0 \right\} $$

Therefore, the projection onto the convex hull of a vector $\mathbf{y}$ can be computed from the solution of $$ \min_{\mathbf{w}} \quad \left \| \sum_{i=1}^m w_i \mathbf{v}_i - \mathbf{y} \right \|_2^2 \quad \text{s.t.} \quad \sum_{i=1}^m w_i=1, w_i \geq 0 $$ The vector $\mathbf{w}$ contains the barycentric coordinates of the projection. I am not aware of a specialized algorithm of computing the optimal $\mathbf{w}$, but any QP solver can do so.

A simple solution method which you can implement yourself can be obtained from a fast projected gradient algorithm ,such as FISTA, and remembering that the constraints are exactly the unit simplex, and there is a simple $O(m \log m)$ algorithm for projecting a point onto the unit simplex. Moreover, you also have the well-known Mirror Descent algorithm for optimizing over the unit simplex, which can be quite fast in practice.

  • $\begingroup$ $+1$ for an elegant answer .. I think the optimal way would be a brute-force of order $O(N^m)$ where $N$ would be the size of the grid, or more generally $O(\prod_{k=1}^m N_k)$ if each $w_k$ is searched on a different grid of size $N_k$. $\endgroup$ – Ahmad Bazzi Feb 12 at 12:10
  • $\begingroup$ @AhmadBazzi, what grid are you referring to? A grid defined on $[0,1]^m$? I believe that this becomes impractical even for a small value of $m$, say, $m \geq 5$. $\endgroup$ – Alex Shtof Feb 12 at 12:24
  • $\begingroup$ True .. that is the optimal and yet an impractical solution. $\endgroup$ – Ahmad Bazzi Feb 12 at 12:33
  • $\begingroup$ @AlexShtof - what's a good reference on mirror descent? $\endgroup$ – ted Feb 12 at 18:52
  • $\begingroup$ This one: users.cecs.anu.edu.au/~xzhang/teaching/bregman.pdf $\endgroup$ – Alex Shtof Feb 12 at 19:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.