Scaled Sparse Linear Regression

Introduction

This paper concerns the simultaneous estimation of the regression coefficients and noise level in a high-dimensional linear model. High-dimensional data analysis is a topic of great current interest due to the growth of applications where the number of unknowns far exceeds the number of data points. Among statistical models arising from such applications, linear regression is one of the best understood. Penalization, convex minimization and thresholding methods have been proposed, tested with real and simulated data, and proved to control errors in prediction, estimation and variable selection under various sets of regularity conditions. These methods typically require an appropriate penalty or threshold level. A larger penalty level may lead to a simple model with large bias, while a smaller penalty level may lead to a complex noisy model due to overfitting. Scale-invariance considerations and existing theory suggest that the penalty level should be proportional to the noise level of the regression model. In the absence of knowledge of the latter level, cross-validation is commonly used to determine the former. However, cross-validation is computationally costly and theoretically poorly understood, especially for the purpose of variable selection and the estimation of regression coefficients. The penalty level selected by cross-validation is called the prediction-oracle in Meinshausen & Bühlmann 2006, who gave an example to show that the prediction-oracle solution does not lead to consistent model selection for the lasso.

Estimation of the noise level in high-dimensional regression is interesting in its own right. Examples include quality control in manufacturing and risk management in finance.

An iterative algorithm

where ${\beta}=(\beta_{1},\ldots,\beta_{p})^{\prime}$ is a vector of regression coefficients. Let the penalty $\rho(t)$ be standardized to ${\dot{\rho}}(0+)=1$, where ${\dot{\rho}}(t)=(d/dt)\rho(t)$. A vector ${\widehat{\beta}}=(\widehat{\beta}_{1},\ldots,\widehat{\beta}_{p})^{\prime}$ is a critical point of the penalized loss (1) if and only if

If the penalized loss (1) is convex in ${\beta}$, then (2) is the Karush–Kuhn–Tucker condition for its minimization.

Given a penalty function $\rho(\cdot)$, one still has to choose a penalty level $\lambda$ to arrive at a solution of (2). Such a choice may depend on the purpose of estimation, since variable selection may require a larger $\lambda$ than does prediction. However, scale-invariance considerations and theoretical results suggest using a penalty level proportional to the noise level $\sigma$. This motivates a scaled penalized least squares estimator as a numerical equilibrium in the following iterative algorithm:

The first step of our implementation is the computation of a solution path $\widehat{\beta}(\lambda)$ of (2) beginning from $\widehat{\beta}(\lambda)=0$ for $\lambda=|X^{\prime}y/n|_{\infty}$. For quadratic spline penalties $\rho(t)$ with $m$ knots, Zhang 2010 developed an algorithm to compute a linear spline path of solutions $\{\lambda^{(t)}\oplus{\widehat{\beta}}^{(t)}:t\geq 0\}$ of (2) to cover the entire range of $\lambda$. This extends the least angle regression solution or lasso path (Osborne et al. 2000a; Osborne et al. 2000b; Efron et al. 2004) from $m=1$ and includes the minimax concave penalty for $m=2$ and the smoothly clipped absolute deviation penalty (Fan & Li 2001) for $m=3$. An R package named plus is available for computing the solution paths for these penalties.

The second step of our implementation is the iteration (3) along the solution path $\beta(\lambda)$ computed in the first step. That is to use the already computed

in (3). For the scaled lasso, we use $a=0$ in (3) and $\rho(t)=t$ in (1) and (2). For the scaled minimax concave penalized selection, we use $a=0$ and the minimax concave penalty $\rho(t)=\int_{0}^{t}(1-x/\gamma)_{+}{\rm d}x$, where $\gamma>0$ regularizes the maximum concavity of the penalty. When $\gamma=\infty$, it becomes the scaled lasso. The algorithm (3) can be easily implemented once a solution path is computed.

Antoniadis 2010 suggested this jointly convex loss function as a way of extending Huber’s robust regression method to high dimensions. For $a=0$ and $\lambda=\widehat{\sigma}\lambda_{0}$ with fixed $\widehat{\sigma}$, $\widehat{\sigma}L_{\lambda_{0}}({\beta},\widehat{\sigma})=L_{\lambda}({\beta})+\widehat{\sigma}^{2}/2$, so that ${\widehat{\beta}}\leftarrow{\widehat{\beta}}(\lambda)$ in (4) minimizes $L_{\lambda_{0}}({\beta},\widehat{\sigma})$ over ${\beta}$. For fixed ${\widehat{\beta}}$, $\widehat{\sigma}^{2}\leftarrow{|}{y}-{X}{\widehat{\beta}}{|}_{2}^{2}/\{(1-a)n\}$ in (3) minimizes $L_{\lambda_{0}}({\widehat{\beta}},\sigma)$ over $\sigma$. During the revision of this paper, we learned that She & Owen 2011 have considered penalizing Huber’s concomitant loss function for outlier detection in linear regression. We summarize some properties of the algorithm (3) with (4) in the following proposition.

Let ${\widehat{\beta}}={\widehat{\beta}}(\lambda)$ be a solution path of (2) with $\rho(t)=t$. The penalized loss function (5) is jointly convex in $({\beta},\sigma)$ and the algorithm (3) with (4) converges to

The resulting estimators ${\widehat{\beta}}={\widehat{\beta}}({X},{y})$ and $\widehat{\sigma}=\widehat{\sigma}({X},{y})$ are scale equivariant in ${y}$ in the sense that ${\widehat{\beta}}({X},c{y})=c{\widehat{\beta}}({X},{y})$ and $\widehat{\sigma}({X},c{y})=|c|\widehat{\sigma}({X},{y})$. Moreover,

Since (5) is not strictly convex, the joint estimator may not be unique for some data $(X,y)$. However, since (5) is strictly convex in $\sigma$, $\widehat{\sigma}$ is always unique in (6) and the uniqueness of $\widehat{\beta}$ follows from that of the lasso estimator $\widehat{\beta}(\lambda)$ at $\lambda=\widehat{\sigma}\lambda_{0}$; ${\widehat{\beta}}(\lambda)$ is unique when the second part of (2) is strict in the sense of not hitting $\pm\lambda$ when $\widehat{\beta}_{j}=0$, which holds almost everywhere in $(X,y)$ for $\lambda>0$. See, for example, Zhang 2010.

Let $\widehat{\sigma}(\lambda)=|y-X\widehat{\beta}(\lambda)|_{2}/\{(1-a)n\}^{1/2}$. For $\lambda_{0}=\{(2/n)\log p\}^{1/2}$, (7) implies that

While the present paper continues our earlier work (Sun & Zhang 2010) by providing further theoretical and numerical justifications for (3) and (4), the estimator has appeared in different forms. In addition to (5) and (6) of Antoniadis 2010, (8) appeared in Zhang 2010. While this paper was in revision, a reviewer called our attention to Belloni et al. 2011, who focused on studying $\widehat{\beta}$ in an equivalent form as square-root lasso. We note that (3) and (4) allow concave penalties and degrees of freedom adjustments as in Zhang 2010.

Theoretical results

Let ${\beta}^{*}$ be a vector of true regression coefficients. An expert with oracular knowledge of ${\beta}^{*}$ would estimate the noise level by the oracle estimator

where $\kappa(\xi,T)$, the compatibility factor (van de Geer & Bühlmann 2009), is defined as

with the cone ${\mathscr{C}}(\xi,T)=\{{u}:{|}{u}_{T^{c}}{|}_{1}\leq\xi{|}{u}_{T}{|}_{1}\}$. Since the prediction error bound $\eta(\lambda,\xi,w,T)$ is valid for all $w$ and $T$, $\tau_{0}$ is related to its minimum over all $w$ and $T$ at the oracle scale $\sigma^{*}$:

In particular, if $\lambda_{0}=A\{(2/n)\log p\}^{1/2}$ with $A>(\xi+1)/(\xi-1)$ and $\eta_{*}(\sigma\lambda_{0},\xi)/\sigma\to 0$, then

Theorem 3.1 extends to the scaled lasso a unification of prediction oracle inequalities for a fixed penalty. With $\lambda=\sigma^{*}\lambda_{0}/(1-\tau_{0})_{+}$, (13) gives $\max\{(\sigma^{*}\tau_{0})^{2},{|}{X}{\widehat{\beta}}-{X}{\beta}^{*}{|}_{2}^{2}/n\}\leq\eta_{*}(\lambda,\xi)$, or

For fixed penalty $\lambda$, the upper bound $\eta(\lambda,\xi,w,T)$ has been previously established for different $w$ and $T$, with possibly different constant factors. Examples include $\eta(\lambda,\xi,\beta^{*},\emptyset)=2\lambda|\beta^{*}|_{1}$ (Greenshtein & Ritov 2004; Greenshtein 2006), $\eta(\lambda,\xi,\beta^{*},S_{\beta^{*}})\lesssim\lambda^{2}|\beta^{*}|_{0}$ with $S_{w}=\{j:w_{j}\neq 0\}$ (van de Geer & Bühlmann 2009), and $\min_{w}\eta(\lambda,\xi,w,S_{w})=\min_{w}\{|X\beta^{*}-Xw|_{2}^{2}/n+O(\lambda^{2}|w|_{0})\}$ (Koltchinskii et al. 2011). In (10), the coefficient for $|Xw-X\beta^{*}|_{2}^{2}/n$ is 1 as in Koltchinskii et al. 2011.

if ${\widetilde{C}}\geq(1+\xi)\max\{2,1/\kappa^{2}(2\xi+1,{\widetilde{T}})\}$ with ${\widetilde{T}}=\{j:|\beta_{j}^{*}|>\lambda\}$. This allows $\beta^{*}$ to have many small elements, as in Zhang & Huang 2008, Zhang 2009 and Ye & Zhang 2010. The bound $\mu(\lambda,\xi)\leq(\xi+1)\lambda|S_{\beta^{*}}|/\{2\kappa^{2}(\xi,S_{\beta^{*}})\}$ improves upon its earlier version in van de Geer & Bühlmann 2009 by a constant factor $4\xi/(\xi+1)\in(2,4)$.

Let ${\widehat{\beta}},\widehat{\sigma},{\beta}^{*},\sigma^{*},z^{*}$ and $\xi$ be as in Theorem 3.1. Set $\tau_{*}=\{\lambda_{0}\mu(\sigma^{*}\lambda_{0},\xi)/\sigma^{*}\}^{1/2}$. (i) The following inequalities hold when $z^{*}\leq(1-\tau_{*}^{2})\lambda_{0}(\xi-1)/(\xi+1)$,

(ii) Let $\lambda_{0}\geq\{(2/n)\log(p/\epsilon)\}^{1/2}(\xi+1)/\{(\xi-1)(1-\tau_{*}^{2})\}$. For all $\epsilon>0$ and $n-2>\log(p/\epsilon)\to\infty$,

If $\lambda_{0}=A\{(2/n)\log p\}^{1/2}$ with $A>(\xi+1)/(\xi-1)$ and $\lambda_{0}\mu(\sigma\lambda_{0},\xi)/\sigma\ll n^{-1/2}$, then

Since $\sigma^{2}\tau_{*}^{2}\approx\mu(\lambda,\xi)\leq 2(\xi+1)\min_{T}\eta(\lambda,2\xi+1,{\beta}^{*},T)$ with $\lambda=\sigma\lambda_{0}$, the rate $\tau_{*}^{2}$ in (18) is essentially the square of that in (13), in view of (12). It follows that the scaled lasso provides a faster convergence rate than does the penalized maximum likelihood estimator for the estimation of the noise level (Städler et al. 2010; Sun & Zhang 2010). In particular, (18) implies that

with $S_{\beta^{*}}=\{j:\beta^{*}_{j}\neq 0\}$, when $\kappa^{2}(\xi,S_{\beta^{*}})$ can be treated as a constant. The bounds in (20) and its general version (18) lead to the asymptotic normality (19) under proper assumptions. Thus, statistical inference about $\sigma$ is justified with the scaled lasso in certain large-$p$-smaller-$n$ cases, for example, when $|\beta^{*}|_{0}(\log p)/\surd n\to 0$ under the compatibility condition (van de Geer & Bühlmann 2009).

The quantities in (22) are used in the uniform uncertainty principle (Candes & Tao 2007) and the sparse Riesz condition (Zhang & Huang 2008). We note that $1-\delta_{a}^{-}$ is the minimum eigenvalue of $X_{A}^{\prime}X_{A}/n$ among $|A|\leq a$, $1+\delta_{a}^{+}$ is the corresponding maximum eigenvalue, and $\theta_{a,b}$ is the maximum operator norm of size $a\times b$ off-diagonal sub-blocks of the Gram matrix $X^{\prime}X/n$.

Suppose ${|}{\beta}^{*}_{S^{c}}{|}_{1}=0$. Then, Theorem 3.2 holds with $\mu(\lambda,\xi)$ replaced by $\lambda|S|(2\xi)/\{(\xi+1)F_{1}(\xi,S)\}$, and for $z^{*}\leq(1-\tau_{*}^{2})\lambda_{0}(\xi-1)/(\xi+1)$,

for all $1\leq q\leq\infty$, where $k=|S|$. In particular, for $\xi=\surd 2$ and $z^{*}\leq(1-\tau_{*}^{2})\lambda_{0}(\surd{2}-1)^{2}$,

The proofs of Theorems 1 and 2 are based on a basic inequality

as a consequence of the Karush–Kuhn–Tucker conditions (2). The version of (25) with ${w}={\beta}^{*}$ is well-known (van de Geer & Bühlmann 2009) and controls $|X\widehat{\beta}(\lambda)-X\beta^{*}|_{2}^{2}$ for sparse $\beta^{*}$. When $|X\widehat{\beta}(\lambda)-X\beta^{*}|_{2}^{2}>|Xw-X\beta^{*}|_{2}^{2}$, (25) controls the excess for sparse $w$ by the same argument. The general ${w}$ is taken in Theorem 1, while ${w}={\beta}^{*}$ is taken in Theorem 2. In both cases, (25) provides the cone condition in (11) and (21). This is used to derive upper and lower bounds for (7), the derivative of the profile loss function $L_{\lambda_{0}}({\widehat{\beta}}(\sigma\lambda_{0}),\sigma)$ with respect to $\sigma$, within a neighborhood of $\sigma/\sigma^{*}=1$. The bounds for the minimizer $\widehat{\sigma}$ then follow from the joint convexity of the penalized loss (5).

2 Estimation after model selection

We have proved that without requiring the knowledge of $\sigma$, the scaled lasso enjoys prediction and estimation properties comparable to the best known theoretical results for known $\sigma$, and the scaled lasso estimate of $\sigma$ enjoys consistency and asymptotic normality properties under proper conditions. However, the lasso estimator may have substantial bias (Fan & Peng 2004; Zhang 2010), and its bias is significant in our own simulation experiments. Although the smoothly clipped absolute deviation and minimax concave penalized selectors were introduced to remove the bias of the lasso (Fan & Li 2001; Zhang 2010), a theoretical study of their scaled version (3) is beyond the scope of this paper. In this subsection, we present theoretical results for another bias removing method: least squares estimation after model selection.

Given an estimator $\widehat{\beta}$ of the coefficient vector $\beta$, the least squares estimator of $\beta$ and the corresponding estimator of the noise level $\sigma$ in the model selected by $\widehat{\beta}$ are

where $\hbox{\rm supp}(\beta)=\{j:\beta_{j}\neq 0\}$. Alternatively, we may use ${\overline{\sigma}}=\big|y-X{\overline{\beta}}\,\big|_{2}\big/\surd(n-|\widehat{\beta}|_{0})$ to estimate the noise level. However, since the effect of this degrees of freedom adjustment is of smaller order than our error bound, we will focus on the simpler (27).

In addition to the compatibility factor $\kappa(\xi,S)$ in (11), we use sparse eigenvalues to study the least squares estimation after the scaled lasso selection. Let $\lambda_{\min}(M)$ be the smallest eigenvalue of a matrix $M$ and $\lambda_{max}(M)$ the largest . For models $T\subset\{1,\ldots,p\}$, define

as the sparse lower eigenvalue of the Gram matrix for models containing $T$ and the sparse upper eigenvalue for models disjoint with $T$. Let $S=\hbox{\rm supp}(\beta^{*})$ and ${\widehat{S}}=\hbox{\rm supp}(\widehat{\beta})$. The following theorem provides prediction and estimation error bounds for (27) after the scaled lasso selection, along with an upper bound for the false positive $|{\widehat{S}}\setminus S|$, a key element in our study.

Let $({\widehat{\beta}},\widehat{\sigma})$ be the scaled lasso estimator in (6) and $({\overline{\beta}},{\overline{\sigma}})$ the least squares estimator (27) in model ${\widehat{S}}$. Let $\beta^{*},\sigma^{*},z^{*},\xi$ and $\tau_{*}$ be as in Theorem 3.2 and $m$ be an integer satisfying $|S|\xi^{2}/\kappa^{2}(\xi,S)<m/\kappa_{+}(m,S)$. If $z^{*}\leq(1-\tau_{*}^{2})\lambda_{0}(\xi-1)/(\xi+1)$, then

with ${\widehat{\lambda}}=\widehat{\sigma}\lambda_{0}\leq\sigma^{*}\lambda_{0}/(1-\tau_{*}^{2})$ and $\sigma^{*}_{m,S}=\max_{B\supseteq S,|B\setminus S|=m}|(y-X\beta^{*})_{B}|_{2}/\surd n$, and

Moreover, in addition to the probability bound for $z^{*}\leq(1-\tau_{*}^{2})\lambda_{0}(\xi-1)/(\xi+1)$ in Theorem 3.2 (ii), for all integers $1\leq m\leq p$,

For Gaussian design matrices, the sparse eigenvalues $\kappa_{-}(m,\emptyset)$ and $\kappa_{+}(m,\emptyset)$ can be treated as constants when $m(\log p)/n$ is small and the eigenvalues of the expected Gram matrix are uniformly bounded away from zero and infinity (Zhang & Huang 2008). Since $\kappa_{-}(m,S)\geq\kappa_{-}(m+|S|,\emptyset)$ and $\kappa_{+}(m,S)\leq\kappa_{+}(m,\emptyset)$, they can be treated as constants in the same sense in Theorem 3.4. Thus, for sufficiently small $|S|(\log p)/n$, we may take an $m$ of the same order as $|S|$. In this case, the difference between $({\overline{\beta}},{\overline{\sigma}})$ and the scaled lasso estimator $(\widehat{\beta},\widehat{\sigma})$ is of no greater order than the difference between $(\widehat{\beta},\widehat{\sigma})$ and the estimation target $(\beta^{*},\sigma^{*})$. Consequently,

As we have mentioned earlier, the key element in our analysis of (27) is the bound $|{\widehat{S}}\setminus S|<m$ in (28). Since this is a weaker assertion than variable consistency ${\widehat{S}}=S$, the conditions of Theorem 3.4 on the design matrix is of a weaker form than the irrepresentability condition for variable selection consistency (Meinshausen & Bühlmann 2006; Zhao & Yu 2006). In Zhang & Huang 2008 and Zhang 2010, upper bounds for the false positive were obtained under a sparse Riesz condition on $\kappa_{-}(m,\emptyset)$ and $\kappa_{+}(m,\emptyset)$.

Numerical results

The top section of Table 3.2 presents our simulation results, while the bottom section includes the simulation results of Fan et al. 2012 for several joint estimators of $({\beta},\sigma)$ using cross-validation, without repeating their experiment. In addition to the bias and the standard error of the ratios $\widehat{\sigma}/\sigma$ for the five original estimators and ${\overline{\sigma}}/\sigma$ for the least squares estimation after model selection, we report the average model size $|\widehat{S}|$ and the relative frequency of sure screening, $\widehat{S}\supseteq S$, as in Fan et al. 2012, where $\widehat{S}=\hbox{\rm supp}(\widehat{\beta})$ is the selected model.

Without post processing, the scaled minimax concave penalized selector with the universal penalty level $\lambda_{2}=\surd\{(2/n)\log p\}$ clearly outperforms other procedures in this example. However, the results of the least squares estimation after model selection at penalty level $\lambda_{2}$ are nearly identical to the top performer for all five methods. In view of the results in average model size and sure screening proportion, the success of post processing at $\lambda_{2}$ is clearly due to the success of model selection. The five methods select too few variables at the larger penalty level $\lambda_{3}$ and too many at the smaller $\lambda_{1}$, both leading to substantial bias in the estimation of $\sigma$ for $r_{0}=0$. For $r_{0}=0.5$, selecting a slightly smaller model does not harm so much since a substantial portion of the effect of the missing variables is explained by the selected variables correlated to them. The minimax concave penalized selector is nearly unbiased in this example, so that it does not need post processing. Cross-validation methods select about 30 variables when the true model size is 3. This over selection is probably the reason for the large bias for most cross-validation methods and large standard error for all of them.

This experiment has the same setting as in the simulation study in Sun & Zhang 2010, where the scaled lasso and the scaled minimax concave penalized selection are called the naive estimators. We provide the description of the simulation settings in Sun & Zhang 2010 in our notation as follows: $(n,p)=(600,3000)$, the ${x}_{j}$ are normalized columns from a Gaussian random matrix with independent and identically distributed rows and correlation ${r_{0}}^{|k-j|}$ between the $j$-th and $k$-th entries within each row, $\gamma=2/(1-\max|{x}_{k}^{\prime}{x}_{j}|/n)$ for the minimax concave penalty and smoothly clipped absolute deviation penalty, the nonzero $\beta^{*}_{j}$ are composed of five blocks of $\beta_{*}(1,2,3,4,3,2,1)^{\prime}$ centered at random multiples $j_{1},\dots,j_{5}$ of 25, $\beta_{*}$ sets ${|}{X}{\beta}^{*}{|}_{2}^{2}=3n$, and ${y}-{X}{\beta}^{*}$ is a vector of independent and identically distributed $N(0,1)$ variables. Thus, the true noise level is $\sigma=1$. We set ${r_{0}}=0.1$ for low correlation between design vectors and ${r_{0}}=0.9$ for high correlation.

2 Real data example

We study a data set containing 18976 probes for 120 rats, which is reported in Scheetz et al. 2006. Our goal is to find probes that are related to that of gene TRIM32, which has been found to cause Bardet–Biedl syndrome, a genetically heterogeneous disease of multiple organ systems including the retina. We consider linear regression with the probe from TRIM32, 1389163_at, as the response variable. As in Huang et al. 2008, we focus on 3000 probes with the largest variances among the 18975 covariates and consider two approaches. The first approach is to regress on these $p=3000$ probes. The second approach is to regress on the 200 probes among the 3000 with the largest marginal correlation coefficients with TRIM32. For the cross-validation lasso, we randomly partition the data 1000 times, each with a training set of size 80 and a validation set of size 40. For each partition, the penalty level $\lambda$ is selected by minimizing the prediction mean squared error in the validation set. Then we compute the lasso estimator with all 120 observations at the penalty level equal to the median of the selected penalty levels with the 1000 random partitions. Since cross-validation tends to choose a larger model, we also consider an adjusted version using the cross-validated error of the least squares estimator after the lasso selection. For the minimax concave penalty, we set $\gamma=2/(1-\sigma_{0.95})=6.37$, where $\sigma_{0.95}$ is the 95% quantile of $|{x}^{\prime}_{k}{x}_{j}|/n$.

Table 4.2 shows the probe sets identified by four methods: the cross-validation lasso, its adjusted version, the scaled lasso at at universal penalty level $\lambda_{2}=\{2(\log p)/n\}^{1/2}$, and the minimax concave penalized selection at the same penalty level. We apply stability selection (Meinshausen & Buhlmann 2010) to check the reliability of selection. Let $W_{1},\dots,W_{p}$ be independent variables with $P(W=0.2)=P(W=1)=1/2$ and

where ${\widehat{\lambda}}$ is the penalty level chosen by individual methods. Stability selection selects variables with nonzero estimated $\widehat{\beta}^{W}_{j}$ over 50 times in 100 replications. We observe that the scaled minimax concave penalized selector produces most sparse and most stable selection, followed by the adjusted cross-validation, the scaled lasso and then the plain cross-validation. The selection results are consistent among the four methods in the sense that the selected models are almost nested. Since the model size is between 6 and 8 by stability selection in all 8 cases and by the scaled minimax concave penalized selection for both $p=200$ and $p=3000$, these two methods provide most consistent results. The scaled lasso and the adjusted cross-validation yield identical lasso and stability selections for $p=200$ and identical stability selection for $p=3000$.

We also compare the prediction performance of the scaled lasso with that of the lasso with the best fixed penalty level. We compute the scaled estimators in 1000 replications. In each replication, the dataset is split at random into a training set with 80 observations and a test set with 40 observations. The prediction mean squared error is computed within the test set, while the scaled estimators and the lasso estimator with fixed penalty level $\lambda$ are computed based on the training set. Figure 1 demonstrates that in prediction, the scaled lasso with $\lambda_{0}$ chosen as $\lambda_{2}=\{2(\log p)/n\}^{1/2}$ performs almost as well as the lasso with the optimal fixed $\lambda$.

In addition, we compare the prediction performance of all the estimators mentioned in this section. In each replication, we compute the penalized maximum likelihood estimator, its bias-correction, and scaled penalization methods based on the training set of 80 observations. For cross-validation, the training set of 80 observations is further partitioned at random 100 times into two groups of sizes 60 and 20, and a penalty level is selected by minimizing the estimated loss in the smaller group for the lasso estimator based on the larger group. This selected penalty level is then used for the lasso with the entire training set. Thus, the cross-validation lasso is also based on the training set with 80 observations. For the penalty level selected by the adjusted cross-validation, two estimators are considered: the lasso estimator and the least squares estimator after the lasso selection. In Table 4.2, we present the medians of the prediction mean squared error and the selected model size in the 200 replications. The scaled lasso has comparable prediction performance as cross-validation. Again, Table 4.2 suggests that original cross-validation tends to choose larger models, while adjusted cross-validation leads to results comparable with the scaled lasso.

Discussion

In the proof of our theoretical results for the scaled lasso, we use oracle inequalities for fixed penalty which unify and somewhat sharpen existing results. We now present this result. Define

as a sharper version of $\eta(\lambda,\xi,{\beta}^{*},T)$ in (10).

with $\eta_{*}(\lambda,\xi)$ in (12). Moreover, in the same event and with $\mu(\lambda,\xi)$ in (16),

The interpretations of (32) and (33) are given in (15) and (17), along with their relationship to several existing results. We note here that the condition $\kappa(\xi,S)\asymp 1$ for (15) and (17), weaker than the parallel condition on the restricted eigenvalue (Bickel et al. 2009), can be slightly weakened by using $F_{1}(\xi,S)$ in (21) (Ye & Zhang 2010).

Acknowledgement

This research was supported by the National Science Foundation and the National Security Agency. We thank Jian Huang for sharing the gene expression data, and reviewers for valuable suggestions.

Appendix

Here we prove Proposition 2.1, Theorem 5.1, Theorem 3.1, Theorem 3.2 and then Theorem 3.4.

(i) Since ${\widehat{\beta}}={\widehat{\beta}}(\sigma\lambda_{0})$ is a solution of (2) at $\lambda=\sigma\lambda_{0}$,

for all $\widehat{\beta}_{j}(\sigma\lambda_{0})\neq 0$. Since $\{j:\widehat{\beta}_{j}(\lambda)\}$ is unchanged in a neighborhood of $\sigma\lambda_{0}$, $[(\partial/\partial\sigma)\{{\widehat{\beta}}(\sigma\lambda_{0})/\sigma\}]_{j}=0$ for $\widehat{\beta}_{j}(\sigma\lambda_{0})=0$. Thus,

(ii) The convergence of (3) and (4) follows from the joint convexity of $L_{\lambda_{0}}(\beta,\sigma)$. The scale invariance follows from $L_{0}(c{\beta},c\sigma;{X},c{y})=cL_{0}({\beta},\sigma;{X},{y})$, where $L_{0}({\beta},\sigma;{X},{y})$ expresses the dependence of (5) on the data $({X},{y})$.

(i) Let $\widehat{\beta}=\widehat{\beta}(\lambda)$. Since $\sigma^{*}z^{*}={|}{X}^{\prime}({y}-{X}{\beta}^{*}){|}_{\infty}/n$ and ${\dot{\rho}}(|\widehat{\beta}_{j}|/\lambda)=1$ for $\widehat{\beta}_{j}\neq 0$, the inner product of $w-\widehat{\beta}$ and the Karush–Kuhn–Tucker condition (2) yield

Since $2({X}{\widehat{\beta}}-{X}{w})^{\prime}({X}\widehat{\beta}-X\beta^{*})={|}{X}{\widehat{\beta}}-{X}{w}{|}_{2}^{2}+{|}{X}{h}{|}_{2}^{2}-{|}{X}{\beta}^{*}-{X}{w}{|}_{2}^{2}$, this gives the basic inequality (25). Let ${h}={\widehat{\beta}}-{\beta}^{*}$. Since $\sigma^{*}z^{*}\leq\lambda(\xi-1)/(\xi+1)$, $\lambda\{{|}{w}{|}_{1}-{|}{\widehat{\beta}}{|}_{1}\}+\sigma^{*}z^{*}{|}{w}-{\widehat{\beta}}{|}_{1}$ is no greater than $b{|}({w}-{\widehat{\beta}})_{T}{|}_{1}+2\lambda{|}{w}_{T^{c}}{|}_{1}-(b/\xi){|}({w}-{\widehat{\beta}})_{T^{c}}{|}_{1}$ with $b=2\xi\lambda/(\xi+1)$. Thus, (25) implies

with $c={|}{X}{\beta}^{*}-{X}{w}{|}_{2}^{2}/(2n)+2\lambda{|}{w}_{T^{c}}{|}_{1}$. For $T=\emptyset$ and $w=\beta^{*}$, (34) directly yields $|Xh|_{2}^{2}/n\leq c=2\lambda|\beta^{*}|_{1}$. For general $\{w,T\}$, we want to prove

It suffices to consider ${|}{X}{h}{|}_{2}^{2}/n\geq 2c$. In this case, ${\widehat{\beta}}-{w}\in{\mathscr{C}}(\xi,T)$ by (34), so that by (11)

Let $x={|}({w}-{\widehat{\beta}})_{T}{|}_{1}$ and $y={|}{X}{h}{|}_{2}^{2}/n$. It follows from (34) and (35) that $ax^{2}+y\leq 2c+2bx$. For such $(x,y)$, $y-2c\leq\max_{x}\{2bx-ax^{2}\}=b^{2}/a$. This gives $y\leq 2c+b^{2}/a=\eta(\lambda,\xi,{w},T)$.

For ${w}={\beta}^{*}$, it suffices to consider the case $y>c=2\lambda{|}{\beta}^{*}_{T^{c}}{|}_{1}$, where the cone condition holds for $\widehat{\beta}-\beta^{*}$. Now, $(x,y)$ satisfies $ax^{2}\leq y\leq c+bx$. The maximum of $y$, attained at $ax^{2}=c+bx$, is

(ii) Let $0<\nu<1$ and $T\subset\{1,\ldots,p\}$. It follows from (34) with ${w}={\beta}^{*}$ that

It suffices to consider $\nu|{h}|_{1}\geq(\xi+1)|{\beta}^{*}_{T^{c}}|_{1}$. In this case

Thus, $(1-\nu){|}{h}_{T^{c}}{|}_{1}\leq(\xi+\nu){|}{h}_{T}{|}_{1}$, or equivalently ${h}\in{\mathscr{C}}\{(\xi+\nu)/(1-\nu),T\}$. It follows from (11) that ${|}{X}{h}{|}_{2}^{2}/n\geq{|}{h}_{T}{|}_{1}^{2}\kappa^{2}\{(\xi+\nu)/(1-\nu),T\}/|T|$, so that

Let $x={|}{h}_{T}{|}_{1}$ and $y={|}{h}_{T^{c}}{|}_{1}$. Write (37) as $ax^{2}+by\leq cx$. Subject to this inequality, the maximum of $x+y$ is $\max_{x\geq 0}\{x+(cx-ax^{2})/b\}$. This maximum, attained at $2ax=b+c$, is $x(b+c)/(2b)=(b+c)^{2}/(4ab)$. Thus,

This gives ${|}{h}{|}_{1}\leq\mu(\lambda,\xi)$ for $\nu|{h}|_{1}\geq(\xi+1)|{\beta}^{*}_{T^{c}}|_{1}$.

Assume $\tau_{0}<1$ without loss of generality. Consider $t\geq\sigma^{*}(1-\tau_{0})$ and the penalty level $\lambda=t\lambda_{0}$ for the lasso. Since $z^{*}\sigma^{*}\leq\sigma^{*}(1-\tau_{0})\lambda_{0}(\xi-1)/(\xi+1)\leq\lambda(\xi-1)/(\xi+1)$ and $\sigma^{*}={|}{y}-{X}{\beta}^{*}{|}_{2}\big/n^{1/2}$, the Cauchy–Schwarz inequality and (32) imply

Since $\eta_{*}^{1/2}(t\lambda_{0},\xi)\leq\sigma^{*}\tau_{0}$ for $t<\sigma^{*}$, the derivative (7) of the loss with $a=0$ satisfies

This implies $\widehat{\sigma}\geq\sigma^{*}(1-\tau_{0})$ by the strict convexity of the profile loss (5) in $\sigma$. For $t>\sigma^{*}$, $\eta_{*}^{1/2}(t\lambda_{0},\xi)\leq t\tau_{0}$ by (10) and (12), so that at $t=\sigma^{*}/(1-\tau_{0})$,

This implies $\sigma^{*}\geq\widehat{\sigma}(1-\tau_{0})$ by the strict convexity of (5) in $\sigma$. Thus, the first part of (13) holds. Moreover,

The proof of Theorem 2 requires the following lemma.

Let $T_{m}$ have the t-distribution with $m$ degrees of freedom. Then, there exists $\epsilon_{m}\to 0$ such that for all $t>0$

Let $x=[m\{e^{2t^{2}/(m-1)}-1\}]^{1/2}$. Since $T_{m}$ has the $t$-distribution,

where $\epsilon_{m}=\{2/(m-1)\}^{1/2}\Gamma\{(m+1)/2\}/\Gamma(m/2)-1\to 0$ as $m\to\infty$.

We need to express $\tau_{*}^{2}$ as a function of $\sigma$ at $\sigma=\sigma^{*}$ in the proof. Define

We have $\tau_{*}^{2}=\phi(\sigma^{*})<1$, $\phi_{-}\leq\phi(\sigma^{*})$ and $\phi_{+}\leq\phi(\sigma^{*})/(1-\phi(\sigma^{*})$.

(i) Consider $z^{*}\leq(1-\phi_{-})\lambda_{0}(\xi-1)/(\xi+1)$. Let $\lambda={t}\lambda_{0}$ and ${h}={\widehat{\beta}}(\lambda)-{\beta}^{*}$. Since ${|}{X}^{\prime}({y}-{X}{\beta}^{*})/n{|}_{\infty}=z^{*}\sigma^{*}$, the Karush–Kuhn–Tucker condition (2) gives

as lower and upper bounds for $(\sigma^{*})^{2}-{|}{y}-{X}{\widehat{\beta}}(\lambda){|}_{2}^{2}/n$. This is a key point in the proof.

For $t\geq\sigma^{*}(1-\phi_{-})$, $z^{*}\sigma^{*}\leq{t}\lambda_{0}(\xi-1)/(\xi+1)=\lambda(\xi-1)/(\xi+1)$, so that (33) in Theorem 5.1 implies ${|}{h}{|}_{1}\leq\mu(t\lambda_{0},\xi)$. It follows (39) that for $t=\sigma^{*}(1-\phi_{-})$,

due to $\phi_{-}=(\xi-1)(\xi+1)^{-1}\phi(\sigma^{*})=(\xi-1)(\xi+1)^{-1}\lambda_{0}\mu(\sigma^{*}\lambda_{0},\xi)/\sigma^{*}$. As in the proof of Theorem 3.1, we find $\widehat{\sigma}/\sigma^{*}\geq 1-\phi_{-}$ by (7) and the strict convexity of (5) in $\sigma$.

Now we prove that $\widehat{\sigma}/\sigma^{*}\leq 1+\phi_{+}$. For $t>\sigma^{*}$, $\mu(t\lambda_{0},\xi)\leq(t/\sigma^{*})\mu(\sigma^{*}\lambda_{0},\xi)$ by (16). Thus, since $(\xi-1)/(\xi+1)+1=2\phi_{+}\{1-\phi(\sigma^{*})\}/\phi(\sigma^{*})$ and $\phi_{+}\leq(1+\phi_{+})\phi(\sigma^{*})$, for $t/\sigma^{*}=1+\phi_{+}$, (39) and (33) imply that

It follows that $\widehat{\sigma}/\sigma^{*}\leq 1+\phi_{+}$ by convexity.

Since $1-\phi_{-}\leq\widehat{\sigma}/\sigma^{*}\leq 1+\phi_{+}$, ${|}{\widehat{\beta}}(\widehat{\sigma}\lambda_{0})-{\beta}^{*}{|}_{1}\leq\mu(\widehat{\sigma}\lambda_{0},\xi)\leq\mu(\sigma^{*}\lambda_{0},\xi)(1+\phi_{+})$. This completes the proof of (18).

Since $e^{a}-1\leq\sum_{k=1}^{\infty}a^{k}/2^{k-1}=a/(1-a/2)$ for any $0<a<2$,

Since $\lambda_{0}\geq\{(2/n)\log(p/\epsilon)\}^{1/2}(\xi+1)/\{(\xi-1)(1-\phi_{-})\}$, this bounds the tail probability of $z^{*}=\max_{j\leq p}|z_{j}|$ by the union bound. Since $n(\sigma^{*}/\sigma)^{2}$ follows the $\chi^{2}_{n}$ distribution, $n^{1/2}(\sigma^{*}/\sigma-1)$ converges to $N(0,1/2)$ in distribution, which then implies (19) by (18) under $\phi(\sigma)=o(n^{-1/2})$.

Let $h=\widehat{\beta}-\beta^{*}$. It follows from the proof of Theorem 3.2 (i) that $z^{*}\leq(1-\phi_{-})\lambda_{0}(\xi-1)/(\xi+1)\leq(\widehat{\sigma}/\sigma^{*})\lambda_{0}(\xi-1)/(\xi+1)$, so that ${\widehat{\lambda}}+z^{*}\sigma^{*}\leq\xi({\widehat{\lambda}}-z^{*}\sigma^{*})$. By (2), $|x_{j}^{\prime}Xh/n|=|x_{j}^{\prime}(y-X\widehat{\beta}-\varepsilon^{*})/n|\geq{\widehat{\lambda}}-z^{*}\sigma^{*}$ for $\widehat{\beta}_{j}\neq 0$. Let $B\subseteq{\widehat{S}}\setminus S$ with $|B|\leq m$. Since $\kappa_{+}(m,S)$ is the upper sparse eigenvalue, $({\widehat{\lambda}}-z^{*}\sigma^{*})^{2}|B|\leq\kappa_{+}(m,S)|Xh|_{2}^{2}/n$. By the basic inequality (25) with $w=\beta^{*}$, $|Xh|_{2}^{2}/n\leq({\widehat{\lambda}}+z^{*}\sigma^{*})|h_{S}|_{1}$ and $h$ is in the cone ${\mathscr{C}}(\xi,S)$. Thus, since $|h_{S}|_{1}^{2}\kappa^{2}(\xi,S)\leq|Xh|_{2}^{2}|S|/n$ by (11), $|Xh|_{2}^{2}/n\leq({\widehat{\lambda}}+z^{*}\sigma^{*})^{2}|S|/\kappa^{2}(\xi,S)$. It follows that $|B|\leq\kappa_{+}(m,S)\xi^{2}|S|/\kappa^{2}(\xi,S)<m$. Since all $B\subseteq{\widehat{S}}\setminus S$ of size $|B|\leq m$ have size $|B|<m$, ${\widehat{S}}\setminus S$ does not have a subset of size $m$. This gives the first inequality in (28).

Let $P_{B}$ be the orthogonal projection to the linear span of $(x_{j},j\in B)$. By the definition of $\sigma^{*}_{m,S}$ and the prediction error bound $|Xh|_{2}^{2}/n\leq\eta_{*}({\widehat{\lambda}},\xi)$ in Theorem 5.1,

This gives the second inequality in (28). The prediction error bound follows from

References