TLDR: Why does ADMM diverge when solving $\ell_1$ regression?
Introduction
I am learning about convex optimisation and wanted to solve a simple exercise that I am having issues with. I want to solve a $\ell_1$ regression using ADMM (Alternating Direction Method of Multiplies, https://stanford.edu/~boyd/admm.html). Not because I have to, but because I want to try it and have a little fun.
Problem set up
The problem is:
$$\min_\beta \|X\beta -y \|_1$$
where $X$ has dimensions $N \times D$, $\beta$ is a $D$-dimensional vector and $y$ has length $N$.
What I wanted to try is using an augmented variable $\alpha$ to make the problem linear. That is, I want to solve
$$\min_{\alpha, \beta} \alpha^\top 1_N \quad \text{s.t.} \quad X\beta-y \leq \alpha \quad \text{and} \quad -X\beta + y \leq \alpha $$
these constraints are quivalent to $|X\beta-y| \leq \alpha$.
I now write down the Augmented Lagrangian:
$$\mathcal{L}=\alpha^\top 1_N + \lambda_1 (\alpha-X\beta+y) + \lambda_2 (\alpha+X\beta-y) + \rho_1/2 \|\alpha-X\beta+y\|_2^2 + \rho_2/2 \|\alpha+X\beta-y\|_2^2$$
with $\lambda_1 \leq 0$ and $\lambda_2 \leq 0$ (Constraint should be positive, so during minimisation the $\lambda$s should be negative. See Bishop, Pattern Recognition Appendix).
I initialise: $$\rho_1 \leftarrow 1$$ $$\rho_2 \leftarrow 1$$ $$\lambda_1 \leftarrow 1_N$$ $$\lambda_2 \leftarrow 1_N$$ $$\alpha \leftarrow |X\beta-y|$$
Now I update in the following order:
We solve $\nabla_\beta \mathcal{L}=0$ and use the pseudo inverse (same as OLS estimator) to get
$$\beta \leftarrow \frac{1}{\lambda_1-\lambda_2 + \rho_1 - \rho_2} (X^\top X)^{-1}X^\top \left[\lambda_1 (\alpha + y) + \lambda_2 (\alpha - y) + \rho_1 (\alpha + y) + \rho_2 (\alpha - y) \right]$$
We solve $\nabla_\alpha \mathcal{L}=0$
$$\alpha \leftarrow - \frac{1}{\rho_1 + \rho_2} \left[1_N + \lambda_1 + \lambda_2 +(\rho_2-\rho_1) (X \beta - y) \right]$$
Next the dual ascent on $\lambda_1$:
$$\lambda_1 \leftarrow \max \{\lambda_1 + \rho_1 (\alpha-X\beta+y), 0\}$$
and $\lambda_2$:
$$\lambda_2 \leftarrow \max \{ \lambda_2 + \rho_2 (\alpha+X\beta-y), 0 \}$$
The Issue
$\alpha$ grows to $-\infty$ quite quickly and the loss goes to $\infty$. I can't seem to make it stop. What am I doing wrong? I tried playing with hyperparameters.
Edit: Google Colab Notebook