Solving the Lasso Problem

Comparing optimization algorithms for Lasso

Authors

Affiliations

University of Washington Autonomous Controls Lab

Published

Oct. 25, 2023

Introduction

Lasso regression is and important regularization technique for linear regression that can also perform variable selection. What this means is that solutions to the lasso problem tend to be sparse (contain zeros) which allows us to rule out certain independent variables in our model. A great resource to familiarize yourself with lasso is this video.

Lasso is of incredible importance in statistics, signal processing, compressed sensing, and image processing. In this post we will look at a variety of optimization technique for solving the lasso problem.

The lasso optimization problem is an unconstrained optimization problem which can be written as follows:

\begin{aligned} \underset{x}{minimize} & \frac{1}{2} ‖ A x - b ‖_{2}^{2} + λ ‖ x ‖_{1} \end{aligned}

where A \in \mathbb{R}^{m \times n}

A first thought to solve this problem might be gradient descent, however the objective function is non-smooth (due to the l1 penalty), so we cannot use gradient descent. A second thought is to use subgradient descent, which is a generalization of gradient descent to non-smooth functions. This would work, but subgradient descent has an extremely slow worst-case convergence rate of \mathcal{O}(1 / \sqrt{t}) (meaning you need four times the iterations to double the accuracy) so we will look at better algorithms.

Reformulating as Quadratic Program

One of the easiest things to do would be reformulating this problem as a Quadratic Program (QP) as follows:

\begin{aligned} \underset{x}{minimize} & \frac{1}{2} ‖ A x - b ‖_{2}^{2} + λ \sum_{i = 1}^{n} t_{i} \\ subject to & - t_{i} \leq x_{i} \leq t_{i} \forall i \in [1, n] \end{aligned}

To see how this was done reference the Mosek Modeling Cookbook

We can then feed this into a QP solver such as OSQP and then get an answer. This works, but it feels wasteful to turn an unconstrained problem into a constrained one and then use a generic QP solver. There should be better algorithms that exploit the structure of our problem where we have a smooth plus a non-smooth term.

ISTA

ISTA or iterative shrinking threshold algorithm is an application of the proximal gradient method to the Lasso problem.

The proximal gradient method solves problems of the following form, where f is differentiable

\begin{aligned} \underset{x}{minimize} & f (x) + g (x) \end{aligned}

The algorithm looks as follows:

x_{k + 1} = {prox}_{η g} (x_{k} - η \nabla f (x_{k}))

where \eta is the gradient descent stepsize for f which will be the inverse of the Lipschitz constant of f.

The proximal operator \boldsymbol{\text{prox}}_{\eta g} is a generalization of the projection operation and is defined as follows

{prox}_{η f} (v) = \underset{x}{argmin} (f (x) + \frac{1}{2 η} ‖ x - v ‖_{2}^{2})

You can think of the proximal operator as returning a point which balances minimizing the function and staying close to the current point. The proximal operator for many function are well known in closed form. For more information on proximal operators and algorithms using proximal operators, refer to .

The idea behind the proximal gradient method is to perform a gradient descent step assuming we are just going to be minimizing the smooth function f and then do an evaluation of the proximal operator for g which can be interpreted as a gradient descent step on the smoothed version of g. More technically we do a gradient descent step on the Moreau envelope of g.

Alternating between these two steps, we eventually minimize our original objective function.

For Lasso we will take,

\begin{aligned} f (x) & = \frac{1}{2} ‖ A x - b ‖_{2}^{2} \\ g (x) & = ‖ x ‖_{1} \end{aligned}

It can be shown that the proximal operator for the l1 norm is the soft threshold operator

{prox}_{η ‖ \cdot ‖_{1}} (v) = S_{η} (v) = sign (v) max (| v | - η, 0)

We can now write out the ISTA iterates as follows

x_{k + 1} = S_{λ / L} (x_{k} - \frac{1}{L} A^{⊤} (A x_{k} - b))

Where L is the maximum eigenvalue of A^TA. This can quickly be computed via power iteration.

It can be shown that this algorithm has a worst-case convergence rate of \mathcal{O}(1 / t) meaning that if we double the number of iterations, we double the accuracy of the solution. This is already better than subgradient method, but is not the best we can do.

FISTA

ISTA was used for a while, but many researchers noticed that it can be painfully slow to converge. In 2009, Beck and Teboulle introduced FISTA (Fast Iterative Shrinking Threshold Algorithm) where they used momentum to accelerate ISTA and were able to achieve the worst-case convergence rate of \mathcal{O}(1 / t^2) , meaning that if we double the number of iterations, we quadruple the accuracy of the solution . FISTA can be thought of as applying ideas from Nesterov’s Accelerated Gradient to ISTA.

This algorithm can be written as follows

\begin{aligned} (1) & x_{k} & = S_{λ / L} (y_{k} - \frac{1}{L} A^{⊤} (A y_{k} - b)) \\ (2) & t_{k + 1} & = \frac{1 + \sqrt{1 + 4 t_{k}^{2}}}{2} \\ (3) & y_{k + 1} & = x_{k} + (\frac{t_{k} - 1}{t_{k + 1}}) (x_{k} - x_{k - 1}) \end{aligned}

ADMM

The final algorithm we will consider for Lasso is the Alternating Direction Method of Multiplers (ADMM). This algorithm was introduced in the mid-1970s, but became popular again after Stephen Boyd et al published their paper in 2011 . This algorithm attempts to solve problems of the following form

\begin{aligned} \underset{x, z}{minimize} & f (x) + g (z) \\ subject to & x = z, \end{aligned}

The iterates of the algorithm are as follows:

\begin{aligned} (4) & x_{k + 1} & = \underset{x}{argmin} L_{ρ} (x, z_{k}, y_{k}) \\ (5) & z_{k + 1} & = \underset{z}{argmin} L_{ρ} (x_{k + 1}, z, y_{k}) \\ (6) & y_{k + 1} & = y_{k} + ρ (x_{k + 1} - z_{k + 1}) \end{aligned}

where

L_{ρ} (x, z, y) = f (x) + g (z) + y^{T} (x - z) + \frac{ρ}{2} ‖ x - z ‖_{2}^{2}

For the case of Lasso, we have

\begin{aligned} (7) & f (x) & = \frac{1}{2} ‖ A x - b ‖_{2}^{2} \\ (8) & g (z) & = ‖ z ‖_{1} \end{aligned}

and the ADMM iterates become

\begin{aligned} (9) & x_{k + 1} & = (A^{T} A + ρ I)^{- 1} (A^{T} b + ρ z_{k} - y_{k}) \\ (10) & z_{k + 1} & = S_{λ / ρ} (x_{k + 1} + y_{k} / ρ) \\ (11) & y_{k + 1} & = y_{k} + ρ (x_{k + 1} - z_{k + 1}) \end{aligned}

Here, \rho > 0 is a stepsize.

Results

Now we will test the QP version of Lasso, ISTA, FISTA, and ADMM to see which is fastest. To generate the data, I generated a random A \in \mathbb{R}^{m \times n} with m < n, then generated a random sparse vector x_{*}, and calculated b=Ax_*.

The stopping criteria for all solver was coming within 0.0001 of the optimal objective function value.

Convergence of ISTA, FISTA, and ADMM with varying stepsizes

Algorithm	Solve Time (sec)
ADMM (rho=50)	0.197
ADMM (rho=100)	0.202
ADMM (rho=10)	1.097
FISTA	1.652
OSQP	2.271
ISTA	8.880

It should be mentioned that the ISTA, FISTA, and ADMM implementations are quite naive and unoptimized, but the OSQP solver is written is pure C. The slowest algorithm by far is ISTA followed by reformulating Lasso as a QP and using OSQP, followed by FISTA, and the fastest algorithm was ADMM. The code to generate the plots can be found here.

References

Proximal Algorithms [link]
Parikh, N. and Boyd, S., 2014. Found. Trends Optim., Vol 1(3), pp. 127–239. Now Publishers Inc. DOI: 10.1561/2400000003
A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems [link]
Beck, A. and Teboulle, M., 2009. SIAM Journal on Imaging Sciences, Vol 2(1), pp. 183-202. DOI: 10.1137/080716542
Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers
Boyd, S., Parikh, N., Chu, E., Peleato, B. and Eckstein, J., 2011. Foundations and Trends\textregistered{} in Machine Learning, Vol 3(1), pp. 1--122. Now Publishers.