I have programmed a logistic regression using the IRLS algorithm. I would like to apply a LASSO penalization in order to automatically select the right features. At each iteration, the following is solved:

$$\mathbf{\left(X^TWX\right) \delta\hat\beta=X^T\left(y-p\right)}$$

Let $\lambda$ be a non-negative real number. I am not penalizing the intercept as suggested in The Elements of. Statistical Learning. Ditto for the already zero coefficients. Otherwise, I subtract a term from the right-hand side:

$$\mathbf{X^T\left(y-p\right)-\lambda\times \mathrm{sign}\left(\hat\beta\right)}$$

However, I am unsure about the modification of the IRLS algorithm. Is it the right way to do?

Edit: Although I was not confident about it, here is one of the solutions I finally came up with. What is interesting is this solution corresponds to what I now understand about LASSO. There are indeed two steps at each iteration instead of merely one:

  • the first step is the same as before : we make an iteration of the algorithm (as if $\lambda=0$ in the formula for the gradient above),
  • the second step is the new one: we apply a soft-thresholding to each component (except for the component $\beta_0$, which corresponds to the intercept) of the vector $\beta$ obtained at the first step. This is called Iterative Soft-Thresholding Algorithm.

$$\forall i \geq 1, \beta_{i}\leftarrow\mathrm{sign}\left(\beta_{i}\right)\times\max\left(0,\,\left|\beta_{i}\right|-\lambda\right)$$

  • $\begingroup$ Still could not get a better convergence by adapting IRLS. :'( $\endgroup$ – Wok Oct 14 '10 at 13:29

This problem is typically solved by fit by coordinate descent (see here). This method is both safer more efficient numerically, algorithmically easier to implement and applicable to a more general array of models (also including Cox regression). An R implementation is available in the R package glmnet. The codes are open source (partly in and in C, partly in R), so you can use them as blueprints.

| cite | improve this answer | |
  • $\begingroup$ @wok Of note, the scikit.learn package also offers efficient implementation in Python for this kind of stuff. $\endgroup$ – chl Oct 12 '10 at 12:03
  • $\begingroup$ The coordinate descent algorithm is interesting. Thanks. Still thinking about it. $\endgroup$ – Wok Oct 12 '10 at 19:53

The LASSO loss function has a discontinuity at zero along each axis, so IRLS is going to have problems with it. I have found a sequential minimal optimisation (SMO) type approach very effective, see e.g.


a version with MATLAB software is


the software is available here:


The Basic idea is to optimise the coefficients one at at time, and test to see if you cross the discontinuity one coefficient at a time, which is easy as you are perfoming a scalar optimisation. It may sound slow, but it is actually pretty efficient (although I expect better algorithms have been developed since - probably by Keerthi or Chih-Jen Lin who are both leading experts in that sort of thing).

| cite | improve this answer | |
  • $\begingroup$ Thanks. I am reading this and thinking about it. However, this would be a huge modification of the current algorithm. $\endgroup$ – Wok Oct 12 '10 at 19:53

You may check the paper: Efficient L1-regularized logistic regression, which is an IRLS-based algorithm for LASSO. Regarding the implementation, the link may be useful for you (http://ai.stanford.edu/~silee/softwares/irlslars.htm).

| cite | improve this answer | |

The IRLS for the LASSO problem is as following:

$$ \arg \min_{x} \frac{1}{2} \left\| A x - b \right\|_{2}^{2} + \lambda \left\| x \right\|_{1} = \arg \min_{x} \frac{1}{2} \left\| A x - b \right\|_{2}^{2} + \lambda {x}^{T} W {x} $$

Where $ W $ is a diagonal matrix - $ {W}_{i, i} = \frac{1}{ \left| {x}_{i} \right| } $.
This comes from $ \left\| x \right\|_{1} = \sum_{i} \left| {x}_{i} \right| = \sum_{i} \frac{ {x}_{i}^{2} } { \left| {x}_{i} \right| } $.

Now, the above is just Tikhonov Regularization.
Yet, since $ W $ depends on $ x $ one must solve it iteratively (Also this cancels the 2 factor in Tikhonov Regularization, As the derivative of $ {x}^{T} W x $ with regard to $ x $ while holding $ x $ as constant is $ \operatorname{diag} \left( \operatorname{sign} \left( x \right) \right) $ which equals to $ W x $):

$$ {x}^{k + 1} = \left( {A}^{T} A + \lambda {W}^{k} \right)^{-1} {A}^{T} b $$

Where $ {W}_{i, i}^{K} = \frac{1}{ \left| {x}^{k}_{i} \right| } $.

Initialization can be by $ W = I $.

Pay attention this doesn't work well for large values of $ \lambda $ and you better use ADMM or Coordinate Descent.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.