Lasso, fractional norm and structured sparse estimation using a Hadamard product parametrization

Abstract

Using a multiplicative reparametrization, it is shown that a subclass of $L_q$ penalties with $q$ less than or equal to one can be expressed as sums of $L_2$ penalties. It follows that the lasso and other norm-penalized regression estimates may be obtained using a very simple and intuitive alternating ridge regression algorithm. As compared to a similarly intuitive EM algorithm for $L_q$ optimization, the proposed algorithm avoids some numerical instability issues and is also competitive in terms of speed. Furthermore, the proposed algorithm can be extended to accommodate sparse high-dimensional scenarios, generalized linear models, and can be used to create structured sparsity via penalties derived from covariance models for the parameters. Such model-based penalties may be useful for sparse estimation of spatially or temporally structured parameters.

Introduction

Consider estimation for the normal linear regression model $y\sim N_n(X\beta ,{\sigma }^{2}I)$, where $X\in {\mathbb{R}}^{n\times p}$ is a matrix of predictor variables and $\beta \in {\mathbb{R}}^p$ is a vector of regression coefficients to be estimated. A least squares estimate is a minimizer of the residual sum of squares $\Vert y-X\beta {\Vert }^2.$ A popular alternative estimate is the lasso estimate (Tibshirani, 1996), which minimizes $\Vert y-X\beta {\Vert }^2+\lambda \Vert \beta {\Vert }_1$, a penalized residual sum of squares that balances fit to the data against the possibility that some or many of the elements of $\beta$ are small or zero. Indeed, minimizers of this penalized sum of squares may have elements that are exactly zero.

There exist a large variety of optimization algorithms for finding lasso estimates (see Schmidt et al. (2007) for a review). However, the details of many of these algorithms are somewhat opaque to data analysts who are not well-versed in the theory of optimization. One exception is the local quadratic approximation (LQA) algorithm of Fan and Li (2001), which proceeds by iteratively computing a series of ridge regressions. Fan and Li (2001) also suggested using LQA for non-convex $L_q$ penalization when $q<1$, and this technique was used by Kabán and Durrant (2008) and Kabán (2013) in their studies of non-convex $L_q$-penalized logistic regression. However, LQA can be numerically unstable for some combinations of models and penalties. To remedy this, Hunter and Li (2005) suggested optimizing a surrogate “perturbed” objective function. This perturbation must be user-specified, and its value can affect the parameter estimate. As an alternative to using local quadratic approximations, Zou and Li (2008) suggest $L_q$-penalized optimization using local linear approximations (LLA). While this approach avoids the instability of LQA, the algorithm is implemented by iteratively solving a series of $L_1$ penalization problems for which an optimization algorithm must be chosen as well.

This article develops a simple alternative technique for obtaining $L_q$-penalized regression estimates for many values of $q\le 1$. The technique is based on a non-identifiable Hadamard product parametrization (HPP) of $\beta$ as $\beta =u\circ v$, where “$\circ$” denotes the Hadamard (element-wise) product of the vectors $u$ and $v$. As shown in Section 2, if $\hat{u}$ and $\hat{v}$ are optimal $L_2$-penalized values of $u$ and $v$, then $\hat{\beta }=\hat{u}\circ \hat{v}$ is an optimal $L_1$-penalized value of $\beta$. An alternating ridge regression algorithm for obtaining $\hat{u}\circ \hat{v}$ is easy to understand and implement, and is competitive with LQA in terms of speed. Furthermore, a modified version of HPP can be adapted to provide fast convergence in sparse, high-dimensional scenarios. In Section 3 we consider extensions of this algorithm for non-convex $L_q$-penalized regression with $q\le 1$. As in the $L_1$ case, $L_q$-penalized linear regression estimates may be found using alternating ridge regression, whereas estimates in generalized linear models can be obtained with a modified version of an iteratively reweighted least squares algorithm. In Section 4 we show how the HPP can facilitate structured sparsity in parameter estimates: The $L_2$ penalty on the vectors $u$ and $v$ can be interpreted as independent Gaussian prior distributions on the elements of $u$ and $v$. If instead we choose a penalty that mimics a dependent Gaussian prior, then we can achieve structured sparsity among the elements of $\hat{\beta }=\hat{u}\circ \hat{v}$. This technique is illustrated with an analysis of brain imaging data, for which a spatially structured HPP penalty is able to identify spatially contiguous regions of differential brain activity. A discussion follows in Section 5.