Solution Algorithms (Coordinate Descent)

Unlike ridge regression, the Lasso does not have a closed-form solution. The L1 norm introduces non-differentiability at zero, preventing the simple matrix inversion approach. We must resort to iterative optimization algorithms.

This limitation is not a setback—it's an opportunity. The structure of the Lasso problem admits specialized algorithms that are remarkably efficient, exploiting the geometry we've developed. The workhorse algorithm for Lasso is coordinate descent, which leverages:

1. Separability: The L1 penalty $\sum_j |\beta_j|$ separates across coordinates
2. Closed-form updates: Each coordinate subproblem has an explicit solution (soft thresholding)
3. Sparsity exploitation: Zero coefficients require no computation
4. Warm starting: Solutions at nearby $\lambda$ values are close, enabling fast path computation

This page develops coordinate descent in detail, then surveys alternative approaches: proximal gradient methods (ISTA/FISTA) and the path algorithm LARS.

Before developing specialized algorithms, let's understand why standard optimization methods struggle with the Lasso.

Gradient Descent Failure:

Standard gradient descent requires computing $ abla f(\boldsymbol{\beta})$. But the L1 norm $|\boldsymbol{\beta}|_1 = \sum_j |\beta_j|$ has no gradient at $\beta_j = 0$:

$$\frac{d}{d\beta_j}|\beta_j| = \begin{cases} +1 & \beta_j > 0 \ -1 & \beta_j < 0 \ \text{undefined} & \beta_j = 0 \end{cases}$$

We can use subgradients (replacing gradients with elements of the subdifferential), but subgradient methods converge slowly—typically $O(1/\sqrt{k})$ to $O(1/k)$ compared to $O(\rho^k)$ for smooth problems.

Newton's Method Failure:

Newton's method requires the Hessian $ abla^2 f(\boldsymbol{\beta})$. The L1 norm has no second derivative—it's not even first-order differentiable everywhere. Attempting Newton's method on Lasso is fundamentally ill-posed.

Why Special Structure Helps:

The Lasso Objective Structure:

Write the Lasso objective as:

$$F(\boldsymbol{\beta}) = \underbrace{\frac{1}{2n}|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|2^2}{f(\boldsymbol{\beta}): \text{smooth, convex}} + \underbrace{\lambda|\boldsymbol{\beta}|1}{g(\boldsymbol{\beta}): \text{non-smooth, convex, separable}}$$

This "$f + g$" structure is the key. We can apply proximal methods that handle $f$ with gradients and $g$ with its proximal operator. Or we can apply coordinate descent, which exploits separability to obtain closed-form coordinate updates.

Coordinate descent is the method of choice for Lasso in practice. It's simple to implement, surprisingly fast, and scales well to high-dimensional problems.

The Core Idea:

Instead of optimizing all coordinates simultaneously, optimize one coordinate at a time while holding others fixed. Cycle through coordinates repeatedly until convergence.

Why This Works for Lasso:

When all coordinates except $\beta_j$ are fixed, the Lasso objective becomes a univariate optimization problem in $\beta_j$:

$$\min_{\beta_j} \left{ \frac{1}{2n}|\mathbf{r}_{-j} - \beta_j \mathbf{x}_j|_2^2 + \lambda|\beta_j| \right}$$

where $\mathbf{r}{-j} = \mathbf{y} - \sum{k eq j} \beta_k \mathbf{x}_k$ is the partial residual (residual excluding feature $j$).

This is a 1D Lasso problem with a closed-form solution: soft thresholding!

The Soft Thresholding Update:

Let $\tilde{\beta}_j = \frac{1}{n}\mathbf{x}j^T\mathbf{r}{-j}$ be the OLS coefficient for feature $j$ regressed against the partial residual. If features are standardized ($\frac{1}{n}|\mathbf{x}_j|_2^2 = 1$):

$$\beta_j^{\text{new}} = S_\lambda(\tilde{\beta}_j) = \text{sign}(\tilde{\beta}_j) \max(|\tilde{\beta}_j| - \lambda, 0)$$

Let's derive the coordinate descent update rigorously, handling the case of non-standardized features.

Setup:

We minimize over $\beta_j$ while holding $\boldsymbol{\beta}_{-j}$ (all other coefficients) fixed:

$$\min_{\beta_j} \underbrace{\frac{1}{2n}|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|2^2}{\text{quadratic in } \beta_j} + \lambda|\beta_j|$$

Expanding the Quadratic:

Let $\mathbf{r}{-j} = \mathbf{y} - \sum{k eq j} \beta_k \mathbf{x}_k$ be the partial residual. Then:

$$|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|2^2 = |\mathbf{r}{-j} - \beta_j \mathbf{x}_j|2^2 = |\mathbf{r}{-j}|_2^2 - 2\beta_j \mathbf{x}j^T \mathbf{r}{-j} + \beta_j^2 |\mathbf{x}_j|_2^2$$

The Univariate Objective:

$$h(\beta_j) = \frac{1}{2n}\left(|\mathbf{r}_{-j}|_2^2 - 2\beta_j \mathbf{x}j^T \mathbf{r}{-j} + \beta_j^2 |\mathbf{x}_j|_2^2\right) + \lambda|\beta_j|$$

Dropping constants and rearranging:

$$h(\beta_j) = \frac{|\mathbf{x}_j|_2^2}{2n}\left(\beta_j - \frac{\mathbf{x}j^T \mathbf{r}{-j}}{|\mathbf{x}_j|_2^2}\right)^2 + \lambda|\beta_j| + \text{const}$$

Defining Key Quantities:

The Update Formula:

The univariate problem $\min_{\beta_j} \frac{c_j}{2}(\beta_j - \tilde{\beta}_j)^2 + \lambda|\beta_j|$ has solution:

$$\beta_j^{\text{new}} = S_{\lambda/c_j}(\tilde{\beta}_j) = \text{sign}(\tilde{\beta}_j) \max\left(|\tilde{\beta}_j| - \frac{\lambda}{c_j}, 0\right)$$

Simplified Form (Standardized Features):

If $\frac{1}{n}|\mathbf{x}_j|_2^2 = 1$ (standardized), then $c_j = 1$ and:

$$\beta_j^{\text{new}} = S_\lambda\left(\frac{1}{n}\mathbf{x}j^T \mathbf{r}{-j}\right)$$

This is the formula implemented in the algorithm above.

General Form (Non-standardized Features):

$$\beta_j^{\text{new}} = \frac{1}{|\mathbf{x}_j|2^2} S{n\lambda/|\mathbf{x}_j|_2^2}\left(\mathbf{x}j^T \mathbf{r}{-j}\right)$$

Or equivalently:

$$\beta_j^{\text{new}} = S_{n\lambda}\left(\mathbf{x}j^T \mathbf{r}{-j}\right) / |\mathbf{x}_j|_2^2$$

A natural question: Does coordinate descent converge? And if so, how fast?

Convergence Guarantee:

For the Lasso objective (convex with separable non-smooth part), coordinate descent converges to the global optimum. This is guaranteed by:

Theorem (Tseng, 2001): For minimizing $f(\boldsymbol{\beta}) + \sum_j g_j(\beta_j)$ where $f$ is convex, differentiable with Lipschitz gradient, and each $g_j$ is convex (not necessarily smooth), coordinate descent converges to the global minimum.

The Lasso satisfies these conditions with $f(\boldsymbol{\beta}) = \frac{1}{2n}|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|_2^2$ and $g_j(\beta_j) = \lambda|\beta_j|$.

Convergence Rate:

Linear Convergence for Lasso:

Although the Lasso objective is not strongly convex (the L1 norm has zero curvature), coordinate descent often exhibits linear convergence in practice. This is because:

1. Restricted Strong Convexity (RSC): The objective is effectively strongly convex in directions that matter, even if not globally.

2. Finite Active Set: After some iterations, the set of non-zero coefficients stabilizes. On this subspace, the objective is strongly convex (just least squares).

3. Exact Zeros: Once $\beta_j = 0$ and the optimality condition is satisfied, it stays zero—no oscillation.

Practical Convergence:

In practice, coordinate descent for Lasso is remarkably fast:

• Often converges in 10-100 iterations for moderate problems
• Scales well to $p > 10^6$ features
• Benefits enormously from warm starting along the regularization path

Stopping Criteria:

Common convergence criteria:

Naive coordinate descent is already good. Optimized coordinate descent is exceptional. Here are key implementation techniques:

1. Warm Starting:

When computing solutions for multiple $\lambda$ values (the Lasso path), initialize each solve with the solution from the previous $\lambda$:

• Compute $\lambda_{\max} = \frac{1}{n}|\mathbf{X}^T\mathbf{y}|_\infty$ (all-zeros solution)
• Create a grid $\lambda_1 = \lambda_{\max} > \lambda_2 > \cdots > \lambda_K$
• Initialize $\boldsymbol{\beta}^{(k)}$ with the solution at $\lambda_{k-1}$

For dense grids with ratio $\lambda_{k+1}/\lambda_k \approx 0.9$, each solve typically requires only 1-3 iterations.

2. Active Set Strategy:

Focus computation on non-zero coefficients:

3. Safe Screening Rules:

Before even starting optimization, eliminate features that are guaranteed to have zero coefficients:

SAFE Rule: If $|\mathbf{x}_j^T\mathbf{y}| \leq n\lambda - |\mathbf{x}_j|_2 \cdot |\mathbf{y}|_2 \cdot \theta$, then $\hat{\beta}_j = 0$.

Here $\theta$ is a geometric bound depending on the current and target $\lambda$ values.

4. Matrix-Free Updates:

For very large $p$, store $\mathbf{X}$ sparsely or use matrix-free products:

• Compute $\mathbf{x}_j^T\mathbf{r}$ via sparse matrix operations
• Update residual as $\mathbf{r} \leftarrow \mathbf{r} + \Delta\beta_j \cdot \mathbf{x}_j$

5. Strong Rules (Heuristic Screening):

STRONG rule: Assume features with $|\mathbf{x}j^T\mathbf{r}| < 2\lambda - \lambda{\text{prev}}$ will have $\beta_j = 0$.

This is not safe (can miss active features), so verify with KKT conditions after convergence. But it's very effective in practice—often eliminates 99% of features.

Proximal gradient methods offer an alternative to coordinate descent, particularly useful when the smooth part $f$ has structure that coordinate descent can't exploit.

The Proximal Gradient Framework:

For objectives $F(\boldsymbol{\beta}) = f(\boldsymbol{\beta}) + g(\boldsymbol{\beta})$ where $f$ is smooth and $g$ is convex:

$$\boldsymbol{\beta}^{(k+1)} = \text{prox}_{t_k g}\left(\boldsymbol{\beta}^{(k)} - t_k abla f(\boldsymbol{\beta}^{(k)})\right)$$

where $t_k > 0$ is the step size and the proximal operator is:

$$\text{prox}{t g}(\mathbf{z}) = \arg\min{\mathbf{u}} \left{ g(\mathbf{u}) + \frac{1}{2t}|\mathbf{u} - \mathbf{z}|_2^2 \right}$$

For Lasso:

With $g(\boldsymbol{\beta}) = \lambda|\boldsymbol{\beta}|_1$, the proximal operator is component-wise soft thresholding:

$$\text{prox}_{t\lambda|\cdot|1}(\mathbf{z}) = S{t\lambda}(\mathbf{z})$$

ISTA (Iterative Shrinkage-Thresholding Algorithm):

$$\boldsymbol{\beta}^{(k+1)} = S_{t\lambda}\left(\boldsymbol{\beta}^{(k)} + \frac{t}{n}\mathbf{X}^T(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}^{(k)})\right)$$

where $t = 1/L$ with $L = |\mathbf{X}|_2^2/n$ being the Lipschitz constant of $ abla f$.

LARS (Least Angle Regression) is an elegant algorithm that computes the entire Lasso path—all solutions as $\lambda$ varies from $\lambda_{\max}$ to 0—in one shot.

Key Insight: Piecewise Linear Path

The Lasso solution path $\hat{\boldsymbol{\beta}}(\lambda)$ is piecewise linear in $\lambda$. Between certain critical values of $\lambda$, coefficients change at constant rates. This means:

1. We don't need to solve Lasso for many $\lambda$ values separately
2. The path can be traced exactly with finite computation
3. Feature entry/exit points can be computed precisely

LARS Algorithm Overview:

1. Initialize: All coefficients zero, residual $\mathbf{r} = \mathbf{y}$
2. Find most correlated feature: $j^* = \arg\max_j |\mathbf{x}_j^T\mathbf{r}|$
3. Move toward OLS: Increase $\beta_{j^}$ in the direction that decreases $|\mathbf{x}_{j^}^T\mathbf{r}|$
4. Maintain equicorrelation: As we move, keep all active features' correlations equal
5. Add new features: When another feature becomes equally correlated, add it to active set
6. (For Lasso): If a coefficient crosses zero, remove it from active set
7. Repeat until all features are added or residual is zero

The Equicorrelation Property:

At any point on the Lasso path, all active features have the same absolute correlation with the current residual:

$$|\mathbf{x}_j^T\mathbf{r}| = \lambda \quad \text{for all } j \in \mathcal{A}$$

where $\mathcal{A}$ is the active set. This is the equicorrelation property.

LARS maintains this property exactly, moving coefficients together keeping correlations equal. When a new feature's correlation "catches up" to the active correlation, it joins the active set.

LARS vs. Lasso Path:

Plain LARS doesn't remove features—coefficients never go to zero. The LARS-Lasso modification handles sign changes:

• If a coefficient crosses zero during the path (sign change), remove it from active set
• This is the "Lasso modification" that produces the actual Lasso path

Computational Complexity:

• LARS-Lasso path computation: $O(p \cdot n \cdot \min(n, p))$
• Same cost as solving one OLS problem with $\min(n, p)$ features
• For sparse solutions, often much faster in practice