A lava attack on the recovery of sums of dense and sparse signals

A lava attack on the recovery of sums of dense and sparse signals

\address{Department of Economics, MIT, Cambridge, MA 02139} \email{[email removed]}

\address{Booth School of Business, University of Chicago, Chicago, IL 60637} \email{[email removed]}

\address{Department of Mathematics, University of Maryland, College Park, MD 20741} \email{[email removed]}

Key words: high-dimensional models, penalization, shrinkage, non-sparse signal recovery

Introduction

Many recently proposed high-dimensional modeling techniques build upon the fundamental assumption of sparsity. Under sparsity, we can approximate a high-dimensional signal or parameter by a sparse vector that has a relatively small number of non-zero components. Various $\ell_1$-based penalization methods, such as the lasso and soft-thresholding, have been proposed for signal recovery, prediction, and parameter estimation within a sparse signal framwork. See frank1993statistical, donoho1995adapting, tibshirani96, fan2001variable, efron2004least, zou2005regularization, zhao2006model, yuan2006model, bunea2007sparsity, candes2007dantzig, fan2008sure, Bickeletal, meinshausen2009lasso, wainwright2009sharp, bunea2010spades, zhang2010nearly, loh2013regularized, and others. By virtue of being based on $\ell_1$-penalized optimization problems, these methods produce sparse solutions in which many estimated model parameters are set exactly to zero.

Another commonly used shrinkage method is ridge estimation. Ridge estimation differs from the aforementioned $\ell_1$-penalized approaches in that it does not produce a sparse solution but instead provides a solution in which all model parameters are estimated to be non-zero. Ridge estimation is thus particularly suitable when the model's parameters or unknown signals contain many very small components, i.e. when the model is dense. See, e.g., hsu2014random. Ridge estimation tends to work better than sparse methods whenever a signal is dense in such a way that it can not be well-approximated by a sparse signal.

In practice, we may face environments that have signals or parameters which are neither dense nor sparse. The main results of this paper provide a model that is appropriate for this environment and a corresponding estimation method with good estimation and prediction properties. Specifically, we consider models where the signal or parameter, $\theta$, is given by the superposition of sparse and dense signals:

Here, $\delta$ is a sparse vector that has a relatively small number of large entries, and $\beta$ is a dense vector having possibly very many small, non-zero entries. Traditional sparse estimation methods, such as lasso, and traditional dense estimation methods, such as ridge, are tailor-made to handle respectively sparse signals and dense signals. However, the model for $\theta$ given above is “sparse+dense" and cannot be well-approximated by either a “dense only" or “sparse only" model. Thus, traditional methods designed for either sparse or dense settings are not optimal within the present context.

Motivated by this signal structure, we propose a new estimation method, called “lava". Let $\ell(\mathrm{data}, \theta)$ be a general statistical loss function that depends on unknown parameter $\theta$, and let $p$ be the dimension of $\theta$. To estimate $\theta$, we propose the “lava" estimator given by

where $\widehat\beta $ and $ \widehat\delta$ solve the following penalized optimization problem:

In the formulation of the problem, $\lambda_2$ and $\lambda_1$ are tuning parameters corresponding to the $\ell_2$- and $\ell_1$- penalties which are respectively applied to the dense part of the parameter, $\beta$, and the sparse part of the parameter, $\delta$. The resulting estimator is then the sum of a dense and a sparse estimator. Note that the separate identification of $\beta$ and $\delta$ is not required in ((ref)), and the lava estimator is designed to automatically recover the combination $\widehat\beta+\widehat\delta$ that leads to the optimal prediction of $\beta + \delta$. Moreover, under standard conditions for $\ell_1$-optimization, the lava solution exists and is unique. In naming the proposal “lava", we emphasize that it is able, or at least aims, to capture or wipe out both sparse and dense signals.

The lava estimator admits the lasso and ridge shrinkage methods as two extreme cases by respectively setting either $\lambda_2=\infty$ or $\lambda_1=\infty$.\footnote{With $\lambda_1 = \infty$ or $\lambda_2 = \infty$, we set $\lambda_1\|\delta\|_1 = 0$ when $\delta = 0$ or $\lambda_2\|\beta\|_2^2 = 0$ when $\beta = 0$ so the problem is well-defined.} In fact, it continuously connects the two shrinkage functions in a way that guarantees it will never produce a sparse solution when $\lambda_2<\infty$. Of course, sparsity is not a requirement for making good predictions. By construction, lava's prediction risk is less than or equal to the prediction risk of the lasso and ridge methods with oracle choice of penalty levels for ridge, lasso, and lava; see Figure (ref). Lava also tends to perform no worse than, and often performs significantly better than, ridge or lasso with penalty levels chosen by cross-validation; see Figures (ref) and (ref).

Note that our proposal is rather different from the elastic net method, which also uses a combination of $\ell_1$ and $\ell_2$ penalization. The elastic net penalty function is $\theta \mapsto \lambda_2 \|\theta\|^2_2 + \lambda_1 \|\theta\|_1$, and thus the elastic net also includes lasso and ridge as extreme cases corresponding to $\lambda_2 = 0$ and $\lambda_1 = 0$ respectively. In sharp contrast to the lava method, the elastic net does not split $\theta$ into a sparse and a dense part and will produce a sparse solution as long as $\lambda_1 > 0$. Consequently, the elastic net method can be thought of as a sparsity-based method with additional shrinkage by ridge. The elastic net processes data very differently from lava (see Figure 2 below) and consequently has very different prediction risk behavior (see Figure 1 below).

We also consider the post-lava estimator which refits the sparse part of the model:

where $ \widetilde \delta$ solves the following penalized optimization problem:

This estimator removes the shrinkage bias induced by using the $\ell_1$ penalty in estimation of the sparse part of the signal. Removing this bias sometimes results in further improvements of lava's risk properties.

We provide several theoretical and computational results about lava in this paper. First, we provide analytic expressions for the finite-sample risk function of the lava estimator as well as for other methods in the Gaussian sequence model and in a fixed design regression model with Gaussian errors. Within this context, we exhibit “sparse+dense" examples where lava significantly outperforms both lasso and ridge. Stein's unbiased risk estimation plays a central role in our theoretical analysis, and we thus derive Stein's unbiased risk estimator (SURE) for lava. We also characterize lava's “Efron's” degrees of freedom (efron2004estimation). Second, we give deviation bounds for the prediction risk of the lava estimator in regression models akin to those derived by Bickeletal for lasso. Third, we illustrate lava's performance relative to lasso, ridge, and elastic net through simulation experiments using penalty levels chosen via either minimizing the SURE or by k-fold cross-validation for all estimators. In our simulations, lava outperforms lasso and ridge in terms of prediction error over a wide range of regression models with coefficients that vary from having a rather sparse structure to having a very dense structure. When the model is very sparse, lava performs as well as lasso and outperforms ridge substantially. As the model becomes more dense in the sense of having the size of the “many small coefficients" increase, lava outperforms lasso and performs just as well as ridge. This is consistent with our theoretical results.

We conclude the introduction by noting that our proposed approach complements other recent approaches to structured sparsity problems such as those considered in fused sparsity estimation (tibshirani2005sparsity and dalalyan2012fused) and structured matrix estimation problems (candes2011robust, chandrasekaran2011rank, POET, and klopp2014robust). The latter line of research studied estimation of matrices that can be written as low rank plus sparse matrices. Our new results are related to but are sharply different from this latter line of work since our focus is on regression problems. Specifically, our chief objects of interest are regression coefficients along with the associated regression function and predictions of the outcome variable. Thus, the target statistical applications of our developed methods include prediction, classification, curve-fitting, and supervised learning. Another noteworthy point is that it is impossible to recover the “dense" and “sparse" components separately within our framework; instead, we recover the sum of the two components. By contrast, it is possible to recover the low-rank component of the matrix separately from the sparse part in some of the structured matrix estimation problems. This distinction serves to highlight the difference between structured matrix estimation problems and the framework discussed in this paper. Due to these differences, the mathematical side of our analysis needs to address a completely different set of issues than are addressed in the aforementioned structured matrix estimation problems.

We organize the rest of the paper as follows. Section 2 defines the lava shrinkage estimator in a canonical Gaussian sequence model, and derives its theoretical risk function. Section 3 defines and analyzes the lava estimator in the regression model. Section 4 provides computational examples, and Section 5 concludes. We give all proofs in the appendix.

Notation. The notation $a_n \lesssim b_n$ means that $a_n \leq C b_n$ for all $n$, for some constant $C$ that does not depend on $n$. The $\ell_{2}$ and $\ell_{1}$ norms are denoted by $\|\cdot\|_2$ (or simply $\|\cdot\|$) and $\| \cdot \|_{1}$, respectively. The $\ell_{0}$-“norm", $\|\cdot\|_0$, denotes the number of non-zero components of a vector, and the $\|.\|_{\infty}$ norm denotes a vector's maximum absolute element. When applied to a matrix, $\|\cdot\|$ denotes the operator norm. We use the notation $a \vee b = \max( a, b)$ and $a \wedge b = \min(a , b)$. We use $x'$ to denote the transpose of a column vector $x$.

The lava estimator in a canonical model

The one dimensional case

Consider the simple problem where a scalar random variable is given by $$Z = \theta+ \epsilon,\quad \epsilon \sim N(0,\sigma^2).$$ We observe a realization $z$ of $Z$ and wish to estimate $\theta$. Estimation will often involve the use of regularization or shrinkage via penalization to process input $z$ into output ${d}(z)$, where the map $z \mapsto {d}(z)$ is commonly referred to as the shrinkage (or decision) function. A generic shrinkage estimator then takes the form $\widehat \theta = {d}(Z)$.

The commonly used lasso method uses $\ell_1$-penalization and gives rise to the lasso or soft-thresholding shrinkage function: $$ {d}_{\mathrm{lasso}}(z)=\ \arg\min_{\theta \in \mathbb{R}} \Big \{ (z-\theta)^2+\lambda_l|\theta| \Big \} = (|z|-\lambda_l/2)_+\text{sign}(z), $$ where $y_+ := \max(y,0)$ and $\lambda_l \geq 0$ is a penalty level. The use of the $\ell_2$-penalty in place of the $\ell_1$ penalty yields the ridge shrinkage function: $$ {d}_{\mathrm{ridge}}(z)=\ \arg\min_{\theta \in \mathbb{R}} \ \Big \{ (z-\theta)^2+\lambda_r |\theta|^2 \Big\} =\frac{z}{1+\lambda_r}, $$ where $\lambda_r \geq 0$ is a penalty level. The lasso and ridge estimators then take the form $$\widehat \theta_{\mathrm{lasso}} = {d}_{\mathrm{lasso}} (Z), \quad \widehat \theta_{\mathrm{ridge}} = {d}_{\mathrm{ridge}} (Z).$$ Other commonly used shrinkage methods include the elastic-net (zou2005regularization), which uses $\theta \mapsto \lambda_2|\theta|^2+\lambda_1|\theta|$ as the penalty function; hard-thresholding; and the SCAD (fan2001variable), which uses a non-concave penalty function.

Motivated by points made in the introduction, we proceed differently. We decompose the signal into two components $$\theta=\beta+\delta,$$ and use the different penalty functions -- the $\ell_2$ and $\ell_1$ -- for each component in order to predict $\theta$ better. We thus consider the penalty function $$(\beta, \delta) \mapsto \lambda_2|\beta|^2+\lambda_1|\delta|,$$ and introduce the “lava" shrinkage function $z \mapsto {d}_{\mathrm{lava}}(z)$ defined by

where ${d}_1(z)$ and ${d}_2(z)$ solve the following penalized prediction problem:

Although the decomposition $\theta=\beta+\delta$ is not unique, the optimization problem ((ref)) has a unique solution for any given $(\lambda_1, \lambda_2)$. The proposal thus defines the lava estimator of $\theta$: $$\widehat\theta_{\mathrm{lava}} = {d}_{\mathrm{lava}}(Z).$$

For large signals such that $|z|>\lambda_1/(2k)$, lava has the same bias as the lasso. This bias can be removed through the use of the post-lava estimator$$ \widehat\theta_{\mathrm{post}-\mathrm{lava}} = {d}_{\mathrm{post}-\mathrm{lava}}(Z),$$ where ${d}_{\mathrm{post}-\mathrm{lava}}(z) := {d}_2(z) + \tilde {d}_1(z),$ and $\tilde {d}_1(z)$ solves the following penalized prediction problem:

The removal of this bias will result in improved risk performance relative to the original estimator in some contexts.

From the Karush-Kuhn-Tucker conditions, we obtain the explicit solution to ((ref)).

Figure (ref) plots the lava shrinkage function along with various alternative shrinkage functions for $z>0$. The top panel of the figure compares lava shrinkage to ridge, lasso, and elastic net shrinkage. It is clear from the figure that lava shrinkage is different from lasso, ridge, and elastic net shrinkage. The figure also illustrates how lava provides a bridge between lasso and ridge, with the lava shrinkage function coinciding with the ridge shrinkage function for small values of the input $z$ and coinciding with the lasso shrinkage function for larger values of the input. Specifically, we see that the lava shrinkage function is a combination of lasso and ridge shrinkage that corresponds to using whichever of the lasso or ridge shrinkage is closer to the 45 degree line.

It is also useful to consider how lava and post-lava compare with the post-lasso or hard-thresholding shrinkage: ${d}_{\text{post-lasso}}(z)=z1\{|z|>\lambda_l/2\}.$ These different shrinkage functions are illustrated in the bottom panel of Figure (ref).

From ((ref)), we observe some key characteristics of the lava shrinkage function:

1) The lava shrinkage admits the lasso and ridge shrinkages as two extreme cases. The lava and lasso shrinkage functions are the same when $\lambda_2=\infty$, and the ridge and lava shrinkage functions coincide if $\lambda_1 = \infty$.

2) The lava shrinkage function ${d}_{\mathrm{lava}}(z)$ is a weighted average of data $z$ and the lasso shrinkage function ${d}_{\mathrm{lasso}}(z)$ with weights given by $1-k$ and $k$.

3) The lava never produces a sparse solution when $\lambda_2<\infty$: If $\lambda_2<\infty$, ${d}_{\mathrm{lava}}(z)=0$ if and only if $z=0$. This behavior is strongly different from elastic net which produces a sparse solution as long as $\lambda_1 > 0$.

4) The lava shrinkage function continuously connects the ridge shrinkage function and the lasso shrinkage function. When $|z|<\lambda_1/(2k)$, lava shrinkage is equal to ridge shrinkage; and when $|z|>\lambda_{1}/(2k)$, lava shrinkage is equal to lasso shrinkage.

5) The lava shrinkage does exactly the opposite of the elastic net shrinkage. The elastic net shrinkage function coincides with the lasso shrinkage function when $|z|<\lambda_1/(2k)$; and when $|z|>\lambda_{1}/(2k)$, the elastic net shrinkage is the same as ridge shrinkage.

The risk function of the lava estimator in the one dimensional case

In the one-dimensional case with $Z\sim N(\theta, \sigma^2)$, a natural measure of the risk of a given estimator $\widehat\theta = {d}(Z)$ is given by

Let ${\mathrm{P}}_{\theta,\sigma}$ denote the probability law of $Z$. Let $\phi_{\theta,\sigma}$ be the density function of $Z$. We provide the risk functions of lava and post-lava in the following theorem. We also present the risk functions of ridge, elastic net, lasso, and post-lasso for comparison.

These results for the one-dimensional case provide a key building block for results in the multidimensional case provided below. In particular, we build from these results to show that the lava estimator performs very favorably relative to, and can substantially dominate, the maximum likelihood estimator, the ridge estimator, and $\ell_1$-based estimators (such as lasso and elastic-net) in interesting multidimensional settings.

Multidimensional case

We consider now the canonical Gaussian model or the Gaussian sequence model. In this case, we have that $$Z\sim N_p(\theta, \sigma^2I_p)$$ is a single observation from a multivariate normal distribution where $\theta=(\theta_1,...,\theta_p)'$ is a $p$-dimensional vector. A fundamental result for this model is that the maximum likelihood estimator $Z$ is inadmissible and can be dominated by the ridge estimator and related shrinkage procedures when $p\geq 3$ (e.g. stein1956).

In this model, the lava estimator is given by $$\widehat{\theta}_{\text{lava}} := (\widehat \theta_{\mathrm{lava},1},..., \widehat \theta_{\mathrm{lava}, p})':=({d}_{\text{lava}}(Z_1),...,{d}_{\text{lava}}(Z_p))',$$ where ${d}_{\text{lava}}(z)$ is the lava shrinkage function as in ((ref)). The estimator is designed to capture the case where $$ \theta=\underbrace{\beta}_{\text{dense part}}+\underbrace{\delta}_{\text{sparse part}} $$ is formed by combining a sparse vector $\delta$ that has a relatively small number of non-zero entries which are all large in magnitude and a dense vector $\beta$ that may contain very many small non-zero entries. This model for $\theta$ is “sparse+dense.” It includes cases that are not well-approximated by “sparse” models - models in which a very small number of parameters are large and the rest are zero - or by “dense” models - models in which very many coefficients are non-zero but all coefficients are of similar magnitude. This structure thus includes cases that pose challenges for estimators such as the lasso and elastic net that are designed for sparse models and for estimators such as ridge that are designed for dense models.

All of the shrinkage estimators discussed in Section 2.1 generalize to the multidimensional case in the same way as lava. Let $z \mapsto {d}_e(z)$ be the shrinkage function associated with estimator $e$ in the one dimensional setting where $e$ can take values in the set $$ \mathcal{E} = \{ \mathrm{lava}, \mathrm{post\text{-}lava}, \mathrm{ridge}, \mathrm{lasso}, \mathrm{post\text{-}lasso}, \mathrm{elastic \ net} \}.$$ We then have a similar estimator in the multidimensional case given by $$\widehat{\theta}_e:= (\widehat \theta_{e,1},..., \widehat \theta_{e, p})':=({d}_{e }(Z_1),...,{d}_{ e}(Z_p))'.$$ The risk calculation from the one dimensional case then caries over to the multidimensional case since $$ \mathrm{R}(\theta, \widehat \theta_e) := {\mathrm{E}} \| \theta - \widehat \theta_e\|_2^2 = \sum_{j=1}^p R(\theta_j, \widehat \theta_{e,j}). $$ Given this fact, we immediately obtain the following result.

These risk functions are illustrated in Figure 1 in a prototypical “sparse+dense" model generated according to the model discussed in detail in Section 2.5. The tuning parameters used in this figure are the best possible (risk minimizing or oracle) choices of the penalty levels found by minimizing the risk expression given in Theorem (ref).

Canonical plug-in choice of penalty levels

We now discuss simple, rule-of-thumb choices for the penalty levels for lasso ($\lambda_l$), ridge ($\lambda_r$) and lava ($\lambda_1, \lambda_2$). In the Gaussian model, a canonical choice of $\lambda_l$ is $$ \lambda_l=2\sigma{\Phi^{-1}(1-c/(2p))}, $$ which satisfies $${\mathrm{P}} \left (\max_{j\leq p}|Z_j-\theta_j|\leq\lambda_l/2 \right)\geq 1-c;$$ see, e.g., donoho1995adapting. Here $\Phi(\cdot)$ denotes the standard normal cumulative distribution function, and $c$ is a pre-determined significance level which is often set to 0.05. The risk function for ridge is simple, and an analytic solution to the risk minimizing choice of ridge tuning parameter is given by $$ \lambda_r=\sigma^2 (p/\|\theta\|_2^2). $$

As for the tuning parameters for lava, recall that the lava estimator in the Gaussian model is

If the dense component $\beta$ were known, then following donoho1995adapting would suggest setting $$ \lambda_1=2\sigma{\Phi^{-1}(1-c/(2p))} $$ as a canonical choice of $\lambda_1$ for estimating $\delta$. If the sparse component $\delta$ were known, one could adopt $$ \lambda_2=\sigma^2(p/\|\beta\|_2^2) $$ as a choice of $\lambda_2$ for estimating $\beta$ following the logic for the standard ridge estimator.

We refer to these choices as the “canonical plug-in" tuning parameters and use them in constructing the risk comparisons in the following subsection. We note that the lasso choice is motivated by a sparse model and does not naturally adapt to or make use of the true structure of $\theta$. The ridge penalty choice is explicitly tied to risk minimization and relies on using knowledge of the true $\theta$. The lava choices for the parameters on the $\ell_1$ and $\ell_2$ penalties are, as noted immediately above, motivated by the respective choices in lasso and ridge. As such, the motivations and feasibility of these canonical choices are not identical across methods, and the risk comparisons in the following subsection should be interpreted within this light.

Some risk comparisons in a canonical Gaussian model

To compare the risk functions of lava, lasso, and ridge estimators, we consider a canonical Gaussian model, where $$\theta_1=3 , \quad \theta_j=0.1q,\quad j=2,...,p,$$ for some $q\geq 0$. We set the noise level to be $\sigma^2=0.1^2$. The parameter $\theta$ can be decomposed as $\theta=\beta+\delta$, where the sparse component is $\delta=(3,0,...,0)'$, and the dense component is $$\beta=(0, 0.1q,...,0.1q)',$$ where $q$ describes the “size of small coefficients." The canonical tuning parameters are $\lambda_l=\lambda_1=2\sigma{\Phi^{-1}(1-c/(2p))}$, $\lambda_r=\sigma^2p/(3+0.1^2q^2(p-1))$ and $\lambda_2=\sigma^2p/(0.1^2q^2(p-1))$.

Figure (ref) (given in the introduction) compares risks of lava, lasso, ridge, elastic net, and the maximum likelihood estimators as functions of the size of the small coefficients $q$, using the ideal (risk minimizing or oracle choices) of the penalty levels. Figure (ref) compares risks of lava, lasso, ridge and the maximum likelihood estimators using the “canonical plug-in" penalty levels discussed above. Theoretical risks are plotted as a function of the size of the small coefficients $q$. We see from these figures that regardless of how we choose the penalty levels -- ideally or via the plug-in rules -- lava strictly dominates the competing methods in this “sparse+dense" model. Compared to lasso, the proposed lava estimator does about as well as lasso when the signal is sparse and does significantly better than lasso when the signal is non-sparse. Compared to ridge, the lava estimator does about as well as ridge when the signal is dense and does significantly better than ridge when the signal is sparse.

In Section 5 we further explore the use of feasible, data-driven choices of penalty levels via cross-validation and SURE minimization; see Figures (ref) and 5. We do so in the context of the Gaussian regression model with fixed regressors. With either cross-validation or SURE minmization, the ranking of the estimators remains unchanged, with lava consistently dominating lasso, ridge, and the elastic net.

stein1956 proved that a ridge estimator strictly dominates maximum likelihood in the Gaussian sequence model once $p \geq 3$. In the comparisons above, we also see that the lava estimator strictly dominates the maximum likelihood estimator; and one wonders whether this domination has a theoretical underpinning similar to Stein's result for ridge. The following result provides some (partial) support for this phenomenon for the lava estimator with the plug-in penalty levels. The result shows that, for a sufficiently large $n$ and $p$, lava does indeed uniformly dominate the maximum likelihood estimator on the compact set $\{\theta=\beta+\delta: \|\beta\|_{\infty}+\|\delta\|_{\infty}<M\}$.

Stein's unbiased risk estimation for lava

stein1981estimation proposed a useful risk estimate based on the integration by parts formula, now commonly referred to as Stein's unbiased risk estimate (SURE). This subsection derives SURE for the lava shrinkage in the multivariate Gaussian model.

Note that

An essential component to understanding the risk is given by applying Stein's formula to calculate $ {\mathrm{E}}[(Z-\theta)'\widehat\theta_{\text{lava}}]$. A closed-form expression for this expression in the one-dimensional case is given in equation ((ref)) in the appendix. The following result provides the SURE for lava in the more general multidimensional case.

Lava in the Regression Model

Definition of Lava in the Regression Model

Consider a fixed design regression model: $$ Y=X\theta_0+U,\quad U\sim N(0,\sigma_u^2I_n), $$ where $Y=(y_1,...,y_n)'$, $X=(X_1,...,X_n)'$, and $\theta_0$ is the true regression coefficient. Following the previous discussion, we assume that $$\theta_0=\beta_0+\delta_0$$ is “sparse+dense” with sparse component $\delta_0$ and dense component $\beta_0$. Again, this coefficient structure includes cases which cannot be well-approximated by traditional sparse models or traditional dense models and will pose challenges for estimation strategies tailored to sparse settings, such as lasso and similar methods, or strategies tailored to dense settings, such as ridge.

In order to define the estimator we shall rely on the normalization condition that

Note that without this normalization, the penalty terms below would have to be modified in order to insure equivariance of the estimator to changes of scale in the columns of $X$.

The lava estimator $\widehat \theta_{\mathrm{lava}}$ of $\theta_0$ solves the following optimization problem:

The lava program splits parameter $\theta$ into the sum of $\beta$ and $\delta$ and penalizes these two parts using the $\ell_2$ and $\ell_1$ penalties. Thus, the $\ell_1$- penalization regularizes the estimator of the sparse part $\delta_0$ of $\theta_0$ and produces a sparse solution $\widehat\delta$. The $\ell_2$-penalization regularizes the estimator of the dense part $\beta_0$ of $\theta_0$ and produces a dense solution $\widehat \beta$. The resulting estimator of $\theta_0$ is then simply the sum of the sparse estimator $\widehat\delta$ and the dense estimator $\widehat\beta$.

A Key Profile Characterization and Some Insights.

The lava estimator can be computed in the following way. For a fixed $\delta$, we minimize $$ \widehat\beta(\delta)=\arg\min_{\beta\in\mathbb{R}^p} \left \{ \frac{1}{n}\|Y-X(\beta+\delta)\|_2^2+\lambda_2\|\beta\|_2^2 \right\}, $$ with respect to $\beta$. This program is simply the well-known ridge regression problem, and the solution is $$ \widehat\beta(\delta)=(X'X+n\lambda_2I_p)^{-1}X'(Y-X\delta). $$ By substituting $\beta=\widehat\beta(\delta)$ into the objective function, we then define an $\ell_1$-penalized quadratic program which we can solve for $\widehat\delta$:

The lava solution is then given by $\widehat\theta=\widehat\beta(\widehat\delta)+\widehat\delta$. The following result provides a useful characterization of the solution.

The theorem shows that solving for the sparse part $\widehat\delta$ of the lava estimator is equivalent to solving for the parameter in a standard lasso problem using transformed data. This result is key to both computation and our theoretical analysis of the estimator.

Degrees of Freedom and SURE

Degrees of freedom is often used to quantify model complexity and to construct adaptive model selection criteria for selecting tuning parameters. In a Gaussian linear regression model $Y\sim N(X\theta_0, \sigma_u^2I_n)$ with a fixed design, we can define the degrees of freedom of the mean fit $X\widehat\theta$ to be $$ \mathrm{df}(\widehat\theta)=\frac{1}{\sigma_u^2}{\mathrm{E}}[(Y-X\theta_0)'X\widehat\theta]; $$ see, e.g., efron2004estimation. Note that this quantity is also an important component of the mean squared prediction risk: $$ {\mathrm{E}}\frac{1}{n}\|X\widehat\theta-X\theta_0\|_2^2=-\sigma_u^2+{\mathrm{E}}\frac{1}{n}\|X\widehat\theta-Y\|_2^2+ \frac{2\sigma_u^2}{n} \mathrm{df}(\widehat\theta). $$

stein1981estimation's SURE theory provides a tractable way of deriving an unbiased estimator of the degrees of freedom, and thus the mean squared prediction risk. Specifically, write $\widehat\theta=d(Y,X)$ as a function of $Y$, conditional on $X$. Suppose $d(\cdot, X):\mathbb{R}^n\to\mathbb{R}^p$ is almost differentiable; see meyer2000degrees and efron2004least). For $f: \mathbb{R}^n \to \mathbb{R}^n$ differentiable at $y$, define $$ \partial_y f(y) := [ \partial f_{ij} (y) ], \quad (i, j) \in\{1,...,n\}^2, \quad \partial f_{ij} (y) : = \frac{\partial}{\partial y_j} f_i(y), $$ $$ \nabla_y\cdot f(y):=\text{tr} (\partial_y f(y)). $$ Let $X_i'$ denote the $i$-th row of $X$, $i=1,...,n$. Then, from stein1981estimation, we have that $$ \frac{1}{\sigma_u^2}{\mathrm{E}}[(Y-X\theta_0)'Xd(Y,X)]={\mathrm{E}}[ \nabla_y\cdot (Xd(Y,X))] = \text{tr} \big (\partial_y[Xd(Y,X)] \big ). $$ An unbiased estimator of the term on the right-hand-side of the display may then be constructed using its sample analog.

In this subsection, we derive the degrees of freedom of the lava, and thus a SURE of its mean squared prediction risk. By Theorem (ref),

where $d_{\mathrm{lava}}(y,X)$ is the lava estimator on the data $(y,X)$ and $d_{\mathrm{lasso}}(\operatorname{\mathsf{K}}_{\lambda_2}^{1/2}y, \operatorname{\mathsf{K}}_{\lambda_2}^{1/2}X)) $ is the lasso estimator on the data $ (\operatorname{\mathsf{K}}_{\lambda_2}^{1/2}y, \operatorname{\mathsf{K}}_{\lambda_2}^{1/2}X)$ with the penalty level $\lambda_1$. The almost differentiability of the map $y \mapsto d_{\mathrm{lasso}}(\operatorname{\mathsf{K}}_{\lambda_2}^{1/2}y, \operatorname{\mathsf{K}}_{\lambda_2}^{1/2}X)$ follows from the almost differentiability of the map $u \mapsto d_{\mathrm{lasso}}(u, \operatorname{\mathsf{K}}_{\lambda_2}^{1/2}X)$, which holds by the results in dossal2011degrees and tibshirani2012degrees.

The following theorem presents the degrees of freedom and SURE for lava. Let $\widehat J=\{j\leq p: \widehat\delta_j\neq 0\}$ be the active set of the sparse component estimator with cardinality denoted by $|\widehat J|$. Recall that $\widetilde X=\operatorname{\mathsf{K}}_{\lambda_2}^{1/2}X$. Let $\widetilde X_{\hat J}$ be an $n\times |\widehat J|$ submatrix of $\widetilde X$ whose columns are those corresponding to the entries in $\widehat J$. Let $A^-$ denote the Moore-Penrose pseudo-inverse of a square matrix $A$.

Post-lava in regression

We can also remove the shrinkage bias in the sparse component introduced by the $\ell_1$-penalization via a post-selection procedure. Specifically, let ($\widehat\beta,\widehat\delta$) respectively denote the lava estimator of the dense and sparse components. Define the post-lava estimator as follows:

Let $X_{\hat J}$ be an $n\times |\widehat J|$ submatrix of $X$ whose columns are selected by $\widehat J$. Then we can partition $\widetilde\delta=(\widetilde\delta_{\hat J}, 0)'$, where $\widetilde \delta_{\hat J}= ( X_{\hat J} ' X_{\hat J} )^- X_{\hat J} '(Y-X\widehat\beta)$. Write $\operatorname{\mathsf{P}}_{\hat J}=X_{\hat J}(X_{\hat J}'X_{\hat J})^-X_{\hat J}'$ and $\operatorname{\mathsf{K}}_{\hat J}=I_n-\operatorname{\mathsf{P}}_{\hat J}$. The post-lava prediction for $X\theta$ is: $$ X\widehat \theta_{\text{post-lava}}= \operatorname{\mathsf{P}}_{\hat J}Y+\operatorname{\mathsf{K}}_{\hat J}X\widehat\beta. $$ In addition, note that the lava estimator satisfies $X\widehat\beta= \operatorname{\mathsf{P}}_{\lambda_2}(Y- X \widehat\delta)$. We then have the following expression of $ X\widehat \theta_{\text{post-lava}}$.

The above lemma reveals that the post-lava corrects the $\ell_1$-shrinkage bias of the original lava fit by adding the projection of the lava residual onto the subspace of the selected regressors. This correction is in the same spirit as the post-lasso correction for shrinkage bias in the standard lasso problem; see belloni2013least.

Deviation Bounds for Prediction Errors

In the following, we develop deviation bounds for the lava prediction error: $\frac{1}{n}\|X\widehat\theta_{\text{lava}}-X\theta_0\|_2^2$. We continue to work with the decomposition $\theta_0=\beta_0+\delta_0$ and will show that lava performs well in terms of rates on the prediction error in this setting. According to the discussion in Section 3.2, there are three sources of prediction error: (i) $ \textsf{D}_{\text{ridge}}(\lambda_2)$, (ii) $ \widetilde X\beta_0$ and (iii) $\operatorname{\mathsf{K}}_{\lambda_2}X(\widehat\delta-\delta_0)$. The behavior of the first two terms is determined by the behavior of the ridge estimator of the dense component $\beta_0$, and the behavior of the third term is determined by the behavior of the lasso estimator on the transformed data.

We assume that $U\sim N(0,\sigma_u^2I_n)$ and that $X$ is fixed. As in the lasso analysis of Bickeletal, a key quantity is the maximal norm of the score: $$ \Lambda= \left \|\frac{2}{n}\widetilde X'\widetilde U \right \|_{\infty}= \left \|\frac{2}{n}X'\operatorname{\mathsf{K}}_{\lambda_2}U \right \|_{\infty}. $$ Following belloni2013least, we set the penalty level for the lasso part of lava in our theoretical development as

and $c>1$ a constant. Note that belloni2013least suggest setting $c = 1.1$ and that $\Lambda_{1-\alpha}$ is easy to approximate by simulation.

Let $S=X'X/n$ and $\bar V_{\lambda_2}$ be the maximum diagonal element of $$V_{\lambda_2}=(S+\lambda_2I_p)^{-1}S(S+\lambda_2I_p)^{-1}\lambda_2^2.$$ Then by the union bound and Mill's inequality:

Thus the choice $\Lambda_{1-\alpha}$ is strictly sharper than the union bound-based, classical choice $\bar \Lambda_{1-\alpha} $. Indeed, $\Lambda_{1-\alpha}$ is strictly smaller than $\bar \Lambda_{1-\alpha}$ even in orthogonal design cases since union bounds are not sharp. In collinear or highly-correlated designs, it is easy to give examples where $\Lambda_{1-\alpha} = o(\bar \Lambda_{1-\alpha})$; see BCW-AOS2014. Thus, the gains from using the more refined choice can be substantial.

We define the following design impact factor: For $\widetilde X = \operatorname{\mathsf{K}}_{\lambda_2}^{1/2} X$, $$\iota (c, \delta_0, \lambda_1,\lambda_2) := \inf_{\Delta \in \mathcal{R}(c, \delta_0, \lambda_1,\lambda_2) } \frac{\|\widetilde X \Delta \|_2/\sqrt{n}}{\| \delta_0\|_1 - \|\delta_0 + \Delta \|_1 + c^{-1} \| \Delta \|_1}, $$ where $\mathcal{R}(c, \delta_0, \lambda_1,\lambda_2) = \{ \Delta \in \mathbb{R}^p\setminus\{0\}: \|\widetilde X \Delta \|^2_2/n \leq 2 \lambda_1 (\| \delta_0\|_1 - \| \delta_0 + \Delta \|_1 + c^{-1} \|\Delta\|_1)\}$ is the restricted set, and where $\iota (c,\delta_0, \lambda_1,\lambda_2):= \infty$ if $\delta_0 = 0$.

The design impact factor generalizes the restricted eigenvalues of Bickeletal and and is tailored for bounding estimation errors in the prediction norm (cf. BCW-AOS2014). Note that in the best case, when the design is well-behaved and $\lambda_2$ is a constant, we have that

where $\kappa>0$ is a constant. Remarks given below provide further discussion.

The following theorem provides the deviation bounds for the lava prediction error.

Simulation Study

The lava and post-lava algorithm can be summarized as follows.

We present a Monte-Carlo analysis based on a Gaussian linear regression model: $Y=X\theta+U$, $U \mid X \sim N(0, I_n)$. The parameter $\theta$ is a $p$-vector defined as $$\theta = (3, 0, \ldots, 0)' + q (0, .1, \ldots, .1)',$$ where $q\geq 0$ denotes the “size of small coefficients". When $q$ is zero or small, $\theta$ can be well-approximated by the sparse vector $(3, 0, \ldots, 0)$. When $q$ is relatively large, $\theta$ cannot be approximated well by a sparse vector. We set $n=100$ and $p=2n$, and compare the performance of $X\widehat\theta$ formed from one of five methods: lasso, ridge, elastic net, lava, and post-lava.\footnote{Results with $p=n/2$, where OLS is also included, are available in supplementary material. The results are qualitatively similar to those given here, with lava and post-lava dominating all other procedures.} The rows of $X$ are generated independently from a mean zero multivariate normal with covariance matrix $\Sigma$. We present results under an independent design, $\Sigma=I$, and a factor covariance structure with $\Sigma=LL'+I$ where the rows of $L$ are independently generated from $N(0,I_3)$. In the latter case, the columns of $X$ depend on three common factors. We focus on a fixed design study, so the design $X$ is generated once and fixed throughout the replications.

To measure performance, we consider the risk measure $\mathrm{R}(\theta, \widehat \theta) = {\mathrm{E}} [\frac{1}{n}\| X\widehat\theta-X\theta\|_2^2]$ where the expectation ${\mathrm{E}}$ is conditioned on $X$. For each estimation procedure, we report the simulation estimate of this risk measure formed by averaging over $B=100$ simulation replications. Figures (ref) and 5 plot the simulation estimate of $\mathrm{R}(\theta, \widehat \theta)$ for each estimation method as a function of $q$, the size of the “small coefficients”. In Figure (ref), all the tuning parameters are chosen via minimizing the SURE as defined in Theorem 3.2; and the tuning parameters are chosen by 5-fold cross-validation in Figure 5. The SURE formula depends on the error variance $\sigma_u^2$, which must be estimated. A conservative preliminary estimator for $\sigma_u^2$ can be obtained from an iterative method based on the regular lasso estimator; see, e.g., belloni2013least. On the other hand, $k$-fold cross-validation does not require a preliminary variance estimator.

The comparisons are similar in both figures with lava and post-lava dominating the other procedures. It is particularly interesting to compare the performance of lava to lasso and ridge. The lava and post-lava estimators perform about as well as lasso when the signal is sparse and perform significantly better than lasso when the signal is non-sparse. The lava and post-lava estimators perform about as well as ridge when the signal is dense and perform much better than ridge when the signal is sparse. When the tuning parameters are selected via cross-validation, the post-lava performs slightly better than the lava when the model is sparse. The gain is somewhat more apparent in the independent design. Additional simulations are presented in our supplementary material (lava).

Discussion

We propose a new method, called “lava", which is designed specifically to achieve good prediction and estimation performance in “sparse+dense” models. In such models, the high-dimensional parameter is represented as the sum of a sparse vector with a few large non-zero entries and a dense vector with many small entries. This structure renders traditional sparse or dense estimation methods, such as lasso or ridge, sub-optimal for prediction and other estimation purposes. The proposed approach thus complements other approaches to structured sparsity problems such as those considered in fused sparsity estimation (tibshirani2005sparsity and dalalyan2012fused) and structured matrix decomposition problems (candes2011robust, chandrasekaran2011rank, POET, and klopp2014robust).

There are a number of interesting research directions that remain to be considered. An immediate extension of the present results would be to consider semi-pivotal estimators akin to the root-lasso/scaled-lasso of BCW-Biometrika and sun2012scaled. For instance, we can define

Thanks to the characterization of Theorem (ref), the method can be implemented by applying root-lasso on appropriately transformed data. The present work could also be extended to accommodate non-Gaussian settings and settings with random designs, and it could also be extended beyond the mean regression problem to more general M- and Z- estimation problems, e.g., along the lines of negahban2009unified.