I have to solve a quadratic problem with simplex constraints, e.g. $$\min \frac{1}{2}x^TQx + qx\\ s.t. \sum\limits_{i=0}^{n} x_i = 1\\ x_i \geq 0$$ with a projected gradient method (Q is positive semidefinite). In particular I have to project specifically the (anti)gradient over the feasible set, then do a stepsize over the direction found.

My problem is that I have to find the active set of the constraint, which is defined as the set of inequality constraints satisfied as equality. But what happens when I'm in a point such as $$(0.5, 0.5)$$? This is clearly feasible, but none of the inequality constraints are satisfied as equality.

Also I don't know how to properly handle the equality constraint during the projection phase.

I also have to do the projection in $$\mathcal{O}(n\log n)$$ time, which hints to sorting, but I don't have a very clear idea on what to do.

EDIT: My main problem is how do I project the gradient onto the active set of constraints? And which is the maximum feasible step when $$d'Qd = 0$$?

• Hint: the constraint $\sum_{i=1}^{n} x_{i}=1$ will always be active. The constraints $x_{i}>=0$ will sometimes be active and sometimes not. – Brian Borchers Jul 17 '19 at 14:43
• @BrianBorchers thank you for the hint, but how can I build the matrix made up by the active constraints? I mean, they have to be inequalities, but if I translate the = into <= I get two linearly dependent rows so the resulting matrix would be singular and I cannot proceed. – Agent 8 Jul 17 '19 at 14:59
• Why don't you write that up and add it to your question? In fact, the system of equations associated with the active constraints will have solutions, but it's not clear where your misunderstanding lies. – Brian Borchers Jul 17 '19 at 15:39
• @BrianBorchers My main problem is that I don't know a clever way to do projection. My professor says it's simple like if we had box constraints, but I can't think of a way to do it. – Agent 8 Jul 18 '19 at 14:12
• Have you tried a simple google search for "projection on the unit simplex"? – Brian Borchers Jul 20 '19 at 15:18

## Projection onto the Unit Simplex

### Problem

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*}

### Solution

The Lagrangian 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} 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$$ and its 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}

Hence it can be solved using Newton Iteration.

I wrote MATLAB code which implements them both at Mathematics StackExchange Question 2338491 - GitHub.
There is a test which compares the result to a reference calculated by CVX.