# Orthogonal Projection onto the Unit Simplex

The Unit Simplex is defined by:

$$\mathcal{S} = \left\{ x \in \mathbb{{R}^{n}} \mid x \succeq 0, \, \boldsymbol{1}^{T} x = 1 \right\}$$

Orthogonal Projection onto the Unit Simplex is defined by:

\begin{alignat*}{3} \arg \min_{x} & \quad & \frac{1}{2} \left\| x - y \right\|_{2}^{2} \\ \text{subject to} & \quad & x \succeq 0 \\ & \quad & \boldsymbol{1}^{T} x = 1 \end{alignat*}

How could one solve this convex optimization problem?

Projection onto the Simplex can be calculated as following.
The Lagrangian in that case is given by:

\begin{align} L \left( x, \mu \right) & = \frac{1}{2} {\left\| x - y \right\|}^{2} + \mu \left( \boldsymbol{1}^{T} x - 1 \right) && \text{} \\ \end{align}

The trick is to leave non negativity constrain implicit.
Hence the Dual Function is given by:

\begin{align} g \left( \mu \right) & = \inf_{x \succeq 0} L \left( x, \mu \right) && \text{} \\ & = \inf_{x \succeq 0} \sum_{i = 1}^{n} \left( \frac{1}{2} { \left( {x}_{i} - {y}_{i} \right) }^{2} + \mu {x}_{i} \right) - \mu && \text{Component wise form} \end{align}

Taking advantage of the Component Wise form the solution is given:

\begin{align} {x}_{i}^{\ast} = { \left( {y}_{i} - \mu \right) }_{+} \end{align}

Where the solution includes the non negativity constrain by Projecting onto $${\mathbb{R}}_{+}$$

The solution is given by finding the $$\mu$$ which holds the constrain (Pay attention, since the above was equality constrain, $$\mu$$ can have any value and it is not limited to non negativity as $$\lambda$$).

The objective function (From the KKT) is given by:

\begin{align} 0 = h \left( \mu \right) = \sum_{i = 1}^{n} {x}_{i}^{\ast} - 1 & = \sum_{i = 1}^{n} { \left( {y}_{i} - \mu \right) }_{+} - 1 \end{align}

The above is a Piece Wise linear function of $$\mu$$.

Since the function is continuous yet it is not differentiable due to its piece wise property theory says we must use derivative free methods for root finding. One could use the Bisection Method for instance.

The function Derivative given by:

\begin{align} \frac{\mathrm{d} }{\mathrm{d} \mu} h \left( \mu \right) & = \frac{\mathrm{d} }{\mathrm{d} \mu} \sum_{i = 1}^{n} { \left( {y}_{i} - \mu \right) }_{+} \\ & = \sum_{i = 1}^{n} -{ \mathbf{1} }_{\left\{ {y}_{i} - \mu > 0 \right\}} \end{align}

In practice, it can be solved using Newton Iteration (Since falling into a joint between 2 sections has almost zero probability).

## Accurate / Exact Solution

If we look at the values of the function $$h \left( \mu \right) = \sum_{i = 1}^{n} { \left( {y}_{i} - \mu \right) }_{+} - 1$$ one could easily infer a method to calculate the accurate solution:

In the above the parameter $$\mu$$ took the values of the vector $${y}_{i}$$ with additional values at the edges (Value greater than the max value of $${y}_{i}$$ and value lower of the min value of $${y}_{i}$$).
By iterating the values one could easily track the 2 values which on each side they have value greater than $$0$$ and lower then $$0$$ (In case one of them is zero, then it is the optimal value of $$\mu$$). Since it is linear function and we have 2 points we can infer all parameters of the model $$y = a x + b$$. Than the optimal value of $$\hat{\mu} = - \frac{b}{a}$$.

I wrote MATLAB code which implements the method with Newton Iteration at Mathematics StackExchange Question 2327504 - GitHub. I extended the method for the case $$\sum {x}_{i} = r, \; r > 0$$ (Pseudo Radius).
There is a test which compares the result to a reference calculated by CVX.

• The constraint regarding $\mu$ is given by $h \left( \mu \right) = 0$. I guess objective function isn't the most clear term :-). But in the case above it is the objective function of the Dual Function. – Royi Apr 3 '18 at 15:47
• Yep. You got it right. – Royi Apr 4 '18 at 13:15
• It seems to me that the minimization problem boils down to find $\mu$ such that $h \left( \mu \right)= 0$ or $\sum_{i = 1}^{n} { \left( {y}_{i} - \mu \right) }_{+} - 1 = 0$, right? – LE Anh Dung Feb 28 '20 at 8:35
• Yep. Indeed this is the case from the above (KKT). – Royi Feb 28 '20 at 8:37
• I'm sorry but your solution seems much more simpler than the very complicated one proposed by Yunmei Chen and Xiaojing Ye. May I miss something? – LE Anh Dung Feb 28 '20 at 8:41

The paper by Condat [1] presents a review and comparison of existing algorithms with a new proposal for projection onto the unit simplex. This paper lists the worst-case complexity and empirical complexity of those algorithms, and presents concise pseudo-code for all algorithms. In particular, the algorithm proposed by Condat takes $$O(n)$$ time in practice, whereas sorting-based methods take $$O(n \log n)$$ time in practice. I have implemented Condat's algorithm in the past, and can vouch for its speed relative to direct sorting-based approaches.

EDIT: Condat has included C and MATLAB implementations of all the algorithms mentioned in his paper here: https://lcondat.github.io/software.html

• Feel like posting your code for Condat's algorithm? I think it could be helpful for a lot of people. – littleO Apr 20 '20 at 4:16
• @littleO Thanks for the suggestion. I've included a link to Condat's webpage which has codes for all the algorithms. My own code was in one of the early versions of Julia, which is almost surely outdated – madnessweasley Apr 20 '20 at 4:24
• My code above is also ‘log(n)’ approach. – Royi Apr 20 '20 at 6:44
• @Royi I don't understand your explanation for the "accurate/exact solution", but each function evaluation in the bisection step takes $\approx n$ operations. Additionally, even evaluating $h(\mu)$ takes $\approx n$ operations, so I'm not sure if $O(\log n)$ is even possible – madnessweasley Apr 20 '20 at 6:52
• @Royi I guess the C code includes the results of the paper, and the MATLAB code only includes a simple sorting-based approach. You are probably correct that he hasn't compared against your exact approach (which might be faster). I'm sorry, I don't have the time to delve deeper into the numerical comparison. Condat seems to have made a fair (for me) numerical comparison with the literature, and because I had a favorable experience with his approach, I just wanted to share it with the community. – madnessweasley Apr 20 '20 at 10:31

The best algorithm to compute the exact solution to this problem can be found in Projection Onto A Simplex.

• would be awesome if the key ideas can be summarized in the post. – Siong Thye Goh Mar 21 '18 at 2:23
• @SiongThyeGoh, It based on the same idea as my solution below. The main difference is that instead of utilizing Newton Method to find the optimal $\mu$ they use some other trick based on the sorted values of the vector $y$. – Royi Mar 21 '18 at 13:15
• This idea is much older than that paper. Late 80's, I believe, a JOTA paper if I remember. – copper.hat Feb 29 '20 at 6:25
• I extracted exact solution from my approach above. – Royi Apr 14 '20 at 9:40
• I think the algorithm from Condat's 2014 paper (linked to in another answer) might be better than this one. – littleO Apr 20 '20 at 4:21