I found the answer in Chapter 3 of the book "Constrained minimization and Lagrange multiplier methods" by Dimitri P. Bertsekas.
Consider the following one-sided constrained minimization problem,
\begin{equation}
\min f(x) \quad \text{subject to} \quad \alpha_j < g_j(x), \qquad j=1,...,r
\end{equation}
The above problem is equivalent to the following problem,
\begin{equation}
\min f(x) \quad \text{subject to} \quad \alpha_j < g_j(x) - u_j, \qquad u_j=0,\quad j=1,...,r
\end{equation}
We consider a multiple method for the above minimization where only the constraints $u_j=0$ are eliminated by means of a quadratic penalty function, viz.
\begin{equation}
\min f(x) + \sum_{j=1}^r \{ \mu_j^k u_j + \frac{1}{2} c_k |u_j|^2 \} \qquad \text{subject to} \quad \alpha_j < g_j(x) -u_j \quad j=1,...,r
\end{equation}
The above minimization can be carried out first with respect to $u_j$ yielding the following equivalent problem
\begin{equation}
\begin{gathered}
\min f(x) + \sum_{j=1}^r p_j[g_j(x),\mu_j^k,c^k] \qquad \text{subject to} \quad x\in R^n \\
\text{where} \quad p_j[g_j(x),\mu_j^k,c^k]=\min_{\alpha_j \leq g_j(x)-u_j} \{ \mu_j^k u_j + \frac{1}{2}c^k |u_j|^2 \}.
\end{gathered}
\end{equation}
Performing the above minimization with respect to $u_j$ we obtain the expression $u_j=-\mu_j^j/c^k$ for $u_j$. We only consider the upper limit for $u_j$, i.e.,
\begin{equation}
g_j(x)-u_j > \alpha_j \longrightarrow g_j(x)-\alpha_j > u_j,
\end{equation}
which means that $g_j(x)-\alpha_j \leq u_j$ hits the boundary of the admissible range. Thus
\begin{equation}
g_j(x)-\alpha_j \leq -\frac{\mu_j^k}{c^k} \longrightarrow \mu_j^k + c^k (g_j(x)-\alpha_j)\leq0
\end{equation}
lies on the surface and outside the surface of admissible range and if the above equation holds then $u_j$ must remain on its admissible surface value, i.e. $u_j \rightarrow g_j(x)-\alpha_j$. Consequently, the value of $u_j$ is calculated as
\begin{equation}
u_j=\begin{cases}
g_j(x)-\alpha_j & \text{if} \quad \mu_j^k + c^k (g_j(x)-\alpha_j)\leq0, \\ -\mu_j^k/c^k & \text{if otherwise}.
\end{cases}
\end{equation}
Substituting the value of $u_j$ into the expression of $p_j$, the augmented Lagrangian function is obtained as
\begin{equation}
p_j[g_j(x),\mu_j^k,c^k]=\begin{cases}
\mu_j^k (g_j(x)-\alpha_j) + \frac{1}{2} c |g_j(x)-\alpha_j|^2 & \text{if} \quad \mu_j^k + c^k (g_j(x) - \alpha_j)\leq0, \\
-\mu_j^{k,2}/2c^k & \text{if otherwise}.
\end{cases}
\end{equation}
In the elastic frictionless contact problem that is in the question, the terms appearing in the expression of $p_j$ are as following,
\begin{equation}
g_j(x)=g_N(x), \quad \mu_j^k=\lambda_j^k,\quad c^k=\varrho,\quad \alpha_j=0,
\end{equation}
which upon substituting in the expression of $p_j$ yields the following augmented Lagrangian function as,
\begin{equation}
p_j[g_j(x),\lambda_j^k]=\begin{cases}
(\lambda_j^k+\frac{\varrho}{2}g_N)g_N, & \text{if} \quad \lambda_j^k + \varrho g_N \leq 0, \\
-\lambda_j^{k,2}/2\varrho & \text{if otherwise}.
\end{cases}
\end{equation}
Note that the first condition in the above equation shows the contact state between the contactor and the target and the second condition shows the separation state between the contactor and the target.
The way that the optimization problem (P1) has been reformulated here, the $x$-update is indeed not easy. An iterative algorithm would be required just to solve that subproblem. So what we should do is reformulate the problem differently so that both updates are easy.
In this case, you can just express the optimization problem (P1) as
$$
\tag{1} \text{minimize} \quad f(x) + \underbrace{\frac12 \| Ax - r \|_2^2}_{g(x)}
$$
and minimize $f+g$ using the Douglas-Rachford method (which is a special case of ADMM).
(Here the optimization variable is $x$, and $f$ is the indicator function of the $\ell_1$-norm ball of radius $b$.)
Evaluating the proximal operator of $g$ reduces to solving a linear system of equations involving the matrix $A$. Evaluating the proximal operator of $f$ is equivalent to projecting a point onto the $\ell_1$-norm ball of radius $b$. (One way to do that is explained on slide 6-15 in Vandenberghe's UCLA 236c notes.)
By the way, because $g$ is differentiable, you can also solve problem (1) using the proximal gradient method or an accelerated proximal gradient method such as FISTA. I'd bet FISTA would be faster than ADMM, and also these methods have two other advantages: 1) There is no need to solve a large linear system of equations at each iteration; 2) Line search procedures are available (so there's no need to laboriously tune your step size, as is often necessary with ADMM).
Best Answer
ADMM steps (from https://web.stanford.edu/~boyd/papers/pdf/admm_distr_stats.pdf) can be
\begin{align*} {{x}} &\leftarrow \arg\min_{{{x}} } f \left( {{x}} \right) + {{u}} ^T \left( {{{x}} } - {{z}} \right) + \frac{\rho}{2} \left\| {{{x}} } - {{z}} \right\|_2^2 \\ &\equiv \arg\min_{{{x}} } f \left( {{x}} \right) + \frac{\rho}{2} \left\| {{{x}} } - {{z}} + {u} \right\|_2^2 \\ {{z}} &\leftarrow \arg\min_{{{z}}} g\left( {{z}} \right) + {{u}} ^T \left( {{{x}} } - {{z}} \right) + \frac{\rho}{2} \left\| {{{y}} } - {{x}} \right\|_2^2 \\ &\equiv \arg\min_{{{{z}}}} g\left( {{z}} \right) + \frac{\rho}{2} \left\| {{{x}} } - {{z}} + {u} \right\|_2^2 \\ {{u}} &\leftarrow {{u}} + \left( {{{x}} } - {{z}} \right) \end{align*}
Let us say $f(x) = \frac{1}{2} \|A x - b \|_2^2$ and $g(z) = \lambda \| z\|_1$. We can exploit proximal operator, that is,
Also, define for brevity (we can use the equivalent scaled-form ), $$F(x) := f \left( {{x}} \right) + \frac{\rho}{2} \left\| {{{x}} } - {{z}} + {u} \right\|_2^2$$
$$G(z) := g \left( {{z}} \right) + \frac{\rho}{2} \left\| {{{x}} } - {{z}} + {u} \right\|_2^2 .$$
Now, just find the gradients and set them to zero, that is,
$$\frac{\partial F(x)}{\partial x} = 0 \Longleftrightarrow \frac{1}{\rho}\partial f(x) + \left(x - z + u \right) = 0 \Longleftrightarrow x = \left(I + \frac{1}{\rho} \partial f \right)^{-1} \left( z - u\right) = \operatorname{prox}_{\frac{1}{\rho} f}\left( z - u\right)$$
and
$$\frac{\partial G(z)}{\partial z} = 0 \Longleftrightarrow z = \left(I + \frac{1}{\rho} \partial g \right)^{-1} \left( x + u\right) = \operatorname{prox}_{\frac{1}{\rho} g}\left( x + u\right).$$
Thus, the ADMM iterative steps are \begin{align*} {{x}^{k+1}} &:= \operatorname{prox}_{\frac{1}{\rho}f}\left( z^{k} - u^{k} \right) \\ {{z}^{k+1}} &:= \operatorname{prox}_{\frac{1}{\rho}g}\left( {{x}^{k+1}} + u^{k} \right) \\ {{u}^{k+1}} &:= {{u}^k} + \left( {{x}^{k+1}} - {{z}^{k+1}} \right) \end{align*}
Now, you can use the prox operators for both affine $f(x)$ and L1 norm $g(z)$.
Appendix
The prox operators for $f(x) = \frac{1}{2} \|A x - b \|_2^2$ and $g(z) = \lambda \| z\|_1$ are given below.
\begin{align} \operatorname{prox}_{\lambda f}\left( x \right) &= \arg\min_{v} \left\{ \frac{1}{2} \|A v - b \|_2^2 + \frac{1}{ 2 \lambda} \left\|x - v \right\|_2^2\right\} \\ \Longrightarrow 0&= A^T\left( Av - b \right) + \left(-\frac{1}{ \lambda} \left( x - v \right) \right) \\ \Longleftrightarrow 0&= \left(A^TA + \frac{1}{ \lambda}I \right)v - \left(A^Tb + \frac{1}{ \lambda} x \right)\\ \Longleftrightarrow v&= \operatorname{prox}_{\lambda f}\left( x \right) = \left(A^TA + \frac{1}{ \lambda} I \right)^{-1}\left(A^Tb + \frac{1}{ \lambda} x \right). \end{align}
\begin{align} \operatorname{prox}_{\lambda g}\left( z \right) &= \arg\min_{v} \left\{ \lambda \| v\|_1 + \frac{1}{ 2} \left\|z - v \right\|_2^2\right\} \\ &= \arg\min_{ \left\{v_i\right\}} \left\{ \sum_i \lambda|v_i| + \frac{1}{ 2} \sum_i \left\|z_i - v_i \right\|_2^2\right\} \end{align} Since the problem is separable, then you can use KKT conditions to obtain so-called soft thresholding operator. Not to make this post too long, I can refer to you for instance this The Proximal Operator of the $ {L}_{1} $ Norm Function which shows the derivation.
I hope this helps you.