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!
