Reformulating the problem with arbitrary boundaries (Preserving $ \lambda $ for the Lagrange Multiplier):
$$\begin{align*}
\arg \min_{x} \quad & \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} \\
\text{subject to} \quad & {\left\| x \right\|}_{1} \leq r \\
& {l}_{i} \leq {x}_{i} \leq {u}_{i} \; \forall i = 1, 2, \ldots, n
\end{align*}$$
The Lagrangian is given by (While the box constraints are implied):
$$ L \left( x, \lambda \right) = \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} + \lambda \left( {\left\| x \right\|}_{1} - r \right) = \sum_{i = 1}^{n} \left[ \frac{1}{2} {\left( {x}_{i} - {y}_{i} \right)}^{2} + \lambda \left| {x}_{i} \right| \right] - \lambda r $$
Namely the problem can be formulated component wise.
The dual function is given by:
$$ q \left( \lambda \right) = \inf_{ {l}_{i} \leq {x}_{i} \leq {u}_{i} } L \left( x, \lambda \right) = \inf_{ {l}_{i} \leq {x}_{i} \leq {u}_{i} } \left\{ \sum_{i = 1}^{n} \left[ \frac{1}{2} {\left( {x}_{i} - {y}_{i} \right)}^{2} + \lambda \left| {x}_{i} \right| \right] - \lambda r \right\} $$
Now, let's think about the solution a little bit.
The solution $ {x}_{i}^{\ast} $ would like to have the same sign of $ {y}_{i} $ as if it would have different sign being $ 0 $ would clearly be a better solution (Lower objective value both for the squared term and for the absolute value term).
This above is the basic intuition for the trick done by @gerw.
Now, imagine we know the sign of each element of the solution - $ {s}_{i} = \operatorname{sign} \left( {x}_{i} \right) $.
Then the problem is given as:
$$\begin{align*}
\arg \min_{x} \quad & \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} \\
\text{subject to} \quad & {s}^{T} x \leq r \\
& {l}_{i} \leq {x}_{i} \leq {u}_{i} \; \forall i = 1, 2, \ldots, n
\end{align*}$$
Then the Lagangian is given by:
$$ L \left( x, \lambda \right) = \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} + \lambda \left( {s}^{T} x - r \right) = \sum_{i = 1}^{n} \left[ \frac{1}{2} {\left( {x}_{i} - {y}_{i} \right)}^{2} + \lambda {s}_{i} {x}_{i} \right] - \lambda r $$
Namely the problem can be formulated component wise.
The dual function is given by:
$$ q \left( \lambda \right) = \inf_{ {l}_{i} \leq {x}_{i} \leq {u}_{i} } L \left( x, \lambda \right) = \inf_{ {l}_{i} \leq {x}_{i} \leq {u}_{i} } \left\{ \sum_{i = 1}^{n} \left[ \frac{1}{2} {\left( {x}_{i} - {y}_{i} \right)}^{2} + \lambda {s}_{i} {x}_{i} \right] - \lambda r \right\} $$
It is easy to show that the solution $ {x}^{\ast} = \operatorname{Proj}_{\left[ l, u \right]} \left( y - \lambda s \right) $ obeys all the KKT conditions of the problem.
One of the KKT conditions is given by:
$$ \lambda \left( {s}^{T} x - r \right) = 0 \Rightarrow \lambda \left( {s}^{T} \operatorname{Proj}_{\left[ l, u \right]} \left( y - \lambda s \right) - r \right) = 0 $$
From the above it is easy to derive the optimal $ \lambda $ and then by it to extract $ {x}^{\ast} $.
The only question is how to have $ s $ without having $ {x}^{\ast} $ first?
Keep in mind the objective function $ f \left( \lambda \right) = \lambda \left( {s}^{T} \operatorname{Proj}_{\left[ l, u \right]} \left( y - \lambda s \right) - r \right) $ which we need to find the root of.
So we can think about 3 exclusive cases:
- $ \operatorname{sign} \left( {l}_{i} \right) = \operatorname{sign} \left( {u}_{i} \right) \Rightarrow {s}_{i} = \operatorname{sign} \left( {l}_{i} \right) = \operatorname{sign} \left( {u}_{i} \right) $.
This also suggest that for next cases $ {u}_{i} \geq 0 $ and $ {l}_{i} \leq 0 $ otherwise we're back into this case.
- $ \operatorname{sign} \left( {y}_{i} \right) = \operatorname{sign} \left( {l}_{i} \right) $ (Implicitly $ {l}_{i} \leq 0 $ and $ {y}_{i} \leq 0 $) hence $ {s}_{i} = \operatorname{sign} \left( {y}_{i} \right) $ hence we must update $ {u}_{i} = 0 $ as in the objective function we can't allow $ {x}_{i}^{\ast} > 0 $ as it means its sign will reverse.
- $ \operatorname{sign} \left( {y}_{i} \right) = \operatorname{sign} \left( {u}_{i} \right) $ (Implicitly $ {u}_{i} \geq 0 $ and $ {y}_{i} \geq 0 $) hence $ {s}_{i} = \operatorname{sign} \left( {y}_{i} \right) $ hence we must update $ {l}_{i} = 0 $ as in the objective function we can't allow $ {x}_{i}^{\ast} < 0 $ as it means its sign will reverse.
So we have set of rules to set $ s $ and update the boundaries (Component wise as this function has component wise property) to validate the objective function.
I wrote a MATLAB code which implements the solution at my Mathematics StackExchange Question 2824418 - GitHub Repository.
The code validates the result to a reference calculated by CVX.