Year: 2021

Abstract:
Lasso is a celebrated method for variable selection in linear models, but it faces challenges when the variables are moderately or strongly correlated. This motivates alternative approaches such as using a non-convex penalty, adding a ridge regularization, or conducting a post-Lasso thresholding. In this paper, we compare Lasso with 5 other methods: Elastic net, SCAD, forward selection, thresholded Lasso, and forward backward selection. We measure their performances theoretically by the expected Hamming error, assuming that the regression coefficients are iid drawn from a two-point mixture and that the Gram matrix is block-wise diagonal. By deriving the rates of convergence of Hamming errors and the phase diagrams, we obtain useful conclusions about the pros and cons of different methods.

heading: 1 INTRODUCTION
text: Variable selection is one of the core problems in high-dimensional data analysis. Consider a linear regression, where the response y ∈ Rn and the design matrix X = [x1, . . . , xp] ∈ Rn×p satisfy that
y = Xβ + z, ‖xj‖ = 1, z ∼ N (0, σ2In). (1) The goal is estimating the support of β (Supp(β)). Lasso (Tibshirani, 1996) is a popular method:
β̂lasso = argminβ { ‖y −Xβ‖2/2 + λ‖β‖1 } . (2)
Lasso has good rates of convergence on the Lq-estimation error or prediction error (Bickel et al., 2009). However, it can be unsatisfactory for variable selection, especially when the columns in the design matrix are moderately or strongly correlated. Zhao & Yu (2006) showed that an irrepresentable condition on X is necessary for Lasso to recover Supp(β) with high probability, and such a condition is restrictive when p is large (Fan & Lv, 2010). Ji & Jin (2012) studied the Hamming error of Lasso and revealed Lasso’s non-optimality by lower-bounding its Hamming error rate. Many alternative strategies were proposed for variable selection, such as using non-convex penalties (Fan & Li, 2001; Zhang, 2010; Shen et al., 2012), adding a ridge regularization (Zou & Hastie, 2005), post-processing on the Lasso estimator (Zou, 2006; Zhou, 2009), and iterative algorithms (Zhang, 2011; Donoho et al., 2012). In this paper, our main interest is to theoretically compare these different strategies.
Existing theoretical studies focused on ‘model selection consistency’ (e.g., Fan & Li (2001); Zhao & Yu (2006); Zou (2006); Meinshausen & Bühlmann (2010); Loh & Wainwright (2017)), which uses P(Supp(β̂) = Supp(β)) to measure the performance of variable selection. However, for many real applications, the study of the Hamming error (i.e., total number of false positives and false negatives) is in urgent need. For example, in genome-wide association studies (GWAS) or Genetic Regulatory Network, the goal is to identify the genes or SNPs that are truly associated with a given phenotype, and we hope to find a multiple testing procedure that simultaneously controls the FDR and maximizes the power (for multiple testing). This problem can be re-cast as minimizing the Hamming error in a special regression setting (Efron, 2004; Jin, 2012; Sun & Cai, 2007). This motivates us to study the Hamming errors of variable selection methods, which were rarely considered in the literature.
We adopt the rare and weak signal model (Donoho & Jin, 2004; Arias-Castro et al., 2011; Jin & Ke, 2016), which is often used in theoretical analysis of sparse linear models. Let p be the asymptotic parameter. Given constants ϑ ∈ (0, 1) and r > 0, we assume that βj’s are iid generated such that
βj = { τp, with probability p, 0, with probability 1− p,
where p = p−ϑ, τp = √ 2r log(p). (3)
As p→∞, ‖β‖0 ≈ p1−ϑ, and a nonzero βj is at the minimax order √
log(p). The two parameters (ϑ, r) capture the sparsity level and signal strength, respectively. We may generalize (3) to let nonzero βj’s take different values in [τp,∞), but the current form is more convenient for presentation. The blockwise covariance structure is frequently observed in real applications. In genetic data, there may exist strong correlations between nearby genetic markers, but the long-range dependence is usually negligible; as a result, the sample covariance matrix is approximately blockwise diagonal (Dehman et al., 2015). In financial data, the sample covariance matrix of stock returns (after common factors are removed) is also approximately blockwise diagonal, where each block corresponds to an industry group (Fan et al., 2015). Motivated by these examples, we consider an idealized setting, where the Gram matrix G = X ′X is block-wise diagonal consisting of 2× 2 blocks
G = diag(B,B, . . . , B,B0), where B = [ 1 ρ ρ 1 ] and B0 = { B, if p is even, 1, if p is odd. (4)
This is an idealization of the blockwise covariance structures in real applications. We may generalize (4) to allow unequal block sizes and unequal off-diagonal entries, but we keep the current form for convenience of presentation. Model (4) is also closely connected to the random designs in compressed sensing (Donoho, 2006). Write X = [X1, X2, . . . , Xn]′. Suppose X1, X2, . . . , Xn are iid generated from N ( 0, n−1Σ), where Σ has the same form as G in (4). In a high-dimensional sparse setting, we have ‖β‖0 n p. Then, G = X ′X ≈ Σ, and due to the blessing of sparsity of β, Gβ ≈ Σβ. As a result, X ′y (sufficient statistic of β) satisfies that X ′y = Gβ +N (0, G) ≈ Σβ +N (0,Σ), and the right hand side reduces to Model (4) (Genovese et al., 2012). In Section 3.5, we formally show that this random design setting is asymptotically equivalent to Model (4), because they share the same rate of convergence for the optimal Hamming error of variable selection.
Now, under model (3) and model (4), we have three parameters (ϑ, r, ρ). They capture the sparsity level, signal strength and design correlations, respectively. Our main results are the explicit convergence rates of Hamming error, as a function of (ϑ, r, ρ), for different methods. We will study six methods: (i) Lasso as in (2); (ii) Elastic net (Zou & Hastie, 2005), which adds an additional L2penalty to (2), (iii) smoothly clipped absolute deviation (SCAD) (Fan & Li, 2001), which replaces the L1-penalty by a non-convex penalty, (iv) thresholded Lasso (Zhou, 2009), which further thresholds the Lasso solution, and two iterative algorithms, (v) forward selection and (vi) forward backward selection (Huang et al., 2016); see Section 3 for a precise description of each method. To our best knowledge, our results are the first that directly compare Hamming errors of these methods.

heading: 2 A PREVIEW OF MAIN RESULTS AND SOME DISCUSSIONS
text: For any β̂, its Hamming error is H(β̂, β) = ∑p j=1 1{β̂j 6= 0, βj = 0}+ ∑p j=1 1{β̂j = 0, βj 6= 0}. We shall show that, for each of the six methods studied here, there exists a function h(ϑ, r, ρ) ∈ [0, 1] such that E[H(β̂, β)] = Lpp1−h(ϑ,r,ρ), where Lp is a multi-log(p) term (a multi-log(p) term is such that Lp · p →∞ and Lp · p− → 0 for any > 0. Noting that the expected number of true relevant variables is p1−ϑ, we are interested in three cases:
• Exact recovery: h(ϑ, r, ρ) > 1. In this case, the expected Hamming error is o(1) as p→∞. It follows that model selection consistency holds.
• Almost full recovery: ϑ < h(ϑ, r, ρ) < 1. In this case, the expected Hamming error does not vanish as p→∞, but it is much smaller than the total number of true relevant variables. Variable selection is still satisfactory (although model selection consistency does not hold).
• No recovery: h(ϑ, r, ρ) ≤ ϑ. In this case, the expected Hamming error is comparable with or much larger than the total number of true relevant variables. Variable selection fails.
We call the two-dimensional space (ϑ, r) the phase space. For each fixed ρ, the phase space is divided into three regions: Region of Exact Recovery (ER), which is the subset {(ϑ, r) : h(ϑ, r, ρ) > 1}, and Region of Almost Full Recovery (AFR) and Region of No Recovery (NR) defined similarly. This gives rise to a phase diagram for each method. We denote the curve separating ER region and AFR region by r = U(ϑ) and the curve separating AFR region and NR region by r = L(ϑ); they are called the upper and lower phase curves, respectively. The phase diagram and phase curves are convenient ways to visualize the rates of convergence of the Hamming error.
Figure 1 shows the phase curves for six methods (explicit expressions are in the theorems in Section 3). These phase curves depend on the correlation parameter ρ. Under our model, for each diagonal block (j, j + 1), it holds that E[x′jy|β] = βj + ρβj+1, where βj , βj+1 ∈ {0, τp}. Therefore, a positive ρ boosts the signal at each individual site (i.e., E[x′jy|β] ≥ βj), while a negative ρ leads to potential ‘signal cancellation’ (i.e., E[x′jy|β] ≤ βj). This is why the phase curves have different shapes for positive and negative ρ. In Figure 1, we plot the phase curves for ρ = 0.5 and ρ = −0.4. For other positive/negative value of ρ, the patterns are similar.
Discussion of SCAD. SCAD is a representative of the non-convex penalization methods. There were inspiring works that demonstrated the advantage of using a non-convex penalty (e.g., Fan & Peng (2004); Loh & Wainwright (2017)). Our results support their insights from a different angle: The phase curve of SCAD is strictly better than that of Lasso, when ϑ < 0.5 and ρ < 0. Furthermore, our results illustrate where the advantage of SCAD comes from — compared with Lasso, it handles ‘signal cancellation’ better. To see this, we recall that ‘signal cancellation’ can only happen for ρ < 0. Moreover, under our model (3), for a diagonal block {j, j + 1}, we get a pair of signals (i.e., βj and βj+1 are both nonzero) with probability 2p = p
1−2ϑ approximately. Therefore, ‘signal cancellation’ becomes problematic only when ϑ < 0.5 and ρ < 0 both hold. This explains why the phase curves of SCAD and Lasso are the same for the other values of ϑ and ρ. We note that in previous studies (e.g., Loh & Wainwright (2017)), the advantage of a non-convex penalty in handling ‘signal cancellation’ is reflected in the conditions of (X,β) required for model selection consistency. Our results support this viewpoint from a different angle.
The performance of SCAD can be further improved by adding an entry-wise thresholding on β̂. We are sure that the phase diagrams of thresholded SCAD will be better than those of SCAD, although we do not have the results in hand (the analysis is too tedious). We are also cautious in comparing SCAD and thresholded Lasso. In our settings, Lasso has no model selection consistency because the signals are too weak (i.e., r is not sufficiently large). In such a case, thresholded Lasso outperforms SCAD in terms of the Hamming error. Alternatively, there are cases where Lasso has no model selection consistency no matter how large the signal strength is (Zhao & Yu, 2006). For those cases, it is possible that SCAD is better than thresholded Lasso (see Wainwright (2009) for a related study).
Discussion of Elastic net. The phase curve of Elastic net is worse than that of Lasso. As we will explain in Section 3.1, Elastic net is a ‘bridge’ between Lasso and marginal regression. Since the phase curve of marginal regression is always worse than that of Lasso for the blockwise diagonal design, we do not benefit from using Elastic net in the current setting. We must note that Elastic net was motivated by genetic applications where several correlated variables are competing as predictors, and where it is implicitly assumed that groups of correlated variables tend to be all relevant or all irrelevant (Zou & Hastie, 2005). This is not captured by our model (3). Therefore, our result has no contradiction to the benefit of Elastic net known in the literature. Our result actually supports that the
advantage of Elastic net comes from ‘group’ appearance of signal variables. We conjecture that for a different signal model (e.g., a group prior on β), Elastic net has a better phase diagram than Lasso.
Discussion of thresholded Lasso. Thresholded Lasso is a representative of improving Lasso by post-processing. There have been inspiring works that demonstrate the advantage of such a postprocessing (van de Geer et al., 2011; Wang et al., 2020; Weinstein et al., 2020). Our results support these insights from a different angle. It is surprising (and very encouraging) that the improvement by post-Lasso thresholding is so significant. We note that Lasso is a 1-stage method, which solves a single optimization to obtain β̂. In comparison, thresholded Lasso is a 2-stage method. Lasso has only one algorithm parameter λ, while thresholded Lasso has two algorithm parameters λ and t (the threshold). In Lasso, we control false positives and false negatives via the same algorithm parameter λ, and it is sometimes hard to find a value of λ that simultaneously controls the two types of errors well. In thresholded Lasso, the two types of errors are controlled separately by two algorithm parameters. This explains why thresholded Lasso has such a big improvement over Lasso. It inspires us to modify other 1-stage methods, such as SCAD, by adding a post-processing step of thresholding. For example, we conjecture that thresholded SCAD also has a strictly better phase diagram than that of Lasso, even for a positive ρ. On the other hand, thresholding is not a free lunch. It yields one more tuning parameter to decide in practice. Our theoretical results are based on ideal tuning parameters. How to properly choose these tuning parameters in a data-driven way is an interesting question. Weinstein et al. (2020) proposed a promising approach, where they used cross-validation to select λ and FDR control by knock-off to select t. We leave it to future work to study the phase diagram with data-driven tuning parameters.
Discussion of the two iterative algorithms. We consider two iterative algorithms, forward selection (‘Forward’) and forward backward selection (FB). The FB algorithm we analyze is a simplified version in Huang et al. (2016), which has only one backward step (after all the forward steps are completed) by thresholding the refitted least-squares solution. Our results show that both methods outperform Lasso, and between these two methods, FB is strictly better than Forward. In the literature, there were very interesting theoretical works showing the advantage of iterative algorithms for variable selection (Donoho et al., 2012; Zhang, 2011). Our results support their insights from a different angle. We discover that, for a wide range of ρ, FB has the best phase diagram among all six methods. This is a very encouraging result. At the same time, we keep in mind that the performance of an iterative algorithm tends to be more sensitive to the form of the design, due to its sequential nature.

heading: 3 MAIN RESULTS
text: Consider model (1), (3), and (4), where we set σ2 = 1 without loss of generality. Let E[H(β̂, β)] be the expected Hamming error, where the expectation is with respect to the randomness of β and z. Let Lp denote a generic multi-log(p) term such that Lpp →∞ and Lpp− → 0 for any > 0. Theorem 1. Under Models (1), (3), and (4), for each of the methods considered in this paper (Lasso, Elastic net, SCAD, thresholded Lasso, forward selection, forward backward selection, as well as marginal regression in Section 3.1), there exists a function h(ϑ, r, ρ) such that E[H(β̂, β)] = Lpp
1−h(ϑ,r,ρ). The explicit expressions of h(ϑ, r, ρ), which may depend on the tuning parameters of a method, are given in Theorems B.1, C.1, D.1-D.3, F.1, G.1, H.1-H.4 of the supplement.
In the main article, to save space, we only present the expressions of the upper phase curve U(ϑ) = U(ϑ; ρ) and the lower phase curve L(ϑ) = L(ϑ; ρ) for each method, which are defined as follows:
U(ϑ; ρ) = inf{r > 0 : h(ϑ, r, ρ) > 1}, L(ϑ; ρ) = inf{r > 0 : h(ϑ, r, ρ) > ϑ}. (5) These two curves describe the phase diagram: The upper phase curve U(ϑ) separates the ER region and AFR region, and the lower phase curve L(ϑ) separates the AFR region and NR region.

heading: 3.1 ELASTIC NET AND LASSO
text: The Elastic net (Zou & Hastie, 2005) is a method that estimates β by β̂EN = argminβ { ‖y −Xβ‖2/2 + λ‖β‖1 + (µ/2)‖β‖ 2} . (6)
Compared with Lasso, it adds an additional L2-penalty to the objective function. Below, we fix µ > 0 and re-parametrize λ = √ 2q log(p), for some constant q > 0. The choice of q affects the exponent,
1− h(ϑ, r, ρ), in the expression of E[H(β̂, β)]. We choose the ideal q that minimizes 1− h(ϑ, r, ρ). The next theorem is proved in the supplement. Theorem 2 (Elastic Net). Under Models (1), (3), and (4), let β̂EN be the Elastic net estimator in (6). Fix µ and write η = ρ/(1 + µ). Let λ = √ 2q log(p) with an ideal choice of q that minimizes the
exponent of E[H(β̂, β)]. The phase curves are given by L(ϑ) = ϑ, and
U(ϑ) = { max {h1(ϑ), h2(ϑ)} , when ρ ≥ 0, max {h1(ϑ), h2(ϑ), h3(ϑ), h4(ϑ)} , when ρ < 0,
where h1(ϑ) = (1 + √ 1− ϑ)2, h2(ϑ) = ( 1−|η| 1−|ρ| + √ 1+η2−2ρη 1−|ρ| )2 (1 − ϑ), h3(ϑ) = 1(1−|ρ|)2 ( 1 + √ 1+η2−2ρη 1+|η| √ 1− 2ϑ )2 , and h4(ϑ) = 1+η 2−2ρη (1−2|ρ|+ρη)2+ (√ 1− ϑ+ 1−|η|1+|η| √ 1− 2ϑ )2 .
Lasso is a special case with µ = 0. By setting µ = 0 (equivalently, η = ρ) in Theorem 2, we obtain the phase curves for Lasso. They agree with the results in Ji & Jin (2012) (but Ji & Jin (2012) does not cover Elastic net).
To see the effect of the L2-penalty, we consider an extreme case where µ→∞. Some elementary algebra shows that (1 + µ)β̂EN converges to the soft-thresholding of X ′y at the threshold λ. In other words, as µ→∞, Elastic net converges to marginal screening (i.e., select variables by thresholding the marginal regression coefficients). At the same time, when µ = 0, (1 + µ)β̂EN equals the Lasso estimate. Hence, Elastic net serves as a bridge between Lasso and marginal regression. In the setting here, the phase diagram of marginal regression is inferior to that of Lasso, and so the phase diagram of Elastic net is also inferior to that of Lasso. See the proposition below and Figure 2: Proposition 1. In Theorem 2, for each fixed ϑ ∈ (0, 1), as µ → 0, U(ϑ) is monotone decreasing and converges to ULasso(ϑ), which is the upper phase curve of Lasso; as µ → ∞, U(ϑ) is monotone increasing and converges to UMR(ϑ), which is the upper phase curve of marginal regression. Furthermore, when ρ ≤ − 12 , UMR(θ) =∞ for all 0 < ϑ ≤ 1 2 (i.e. exact recovery is impossible to achieve no matter how large r is).

heading: 3.2 SMOOTHLY CLIPPED ABSOLUTE DEVIATION PENALTY (SCAD)
text: SCAD (Fan & Li, 2001) is a non-convex penalization method. For any a > 2, it defines a penalty function qλ(θ) on (0,∞) by qλ(θ) = ∫ θ 0 q′λ(t)dt, where q ′ λ(θ) = λ { I(θ ≤ λ) + (aλ−θ)+(a−1)λ I(θ > λ) } . The resulting penalty function qλ(·) coincides with the L1-penalty in (0, λ] and becomes a constant in [aλ,∞). Let Qλ(β) = ∑p j=1 qλ(|βj |). Then, SCAD estimates β by
β̂SCAD = argminβ { ‖y −Xβ‖2/2 +Qλ(β) } . (7)
The following theorem is proved in the supplemental material (see Figure 3, left panel): Theorem 3 (SCAD). Under Models (1), (3), and (4), let β̂SCAD be the SCAD estimator in (7). Fix a ∈ (2, 21−|ρ| ). Let λ = √ 2q log(p) with an ideal choice of q that minimizes the rates of convergence of the expected Hamming error. The phase curves are given by L(ϑ) = ϑ, and
U(ϑ) = { max {h1(ϑ), h2(ϑ), h3(ϑ)} , when ρ ≥ 0, max {h1(ϑ), h2(ϑ), h4(ϑ), h5(ϑ)} , when ρ < 0,
where h1(ϑ) = (1 + √ 1− ϑ)2, and h2(ϑ) = ( 1 + √ 1+|ρ| 1−|ρ| )2 (1−ϑ), h4(ϑ) = (√ 1−2ϑ 1−|ρ|2 + 1 1−|ρ| )2 ,
h3(ϑ) = ( 3+ρ 2(1−ρ2) √ 1+ρ 1−ρ √ 1− ϑ+ 12 √ 2(1−2ϑ) 1+ρ − (1−ϑ) (1−ρ)2 )2 , and
h5(ϑ) =  ( 5+3|ρ| 1−|ρ| ) (1− ϑ), if √ 1−2ϑ 1−ϑ ≥ 3−4|ρ|−3ρ2 (1−|ρ|) √ 1+|ρ| 5+3|ρ| , 1 (1−|ρ|)2 (√ 1+|ρ| 1−|ρ| √ 1− ϑ+ √ 1−|ρ| 1+|ρ| √ 1− 2ϑ )2 , if √ 1−2ϑ 1−ϑ ≤ (1+|ρ|)(1−2|ρ|) 1−|ρ| ,
h6(ϑ) other wise,
with
h6(ϑ) =  √ 1− 2ϑ 1− ρ2 + 1−2|ρ| 1−|ρ| √ 1−2ϑ 1−ρ2 + √[( 1−2|ρ| 1−|ρ| )2 + 1−|ρ|1+|ρ| ] (1− ϑ)− 1−2ϑ(1+|ρ|)2 (1− |ρ|) [( 1−2|ρ| 1−|ρ| )2 + 1−|ρ|1+|ρ| ]  2 .
Note that the phase curves of Lasso are given in Theorem 2 by setting η = ρ. We compare SCAD with Lasso. When ρ < 0, the upper phase curve in Theorem 3 is strictly lower than that of Lasso (see Figure 3, middle and right panels). When ρ ≥ 0, the upper phase curve in Theorem 3 is sometimes higher than that of Lasso. Note that we restrict a < 21−|ρ| in Theorem 3. In fact, a larger a may be preferred for ρ ≥ 0. The next proposition is about using an optimal a. Proposition 2. In the SCAD estimator, we choose a = a∗ and λ = √ 2q∗ log(p) such that (a∗, q∗) = (a∗(ϑ, r, ρ), q∗(ϑ, r, ρ)) minimize the rates of convergence of the expected Hamming error among all choices of (a, q). Let U∗(ϑ) be the resulting upper phase curve for SCAD. Then, U(ϑ) = ULasso(ϑ) when ρ ≥ 0, and U(ϑ) < ULasso(ϑ) when ρ < 0.
The phase curves of SCAD are insensitive to the choice of a. When a < 0, the optimal a∗ can be any value in (2, 21−|ρ| ). When ρ ≥ 0, there exists a constant c = c(ϑ, ρ) such that the optimal a
∗ is any value in (c,∞). As a → ∞, the SCAD penalty reduces to the L1-penalty. This explains why the phase curve of SCAD is the same as that of Lasso when ρ ≥ 0.

heading: 3.3 THRESHOLDED LASSO
text: Let β̂Lasso be the Lasso estimator in (2). The thresholded Lasso estimator β̂TL is obtained by applying coordinate-wise hard-thresholding to the Lasso estimator:
β̂TLj = β̂ Lasso · 1{|β̂Lasso| > t}, 1 ≤ j ≤ p. (8)
Theorem 4 (Thresholded Lasso). Under Models (1), (3), and (4), let β̂TL be the thresholded Lasso estimator in (8). Let λ = √ 2q log(p) and t = √ 2w log(p) with the ideal (q, w) that minimize the exponent of the expected Hamming error. The phase curves are given by L(ϑ) = ϑ, and
U(ϑ) = { max {h1(ϑ), h2(ϑ)} , when ρ ≥ 0, max {h1(ϑ), h2(ϑ), h3(ϑ)} , when ρ < 0,
where h1(ϑ) = (1 + √ 1− ϑ)2, h2(ϑ) = 4(1−ϑ)1−ρ2 , and h3(ϑ) = ( 1 + 1+|ρ|2 √ 1−ϑ 1−ρ2 + 1−|ρ| 2 √ 1−2ϑ 1−ρ2 )2 .
See Figure 4 for a comparison with Lasso (a special case of t = 0). With the flexibility of using an optimal t, the phase diagram of thresholded Lasso is always better than that of Lasso.
Theorem 4 also gives other interesting facts about thresholded Lasso. First, the shape of phase curves is much less affected by the sign of ρ. This differs from Lasso, Elastic net, and SCAD, for which the shape of phase curves is significantly different for positive and negative ρ. Second, the optimal λ in thresholded Lasso is considerably smaller than the optimal λ in Lasso (it can be seen from the proofs of Theorem 4 and Theorem 2). This is because the λ in thresholded Lasso only serves to control false negatives, but the λ in Lasso is used to simultaneously control false positives and false negatives, hence, cannot be too small. We observe the same phenomenon in simulations; see Section 4.

heading: 3.4 FORWARD SELECTION AND FORWARD BACKWARD SELECTION
text: Forward selection is a classical textbook method for variable selection. Write X = [x1, x2, . . . , xp], where xi ∈ Rn for 1 ≤ i ≤ p. For any subset A ⊂ {1, 2, . . . , p}, let P⊥A be the projection onto the orthogonal complement of the linear space spanned by {xi : i ∈ A}. Given a threshold t > 0, the forward selection algorithm initializes with S0 = ∅ and r̂0 = y. At the kth iteration, compute
i∗ = argmaxi/∈Sk−1 |x ′ ir̂k−1|, δ = |x′i∗ r̂k−1|/‖P⊥Sk−1xi∗‖.
If δ > t, compute Sk = Sk−1 ∪ {i∗} and r̂k = P⊥Sky; otherwise, output β̂ forward as the least-squares estimator restricted to Sk−1. The stopping rule of δ ≤ t is equivalent to measuring the decrease of the residual sum of squares. The following theorem is proved in the supplemental material:
Theorem 5 (Forward Selection). Under Models (1), (3), and (4), let β̂forward be the estimator from forward selection. Let t = √ 2q log(p) with the ideal q that minimizes the exponent of the expected Hamming error. The phase curves are given by L(ϑ) = ϑ, and
U(ϑ) = { max {h1(ϑ), h2(ϑ), h3(ϑ)} , when ρ ≥ 0, max {h1(ϑ), h2(ϑ), h3(ϑ), h4(ϑ)} , when ρ < 0,
with h1(ϑ)=(1+ √ 1− ϑ)2, h2(ϑ) = 2(1−ϑ)1−|ρ| , h3(ϑ) = (1+ √ 1−2ϑ)2 1−ρ2 , h4(ϑ) = (√ 1−2ϑ 2(1−|ρ|) + 1 1−|ρ| )2 .
Forward backward selection (FB) modifies forward selection by allowing to drop variables. We use the FB algorithm in Huang et al. (2016), where the backward step is conducted after all the forward steps are finished. For a threshold v > 0, it applies entry-wise thresholding on β̂forward:
β̂FBj = β̂ forward j · 1{|β̂forwardj | > v}, 1 ≤ j ≤ p. (9)
Theorem 6 (Forward Backward Selection). Under Models (1), (3), and (4), let β̂FB be the estimator from forward selection. Let t = √ 2q log(p) and v = √ 2u log(p) with the ideal (q, u) that minimize the exponent of the expected Hamming error. When ρ ≥ 0, the phase curves are given by L(ϑ) = ϑ, and
U(ϑ) = max {h1(ϑ), h2(ϑ), h∗3(ϑ)} ,
where h1(ϑ) and h2(ϑ) are the same as in Theorem 5 and h∗3(ϑ) = ( √ 1−ϑ+ √ 1−2ϑ)2 1−ρ2 . When ρ < 0,
U(ϑ) ≤ max {g1(ϑ), g2(ϑ), g3(ϑ), g4(ϑ)} ,
where g1(ϑ) = (vmin(ϑ) + √ 1− ϑ)2, g2(ϑ) = 2(1−ϑ)1−|ρ| , g3(ϑ) = (√ 1−2ϑ 1−ρ2 + vmin(ϑ) )2 , g4(ϑ) =(√
1−2ϑ 2(1−|ρ|) + tmin(ϑ) 1−|ρ| )2 , vmin(ϑ) = max { 1, √ 1−ϑ 1−ρ2 } , and tmin(ϑ) = max {√ 2 2 , vmin(ϑ) 1+|ρ|/ √ 1−ρ2 } .
Theorem 6 gives U(ϑ) for ρ ≥ 0 and an upper bound of it for ρ < 0. Combining it with Theorems 2 and 5, we conclude that the upper phase curve of FB is always better than those of Lasso and forward selection (for ρ < 0, the upper bound here is already better than U(ϑ) for the other two methods).
We remark that we did obtain the exact phase curve for ρ < 0 in the proof of Theorem 6. It is just too complicated and space-consuming to present it in the main text. However, given specific values of (ϑ, ρ), we can always plot the exact phase curve using the (complicated) formulas in the supplement. In Figures 1 and 5, the phase curves of FB are indeed the exact ones.

heading: 3.5 CONNECTION TO THE RANDOM DESIGN MODEL
text: Consider the random design as mentioned in Section 1. The minimax Hamming error isH∗(ϑ, r, ρ) = inf β̂ E[H(β̂, β)], where the infimum is taken over all methods β̂ and the expectation is with respect to the randomness of (X,β, z). We can define H∗(ϑ, r, ρ) in the same way for our current model (4). The minimax Hamming error is related to the statistical limit of the model setting, but not any specific method. The next theorem shows that, when n sp = p1−ϑ (we allow both p ≤ n and p > n), the convergence rate of the minimax Hamming error is the same under two models.
Theorem 7. Under Models (1) and (3), suppose X is independent of (β, z) and its rows are iid generated from N (0, n−1Σ), with Σ having the same form as G in (4). Suppose n = pω, with ω > 1−ϑ (note: this allows ω < 1, which corresponds to n p). There exists a number h∗∗(ϑ, r, ρ) such that the minimax Hamming error satisfies that H∗(ϑ, r, ρ) = Lpp1−h
∗∗(ϑ,r,ρ). Furthermore, if we instead have X ′X = Σ (i.e., model (4)), then it also holds that H∗(ϑ, r, ρ) = Lpp1−h ∗∗(ϑ,r,ρ).

heading: 4 SIMULATIONS
text: In Experiments 1-3, (n, p) = (1000, 300). In Experiment 4, (n, p) = (500, 1000).
Experiment 1 (block-wise diagonal designs). We generate (X,β) as in (3)-(4). For each method, we select the ideal tuning parameters that minimize the average Hamming error over 50 repetitions. The averaged Hamming errors and its standard deviations under the ideal tuning parameters over 500 repetitions are reported below. The results are consistent with the theoretical phase diagrams (see Figure 1). E.g., thresholded Lasso and forward backward selection are the two methods that perform the best; Lasso is more unsatisfactory when ρ < 0; SCAD improves Lasso when ρ < 0.
Experiment 2 (general designs). In the Toeplitz design, we let (X ′X)i,j = 0.7|i−j| and set (ϑ, r) = (0.1, 2.5). In the factor model design, we let X ′X = BB′− diag(BB′) + Ip, where entries of B ∈ Rp×2 are iid from Unif(0, 0.6), and set (ϑ, r) = (0.1, 1.5). Same as in Experiment 1, we use the ideal tuning parameters. The averaged Hamming errors and its standard deviations are reported below. The Toeplitz design is a setting where each variable is only highly correlated with a few other variables. The factor model design is a setting where a variable is (weakly) correlated with all the other variables. The results are quite similar to those in Experiment 1. This confirms that the insight gained in the study of the block-wise diagonal design continues to apply to more general designs.
Experiment 3 (tuning parameters). Fix (ϑ, r) = (0.1, 1.5) and ρ ∈ {±0.5} in the block-wise diagonal design. We study the effect of tuning parameters in Lasso, thresholded Lasso (ThreshLasso), forward selection (ForwardSelect), and forward backward selection (FB). In (a)-(b), we show the heatmap of averaged Hamming error (over 50 repetitions) of ThreshLasso for a grid of (t, λ); when t = 0, it reduces to Lasso. In (c)-(d), we show the Hamming error of FB for a grid of (v, t); when v = 0, it reduces to ForwardSelect. Cyan points are theoretically optimal tuning parameters (formulas are in proofs of theorems). Red points are empirically optimal tuning parameters that minimize the averaged Hamming error. The theoretical tuning parameter values are quite close to the empirically optimal ones. Moreover, the optimal λ in ThreshLasso is smaller than the optimal λ in Lasso.
Experiment 4 (p > n and random designs). Fix (n, p, ϑ, r) = (500, 1000, 0.5, 1.5). We simulate data from the random design setting in Theorem 7. We study the average Hamming error over 500 repetitions (tuning parameters are set in the same way as in Experiment 1). See Table 3. We have some similar observations as before: e.g., ThreshLasso and FoBackward are still the best two,

Review summaries:
This paper investigates the theoretical performance of a number of variables selection methods when the columns of the design matrix are correlated. They focus specifically on the case where the Gram matrix for the design is block-wise diagonal (2 x 2 blocks) with the off-diagonals within each block given by a correlation parameter and diagonals set to 1. The quality of each estimator is evaluated by using its expected Hamming distance. Extending ideas presented in _UPS delivers optimal phase diagram in high-dimensional variable selection [Ji and Jin]_ and others, phase diagrams are constructed with phase curves (derived in the appendices) establishing the boundary between each region: exact recovery, almost full recovery, and no recovery. A discussion of each variable selection method provides insight and a comparison between them is included. 

The goal of the paper is to theoretically compare several variable selection schemes related to Lasso. To that end, the authors propose a comprehensive and visual methodology: phase diagrams. Under an idealised setting, they show that, depending on the strength and sparsity of the signal, variants of Lasso fall in one of the three following regimes, Full recovery, Almost Full Recovery, No Recovery. Their main results consist in the exact and non-asymptotic characterisation of those regimes. Once the phase diagrams are computed, the efficiency of algorithms can easily be compared. Most demonstrations rely on the fact that, under the “blockwise correlation” setting, methods can be  reduced to bivariate (Lasso-type) sub-problems.

The paper compares several variable selection methods over block-wise diagonal designs. The performance is theoretically evaluated by the expected Hamming error. Some simulations are implemented to justify the empirical performance of the methods in the comparison.

This paper performs theoretical analyses of several variable selection approaches such as Lasso, Elastic net, SCAD, thresholded Lasso, forward selection, and forward backward selection. While the previous studies focused on the model selection consistency, this paper focused on the impassibleness of the model selection consistency based on Hamming errors. Its theoretical analyses show that SCAD has better performance than Lasso, Elastic net is worse than Lasso, and other approaches are superior to Lasso.