Let $\Delta_n := \{x \in \mathbb{R}^n | x \ge 0, \sum_{1 \le i \le n}x_i = 1\}$ be the $n$-simplex. For $a, b \in \mathbb R^n$, with $\Delta_n \not \ni a$, consider the problem of computing the following value exactly $$\alpha := \min_{x \in \Delta_n}\|x-a\| + \langle b,x\rangle.$$ The case when $b = 0$ corresponds to the problem of projecting the point $a$ unto $\Delta_n$. This special case is well-known to be exactly solvable in $\mathcal{O}(n)$ flops.

Now, for the general problem, it's not hard to rewrite

\begin{equation} \begin{split} \alpha &:= \min_{x \in \Delta_n}\max_{\|y\| \le 1} \langle y, x - a \rangle + \langle b, x\rangle = \max_{\|y\| \le 1}\left(\min_{x \in \Delta_n}\langle y + b, x\rangle\right) - \langle y, a\rangle\\ &= \max_{\|y\| \le 1}\left(\min_{1 \le i \le n}y_i + b_i\right) - \langle y, a \rangle = \max_{t \in \mathbb R}t - \min_{\|y\|^2 \le 1,\;y_i \ge t - b_i \forall i}\langle y,a\rangle. \end{split} \end{equation}

Question: Given $c \in \mathbb R^n$ ($c_i \equiv b_i - t$), can one compute the value of $$\min_{\|y\|^2 \le 1,\;y \le c}-\langle a,y\rangle$$ analytically ? For which values of $c$ does the latter problem have a solution ?

No polynomial time algorithm for exact solution ?

Observation: In low dimensional cases ($n = 1, 2, 3$), when non-degenerate, it's not hard to sketch that this problem has solutions which are piece-wise polynomials (or square roots of such) in the $a_i$'s and $c_i$'s, the number of pieces being in the order of $2^n$.

Update: Algorithm based on @fedja's answer + proof of Q-linear convergence in the "small perturbation" regime

For generality, let $a \in \mathbb R^n$ and $C$ be a "simple" (to be clarified) closed convex subset of $\mathbb R^n$ not containing the point $a$. Consider the problem \begin{equation} \text{minimize } \|x - a\| + \langle b, x\rangle\text{ subject to }x \in C. \end{equation}

Fedja's idea. The idea is to introduce a radial variable $r := \|x-a\|^{-1}$. Indeed, using the well-known elementary inequality \begin{equation} t + t^{-1} \ge 2\; \forall t > 0,\text{ with equality iff } t = 1, \end{equation} it follows that $\forall x \in \mathbb{R}^n$ and $\forall r > 0$, we have $$\|x-a\| \le \frac{1}{2}(r\|x-a\|^2 + r^{-1}),$$ and this bound is attained at $r = \|x-a\|^{-1}.$ Thus completing the square in $x$, the optimal value $\alpha$ for the problem can be rewritten in the form \begin{equation} 2\alpha- b^Ta = \min_{r > 0,x \in C}r\|x - (a + r^{-1}b)\|^2 + (1-b^Tb)r^{-1}. \end{equation}

The algorithm. Based on Fedja's idea of introducing the radial variable $r = \|x-a\|^{-1}$, the following alternating iterative scheme for solving the above problem exactly, is natural \begin{equation} x^{(k + 1)} = \mathrm{proj}_C(a + b / r^{(k)}),\; r^{(k+1)} := \|x^{(k + 1)} - a\|^{-1}, \end{equation} with $r^{(0)} > 0$. Next, we will proof some interesting things regarding the numerics for the algorithm so-obtained.

Q-linear convergence in the "small perturbation" regime: $\|b\| < 1$. Eliminating the $x$ variable from the above iterates, we see that the $r$-update can be rewritten as a Picard process \begin{equation}r^{(k+1)} = \|\mathrm{proj}_C(a + b/r^{(k)}) - a\|^{-1}. \end{equation} Let's proof that this process converges after essentially $\mathcal O(1)$ rounds. This would mean that the cost of solving the "$b \ne 0$" problem is essentially the cost of projecting onto $C$, i.e the cost of solving the "$b=0$" problem.

Assumption. Let's make the following minimal assumption: There is an oracle for exactly projecting onto $C$.

For example, this assumption holds for polyhedra, $\ell_p$ balls, etc. In order to only concentrate on the "heart of the issue'', let's take for granted that the above process has a fixed-point $r^{(*)} > 0$. Now,

\begin{eqnarray*} \begin{split} \left|\frac{1}{r^{(k+1)}} - \frac{1}{r^{(*)}}\right| &= \Big|\|\mathrm{proj}_{ C}(a + b/r^{(k)}) - a\| -\|\mathrm{proj}_{C}(a + b/r^{(*)}) - a\|\big|\\ &\le \|\mathrm{proj}_{C}(a + b/r^{(k)}) - \mathrm{proj}_{C}(a + b/r^{(*)})\| \\ &\le \|b/r^{(k)} - b/r^{(*)}\| = \|b\|\Big|\frac{1}{r^{(k)}} - \frac{1}{r^{(*)}}\big|, \end{split} \end{eqnarray*} where the first inequality is the "reverse triangle inequality'' and the second follows from the nonexpansivity of the projection operator onto a closed convex set. Thus if $\|b\| < 1$, then the residuals $\Big|\frac{1}{r^{(k)}} - \frac{1}{r^{(*)}}\Big|$ decay exponentially fast (i.e Q-linearly) at rate $\|b\|$.

  • 1
    $\begingroup$ btw, projection onto simplex is solvable in $O(n)$ also, no need of the extra $\log n$ in there.... $\endgroup$
    – Suvrit
    Dec 6, 2015 at 3:08
  • $\begingroup$ @Suvrit: For the problem of projecting onto a simplex, there are indeed randomized algorithms (for example based on randomized median-finding) with linear average running time, but no algorithm with linear worst-case running time is known, nor is likely to exist... $\endgroup$
    – dohmatob
    Dec 6, 2015 at 10:44
  • 2
    $\begingroup$ But there exist methods for worst-case linear time median finding (they might be empirically a bit slower, because of constants in the big-Oh) --- so the true complexity of this problem is $O(n)$ (using pivoting methods with worst-case linear time median finding routines). And to disprove your claim on "no method is likely to exist" have a look at: gipsa-lab.grenoble-inp.fr/~laurent.condat/publis/… $\endgroup$
    – Suvrit
    Dec 6, 2015 at 15:08
  • $\begingroup$ Please re-read the Condat paper carefully. The algorithm proposed in the paper (Algorithm 3) has linear expected running time, not linear worst-case running time. Even the authors of that paper admit this fact cautiously. I'll offer you a dollar if you can exhibit (of course, with proofs!) an algorithm for projecting onto the simplex with linear worst-case running time :) $\endgroup$
    – dohmatob
    Dec 6, 2015 at 15:34
  • 1
    $\begingroup$ Thank you for writing up a verbose version of fedja's answer to make it easier to follow for folks like me $\endgroup$
    – Yi Liu
    Dec 21, 2015 at 23:39

1 Answer 1


I would try to approach your original problem a bit differently. Note that $|x-a|\le \frac 12(r|x-a|^2+r^{-1})$ and the equality is attained for $r=|x-a|^{-1}$. Thus, $$ \min_x[|x-a|+\langle b,x\rangle]=\frac 12\min_r\min_x[r|x-a+r^{-1}b|^2+r^{-1}(1-|b|^2)]+\langle a,b\rangle $$ However, the inner minimum can be now found rather quickly and finding the outer one is a one-dimensional problem. Moreover, I suspect that even the naive algorithm of choosing $r$ arbitrarily to start with, finding the corresponding $x$, and then changing $r$ to something stupid like $(|x-a|^{-1}+r)/2$ has a decent chance to work but checking that would require some accurate analysis of the properties of the inner minimum as a function of $r$. What is true, however, is that if the minimizer $x$ for some $r$ satisfies $|x-a|=r^{-1}$, then you have the true minimum for the original expression, so, at least, you can recognize the solution when you see it this way.


I thought it might make sense to add a few more remarks to what I wrote already.

Let's look at what we have so far and see if we can convert it into something that is guaranteed to work. Put $\rho=r^{-1}$. Then the main idea was to observe that instead of the original minimization problem, we can solve the problem $$ \Psi(x,\rho)=\frac 12\left[\frac{|x-a|^2}\rho+\rho\right]+\langle b,x\rangle\to\min $$ and that in this new problem the minimization in $x$ is equivalent to finding the closest point to $a-\rho b$.

Since $\Psi$ is jointly convex in $x,\rho$ and the best $\rho$ has physical meaning of the distance from $a$ to $x$, the outer minimization problem is just the problem of finding the minimum of a convex function over a finite interval of length about the diameter of the convex body $C$ we project upon, so, even if the iterations fail, we can always resort to the good old bisection, which is guaranteed to give us the answer with the precision about $2^{-m/2}|b|\operatorname{diam} C$ where $m$ is the number of steps.

This is not an "exact" solution, of course. However, if you implement an exact algorithm on any real machine, you'll still never be able to get anything better than the machine precision. So, I'm not really sure how much sense the request for an exact solution makes here from the purely practical viewpoint. Still, notice that the closest point to $a-\rho b$ moves along a straight line as long as we stay in the relative interior of the same $k$-dimensional face of the simplex. Unfortunately, this gives us exponentially many lines, so just tracing all of them is impractical. However, we can always first run the bisection for a while until we reduce the situation to a sufficiently small interval, after which we can use only the lines that correspond to the points on that interval and, in principle, we can try to trace them. The ideal scenario is that after a few iterations in the bisection, we get the values of $x$ lying in the relative interior of the same $k$-dimensional face for both endpoints of the remaining interval, in which case we have just one line left and we can finish in a single step. It is not guaranteed to always happen, of course, but I suspect that it may work more often than not.

  • $\begingroup$ The introduction of the radial parameter r is brilliant, and didn't occur to me. Indeed your proposed scheme converges in essentially 1 iteration. I've augmented my question with a rigorous proof of this fact (with a slight modification of the r update). You have my upvote. I leaving the bounty open for a while to see if I can attract answers of similar quality. $\endgroup$
    – dohmatob
    Dec 20, 2015 at 17:57
  • $\begingroup$ @dohmatob Thanks. However, are you sure that it is $\gamma^{-1}b$ and not just $b$ in your scaling? It looks to me that when you change $C$ to $\gamma C$ and $a$ to $\gamma a$, you would like the whole expression, including the scalar product, to be multiplied by $\gamma$ to get the desired homogeneity. Of course, the case of small $|b|$ is very nice because the change in the distance always dominates the change in the scalar product but it would be interesting to see what happens with moderate size $|b|$ as well. I would just run the iterations and see before trying to prove anything. $\endgroup$
    – fedja
    Dec 21, 2015 at 4:14
  • $\begingroup$ Yes you're right, it was too good to be true. The $\gamma^{-1}$ factor next to the $b$ iwas an error (fixed). I've scaled down my claim to the statement: If $\|b\| < 1$, then Fedja's algorithm convergences Q-linearly at rate $\|b\|$. The proof should contain no gaps now. About why, I'm interested in making this rigorous, I stated in the original problem that I wanted the solution method to be exact. Since all we have at the moment is an iterative scheme, the least of things would be to show that it converges in finitely many steps (or a slightly weaker statement). Thus the proof above. $\endgroup$
    – dohmatob
    Dec 21, 2015 at 14:49
  • $\begingroup$ BTW, indeed numerically, the algorithm is observed to converge in finitely-many steps (typically less than $5$). This is what motivated the proof. $\endgroup$
    – dohmatob
    Dec 21, 2015 at 14:53
  • 1
    $\begingroup$ @dohmatob Yeah, but the same degeneracy arises also when $|b|=1$ and the ray $a-\rho b$ intersects the interior of the simplex (I doubt you can exclude such configurations as well). So, sometimes the iterations work poorly while the bisection never fails. What I really wonder about at this moment is how complicated the projection of a line to the simplex can be. If we could show that the corresponding broken line never has too many very short segments, we could make something provably finite out of the current scheme. $\endgroup$
    – fedja
    Dec 21, 2015 at 21:19

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.