Feature Selection in Generalized Linear models via the Lasso: To Scale or Not to Scale?

We begin with providing the background notation and facts for the linear regression model and its extensions to generalized linear models. With this requisite notation and background material, we will then be able to explain the main contribution of our paper without unnecessary verbosity.

We denote the $n\times p$ regression matrix (or feature matrix) containing the $p$ features $\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{p}$ as

and the corresponding regression response vector as $\boldsymbol{Y}\in\mathbb{R}^{n}$ (with $\boldsymbol{y}$ being the realization of the random vector $\boldsymbol{Y}$). We make the standard assumptions (see, for example, [6]) that the regression function, or the mean of $\boldsymbol{Y}$ conditional on ${\color[rgb]{0,0,1}\mathbf{X}}$, satisfies $\boldsymbol{Y}=\mathbb{E}_{{\color[rgb]{0,0,1}\mathbf{X}}}[\boldsymbol{Y}]+\boldsymbol{\epsilon}$ for some noise vector $\boldsymbol{\epsilon}$ with conditional expectation $\mathbb{E}_{\color[rgb]{0,0,1}\mathbf{X}}[\boldsymbol{\epsilon}]=\boldsymbol{0}$ and conditional variance $\mathbb{V}\mathrm{ar}_{\color[rgb]{0,0,1}\mathbf{X}}(\boldsymbol{\epsilon})=\sigma^{2}\mathbf{I}_{n}$, where $\sigma$ is an unknown parameter. Recall that in a linear model, we assume that

is a linear function of some model coefficients $\boldsymbol{\beta}\in\mathbb{R}^{p}$ and $\beta_{0}\in\mathbb{R}$, the last one being called the intercept and corresponding to the constant feature $\boldsymbol{1}\in\mathbb{R}^{n}$. Thus, the model for $\boldsymbol{Y}$ is

Define the projection (idempotent) matrix $\mathbf{C}:=\mathbf{I}_{n}-\boldsymbol{1}\boldsymbol{1}^{\top}/n$ and let

be the empirical variance of the components of the $i$-th feature vector. We define the Lasso estimate [7] of $\boldsymbol{\beta}$ as the solution to the penalized least squares:

where $\lambda>0$ is a suitably chosen regularization parameter and the intercept $b_{0}$ is not penalized. It will shortly become clear why the definition (1) of the lasso estimator appears prima facie to be different from the one used in textbooks [5] and statistical packages [2]. We remind the reader that the main advantage of the Lasso regularization is that many of the components of $\hat{\boldsymbol{\beta}}_{\lambda}$ are estimated as zeros, making it very easy to tell which features are important and which are not. This phenomenon is behind the well-known ability of the Lasso estimator to perform simultaneous shrinkage and model selection [5, 7].

Since the true $\boldsymbol{\beta}$ is not known apriori, in this section we explain how to approximately compute (5) using an iterative reweighted least squares (IRLS) method. We begin by reviewing the well-known IRLS for computing the estimator (4) and then explain how it is easily modified to approximately compute our proposed estimator (5). To this end, we introduce the notation $\breve{\boldsymbol{b}}:=[b_{0},\boldsymbol{b}^{\top}]^{\top}$ and $\breve{{\color[rgb]{0,0,1}\mathbf{X}}}:=[\boldsymbol{1},{\color[rgb]{0,0,1}\mathbf{X}}]$, so that $\breve{{\color[rgb]{0,0,1}\mathbf{X}}}\breve{\boldsymbol{b}}=\boldsymbol{1}b_{0}+{\color[rgb]{0,0,1}\mathbf{X}}\boldsymbol{b}.$ Then, computing (4) is equivalent to minimizing

with respect to $\breve{\boldsymbol{b}}$. This problem is nonlinear, but it can be solved by successive and repeated linearizations of $r(\breve{\boldsymbol{b}})$. Suppose that at iteration $t$, we have a current best guess $\breve{\boldsymbol{b}}_{t}$ for the minimizer of (4). Given this $\breve{\boldsymbol{b}}_{t}$, we consider the quadratic multivariate Taylor approximation to the cross-entropy loss at the point $\breve{\boldsymbol{b}}_{t}$:

Then, we update $\breve{\boldsymbol{b}}_{t}$ to $\breve{\boldsymbol{b}}_{t+1}$ by computing the linear Lasso estimator:

This computation is then iterated until convergence [2].

To keep the mathematical detail simple, we henceforth use the Logit model to illustrate the computations that are typically required for all GLMs. In the Logit model the binary response $Y_{1},\ldots,Y_{n}$ are assumed to be independent Bernoulli random variables with conditional mean

yielding the cross-entropy loss:

Then, the gradient $\nabla r(\breve{\boldsymbol{b}}_{t})$ and Hessian $\mathbf{H}(\breve{\boldsymbol{b}}_{t})$ are given by the formulas:

If we define the quantities:

then straightforward algebra shows that the estimator in (6) can also be written in terms of the following (weighed) regularized least-squares:

We can eliminate the intercept term $b_{0}$ from the optimization by applying the centering (projection) matrix

to both ${\color[rgb]{0,0,1}\mathbf{X}}_{t}$ and $\boldsymbol{y}_{t}$. Once the $\boldsymbol{b}_{t+1}$ is computed, then we recover $b_{0,t+1}=\boldsymbol{w}_{t}^{\top}(\boldsymbol{y}_{t}-{\color[rgb]{0,0,1}\mathbf{X}}_{t}\boldsymbol{b}_{t+1})/\|\boldsymbol{w}_{t}\|^{2}$. This gives the following formulas for updating $(b_{0,t},\boldsymbol{b}_{t})$ to $(b_{0,t+1},\boldsymbol{b}_{t+1})$:

We iterate for $t=1,2,\ldots$ until a convergence criterion is met. This iterative reweighed penalized least squares is summarized in the following pseudo-algorithm [2], which assumes that we compute the Lasso estimate $(\hat{\beta}_{0,\lambda},\hat{\boldsymbol{\beta}}_{\lambda})$ on a grid of $m$ values (with $m=200$ being typical values):

where $\lambda_{m}$ is large enough so that $\hat{\boldsymbol{\beta}}_{\lambda_{m}}=\boldsymbol{0}$; see [2].

Line 1 in Algorithm 1 calls the subroutine $\mathrm{lassoCD}(\boldsymbol{b},\boldsymbol{y},{\color[rgb]{0,0,1}\mathbf{X}},\epsilon,\lambda,\boldsymbol{\eta})$ to compute the Lasso estimate:

to within an error tolerance of $\epsilon$ via the method of coordinate descent. For completeness, we include the pseudocode of this subroutine in the Appendix and refer the reader to [2] for more details of this well-known algorithm.

We now describe our proposed method for approximately computing the Lasso estimate (5). The basic idea is to modify the standard linearization (6) by replacing the penalty term $\sum_{i=1}^{p}\eta_{i}\times|b_{i}|$ with $\sum_{i=1}^{p}\eta_{t,i}\times|b_{i}|$, where $\boldsymbol{\eta}_{t}$ is determined from the columns of the Hessian of the cross-entropy loss evaluated at the current estimate $\breve{\boldsymbol{b}}_{t}$. This gives an iterative reweighed least squares in which the scaling in the lasso penalty term is modified from one iteration to the next:

We implement this iterative rescaling with only one minor modification of Algorithm 1. In particular, we replace line 1 in Algorithm 1 with the variances of the features in ${\color[rgb]{0,0,1}\mathbf{X}}_{t}$ weighted by $\sqrt{\mu_{t,i}(1-\mu_{t,i})},\;\forall i$:

In other words, rather than computing the square root of the diagonal elements of the matrix ${\color[rgb]{0,0,1}\mathbf{X}}^{\top}\mathbf{C}{\color[rgb]{0,0,1}\mathbf{X}}/n$ (as currently done in line 1 of Algorithm 1), we instead compute the square root of the diagonal elements of the matrix ${\color[rgb]{0,0,1}\mathbf{X}}_{t}^{\top}\mathbf{C}_{t}{\color[rgb]{0,0,1}\mathbf{X}}_{t}/n$. This is the only significant difference between the implementations of the current widely-used constant scaling and our proposed iteratively rescaled Lasso (IRL) algorithm.

In this section, we provide empirical comparisons between the IRL method and the state-of-the-art Glmnet implementation [2] for Logistic and Poisson regression. We provide numerical results for experiments conducted on both simulated and real data.

In this paper, we proposed a novel scaling of the features for Lasso penalized nonlinear regression. The current widely-used standardization of the feature matrix is performed only once, typically as a preprocessing step, and is independent of the true regression coefficient $\boldsymbol{\beta}$. In contrast, our rescaling proposal uses a standardization which depends on the current best estimate of $\boldsymbol{\beta}$ and is thus updated from one iteration to the next during the course of the nonlinear optimization.

Numerical experiments suggest that the proposed iteratively rescaled Lasso estimator, while never harming the model-selection performance of the estimator in any of the examples considered, has the potential to significantly improve the bias and model selection performance of the estimator when the features matrix is sparse and the signal is strong. In addition, the proposed estimator is just as easy to compute as the vanilla Lasso estimator and only requires a minor modification of one line of the original algorithm.

Future areas of research include providing a theoretical framework of precise conditions under which a significant improvement in the statistical accuracy of the IRL estimator is guaranteed.