Constrained Optimization for Probability Masses

Question

It is well known that the weighted Pool Adjacent Violators Algorithm can be used to solve the following quadratic program:

$\arg \min_x \sum_{i = 1}^k w_i(x_i - y_i)^2$

$\text{such that } x_j \leq x_{j+1}$

This can be used in optimization problems with monotonic constraints, where the off diagonal of the Hessian is ignored. Very relevant to my interest, this problem is solved in $O(k)$ time.

My question is where there is some sort of similar solver for probability constraints? In other words, I am interested in solving the quadratic program

$\arg \min_p \sum_{i = 1}^k w_i(p_i - y_i)^2$

$\text{such that } p_j \geq 0 \text{ and } \sum_{i = 1}^n p_i = 1$

Of course, you could use a generic quadratic solver for this problem. However, for my problem, the complexity of the solver is vital, as in the statistical problem of interest the length of $p$ can grow linearly with the sample size.

Answer 1

In case anyone is curious about this, here's the alternative I've come up with in the mean time.

Ideally, I would have liked a solution to

$\arg \min_p \sum_{i = 1}^k w_i(p_i - y_i)^2$

$\text{such that } p_j \geq 0 \text{ and } \sum_{i = 1}^n p_i = 1$

which could be solved in $O(k)$. With this solver, we could have approximated the negative log likelihood function as a quadratic function of $p$ (ignoring the off diagonal of the Hessian) and then used our solver for Sequential Quadratic Programming.

Without having a solver which could be used in linear time, I switched to a gradient ascent step. The probability constraints were respected by manipulating the stepping direction. To do this, I start with the standard step direction of the gradient. Then, the step direction of inactive points (i.e. parameters such that $p_j = 0$) were set to 0 (there's another step of the algorithm that will "resurrect" these point should they be necessary). Second, for each active point, I subtracted off the average gradient of the active points. This insures that the stepping direction sums to 0, insuring that if $\sum p^{(t)} = 1$, then $\sum p^{(t+1)} = 1$ as well. It is trivial to show that this still results in a directional step that will increase the log likelihood function. Finally, if the proposed step is too far (i.e. $p^{(t+1)} < 0$), the entire step is scaled back until the non-negative constraint is met.

How does it work? About as expected. That is to say, it is very helpful in my situation: it gets my algorithm (which includes several other updating steps, including a reparameterization of the probability parameters such that the new parameters are increasing on $\mathbb{R}$ which then uses the pool adjacent violators step mentioned earlier) out of the special cases in would get stuck in before. But were this algorithm to use the gradient ascent step alone to update the baseline probabilities, it would be much too slow.

If there's a more clever way to do this, I'm still very interested!