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\|$.

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\|$.

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\| + b^Tx\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\|$.

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\|$.

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 ?

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\| + b^Tx\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\|$.

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\| + b^Tx\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\|$.