Iteratively Reweighted $\ell_1$-Penalized Robust Regression

TL;DR

This paper develops an iteratively reweighted $ ext{l}_1$-penalized robust regression method that achieves strong statistical guarantees and efficiency in high-dimensional settings with heavy-tailed errors.

Contribution

It introduces a novel adaptive Huber regression estimator with theoretical guarantees under weak conditions and analyzes its computational complexity.

Findings

Estimator satisfies exponential deviation bounds.

Achieves oracle convergence rate and variable selection consistency.

Requires $O( ext{log } s + ext{log log } d)$ iterations.

Abstract

This paper investigates tradeoffs among optimization errors, statistical rates of convergence and the effect of heavy-tailed errors for high-dimensional robust regression with nonconvex regularization. When the additive errors in linear models have only bounded second moment, we show that iteratively reweighted $\ell_1$-penalized adaptive Huber regression estimator satisfies exponential deviation bounds and oracle properties, including the oracle convergence rate and variable selection consistency, under a weak beta-min condition. Computationally, we need as many as $O(\log s + \log\log d)$ iterations to reach such an oracle estimator, where $s$ and $d$ denote the sparsity and ambient dimension, respectively. Extension to a general class of robust loss functions is also considered. Numerical studies lend strong support to our methodology and theory.

Tables

Table 1. Table 1: Simulation results for the Lasso, SCAD, Huber-SCAD, MC+ and Huber-MC+ estimators under the homoscedastic model.

Error dist.		Lasso	SCAD	Huber-SCAD	MC+	Huber-MC+
Normal	TP	6.00(0)	6.00(0)	6.00(0)	6.00(0)	6.00(0)
	FP	24.44(14.25)	3.11(4.53)	2.19(3.87)	0.84(1.91)	0.53(1.27)
	RE1	1.00	0.23(0.12)	0.22(0.11)	0.19(0.09)	0.19(0.09)
	RE2	1.00	0.32(0.13)	0.33(0.13)	0.30(0.12)	0.30(0.12)
Skewed $t$	TP	4.74(1.37)	4.87(1.35)	4.74(1.39)	3.97(1.67)	3.97(1.62)
	FP	20.78(17.10)	18.49(9.65)	11.48(8.82)	4.28(4.24)	2.76(3.23)
	RE1	1.00	0.88(0.22)	0.73(0.23)	0.73(0.21)	0.65(0.22)
	RE2	1.00	0.91(0.17)	0.86(0.19)	0.94(0.20)	0.88(0.23)
Lognormal	TP	5.68(0.87)	5.71(0.84)	6.00(0.07)	5.49(1.14)	5.97(0.37)
	FP	29.70(16.66)	16.75(8.70)	3.80(4.52)	4.32(4.62)	0.91(1.95)
	RE1	1.00	0.54(0.26)	0.15(0.12)	0.42(0.32)	0.13(0.11)
	RE2	1.00	0.62(0.26)	0.23(0.14)	0.60(0.30)	0.22(0.14)
Pareto	TP	5.64(1.09)	5.67(1.01)	6.00(0)	5.44(1.35)	5.98(0.35)
	FP	28.30(16.21)	14.69(8.97)	2.91(4.34)	3.48(3.39)	0.71(1.71)
	RE1	1.00	0.51(0.30)	0.14(0.08)	0.40(0.25)	0.13(0.17)
	RE2	1.00	0.58(0.26)	0.21(0.11)	0.57(0.28)	0.22(0.22)

Table 2. Table 2: Simulation results for the Lasso, SCAD, Huber-SCAD, MC+ and Huber-MC+ estimators under the heteroscedastic model.

Error dist.		Lasso	SCAD	Huber-SCAD	MC+	Huber-MC+
Normal	TP	6.00(0)	6.00(0)	5.96(0.40)	6.00(0)	5.98(0.28)
	FP	22.71(16.51)	3.29(5.76)	0.31(1.68)	0.88(2.03)	0.13(0.70)
	RE1	1.00	0.28(0.17)	0.21(0.16)	0.24(0.14)	0.16(0.18)
	RE2	1.00	0.38(0.19)	0.31(0.16)	0.36(0.18)	0.25(0.14)
Skewed $t$	TP	4.93(1.59)	5.04(1.53)	5.83(0.65)	4.58(1.76)	5.52(1.17)
	FP	22.99(18.62)	18.21(10.83)	2.71(4.29)	4.99(5.14)	0.92(2.42)
	RE1	1.00	0.83(0.30)	0.26(0.26)	0.69(0.29)	0.27(0.27)
	RE2	1.00	0.87(0.28)	0.34(0.26)	0.87(0.29)	0.35(0.28)
Lognormal	TP	5.74(0.96)	5.77(0.91)	6.00(0)	5.65(1.14)	6.00(0)
	FP	26.61(16.51)	11.28(9.50)	1.23(3.55)	2.62(3.45)	0.30(0.86)
	RE1	1.00	0.45(0.28)	0.14(0.13)	0.35(0.23)	0.12(0.09)
	RE2	1.00	0.53(0.26)	0.21(0.15)	0.50(0.26)	0.19(0.12)
Pareto	TP	5.67(1.19)	5.67(1.18)	5.97(0.42)	5.59(1.29)	5.95(0.55)
	FP	25.56(16.04)	10.13(10.31)	0.61(2.03)	2.80(4.06)	0.23(0.91)
	RE1	1.00	0.46(0.29)	0.14(0.12)	0.39(0.28)	0.15(0.23)
	RE2	1.00	0.55((0.29)	0.22(0.16)	0.54(0.31)	0.23(0.35)

Equations

$$y={\mathbf{x}}^{\mathrm{T}}{\mathbf{\beta }}^{\ast }+\epsilon =\sum _{j=1}^d{\beta }_j^{\ast }{x}_j+\epsilon ,$$

$$\displaystyle\ell_{\tau}(x)=\left\{\begin{array}[]{ll}x^{2}/2{}{}&\mbox{if }|x|\leq\tau,\\ \tau|x|-\tau^{2}/2{}{}&\mbox{if }|x|>\tau,\end{array}\right.$$

$$\displaystyle\widehat{}\bm{\beta}\in\mathop{\mathrm{argmin}}_{\bm{\beta}\in\mathbb{R}^{d}}\bigg{\{}\widehat{\mathcal{L}}_{\tau}(\bm{\beta})+\sum_{j=1}^{d}p_{\lambda}(\beta_{j})\bigg{\}},$$

$$\hat{}{\mathbf{\beta }}^{\text{mbox}H-Lasso}\in {\mathrm{argmin}}_{\mathbf{\beta }\in {\mathbb{R}}^d}\{{\hat{\mathcal{L}}}_{\tau }(\mathbf{\beta })+\lambda \Vert \mathbf{\beta }{\Vert }_1\},$$

$$\Vert \hat{}{\mathbf{\beta }}^{\text{mbox}H-Lasso}-{\mathbf{\beta }}^{\ast }{\Vert }_2\le c_1{\sigma }_2(\log d+t{)}^{1/2}\sqrt{\frac{s}{n}}\text{mbox}and$$

$$\Vert \hat{}{\mathbf{\beta }}^{\text{mbox}H-Lasso}-{\mathbf{\beta }}^{\ast }{\Vert }_1\le c_2{\sigma }_2(\log d+t{)}^{1/2}\frac{s}{\sqrt{n}}$$

$$\displaystyle\min_{\bm{\beta}=(\beta_{1},\ldots,\beta_{d})^{{\rm T}}}\bigg{\{}\widehat{\mathcal{L}}_{\tau}(\bm{\beta})+\sum_{j=1}^{d}p^{\prime}_{\lambda}(|\widehat{\beta}_{j}^{(\ell-1)}|)|\beta_{j}|\bigg{\}}~{}~{}~{}~{}~{}~{}~{}({\rm P}_{\ell})$$

$$\underset{\mathbf{\beta }\in {\mathbb{R}}^d}{\min }\{{\hat{\mathcal{L}}}_{\tau }(\mathbf{\beta })+\Vert \mathbf{\lambda }\circ \mathbf{\beta }{\Vert }_1\},$$

$$\nabla {\hat{\mathcal{L}}}_{\tau }(\hat{}\mathbf{\beta })+\mathbf{\lambda }\circ \mathbf{\xi }={\text{0}}_d\text{mbox}forsome\mathbf{\xi }=({\xi }_1,\dots ,{\xi }_d{)}^{\mathrm{T}}\in \partial \Vert \hat{}\mathbf{\beta }{\Vert }_1\subseteq [-1,1{]}^d,$$

$${\omega }_{\mathbf{\lambda }}(\mathbf{\beta }):=\underset{\mathbf{\xi }\in \partial \Vert \mathbf{\beta }{\Vert }_1}{\min }\Vert \nabla {\mathcal{L}}_{\tau }(\mathbf{\beta })+\mathbf{\lambda }\circ \mathbf{\xi }{\Vert }_{\infty },\mathbf{\beta }\in {\mathbb{R}}^d.$$

$$\underset{\mathbf{\beta }\in {\mathbb{R}}^d}{\min }\{{\hat{\mathcal{L}}}_{\tau }(\mathbf{\beta })+\Vert {\mathbf{\lambda }}^{(\ell -1)}\circ \mathbf{\beta }{\Vert }_1\},$$

$${\mathcal{E}}_1(r,l,\kappa )=\left\{\underset{\mathbf{\beta }\in {\mathbf{\beta }}^{\ast }+\mathbb{B}(r)\cap \mathbb{C}(l)}{\inf }\frac{⟨\nabla {\hat{\mathcal{L}}}_{\tau }(\mathbf{\beta })-\nabla {\hat{\mathcal{L}}}_{\tau }({\mathbf{\beta }}^{\ast }),\mathbf{\beta }-{\mathbf{\beta }}^{\ast }⟩}{\Vert \mathbf{\beta }-{\mathbf{\beta }}^{\ast }{\Vert }_2^2}\ge \kappa \right\},$$

$${\mathbf{w}}^{\ast }=\nabla {\hat{\mathcal{L}}}_{\tau }({\mathbf{\beta }}^{\ast })-\nabla {\mathcal{L}}_{\tau }({\mathbf{\beta }}^{\ast })\text{mbox}and{b}_{\tau }^{\ast }=\Vert {\mathcal{L}}_{\tau }({\mathbf{\beta }}^{\ast }){\Vert }_2,$$

$$\lambda \ge s^{-1/2}{b}_{\tau }^{\ast },r>2.5{\kappa }^{-1}{s}^{1/2}\lambda .$$

$$\Vert \tilde{}{\mathbf{\beta }}^{(1)}-{\mathbf{\beta }}^{\ast }{\Vert }_2\le 2.5{\kappa }^{-1}{s}^{1/2}\lambda .$$

$$\Vert \tilde{}{\mathbf{\beta }}^{(1)}-{\mathbf{\beta }}^{\ast }{\Vert }_2\lesssim {\sigma }_2\sqrt{\frac{s\log (d)}{n}}\text{mbox}withhighprobabilityaslongasn\gtrsim s\log (d).$$

$$\lambda \ge s^{-1/2}{b}_{\tau }^{\ast },\kappa \gamma >0.5p^{\prime }(\gamma ),$$

$$0.5p^{\prime }(\gamma )(c^2+1{)}^{1/2}+2=c\kappa \gamma .$$

$$r^{\mathrm{crude}}:=c\gamma s^{1/2}\lambda \le r.$$

$$\Vert \tilde{}{\mathbf{\beta }}^{(\ell )}-{\mathbf{\beta }}^{\ast }{\Vert }_2$$

$$\displaystyle\leq\delta\cdot\|\widetilde{}\bm{\beta}^{(\ell-1)}-\bm{\beta}^{*}\|_{2}+\underbrace{\kappa^{-1}\big{\{}\|p_{\lambda}^{\prime}(|\bm{\beta}^{*}_{\mathcal{S}}|-\gamma\lambda)\|_{2}+\|\bm{w}^{*}_{{\mathcal{S}}}\|_{2}+s^{1/2}\epsilon_{\ell}+b^{*}_{\tau}\big{\}}}_{=:r^{{\rm ora}}}.$$

$$\Vert \tilde{}{\mathbf{\beta }}^{(\ell )}-{\mathbf{\beta }}^{\ast }{\Vert }_2\le {\delta }^{\ell -1}{r}^{\mathrm{crude}}+(1-\delta {)}^{-1}{r}^{\mathrm{ora}}\text{mbox}forany\ell \ge 2.$$

$$\underset{\mathrm{shrinkage}\mathrm{bias}}{\underbrace{\Vert p_{\lambda }^{\prime }(|{\mathbf{\beta }}_{\mathcal{S}}^{\ast }|-\gamma \lambda ){\Vert }_2}},\underset{\mathrm{oracle}\mathrm{rate}\mathrm{plus}approx.\mathrm{bias}}{\underbrace{\Vert {\mathbf{w}}_{\mathcal{S}}^{\ast }{\Vert }_2+b_{\tau }^{\ast }}}\text{mbox}and\underset{\mathrm{optimization}\mathrm{error}}{\underbrace{s^{1/2}{ϵ}_{\ell }}}.$$

$$\hat{}{\mathbf{\beta }}^{\mathrm{ora}}={\mathrm{argmin}}_{\mathbf{\beta }:{\mathbf{\beta }}_{{\mathcal{S}}^{\mathrm{c}}}=\text{0}}{\hat{\mathcal{L}}}_{\tau }(\mathbf{\beta })={\mathrm{argmin}}_{\mathbf{\beta }:{\mathbf{\beta }}_{{\mathcal{S}}^{\mathrm{c}}}=\text{0}}\frac{1}{n}\sum _{i=1}^n{\ell }_{\tau }(y_i-{\mathbf{x}}_{i,\mathcal{S}}^{\mathrm{T}}{\mathbf{\beta }}_{\mathcal{S}}).$$

$$\Vert \hat{}{\mathbf{\beta }}^{\mathrm{ora}}-{\mathbf{\beta }}^{\ast }{\Vert }_2\lesssim \Vert {\mathbf{w}}_{\mathcal{S}}^{\ast }{\Vert }_2+b_{\tau }^{\ast }.$$

$${\mathcal{E}}_2(r,l,\kappa ):=\left\{\underset{({\mathbf{\beta }}^{\prime },{\mathbf{\beta }}^{\prime \prime })\in \mathcal{C}(r,l)}{\inf }\frac{⟨\nabla {\hat{\mathcal{L}}}_{\tau }({\mathbf{\beta }}^{\prime })-\nabla {\hat{\mathcal{L}}}_{\tau }({\mathbf{\beta }}^{\prime \prime }),{\mathbf{\beta }}^{\prime }-{\mathbf{\beta }}^{\prime \prime }⟩}{\Vert {\mathbf{\beta }}^{\prime }-{\mathbf{\beta }}^{\prime \prime }{\Vert }_2^2}\ge \kappa \right\}.$$

$${\mathbf{w}}^{\mathrm{ora}}=\nabla {\hat{\mathcal{L}}}_{\tau }(\hat{}{\mathbf{\beta }}^{\mathrm{ora}}),$$

$$1+0.5p^{\prime }({\gamma }_0)(c_0^2+1{)}^{1/2}=c_0\kappa {\gamma }_0.$$

$$\displaystyle\bigg{\{}\|\bm{w}^{{\rm ora}}\|_{\infty}\leq\frac{p^{\prime}(\gamma_{0})}{2}\lambda\bigg{\}}\cap\bigg{\{}\|\widehat{}\bm{\beta}^{{\rm ora}}-\bm{\beta}^{*}\|_{\infty}\leq\frac{\lambda}{5\delta\kappa}\bigg{\}}$$

$$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}\cap\big{\{}\|\widehat{}\bm{\beta}^{{\rm ora}}-\bm{\beta}^{*}\|_{\Sigma}\leq r\big{\}}\cap\mathcal{E}_{2}(r,l,\kappa),$$

$$\displaystyle{\rm S}=\mathbb{E}(\bm{x}_{{\mathcal{S}}}\bm{x}_{{\mathcal{S}}}^{{\rm T}}),~{}\mbox{ where }\bm{x}_{{\mathcal{S}}}\in\mathbb{R}^{s}\mbox{ is the subvector of $\bm{x}$ indexed by ${\mathcal{S}}$}.$$

$${\rho }_{\mathbf{x}}$$

Full text

62A01,

62J07,

Adaptive Huber regression,

convex relaxation,

heavy-tailed noise,

nonconvex regularization,

optimization error,

oracle property,

oracle rate,

\startlocaldefs\endlocaldefs

and

1. 1 Introduction

2. 1.1 Related literature

3. 1.2 Notation

4. 2 Regularized Huber $M$-estimation

5. 2.1 $\ell_{1}$-penalized Huber regression

6. 2.2 Iteratively reweighted $\ell_{1}$-penalized Huber regression

7. 2.3 Deterministic analysis

8. 2.4 Random analysis

9. 3 Optimization Algorithm

10. 3.1 LAMM algorithm

11. 3.2 Complexity theory

12. 4 Extension to General Robust Losses

13. 5 Numerical Study

14. 6 Preliminaries

15. 7 Proofs of Propositions

16. 7.1 Proof of Proposition 2.1

17. 7.2 Proof of Proposition 2.2

18. 7.3 Proof of Proposition 2.3

19. 7.4 Proof of Proposition 2.4

20. 7.5 Proof of Proposition 2.5

21. 8 Proofs of Theorems

22. 8.1 Proof of Theorem 2.1

23. 8.2 Proof of Theorem 2.2

24. 8.3 Proof of Theorem 2.3

25. 8.3.1 Proof of Lemma 8.1

26. 8.3.2 Proof of Lemma 8.2

27. 8.4 Proof of Theorem 3.1

28. 8.4.1 Technical lemmas

29. 8.4.2 Proof of the theorem

30. 8.5 Proof of Theorem 3.2

31. 8.5.1 Technical lemmas

32. 8.5.2 Proof of the theorem

1 Introduction

Suppose we observe independent and identically distributed (i.i.d.) data vectors $\{(y_{i},\bm{x}_{i}):1\leq i\leq n\}$ from $(y,\bm{x})$ that follows the linear model

[TABLE]

where $\bm{x}=(x_{1},\ldots,x_{d})^{{\rm T}}\in\mathbb{R}^{d}$ with $x_{1}\equiv 1$ is the predictor, $\bm{\beta}^{*}=(\beta^{*}_{1},\ldots,\beta_{d}^{*})^{{\rm T}}\in\mathbb{R}^{d}$ is the vector of regression coefficients with $\beta^{*}_{1}$ denoting the intercept, and $\varepsilon$ is an error term satisfying $\mathbb{E}(\varepsilon|\bm{x})=0$. This setting includes the location-scale model in which $\varepsilon=\sigma(\bm{x})e$, $\sigma(\cdot):\mathbb{R}^{d}\mapsto\mathbb{R}$ is an unknown function, and $e$ is independent of $\bm{x}$ and satisfies $\mathbb{E}(e)=0$. We considers the high-dimensional regime, where the number of features $d$ exceeds the sample size $n$ and $\bm{\beta}^{*}$ is $s$-sparse. Of particular interest is the case where the error variable is asymmetric and heavy-tailed with only bounded second moment.

Since the invention of Lasso two decades ago (Tibshirani, 1996), a variety of variable selection methods have been developed for finding a small group of covariates that are associated with the response from a large pool. The Lasso estimator $\widehat{}\bm{\beta}^{\mbox{{\tiny Lasso}}}$ solves the convex optimization problem $\min_{\bm{\beta}\in\mathbb{R}^{d}}\,(2n)^{-1}\sum_{i=1}^{n}(y_{i}-\bm{x}_{i}^{{\rm T}}\bm{\beta})^{2}+\lambda\|\bm{\beta}\|_{1}$, where $\lambda>0$ is the regularization parameter. The Lasso is an $\ell_{1}$-penalized least squares method in nature: the quadratic loss is used as a goodness of fit measure and the $\ell_{1}$-norm induces sparsity. To achieve better performance under different circumstances, several Lasso variants have been proposed and studied; see, Fan and Li (2001), Zou and Hastie (2005), Zou (2006), Yuan and Lin (2006), Belloni, Chernozhukov and Wang (2011), Sun and Zhang (2012) and Bogdan et al. (2015), to name a few. We refer to Bühlmann and van de Geer (2011), Hastie, Tibshirani and Wainwright (2015) and Wainwright (2019) for comprehensive and systematic introductions of high-dimensional statistical methods and theory.

As a general regression analysis method, the Lasso, along with many of its variants, has two potential downsides. First, the regularized least squares methods are sensitive to the tails of error distributions, even though various alternative penalties have been proposed to achieve better model selection performance. Consider a Lasso-type estimator that solves the penalized empirical risk minimization $\min_{\bm{\beta}\in\mathbb{R}^{p}}\{(1/n)\sum_{i=1}^{n}\ell(y_{i}-\bm{x}_{i}^{{\rm T}}\bm{\beta})+\lambda\|\bm{\beta}\|_{1}\}$, where $\ell(x):\mathbb{R}\mapsto[0,\infty)$ is a general loss function. The effects of the loss and noise on estimation error are coded in the vector $\{\ell^{\prime}(\varepsilon_{i})\}_{i=1}^{n}$. If $\ell$ is the quadratic loss, this vector is likely to have relatively many large coordinates when $\varepsilon$ is heavy-tailed. As a result, the combination of the rapid growth of $\ell$ with heavy-tailed sampling distribution inevitably leads to outliers, which will eventually be translated into spurious discoveries. Secondly, the $\ell_{1}$-penalty introduces nonnegligible estimation bias (Fan and Li, 2001; Zou, 2006). For correlated designs, the bias of the Lasso may offset true signals and creates spurious effects, leading to inconsistency in support recovery. Technically, this is expressed as the irrepresentable condition for the selection consistency of the Lasso (Zhao and Yu, 2006). Under restricted eigenvalue type conditions, the Lasso and its sorted variant Slope (Bellec, Lecué and Tsybakov, 2018; Alquier, Cottet and Lecué, 2019) do achieve rate optimality for prediction and coefficient estimation. However, they do not benefit much from strong signals because the bias does not diminish as signal strengthens. Under the restricted isometry property (on the design) and Gaussian errors, Ndaoud (2019) derived the lower bound for the minimax risk: $\inf_{\widehat{}\bm{\beta}}\sup_{\bm{\beta}^{*}\in\Omega(s,a)}\mathbb{E}\|\widehat{}\bm{\beta}-\bm{\beta}^{*}\|_{2}^{2}\gtrsim\sigma^{2}s/n$ when $a\gtrsim\sigma\sqrt{\log(ed/s)/n}$, where $\Omega(s,a)=\{\bm{\beta}\in\mathbb{R}^{d}:\|\bm{\beta}\|_{0}\leq s,\min_{j:\beta_{j}\neq 0}|\beta_{j}|\geq a\}$. For estimating such sparse vectors with sufficiently strong signals, Lasso can not achieve the oracle rate without the strong irrepresentable condition (Meinshausen and Bühlmann, 2006; Zhao and Yu, 2006), which is a condition on how strongly the important and unimportant variables can be correlated. This condition, however, is in general very restrictive; see Zou (2006) for counterexamples and numerical demonstrations.

In the presence of heavy-tailed noise, outliers occur more frequently and may have a significant impact on (regularized) empirical risk minimization when the loss grows quickly. When the regression error $\varepsilon$ only has finite second moment, the Lasso still achieves the minimax rate $\sqrt{s\log(d)/n}$ (under $\ell_{2}$-norm) but with worse deviations (Lecué and Mendelson, 2018). To reduce the ill-effects of outliers, a widely recognized strategy is to use a robust loss function that is globally Lipschitz continuous and locally quadratic. A prototypical example is the Huber loss (Huber, 1964):

[TABLE]

where $\tau>0$ is a robustification parameter that controls the tradeoff between the robustness and bias. The second important issue is the choice of sparsity-inducing penalty. In order to eliminate the nonnegligible estimation bias introduced by convex regularization, Fan and Li (2001) introduced a family of folded-concave penalties, including the smoothly clipped absolute deviation (SCAD) penalty (Fan and Li, 2001), minimax concave (MC+) penalty (Zhang, 2010a), and the capped $\ell_{1}$-penalty (Zhang, 2010b; Shen, Pan and Zhu, 2012). These ideas motivate the following nonconvex (folded concave) regularized $M$-estimator

[TABLE]

where $\widehat{\mathcal{L}}_{\tau}(\bm{\beta}):=(1/n)\sum_{i=1}^{n}\ell_{\tau}(y_{i}-\bm{x}_{i}^{{\rm T}}\bm{\beta})$ is the empirical loss, $\tau>0$ is a robustification parameter, and $p_{\lambda}:\mathbb{R}\mapsto[0,\infty)$ is a concave penalty function with a regularization parameter $\lambda>0$. We refer to Zhang and Zhang (2012) for a comprehensive survey of folded concave regularized methods.

In practice, it is inherently difficult to solve the nonconvex optimization problem (5) directly. Statistical properties, such as the rate of convergence under various norms and oracle properties, are established for either the hypothetical global optimum that is unobtainable by any practical algorithm in polynomial time, or some local optimum that exists somewhere like a needle in a haystack. To mitigate the gap between statistical theory and algorithmic complexity, we apply an iteratively reweighted $\ell_{1}$-penalized algorithm, which originates from Zou and Li (2008), to adaptive Huber regression. This multi-step regularized robust regression procedure (provably) yields an estimator with desired oracle properties, and is computationally efficient because it only involves solving a sequence of (unconstrained) convex programs. Our theoretical analysis is based on and improves upon Fan et al. (2018), who established the statistical and algorithmic theory for the iteratively reweighted $\ell_{1}$-penalized least squares regression estimator. The aim of this paper is to explore a general class of robust loss functions, typified by the Huber loss, not merely for the purpose of generality but owing to a real downside of the quadratic loss. Typified by the Huber loss, our general principle applies to a class of robust loss functions as will be discussed in Section 4. Software implementing the proposed procedure and reproducing our computational results is available at https://github.com/XiaoouPan/ILAMM.

1.1 Related literature

Nonasymptotic or finite-sample analysis of regularized regression methods beyond least squares, such as regularized empirical risk minimization (ERM) or $M$-estimation with a non-quadratic loss, can be divided into three categories depending on the form of the regularizer/penalty.

$\ell_{1}$-regularization: For high-dimensional sparse linear models, Minsker (2015) and Fan, Li and Wang (2017), respectively, proposed a robust version of Lasso based on geometric median and $\ell_{1}$-penalized Huber’s $M$-estimator. Both estimators achieve sub-Gaussian deviation bounds when the regression error only has finite variance. In such a heavy-tailed case, Lecué and Mendelson (2018) showed that the Lasso still achieves the minimax rate under expectation but with much worse deviations. When the regression error only has finite $(1+\delta)$-th absolute moment for some $\delta\in(0,1)$, Sun, Zhou and Fan (2020) established exponential deviation bounds for $\ell_{1}$-penalized adaptive Huber regression estimator with a more delicate choice of the robustification parameter. For more general penalized $M$-estimators with a convex and Lipschitz continuous loss, Alquier, Cottet and Lecué (2019) established both estimation bounds and sharp oracle inequalities. Their results do not require a local strong convexity on the loss function, thus also including the hinge loss and quantile regression loss. For nonconvex loss functions with a redescending derivative, typified by Tukey’s bisquare loss, Mei, Bai and Montanari (2018) proved the statistical consistency of the $\ell_{1}$-penalized estimator subject to an $\ell_{2}$-constraint to stationary points.

Folded concave regularization: For folded concave penalized $M$-estimators subject to a convex side constraint, Loh and Wainwright (2015) and Loh and Wainwright (2017) were among the first to provide rigorous statistical and algorithmic theory for local optima. They quantified statistical accuracy by providing bounds on $\ell_{1}/\ell_{2}$- and prediction errors between stationary points and the population-level optimum. They also provided conditions under which the stationary point is unique, and proposed a composite gradient algorithm for provably solving the constrained optimization problem efficiently. In the context of generalized linear models with a sufficiently smooth link function and bounded covariates, Loh and Wainwright (2017) proved under the scaling $n\gtrsim s^{3}\log(d)$ that the nonconvex regularized program subject to an $\ell_{1}$-ball constraint has a unique stationary point given by the oracle estimator with high probability. For linear regression with symmetric heavy-tailed errors, Loh (2017) studied statistical consistency and asymptotic normality of nonconvex regularized robust $M$-estimators (also subject to an $\ell_{1}$-ball constraint) with a locally strongly convex loss. For sub-Gaussian covariates, Loh (2017) proved the uniqueness of a stationary point which has $\ell_{2}$- and $\ell_{1}$-error bounds in the order of $\sqrt{s\log(d)/n}$ and $s\sqrt{\log(d)/n}$, respectively. Furthermore, under the scaling $n\gtrsim\max\{s^{2},s\log(d)\}$ and the beta-min condition $\|\bm{\beta}^{*}_{{\mathcal{S}}}\|_{\min}\gtrsim\lambda+\sqrt{\log(s)/n}$, this stationary point coincides with the oracle estimator.

In this paper, we address nonconvex regularized robust regression from a different angle. Motivated by the local linear approximation (LLA) algorithm proposed by Zou and Li (2008), we apply an iteratively reweighted $\ell_{1}$-penalized algorithm to adaptive Huber regression, which involves solving a sequence of (unconstrained) convex programs. We simultaneously analyze the statistical property and algorithmic complexity of the solutions produced by such an iterative procedure. For sub-exponential covariates and asymmetric error with finite variance, we show that the multi-step penalized estimator, after $\mathcal{O}(\log s+\log(\log d))$ iterations, achieves exponential deviation bounds with $\ell_{2}$- and $\ell_{1}$-errors in the order of $\sqrt{s/n}$ and $s/\sqrt{n}$, respectively, under the scaling $n\gtrsim s\log(d)$ and the above beta-min condition. The strong oracle property can be obtained under slightly stronger moment condition and the scaling $n\gtrsim\max\{s^{2},s\log(d)\}$.

$\ell_{0}$-regularization: Another popular class of sparse recovery algorithms is based on directly solving $\ell_{0}$-constrained or $\ell_{0}$-penalized empirical risk minimizations, which naturally produces sparse solutions. Such a formation is NP-hard, and believed to be intractable in practice. Despite its computational hardness, many practically useful algorithms have been proposed to solve $\ell_{0}$-regularized ERM, while the statistical properties beyond least squares regression are much less studied. We refer to Bertsimas, Pauphilet and Van Parys (2020) and Hastie, Tibshirani and Tibshirani (2020) for two comprehensive survey articles on $\ell_{0}$-regularized regression methods.

The idea of having the robustification parameter grow with the sample size in order to achieve exponential deviations even when the sampling distribution only has finite variance dates back to Catoni (2012) in the context of mean estimation. Therefore, the robustness considered in this paper is primarily about nonasymptotic exponential deviation of the estimator versus polynomial tail of the error distribution. The resulting procedure does sacrifice a fair amount of robustness to adversarial contamination of the data. To echo the message in Catoni (2012), the motivation of this work is different from and should not be confused with the classical notion of robust statistics.

From a variable selection perspective, this paper focuses on oracle properties of multi-step penalized robust regression estimators when the signal is sufficiently strong. While allowing for heavy-tailed noise, the high-dimensional feature vector $\bm{x}\in\mathbb{R}^{d}$ is assumed to have either sub-exponential or sub-Gaussian tails. For more complex problems in which both the covariates and noise can be (i) heavy-tailed and/or (ii) adversarially contaminated, the estimator obtained by minimizing a robust loss function is still sensitive to outliers in the feature space. To achieve robustness in both feature and response spaces, recent years have witnessed a rapid development of the “median-of-means” (MOM) principle, which dates back to Nemirovsky and Yudin (1983) and Jerrum, Valiant and Vazirani (1986), and a variety of MOM-based procedures for regression and classification in both low- an high-dimensional settings (Lecué and Lerasle, 2018; Lugosi and Mendelson, 2019a; Chinot, Lecué and Lerasle, 2019, 2020; Lugosi and Mendelson, 2020; Lecué and Lerasle, 2020). We refer to Lugosi and Mendelson (2019b) for a recent survey. An interesting open problem is how to efficiently incorporate the MOM principle with nonconvex regularization or iteratively reweighted $\ell_{1}$-regularization so as to achieve high degree of robustness and variable selection consistency simultaneously.

1.2 Notation

Let us summarize our notation. For every integer $k\geq 1$, we use $\mathbb{R}^{k}$ to denote the $k$-dimensional Euclidean space. The inner and Hadamard products of any two vectors $\bm{u}=(u_{1},\ldots,u_{k})^{{\rm T}},\bm{v}=(v_{1},\ldots,v_{k})^{{\rm T}}\in\mathbb{R}^{k}$ are defined by $\bm{u}^{{\rm T}}\bm{v}=\langle\bm{u},\bm{v}\rangle=\sum_{i=1}^{k}u_{i}v_{i}$ and $\bm{u}\circ\bm{v}=(u_{1}v_{1},\ldots,u_{k}v_{k})^{{\rm T}}$, respectively. We use $\|\cdot\|_{p}$ $(1\leq p\leq\infty)$ to denote the $\ell_{p}$-norm in $\mathbb{R}^{k}$: $\|\bm{u}\|_{p}=(\sum_{i=1}^{k}|u_{i}|^{p})^{1/p}$ and $\|\bm{u}\|_{\infty}=\max_{1\leq i\leq k}|u_{i}|$. Moreover, we write $\|\bm{u}\|_{\min}=\min_{1\leq i\leq k}|u_{i}|$. For $k\geq 2$, $\mathbb{S}^{k-1}=\{\bm{u}\in\mathbb{R}^{k}:\|\bm{u}\|_{2}=1\}$ denotes the unit sphere in $\mathbb{R}^{k}$. For any function $f:\mathbb{R}\mapsto\mathbb{R}$ and vector $\bm{u}=(u_{1},\ldots,u_{k})^{{\rm T}}\in\mathbb{R}^{k}$, we write $f(\bm{u})=(f(u_{1}),\ldots,f(u_{k}))^{{\rm T}}\in\mathbb{R}^{k}$.

For $k\geq 2$, ${\rm I}_{k}$ represents the identity/unit matrix of size $k$. For any $k\times k$ symmetric matrix $\Sigma\in\mathbb{R}^{k\times k}$, $\|\Sigma\|_{2}$ is the operator norm of $\Sigma$, and we use $\lambda_{\min}(\Sigma)$ and $\lambda_{\max}(\Sigma)$ to denote the minimal and maximal eigenvalues of $\Sigma$, respectively. For a positive semidefinite matrix $\Sigma\in\mathbb{R}^{k\times k}$, $\|\cdot\|_{\Sigma}$ denotes the norm linked to $\Sigma$ given by $\|\bm{u}\|_{\Sigma}=\|\Sigma^{1/2}\bm{u}\|_{2}$, $\bm{u}\in\mathbb{R}^{k}$. For any two real numbers $u$ and $v$, we write $u\vee v=\max(u,v)$ and $u\wedge v=\min(u,v)$. For any integer $d\geq 1$, we write $[d]=\{1,\ldots,d\}$. For any set $\mathcal{S}$, we use $|\mathcal{S}|$ to denote its cardinality, i.e., the number of elements in $\mathcal{S}$.

2 Regularized Huber $M$-estimation

We first revisit the $\ell_{1}$-penalized Huber regression estimator in Section 2.1, and point out two different regimes for the robustification parameter $\tau$. In Section 2.2, we propose a multi-step procedure, which is closely related to folded concave regularized Huber regression, for fitting high-dimensional sparse models with heavy-tailed noise. This multi-step penalized robust regression method not only is computationally efficient, but also achieves optimal rate of convergence and oracle properties, as will be studied in Sections 2.3 and 2.4. Throughout, ${\mathcal{S}}={\rm supp}(\bm{\beta}^{*})=\{1\leq j\leq d:\beta_{j}^{*}\neq 0\}\subseteq[d]$ denotes the active set and $s=|{\mathcal{S}}|$ is the sparsity.

2.1 $\ell_{1}$-penalized Huber regression

Given i.i.d. observations $\{(y_{i},\bm{x}_{i})\}_{i=1}^{n}$ from the linear model (1), consider the $\ell_{1}$-regularized Huber $M$-estimator, which we refer to as the Huber-Lasso,

[TABLE]

where $\widehat{\mathcal{L}}_{\tau}(\cdot)$ is the emprical loss function defined in (5). Statistical properties of the penalized Huber $M$-estimator have been studied by Lambert-Lacroix and Zwald (2011), Fan, Li and Wang (2017), Loh (2017) and Alquier, Cottet and Lecué (2019) under different assumptions. A less-noticed problem is the connection between the robustification parameter and the error distribution, which in turn quantifies the tradeoff between robustness and unbiasedness. Recent studies by Sun, Zhou and Fan (2020) reveal that the use of Huber loss is particularly suited for heavy-tailed problems in both low and high dimensions. With a properly chosen robustification parameter, calibrated by the noise level, sample size and parametric dimension, the effects of the heavy-tailed noise can be removed or dampened.

Remark 2.1**.**

In practice, it is natural to leave the intercept or a given subset of the parameters unpenalized in the penalized $M$-estimation framework. Denote by $\mathcal{R}$ be a user-specified index set of unpenalized parameters, which contains at least index 1. A modified Huber-Lasso estimator is then defined as the solution to $\min_{\bm{\beta}\in\mathbb{R}^{d}}\{\widehat{\mathcal{L}}_{\tau}(\bm{\beta})+\lambda\|\bm{\beta}_{\mathcal{R}^{{\rm c}}}\|_{1}\}$, where $\|\bm{\beta}_{\mathcal{R}^{{\rm c}}}\|_{1}=\sum_{j\in\mathcal{R}^{{\rm c}}}|\beta_{j}|$. Similar theoretical analysis can be carried out with slight modifications, and thus will be omitted for ease of exposition.

We first impose the following assumptions on the data generating process. The (random) covaraite vectors are assumed to be sub-exponential/sub-gamma (Boucheron, Lugosi and Massart, 2013), and we allow the regression errors to be heavy-tailed and asymmetric.

Condition 2.1**.**

There exist some constant $\sigma_{\bm{x}},c_{0}\geq 1$ such that $\mathbb{P}(|\bm{u}^{{\rm T}}\bm{x}|\geq\sigma_{\bm{x}}t)\leq c_{0}e^{-t}$ for all $\bm{u}\in\mathbb{S}^{d-1}$ and $t\geq 0$. For simplicity, we set $c_{0}=1$. Moreover, $\Sigma=\mathbb{E}(\bm{x}\bm{x}^{{\rm T}})$ is positive definite with $\rho_{l}=\lambda_{\min}(\Sigma)>0$. The regression error $\varepsilon$ satisfies $\mathbb{E}(\varepsilon|\bm{x})=0$ and $\mathbb{E}(\varepsilon^{2}|\bm{x})\leq\sigma_{2}^{2}$ almost surely.

Theorem 2.1**.**

Assume that Condition 2.1 holds for model (1). For every $t>0$, any optimal solution $\widehat{\bm{\beta}}^{\mbox{{\tiny H-Lasso}}}$ to the convex program (6) with $\tau\asymp\sigma_{2}\sqrt{n/(\log d+t)}$ and $\lambda\asymp\sigma_{2}\sqrt{(\log d+t)/n}$ satisfies

[TABLE]

with probability at least $1-3e^{-t}$ as long as $n\geq c_{3}(s\log d+t)$, where $c_{1}$–$c_{3}$ are constants that are independent of $(n,d,s)$ and $t$.

Theorem 2.1 provides the error bounds for the one-step penalized estimator, and paves the way for our subsequent analysis for the multi-step procedure. Theorem 2.1 is a modified version of Theorem B.2 in Sun, Zhou and Fan (2020) (when $\delta=1$) with an explicit relation between deviation bound and confidence level under slightly relaxed moment condition on the design. When the (conditional) distribution of $\varepsilon$ is symmetric, $\bm{\beta}^{*}$ can be identified as $\bm{\beta}^{*}\in\mathop{\mathrm{argmin}}_{\bm{\beta}\in\mathbb{R}^{d}}\mathbb{E}\widehat{\mathcal{L}}_{\tau}(\bm{\beta})$. Then, with a fixed $\tau$ (e.g. $\tau\asymp\sigma_{2}$), Theorem 2.1 can also be obtained as a special case of Theorem 2.1 in Alquier, Cottet and Lecué (2019) when the feature vector $\bm{x}$ is sub-Gaussian.

2.2 Iteratively reweighted $\ell_{1}$-penalized Huber regression

For fitting sparse regression models, the Lasso-type estimators typically exhibit a suboptimal rate of convergence, as compared to the oracle rate achieved by nonconvex regularization methods, under a minimum signal strength condition (Zhang and Zhang, 2012; Ndaoud, 2019), also known as the beta-min condition (Bühlmann and van de Geer, 2011, Section 7.4). However, as noted previously, directly solving the nonconvex optimization problem (5) is computationally challenging. Moreover, statistical properties are only established for the hypothetical global optimum (or some stationary point), which is typically unobtainable by any polynomial time algorithm.

Inspired by the local linear approximation to folded concave penalties (Zou and Li, 2008), we consider a multi-stage procedure that solves a sequence of convex programs up to a prespecified optimization precision. This is an iteratively reweighted $\ell_{1}$-penalized algorithm, which is similar in spirit to the iteratively reweighted basis-pursuit algorithms studied in Gaïffas and Lecué (2011). Let $p_{\lambda}(\cdot)$ be a differentiable penalty function as in (5) and recall that $\widehat{\mathcal{L}}_{\tau}(\cdot)$ is the empirical loss function. Starting with an initial estimate $\widehat{}\bm{\beta}^{(0)}=(\widehat{\beta}_{1}^{(0)},\ldots,\widehat{\beta}_{d}^{(0)})^{{\rm T}}$, consider a sequence of convex optimization programs $\{({\rm P}_{\ell})\}_{\ell\geq 1}$:

[TABLE]

for $\ell=1,2,\ldots$, where $\widehat{}\bm{\beta}^{(\ell)}=(\widehat{\beta}_{1}^{(\ell)},\ldots,\widehat{\beta}_{d}^{(\ell)})^{{\rm T}}$ is the optimal solution to program $({\rm P}_{\ell})$. Following Zhang and Zhang (2012), we assume the following conditions on the penalty function $p_{\lambda}$.

Condition 2.2**.**

The penalty function $p_{\lambda}$ is of the form $p_{\lambda}(t)=\lambda^{2}p(t/\lambda)$ for $t\in\mathbb{R}$, where $p:\mathbb{R}\mapsto[0,\infty)$ satisfies: (i) $p(t)=p(-t)$ for all $t$ and $p(0)=0$; (ii) $p$ is nondecreasing on $[0,\infty)$; (iii) $p$ is differentiable almost everywhere on $(0,\infty)$ and $\lim_{t\to 0^{+}}p^{\prime}(t)=1$; (iv) $p^{\prime}(t_{1})\leq p^{\prime}(t_{2})$ for all $t_{1}\geq t_{2}>0$.

Prototypical examples of the penalty function $p(\cdot)$ in Condition 2.2 include the $\ell_{1}$-function, the SCAD penalty (Fan and Li, 2001), MC+ penalty (Zhang, 2010a), and capped-$\ell_{1}$ function (Zhang, 2010b).

(SCAD) $p(t)=\int_{0}^{|t|}\min\{1,(1-\frac{u-1}{a-1})_{+}\}{\rm d}u$ and $p^{\prime}(t)=\mathop{\mathrm{sign}}(t)\min\{1,(1-\frac{|t|-1}{a-1})_{+}\}$ for some $a>2$. By a Bayesian argument, Fan and Li (2001) suggested the choice of $a=3.7$. 2. 2.

(MC+) $p(t)=\int_{0}^{|t|}(1-u/a)_{+}{\rm d}u$ and $p^{\prime}(t)=\mathop{\mathrm{sign}}(t)(1-|t|/a)_{+}$ for some $a>1$. 3. 3.

(Capped-$\ell_{1}$) $p(t)=\min(1,|t|)$ and $p^{\prime}(t)=I(|t|\leq 1)$.

For each $\ell\geq 1$, program $({\rm P}_{\ell})$ corresponds to a weighted $\ell_{1}$-penalized empirical Huber loss minimization of the form

[TABLE]

where $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{d})^{{\rm T}}$ is a $d$-vector of regularization parameters with $\lambda_{j}\geq 0$. By convex optimization theory, any optimal solution $\widehat{}\bm{\beta}$ to the convex program (9) satisfies the first-order condition

[TABLE]

where $\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta})=(-1/n)\sum_{i=1}^{n}\ell_{\tau}^{\prime}(y_{i}-\bm{x}_{i}^{{\rm T}}\bm{\beta})\bm{x}_{i}$.

Definition 2.1**.**

Following the terminology in Fan et al. (2018), for a prespecified tolerance level $\epsilon>0$, we say $\widetilde{}\bm{\beta}$ is an $\epsilon$-optimal solution to (9) if $\omega_{\bm{\lambda}}(\widetilde{}\bm{\beta})\leq\epsilon$, where

[TABLE]

In view of Definition 2.1, for a prespecified sequence of tolerance levels $\{\epsilon_{\ell}\}_{\ell\geq 1}$, we use $\widetilde{}\bm{\beta}^{(\ell)}=(\widetilde{\beta}_{1}^{(\ell)},\ldots,\widetilde{\beta}_{d}^{(\ell)})^{{\rm T}}$ to denote an $\epsilon_{\ell}$-optimal solution to program $({\rm P}_{\ell})$, that is,

[TABLE]

where $\bm{\lambda}^{(\ell-1)}:=p_{\lambda}^{\prime}(\widetilde{}\bm{\beta}^{(\ell-1)})$. For simplicity, we consider a trivial initial estimator $\widetilde{}\bm{\beta}^{(0)}=\textbf{0}$. Since $p_{\lambda}^{\prime}(|\widetilde{\beta}_{j}^{(0)}|)=p^{\prime}_{\lambda}(0)=\lambda$ for $j=1,\ldots,d$, the program $({\rm P}_{1})$ coincides with that in (6). In Section 3, we will describe an iterative local adaptive majorize-minimization (I-LAMM) algorithm which produces an $\epsilon$-optimal solution to (9) after a few iterations.

The above procedure is sequential, and can be categorized into two stages: contraction ($\ell=1$) and tightening ($\ell\geq 2$). As we will see in the next subsection, even starting with a trivial initial estimator that is fairly remote from the true parameter, the contraction stage will produce a reasonably good estimator whose statistical error is of the order $\sqrt{\log(d)\cdot s/n}$. Essentially, the contraction stage is equivalent to the $\ell_{1}$-penalized Huber regression in (6). A tightening stage further refines this coarse contraction estimator consecutively, and eventually gives rise to an estimator that achieves the oracle rate $\sqrt{s/n}$ under a weak beta-min condition.

2.3 Deterministic analysis

To analyze the statistical properties of $\{\widetilde{}\bm{\beta}^{(\ell)}\}_{\ell\geq 1}$, we first define a “good” event regarding the restricted strong convexity (RSC) property of the empirical Huber loss over a local $\ell_{1}$-cone.

Definition 2.2**.**

For some $r,l,\kappa>0$, define the event

[TABLE]

where $\mathbb{B}(r)=\mathbb{B}^{d}(r)=\{\bm{\delta}\in\mathbb{R}^{d}:\|\bm{\delta}\|_{2}\leq r\}$ is an $\ell_{2}$-ball and $\mathbb{C}(l):=\{\bm{\delta}\in\mathbb{R}^{d}:\|\bm{\delta}\|_{1}\leq l\|\bm{\delta}\|_{2}\}$ is an $\ell_{1}$-cone. Here $\bm{\beta}+\mathbb{B}(r)\cap\mathbb{C}(l)=\{\bm{\beta}+\bm{\delta}:\bm{\delta}\in\mathbb{B}(r)\cap\mathbb{C}(l)\}$.

Throughout the following, we assume that the penalty function $p_{\lambda}(\cdot)$ satisfies Condition 2.2. Moreover, define

[TABLE]

where $\mathcal{L}_{\tau}(\bm{\beta})=\mathbb{E}\widehat{\mathcal{L}}_{\tau}(\bm{\beta})$ is the population loss. Here $\bm{w}^{*}\in\mathbb{R}^{d}$ is the centered gradient vector which corresponds to the stochastic error, and $b^{*}_{\tau}$ denotes the (deterministic) approximation bias induced by the Huber loss. See Lemma 6.1 in the Supplementary Material for an upper bound on the bias.

Remark 2.2**.**

In this paper, we introduce the bias term $b^{*}_{\tau}$ into the results primarily because the error distribution, if not specified, can be asymmetric. This term is typically nullified in the literature due to two reasons. First, under the symmetry assumption that $\varepsilon$ (conditional on $\bm{x}$) is symmetric around zero, then for any given $\tau>0$, $\mathbb{E}\{\ell^{\prime}_{\tau}(\varepsilon)\bm{x}\}=0$. Secondly, it is sometimes assumed that $\bm{x}_{i}$ and $\varepsilon_{i}$ are independent, and both have zero means. Again, for any $\tau>0$, it follows that $\mathbb{E}\{\ell^{\prime}_{\tau}(\varepsilon)\bm{x}\}=\mathbb{E}\{\ell^{\prime}_{\tau}(\varepsilon)\}\cdot\mathbb{E}(\bm{x})=0$. In these two scenarios, the bias $b^{*}_{\tau}$ vanishes for any given $\tau$.

Proposition 2.1**.**

Let $\lambda,r,\kappa>0$ satisfy

[TABLE]

Then, conditioned on the event $\mathcal{E}_{1}(r,l,\kappa)\cap\{\lambda\geq 2(\|\bm{w}^{*}\|_{\infty}+\epsilon_{1})\}$ with $l=6s^{1/2}$, any $\epsilon_{1}$-optimal solution $\widetilde{}\bm{\beta}^{(1)}$ of program $({\rm P}_{1})$ satisfies

[TABLE]

Proposition 2.1 is deterministic in the sense that the error bound (14) holds conditioning on the event $\mathcal{E}_{1}(r,l,\kappa)\cap\{\lambda\geq 2(\|\bm{w}^{*}\|_{\infty}+\epsilon_{1})\}$. Under Condition 2.1 (sub-exponential design and heavy-tailed error with finite variance), we will establish the delicate choices of $\lambda,\epsilon_{1},r$ and sample size requirement in order that this event occurs with high probability. Specifically, we will show that

[TABLE]

Next, we investigate the statistical properties of $\{\widetilde{}\bm{\beta}^{(\ell)}\}_{\ell\geq 2}$ in the tightening stage. We impose a minimum signal strength condition on $\|\bm{\beta}^{*}_{{\mathcal{S}}}\|_{\min}=\min_{j\in{\mathcal{S}}}|\beta^{*}_{j}|$, so that the error rate obtained in Proposition 2.1 is improvable (Zhang and Zhang, 2012; Ndaoud, 2019). Recall that $s=|{\mathcal{S}}|$.

Proposition 2.2**.**

Assume there exists some $\gamma>0$ such that $p^{\prime}(\gamma)>0$. Let

[TABLE]

and choose $c>0$ so that

[TABLE]

Set $l=\{2+\frac{2}{p^{\prime}(\gamma)}\}(c^{2}+1)^{1/2}s^{1/2}+\frac{2}{p^{\prime}(\gamma)}s^{1/2}$ and let $r>0$ satisfy

[TABLE]

Under the minimum signal strength condition $\|\bm{\beta}^{*}_{{\mathcal{S}}}\|_{\min}\geq\gamma\lambda$, and conditioned on event $\mathcal{E}_{1}(r,l,\kappa)\cap\{\lambda\geq\frac{2}{p^{\prime}(\gamma)}(\|\bm{w}^{*}\|_{\infty}+\max_{\ell\geq 1}\epsilon_{\ell})\}$, the $\epsilon_{\ell}$-optimal solutions $\widetilde{}\bm{\beta}^{(\ell)}$ $(\ell\geq 2)$ satisfy

[TABLE]

where $\delta=0.5p^{\prime}(\gamma)/(\gamma\kappa)\in(0,1)$. Furthermore, it holds

[TABLE]

Proposition 2.2 unveils how the tightening stage improves the statistical rate: every tightening step shrinks the estimation error from the previous step by a $\delta$-fraction. The second term on the right-hand side of (2.2) or (19) dominates the $\ell_{2}$-error, and up to constant factors, consists of three components,

[TABLE]

We identify $\|p_{\lambda}^{\prime}(|\bm{\beta}^{*}_{\mathcal{S}}|-\gamma\lambda)\|_{2}$ as the shrinkage bias induced by the penalty function. This explains the limitation of the $\ell_{1}$-penalty $p_{\lambda}(t)=\lambda|t|$ whose derivative $p_{\lambda}^{\prime}(t)=\lambda\mathop{\mathrm{sign}}(t)$ ($t\neq 0$) does not vanish regardless of the signal strength. Intuitively, choosing a proper penalty function $p_{\lambda}(\cdot)$ with a descending derivative reduces the bias as signal strengthens. The second term, $\|\bm{w}^{*}_{{\mathcal{S}}}\|_{2}+b^{*}_{\tau}$, reveals the oracle property. To see this, consider the oracle estimator defined as

[TABLE]

Since $s=|{\mathcal{S}}|\ll n$, the finite sample theory for Huber’s $M$-estimation in low dimensions (Sun, Zhou and Fan, 2020) applies to $\widehat{}\bm{\beta}^{{\rm ora}}$, indicating that with high probability,

[TABLE]

According to Definition 2.1, the last term $s^{1/2}\epsilon_{\ell}$ demonstrates the optimization error, which will be discussed in Section 3.

The above results provide conditions under which the sequence of estimators $\{\widetilde{}\bm{\beta}^{(\ell)}\}_{\ell\geq 1}$ satisfy the contraction property and, meanwhile, fall in a local neighborhood of $\bm{\beta}^{*}$. Another important feature of the proposed procedure is that the resulting estimator satisfies the strong oracle property, as demonstrated by the following result. Let $\{\widehat{}\bm{\beta}^{(\ell)}\}_{\ell\geq 1}$ be any optimal solutions to the convex programs $\{({\rm P}_{\ell})\}_{\ell\geq 1}$ in (8) with $\widehat{}\bm{\beta}^{(0)}=\textbf{0}$. Similarly to Definition 2.2, we define the following event in regard of the restricted strong convexity of the empirical Huber loss. For some $r,l,\kappa>0$,

[TABLE]

where $\mathcal{C}(r,l)=\{(\bm{\beta}_{1},\bm{\beta}_{2}):\bm{\beta}_{1}\in\bm{\beta}_{2}+\mathbb{B}(r)\cap\mathbb{C}(l),\bm{\beta}_{2}\in\bm{\beta}^{*}+\mathbb{B}_{\Sigma}(r),{\rm supp}(\bm{\beta}_{2})\subseteq{\mathcal{S}}\}$ and $\mathbb{B}_{\Sigma}(r)=\{\bm{\beta}\in\mathbb{R}^{d}:\|\bm{\beta}\|_{\Sigma}\leq r\}$. Moreover, define the “oracle” score $\bm{w}^{{\rm ora}}\in\mathbb{R}^{d}$ as

[TABLE]

which satisfies $\bm{w}^{{\rm ora}}_{{\mathcal{S}}^{{\rm c}}}=\textbf{0}$.

Proposition 2.3**.**

Suppose there exist constants $\gamma_{1}>\gamma_{0}>0$ such that $p^{\prime}(\gamma_{0})\in(0,1/2]$, $p^{\prime}(\gamma_{1})=0$. For a prespecified $\delta\in(0,1)$, let $\kappa\geq 1.25/(\delta\gamma_{0})$ and choose $c_{0}>0$ so that

[TABLE]

Moreover, set $l=\{2+\frac{2}{p^{\prime}(\gamma_{0})}\}(c_{0}^{2}+1)^{1/2}s^{1/2}$ and let $r\geq c_{0}\gamma_{0}s^{1/2}\lambda$. Then, conditioned on the event

[TABLE]

the strong oracle property holds under the minimum signal strength condition $\|\bm{\beta}^{*}_{{\mathcal{S}}}\|_{\min}\geq(\gamma_{0}+\gamma_{1})\lambda$: $\widehat{}\bm{\beta}^{(\ell)}=\widehat{}\bm{\beta}^{{\rm ora}}$ for all $\ell\geq\lceil\log(s^{1/2}/\delta)/\log(1/\delta)\rceil$.

The proofs of Propositions 2.1, 2.2 and 2.3 are provided in the Supplementary Material.

2.4 Random analysis

In this section, we complement the previous deterministic results with probabilistic bounds on the random events of interest. To be more specific, events $\mathcal{E}_{1}(r,l,\kappa)$ and $\mathcal{E}_{2}(r,l,\kappa)$ correspond to the RSC properties of $\widehat{\mathcal{L}}_{\tau}(\cdot)$. The order of the regularization parameter $\lambda$ depends on $\|\bm{w}^{*}\|_{\infty}$, where $\bm{w}^{*}=\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta}^{*})-\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})$ is the centered score function evaluated at $\bm{\beta}^{*}$. The oracle convergence rate depends on the $\ell_{2}$-norm of $\bm{w}^{*}_{{\mathcal{S}}}\in\mathbb{R}^{s}$, the subvector of $\bm{w}^{*}$ indexed by ${\mathcal{S}}$.

Under Condition 2.1, $\bm{x}=(x_{1},\ldots,x_{d})^{{\rm T}}$ is sub-exponential and $\Sigma=\mathbb{E}(\bm{x}\bm{x}^{{\rm T}})=(\sigma_{jk})_{1\leq j,k\leq d}$ is positive definite. Here we do not require the components of $\bm{x}$ to have zero means. Moreover, given the true active set ${\mathcal{S}}\subseteq[d]$ of $\bm{\beta}^{*}$, we define the following $s\times s$ principal submatrix of $\Sigma$:

[TABLE]

Throughout, “$\lesssim$” and “$\gtrsim$” stand for “$\leq$” and “$\geq$”, respectively, up to constants that are independent of $(n,d,s)$ but might depend on those in Condition 2.1. In particular, define

[TABLE]

which is a constant depending only on $\sigma_{\bm{x}}$.

Proposition 2.4**.**

Assume Condition 2.1 holds, and let $\rho_{u}=\lambda_{\max}(\bm{\Sigma})$. Then, for any $t\geq 0$,

[TABLE]

holds with probability at least $1-e^{-t}$ as long as $\tau\geq\max\{C_{1}\sigma_{2},C_{2}r\}$ and $n\geq\log(2d)$, where $C_{0},C_{1}$ are absolute constants and $C_{2}$ depends only on $\sigma_{\bm{x}}$.

The next proposition provides high probability bounds on $\|\bm{w}^{*}\|_{\infty}$ and $\|\bm{w}^{*}_{{\mathcal{S}}}\|_{2}$, where $\bm{w}^{*}=\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta}^{*})-\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})$.

Proposition 2.5**.**

Assume Condition 2.1 holds. For any $t>0$, the centered score $\bm{w}^{*}\in\mathbb{R}^{d}$ satisfies

[TABLE]

with probability at least $1-e^{-t}$, and

[TABLE]

with probability at least $1-e^{-t}$.

Similarly to Theorem 2.1, Propositions 2.4 and 2.5 are also modified versions of Lemmas C.4 and C.6 in Sun, Zhou and Fan (2020) under a weaker sub-exponential condition on the feature vector $\bm{x}$. Therefore in the proofs, we only provide the necessary steps that help improve upon the existing results. Together, Propositions 2.4 and 2.5 reveal the impact of the robustification parameter on the statistical properties of the resulting estimator. As discussed in Section 2.3 above, the order of $\|\bm{w}(\bm{\beta}^{*})_{{\mathcal{S}}}\|_{2}$ determines the oracle rate of convergence. In Theorem 2.2, we show that after only a small number of iterations, the proposed procedure leads to an estimator that achieves the oracle rate of convergence. Recall from Section 2.2 that $\{\widetilde{}\bm{\beta}^{(\ell)}\}_{\ell=1,2,\ldots}$ is a sequence of $\epsilon_{\ell}$-optimal solutions of the convex programs (8), initialized at $\widetilde{}\bm{\beta}^{(0)}=\textbf{0}$.

Theorem 2.2**.**

Assume Conditions 2.1 and 2.2 hold, and there exist some $\gamma_{1}>\gamma_{0}>0$ such that

[TABLE]

Given $t\geq 0$, suppose the sample size satisfies $n\gtrsim s\log d+t$, and $\epsilon_{\ell}\leq\sqrt{1/n}$ for all $\ell\geq 1$. Moreover, suppose that we choose a regularization parameter $\lambda\asymp\sigma_{2}\sqrt{(\log d+t)/n}$, and let $\tau$ satisfy

[TABLE]

Then, under the minimum signal strength condition $\|\bm{\beta}^{*}_{{\mathcal{S}}}\|_{\min}\geq(\gamma_{0}+\gamma_{1})\lambda$, the multi-stage estimator $\widetilde{}\bm{\beta}^{(T)}$ with $T\gtrsim\lceil\frac{\log(\log d+t)}{\log(1/\delta)}\rceil$ satisfies the bounds

[TABLE]

with probability at least $1-3e^{-t}$.

We refer to the conclusion of Theorem 2.2 as the weak oracle property in the sense that the proposed estimator achieves the same rate of convergence as the oracle $\widehat{}\bm{\beta}^{{\rm ora}}$ which knows a priori the support ${\mathcal{S}}$ of $\bm{\beta}^{*}$. We keep the two terms $\tau(s+t)/n$ and $b^{*}_{\tau}$ in the upper bounds of (32) to keep track the impact of $\tau$ on the estimator error: the former is part of the stochastic error and the latter characterizes the bias. Below are two cases that are of general interests.

(Symmetry) As discussed in Remark 2.2, if $\varepsilon$ (conditional on $\bm{x}$) is symmetric around zero, then $b^{*}_{\tau}=0$ for any $\tau>0$. To certify (2.2), $\tau$ can be taken as a constant-multiple of $\sigma_{2}$, and the resulting error bounds become

[TABLE]

with probability at least $1-3e^{-t}$. 2. 2.

(Asymmetry) When the conditional distribution of $\varepsilon_{i}$ is asymmetric, there will be a bias-robustness tradeoff. If $\varepsilon$ only has bounded second moment, by Lemma 6.1 in the Supplementary Material we have $b^{*}_{\tau}\lesssim\sigma_{2}^{2}\tau^{-1}$ although $\tau b^{*}_{\tau}\to 0$ as $\tau\to\infty$. Then, the multi-step iterative estimator $\widetilde{}\bm{\beta}^{(T)}$ with $\tau\asymp\sigma_{2}\sqrt{n/(s+\log d+t)}$ satisfies, under the scaling $n\gtrsim s\log d+t$, that

[TABLE]

with probability at least $1-3e^{-t}$.

A more intriguing result, as revealed by the following theorem, is that our estimator achieves the strong oracle property, namely, it coincides with the oracle with high probability. Here we need slightly stronger moment conditions than those in Condition 2.1, that is, the random predictor $\bm{x}$ is sub-Gaussian and the noise variable $\varepsilon$ satisfies an $L_{2+\eta}$-$L_{2}$ norm equivalence for some $\eta\in(0,1]$.

Condition 2.3**.**

There exists $\sigma_{\bm{x}}\geq 1$ such that $\mathbb{P}(|\bm{u}^{{\rm T}}\bm{x}|\geq\sigma_{\bm{x}}t)\leq 2e^{-t^{2}/2}$ for all $\bm{u}\in\mathbb{S}^{d-1}$ and $t\geq 0$. Moreover, $\Sigma=(\sigma_{jk})=\mathbb{E}(\bm{x}\bm{x}^{{\rm T}})$ satisfies $\rho_{l}=\lambda_{\min}(\Sigma)>0$ and

[TABLE]

for some $A_{0}>0$. The random error $\varepsilon$ satisfies $\mathbb{E}(\varepsilon|\bm{x})=0$ and $\mathbb{E}(\varepsilon^{2}|\bm{x})\leq\sigma_{2}^{2}$, and $\{\mathbb{E}(|\varepsilon|^{2+\eta}|\bm{x})\}^{1/(2+\eta)}\leq a_{\eta}\{\mathbb{E}(\varepsilon^{2}|\bm{x})\}^{1/2}$ (almost surely) for some $\eta\in(0,1]$ and $a_{\eta}>1$. Moreover, $\varepsilon$ satisfies the anti-concentration property: there exists a constant $a_{0}>0$ such that

[TABLE]

Theorem 2.3**.**

Assume Conditions 2.2 and 2.3 hold, and there exist some $\gamma_{1}>\gamma_{0}>2.5/\rho_{l}$ such that $p^{\prime}(\gamma_{0})>0$, $p^{\prime}(t)=0$ for all $t\geq\gamma_{1}$. For any $t\geq 0$ and $q\geq\max(s,\log d)$, let $\lambda\asymp\sigma_{2}\sqrt{(\log d+t)/n}$ and $\tau\asymp\sigma_{2}\sqrt{n/(q+t)}$. Moreover, assume the sample size satisfies $n\gtrsim\max\{s\log d+t,(q+t)^{1+1/\eta}(\log d)^{-1/\eta}\}$, and the beta-min condition $\|\bm{\beta}^{*}_{{\mathcal{S}}}\|_{\min}\geq(\gamma_{0}+\gamma_{1})\lambda$. Then, with probability at least $1-8e^{-t}$, $\widehat{}\bm{\beta}^{(\ell)}=\widehat{}\bm{\beta}^{{\rm ora}}$ provided $\ell\geq\lceil\log(s^{1/2}/\delta)/\log(1/\delta)\rceil$, where $\delta:=2.5/(\rho_{l}\gamma_{0})\in(0,1)$.

Theorem 2.3 provides a useful complement to Theorem 2 in Loh (2017), and differs from it in two aspects. First, the latter studies the estimator obtained by solving the folded concave penalized optimization program in (5) subject to an $\ell_{1}$-ball constraint in order to ensure the existence of local/global optima. Secondly, Theorem 2 in Loh (2017) establishes the strong oracle property for any stationary point $\widetilde{}\bm{\beta}$ of the program (5) (with an $\ell_{1}$-ball constraint) that falls inside a local neighborhood of $\bm{\beta}^{*}$. In contrast, Theorem 2.3 concerns the strong oracle property of the proposed iteratively reweighted $\ell_{1}$-penalized estimator obtained by solving a sequence of (unconstrained) convex programs (8).

Remark 2.3**.**

A direct consequence of the strong oracle property is variable selection consistency, saying that

[TABLE]

In particular, assume Condition 2.3 holds with $\eta=1$, implying that $\varepsilon$ satisfies an $L_{3}$-$L_{2}$ norm equivalence. Then, Theorem 2.3 implies that the multi-step estimator $\widehat{}\bm{\beta}^{(\ell)}$ with $\lambda\asymp\sigma_{2}\sqrt{(\log d)/n}$, $\tau\asymp\sigma_{2}\sqrt{n/(s+\log d)}$ and $\ell\asymp\log s$ achieves variable selection consistency as $n,d\to\infty$ under the scaling $n\gtrsim\max(s\log d,s^{2})$ and the necessary beta-min condition $\|\bm{\beta}^{*}_{{\mathcal{S}}}\|_{\min}\gtrsim\sigma_{2}\sqrt{(\log d)/n}$ (Ndaoud, 2019).

As discussed earlier, Lasso (Tibshirani, 1996) achieves desirable risk properties, in terms of both estimation and prediction, under mild conditions, yet its variable selection consistency requires much stronger assumptions Meinshausen and Bühlmann (2006); Zhao and Yu (2006); Wainwright (2009). In addition to sub-Gaussian errors, it requires a stronger beta-min condition—$\|\bm{\beta}^{*}_{{\mathcal{S}}}\|_{\min}\gtrsim\sigma_{2}\sqrt{(s\log d)/n}$, and the irrepresentable condition

[TABLE]

See, for example, Chapter 7 in Bühlmann and van de Geer (2011) and Section 7.5 in Wainwright (2019).

3 Optimization Algorithm

In this section, we use the local adaptive majorize-minimize (LAMM) principal (Fan et al., 2018) to derive an iterative algorithm for solving each subproblem $({\rm P}_{\ell})$ in (8):

[TABLE]

where $\bm{\lambda}^{(\ell-1)}=(\lambda^{(\ell-1)}_{1},\ldots,\lambda^{(\ell-1)}_{d})^{{\rm T}}\in\mathbb{R}^{d}$ with $\lambda^{(\ell-1)}_{j}\geq 0$. Specifically, $\lambda^{(\ell-1)}_{j}=0$ for some $j$ means that the $j$-th coefficient is not penalized.

3.1 LAMM algorithm

To minimize a nonlinear function $f(\cdot)$ on $\mathbb{R}^{d}$, at a given point $\bm{\beta}^{(k)}$, the majorize-minimize (MM) algorithm first majorizes it by another function $g(\cdot\,|\bm{\beta}^{(k)})$, which satisfies

[TABLE]

and then compute $\bm{\beta}^{(k+1)}:=\mathop{\mathrm{argmin}}_{\bm{\beta}\in\mathbb{R}^{d}}g(\bm{\beta}|\bm{\beta}^{(k)})$ (Lange, Hunter and Yang, 2000). The objective value of such an algorithm is non-increasing in each step, because

[TABLE]

where inequality (i) is due to the marization property of $g(\cdot|\bm{\beta}^{(k)})$ and inequality (ii) follows from the definition $\bm{\beta}^{(k+1)}$. Fan et al. (2018) observed that the global majorization requirement is not necessary. It only requires the local properties

[TABLE]

for the inequalities in (36) to hold.

Using the above principle, it suffices to locally majorize the objective function $\widehat{\mathcal{L}}_{\tau}(\bm{\beta})$ in the penalized optimization problem. At the $k$-th step with working parameter vector $\bm{\beta}^{(\ell,k-1)}$, we use an isotropic quadratic function, that is,

[TABLE]

to locally majorize $\widehat{\mathcal{L}}_{\tau}(\bm{\beta})$ such that

[TABLE]

where $\phi^{(\ell,k)}$ is a proper quadratic coefficient at the $k$-th update, and $\bm{\beta}^{(\ell,k)}$ is the solution to

[TABLE]

It is easy to see that $\bm{\beta}^{(\ell,k)}$ takes a simple explicit form

[TABLE]

where $S_{\textnormal{soft}}(\bm{\beta},\bm{\lambda}):=(\mathop{\mathrm{sign}}(\beta_{j})\max\{|\beta_{j}|-\lambda_{j},0\})_{j=1,\ldots,d}$ is the soft-thresholding operator. For simplicity, we summarize and define the above update as $\bm{\beta}^{(\ell,k)}=T_{\bm{\lambda}^{(\ell-1)},\phi^{(\ell,k)}}(\bm{\beta}^{(\ell,k-1)})$. Using this simple update formula of $\bm{\beta}$, we iteratively search for the pair $(\phi^{(\ell,k)},\bm{\beta}^{(\ell,k)})$ that ensures the local majorization (39). Starting with an initial quadratic coefficient $\phi=\phi_{0}$, say $10^{-4}$, we iteratively increase $\phi$ by a factor of $\gamma_{u}>1$ and compute

[TABLE]

until the local property (39) holds. This routine is summarized in Algorithm 1.

3.2 Complexity theory

To investigate the complexity theory of the proposed algorithm, we first impose the following standard regularity conditions on the objective function.

Condition 3.1**.**

$\nabla\widehat{\mathcal{L}}_{\tau}(\cdot)$ is $L$-Lipschitz continuous for some $L>0$, that is, $\|\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta}_{1})-\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta}_{2})\|_{\infty}\leq L\|\bm{\beta}_{1}-\bm{\beta}_{2}\|_{2}$ for any $\bm{\beta}_{1},\bm{\beta}_{2}\in\mathbb{R}^{d}$.

Our next theorem characterizes the computational complexity in the contraction stage. Recall that $\bm{\lambda}^{(0)}=(\lambda,\ldots,\lambda)^{{\rm T}}\in\mathbb{R}^{d}$.

Theorem 3.1**.**

Assume Condition 3.1 holds and the optimal solution $\widehat{}\bm{\beta}^{(1)}$ satisfies $\|\widehat{}\bm{\beta}^{(1)}-\bm{\beta}^{*}\|_{2}\lesssim s^{1/2}\lambda$. Then, to attain an $\epsilon_{{\rm c}}$-optimal solution $\widetilde{}\bm{\beta}^{(1)}$, i.e. $\omega_{\bm{\lambda}^{(0)}}(\widetilde{}\bm{\beta}^{(1)})\leq\epsilon_{{\rm c}}$, in the contraction stage, we need as many as $C_{1}L^{2}(1+\gamma_{u})^{2}(\|\bm{\beta}^{*}\|_{2}+s^{1/2}\lambda)^{2}/{\epsilon^{2}_{{\rm c}}}$ LAMM iterations in (40), where $C_{1}>0$ is a constant independent of $(n,d,s)$.

The sublinear rate in the contraction stage is due to the lack of global strong convexity of the loss function in this stage, because we start with a naive initial value ${\bf 0}$. Once we enter the contracting region where the estimator is relatively closer to the underlying true parameter vector, the problem becomes strongly convex (at least with high probability). This endows the algorithm a linear rate of convergence. Our next theorem provides a formal statement on the geometric convergence rate of LAMM for solving each subproblem in the tightening stage. To this end, we describe a variant of the sparse eigenvalue condition.

Definition 3.1** (LSE—Localized Sparse Eigenvalue).**

Given $r,\tau>0$ and an integer $m\geq|{\mathcal{S}}|$, the localized sparse eigenvalues are defined as

[TABLE]

where $\mathbb{C}_{0}(m):=\{\bm{u}\in\mathbb{S}^{d-1}:{\mathcal{S}}\subseteq{\rm supp}(\bm{u}),\,|{\rm supp}(\bm{u})|\leq m\}$ denotes a sparse cone.

Condition 3.2**.**

We say an LSE$(C_{0})$ condition holds for some $C_{0}\geq 1$ if there exist an integer $s^{\prime}\lesssim s$ and constants $\kappa^{*},\kappa_{*},C_{1}>0$ such that

[TABLE]

Note that if a vector $\bm{u}\in\mathbb{R}^{d}$ belongs to the sparse cone $\mathbb{C}_{0}(m)$ for some $m\geq 1$, by Cauchy-Schwarz inequality we have $\|\bm{u}\|_{1}\leq m^{1/2}\|\bm{u}\|_{2}$. This implies that $\bm{u}$ also falls into the $\ell_{1}$-cone $\mathbb{C}(m^{1/2})$ defined in (11). Proposition 2.4 will remain valid, possibly with different constants, if the $\ell_{1}$-cone $\mathbb{C}(l)$ therein is replaced by a sparse cone. Note that Proposition 2.4 controls the minimum LSE. Similar results can be obtained to bound the maximum LSE from above.

Theorem 3.2**.**

Assume LSE$(C_{0})$ condition holds for a sufficiently large $C_{0}>1$ and $\tau\gtrsim r\gtrsim\lambda\sqrt{s}$. To obtain an $\epsilon_{{\rm t}}$-optimal solution $\widetilde{\bm{\beta}}^{(\ell)}$, i.e. $\omega_{\bm{\lambda}^{(\ell-1)}}(\widetilde{\bm{\beta}}^{(\ell)})\leq\epsilon_{{\rm t}}$, in the $\ell$-th subproblem for $\ell\geq 2$, we need as many as $C_{1}\log(C_{2}s^{1/2}{\lambda}/{\epsilon_{{\rm t}}})$ LAMM iterations in (40), where $C_{1}$ and $C_{2}$ are positive constants.

We summarize the above two theorems in the following result, which characterizes the computational complexity of the whole algorithm.

Corollary 3.1**.**

Assume that the conditions in Theorems 3.1 and 3.2 hold. To achieve a sequence of approximate solutions $\{\widetilde{\bm{\beta}}^{(\ell)}\}_{\ell=1}^{T}$ such that $\omega_{\bm{\lambda}^{(0)}}(\widetilde{\bm{\beta}}^{(1)})\leq\epsilon_{{\rm c}}\lesssim\lambda$ and $\omega_{\bm{\lambda}^{(\ell-1)}}(\widetilde{\bm{\beta}}^{(\ell)})\leq\epsilon_{{\rm t}}\lesssim\sqrt{1/n}$ for $2\leq\ell\leq T$, the required number of LAMM iterations is of the order $C_{1}\epsilon_{{\rm c}}^{-2}+C_{2}(T-1)\log(\epsilon_{{\rm t}}^{-1})$, where $C_{1}$ and $C_{2}$ are positive constants independent of $(n,d,s)$.

4 Extension to General Robust Losses

Thus far, we have restricted our attention to the Huber loss. As a representative robust loss function, the Huber loss has the merit of being (i) globally $\tau$-Lipschitz continuous, and (ii) locally quadratic. A natural question arises that whether similar results, both statistical and computational, remain valid for more general loss functions that possess the above two features. In this section, we introduce a class of loss functions which, combined with folded concave regularization, leads to statistically optimal estimators that are robust against heavy-tailed errors.

Condition 4.1** (Globally Lipschitz and locally quadratic loss functions).**

Consider a general loss function $\ell_{\tau}(\cdot)$ that is of the form $\ell_{\tau}(x)=\tau^{2}\ell(x/\tau)$ for $x\in\mathbb{R}$, where $\ell:\mathbb{R}\mapsto[0,\infty)$ is convex and satisfies: (i) $\ell^{\prime}(0)=0$ and $|\ell^{\prime}(x)|\leq c_{1}$ for all $x\in\mathbb{R}$; (ii) $\ell^{\prime\prime}(0)=1$ and $\ell^{\prime\prime}(x)\geq c_{2}$ for all $|x|\leq c_{3}$; and (iii) $|\ell^{\prime}(x)-x|\leq c_{4}x^{2}$ for all $x\in\mathbb{R}$, where $c_{1}$–$c_{4}$ are positive constants.

Note that Condition 4.1 excludes some important Lipschitz continuous functions, such as the check function for quantile regression and the hinge loss for classification, which do not have a local strong convexity. The recent works Alquier, Cottet and Lecué (2019), Chinot, Lecué and Lerasle (2019) and Chinot, Lecué and Lerasle (2020) established optimal estimation and excess risk bounds for (regularized) empirical risk minimizers and MOM-type estimators based on general convex and Lipschitz loss functions even without a local quadratic behavior. Our work complements the existing results on $\ell_{1}$-regularized ERM by showing oracle properties of nonconvex regularized methods under stronger signals. For this reason, we need an additional local strong convexity condition on the loss. It remains unclear whether the oracle rates or variable selection consistency can still be achieved without such a local curvature of the loss function.

We now discuss the implications of the three properties in Condition 4.1. First, since $\ell_{\tau}^{\prime}(x)=\tau\ell^{\prime}(x/\tau)$, it follows from property (i) that $\sup_{x\in\mathbb{R}}|\ell_{\tau}^{\prime}(x)|\leq c_{1}\tau$. The boundedness of $|\ell_{\tau}^{\prime}|$ facilitates the use of Bernstein’s inequality on deriving upper bounds for the random quantities $\|\bm{w}^{*}\|_{\infty}$ and $\|\bm{w}^{*}_{{\mathcal{S}}}\|_{2}$ as in Proposition 2.5, where $\bm{w}^{*}=\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta}^{*})-\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})$ with $\widehat{\mathcal{L}}_{\tau}(\bm{\beta})=(1/n)\sum_{i=1}^{n}\ell_{\tau}(y_{i}-\bm{x}_{i}^{{\rm T}}\bm{\beta})$. Next, note that $\ell_{\tau}^{\prime\prime}(x)=\ell^{\prime\prime}(x/\tau)$. Property (ii) indicates that $\ell_{\tau}$ is strongly convex on $[-c_{3}\tau,c_{3}\tau]$, which turns out to be the key factor in establishing the restricted strong convexity condition on $\widehat{\mathcal{L}}_{\tau}$. See Proposition 2.4 and Lemma 8.1. Lastly, property (iii) is particularly useful when the error distribution is asymmetric. Even though it can be shown under property (i) that $\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta}^{*})$ is concentrated around its expected value $\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})$ with high probability, $\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})=-\mathbb{E}\{\ell^{\prime}_{\tau}(\varepsilon)\bm{x}\}$ is typically nonzero when the conditional distribution of $\varepsilon$ is asymmetric. However, since $\mathbb{E}(\varepsilon|\bm{x})=0$, we have $\mathbb{E}\{\ell_{\tau}^{\prime}(\varepsilon)|\bm{x}\}=\mathbb{E}\{\ell_{\tau}^{\prime}(\varepsilon)-\varepsilon|\bm{x}\}=\tau\mathbb{E}\{\ell^{\prime}(\varepsilon/\tau)-\varepsilon/\tau|\bm{x}\}$. Together with property (iii), this implies

[TABLE]

We thus use $b^{*}_{\tau}=\|\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})\|_{2}$ to quantify the bias; see Lemma 6.1 and Theorem 2.2.

Below we list five examples of $\ell(\cdot)$ (including the Huber loss) that satisfy Condition 4.1.

(Huber loss): $\ell(x)=x^{2}/2\cdot I(|x|\leq 1)+(|x|-1/2)\cdot I(|x|>1)$ with $\ell^{\prime}(x)=xI(|x|\leq 1)+\mathop{\mathrm{sign}}(x)I(|x|>1)$ and $\ell^{\prime\prime}(x)=I(|x|\leq 1)$. Moreover,

[TABLE] 2. 2.

(Pseudo-Huber loss I): $\ell(x)=\sqrt{1+x^{2}}-1$, whose first and second derivatives are

[TABLE]

respectively. It is easy to see that $\sup_{x\in\mathbb{R}}|\ell^{\prime}(x)|\leq 1$ and $\ell^{\prime\prime}(x)\geq(1+c^{2})^{-3/2}$ for all $|x|\leq c$ and $c>0$. Moreover, since $\ell^{\prime\prime\prime}(x)=-3x(1+x^{2})^{-5/2}$ satisfies $|\ell^{\prime\prime\prime}(x)|<0.9$ for all $x$, it follows from Taylor’s theorem and Lagrange error bound that $|\ell^{\prime}(x)-x|=|\ell^{\prime}(x)-\ell^{\prime}(0)-\ell^{\prime\prime}(0)x|\leq 0.45x^{2}$. 3. 3.

(Pseudo-Huber loss II): $\ell(x)=\log\{(e^{x}+e^{-x})/2\}$, whose first and second derivatives are, respectively,

[TABLE]

It follows that $\sup_{x\in\mathbb{R}}|\ell^{\prime}(x)|\leq 1$ and $\ell^{\prime\prime}(x)\geq 4(e^{c}+e^{-c})^{-2}$ for all $|x|\leq c$ and $c>0$. Moreover, we calculate the third derivative $\ell^{\prime\prime\prime}(x)=-8(e^{x}-e^{-x})(e^{x}+e^{-x})^{-4}$ that satisfies $|\ell^{\prime\prime\prime}(x)|<0.4$. Again, by Taylor’s theorem and Lagrange error bound, $|\ell^{\prime}(x)-x|\leq 0.2x^{2}$. 4. 4.

(Smoothed Huber loss I): The Huber loss is twice differentiable in $\mathbb{R}$, except at $\pm 1$. Modifying the Huber loss gives rise to the following function that is twice differentiable everywhere:

[TABLE]

whose first and second derivatives are

[TABLE]

Direct calculations show that $\sup_{x\in\mathbb{R}}|\ell^{\prime}(x)|\leq 1/2$ and $\ell^{\prime\prime}(x)\geq 1-c$ for all $|x|\leq c$ and $0<c<1$. Since $\ell^{\prime\prime}$ is 1-Lipschitz continuous, we have $|\ell^{\prime}(x)-x|\leq x^{2}/2$. 5. 5.

(Smoothed Huber loss II): Another smoothed version of the Huber loss function is

[TABLE]

The derivative of this function is used in Catoni and Giulini (2017) for mean vector estimation. We compute

[TABLE]

It is easy to see that $\sup_{x\in\mathbb{R}}|\ell^{\prime}(x)|\leq 2\sqrt{2}/3$ and $\ell^{\prime\prime}(x)\geq 1-c^{2}/2$ for all $|x|\leq c$ and $0<c<\sqrt{2}$. Noting that $\ell^{\prime\prime}$ is $\sqrt{2}$-Lipschitz continuous, it holds $|\ell^{\prime}(x)-x|\leq x^{2}/\sqrt{2}$.

The loss functions discussed above, along with their derivatives up to order three, are plotted in Figure 1 except for the Huber loss. Provided that the loss function $\ell_{\tau}(\cdot)$ satisfies Condition 4.1, all the theoretical results in Sections 2.3 and 2.4 remain valid only with different constants. It is worth noticing that the four loss functions described in examples 2–5 also have Lipschitz continuous second derivatives; see Figure 1. In fact, if the function $\ell$ satisfies $\ell^{\prime}(0)=0$, $\ell^{\prime\prime}(0)=1$ and has $L_{2}$-Lipschitz second derivative, then property (iii) in Condition 4.1 holds with $c_{2}=L_{2}/2$. The Lipschitz continuity of $\ell^{\prime\prime}(\cdot)$ also helps remove the anti-concentration condition (34) on $\varepsilon$.

5 Numerical Study

In this section, we compare the empirical performance of the proposed multi-step penalized robust regression estimator with several benchmark methods, such as the Lasso (Tibshirani, 1996), the SCAD and MC+ penalized least squares (Fan and Li, 2001; Zhang, 2010a). All the computational results presented below are reproducible using software available at https://github.com/XiaoouPan/ILAMM.

We generate data vectors $\{(y_{i},\bm{x}_{i})\}_{i=1}^{n}$ from two types of linear models:

(Homoscedastic model): $y_{i}=\bm{x}_{i}^{{\rm T}}\bm{\beta}^{*}+\varepsilon_{i}$ with $\bm{x}_{i}\sim\mathsf{N}(0,{\rm I}_{d})$, $i=1,\ldots,n$. 2. 2.

(Heteroscedastic model): $y_{i}=\bm{x}_{i}^{{\rm T}}\bm{\beta}^{*}+c^{-1}(\bm{x}_{i}^{{\rm T}}\bm{\beta}^{*})^{2}\varepsilon_{i}$ with $\bm{x}_{i}\sim\mathsf{N}(0,{\rm I}_{d})$ for $i=1,\ldots,n$, where the constant $c$ is chosen as $c=\sqrt{3}\,\|\bm{\beta}^{*}\|_{2}^{2}$ such that $\mathbb{E}\{c^{-1}(\bm{x}_{i}^{{\rm T}}\bm{\beta}^{*})^{2}\}^{2}=1$, and therefore the variance of the noise is the same as that of $\varepsilon_{i}$.

In addition, we consider the following four error distributions:

Normal distribution $\mathsf{N}(\mu,\sigma^{2})$ with mean $\mu=0$ and standard deviation $\sigma=1.5$; 2. 2.

Skewed generalized $t$ distribution ${\sf sgt}(0,5,0.75,2,2.5)$ (Theodossiou, 1998) with mean $\mu=0$, variance $\sigma^{2}=q/(q-2)=5$, $q=2.5$, skewness parameter $\lambda=0.75$ and shape parameter $p=2$; 3. 3.

Lognormal distribution $\mathsf{LN}(\mu,\sigma^{2})$ with log location parameter $\mu=0$ and log shape parameter $\sigma=1.2$; 4. 4.

Pareto distribution ${\sf Par}(x_{m},\alpha)$ with scale parameter $x_{m}=2$ and shape parameter $\alpha=2.2$.

Except for the normal distribution, all the other three are skewed and heavy-tailed. To meet our model assumption, we subtract the mean from the lognormal and Pareto distributions.

In both homoscedastic and heteroscedastic models, the sample size $n=100$, the ambient dimension $d=1000$ and the sparsity parameter $s=6$. The true vector of regression coefficients is $\bm{\beta}^{*}=(4,3,2,-2,-2,2,0,\dots,0)^{{\rm T}}$, where the first $6$ elements are non-zero and the rest are all equal to [math]. We apply the proposed TAC (Tightening After Contraction) algorithm to compute all the estimators with tuning parameters $\lambda$ and $\tau$ chosen via three-fold cross-validation. To be more specific, we first choose a sequence of $\lambda$ values the same way as in the glmnet algorithm (Friedman, Hastie and Tibshirani, 2010). Guided by its theoretically “optimal” magnitude, the candidate set for $\tau$ is taken to be $\{2^{j}\widehat{\sigma}^{\mbox{{\tiny MAD}}}\sqrt{n/\mathrm{log}(nd)}:j=-2,-1,0,1,2\}$, where $\widehat{\sigma}^{\mbox{{\tiny MAD}}}:=\mathrm{median}\{|\widehat{\bm{R}}-\mathrm{median}(\widehat{\bm{R}})|\}/\Phi^{-1}(3/4)$ is the median absolute deviation (MAD) estimator using the residuals $\widehat{\bm{R}}=(\widehat{r}_{1},\ldots,\widehat{r}_{n})^{{\rm T}}$ obtained from the Lasso.

To highlight the tail robustness and oracle property of our algorithm, we consider the following four measurements to assess the empirical performance:

True positive, TP, which is the number of signal variables that are selected; 2. 2.

False positive, FP, which is the number of noise variables that are selected; 3. 3.

Relative error, RE1 and RE2, which is the relative error of an estimator $\widehat{}\bm{\beta}$ with respect to the Lasso under $\ell_{1}$- and $\ell_{2}$-norms:

[TABLE]

Tables 1 and 2 summarize the averages of each measurement, TP, FP, RE1, and RE2 with their standard deviations in brackets, over 200 replications under both homoscedastic and heteroscadastic models. RE1 and RE2 for Lasso are defined to be one, so we omit their standard deviations. Here, Huber-SCAD and Huber-MC+ signify the proposed two-stage algorithm using the SCAD and MC+ penalties, respectively. When the noise distributions are heavy-tailed and/or skewed, we see that Huber-SCAD and Huber-MC+ outperform SCAD and MC+, respectively, with fewer spurious discoveries (false positives), smaller estimation errors and less variability. Under the homoscedastic normal model, Huber-SCAD and Huber-MC+ perform similarly to their least squares counterparts; while under heteroscedasticity, the proposed algorithm exhibits a notable advantage over existing methods on selection consistency even though the error is normally distributed. In summary, these numerical studies validate our expectations that the proposed robust regression algorithm improves the Lasso as a general regression analysis method on two aspects: robustness against heavy-tailed (and even heteroscedastic) noise and selection consistency.

To further visualize the advantage of the multi-step penalized robust regression methods over the existing ones (e.g. Lasso, SCAD and MC+), we draw the receiver operating characteristic (ROC) curve, which is the plot of true positive rate (TPR) against false positive rate (FPR) at various regularization parameters. Specifically, TPR and FPR are defined, respectively, as the ratio of true positive to $s$ and the ratio of false positive to $d-s$. We generate data vectors $\{(y_{i},\bm{x}_{i})\}_{i=1}^{n}$ from both homoscedastic and heteroscedastic models with sample size $n=100$, dimension $d=1000$ and sparsity $s=10$. The true vector of regression coefficients is $\bm{\beta}^{*}=(1.5,1.5,\dots,1.5,0,\dots,0)^{{\rm T}}$, where the first $10$ elements are non-zero with weaker signals than the previous experiment, and the rest are all equal to [math]. We apply the proposed TAC algorithm to implement all the five methods, Lasso, SCAD, MC+, Huber-SCAD and Huber-MC+, with a sequence of $\lambda$ values chosen as before and $\tau$ as $\widehat{\sigma}^{\mbox{{\tiny MAD}}}\sqrt{n/\mathrm{log}(nd)}$. For each combination of $\lambda$ and $\tau$, the empirical FPR and TPR are computed based on 200 simulations.

Figures 2 and 3 indicate evident advantage of Huber-SCAD and Huber-MC+ over their least squares counterparts: the robust methods have a greater area under the curve (AUC) when the noise distribution is heavy-tailed and/or skewed in both homoscedastic and heteroscedastic models. Surprisingly, even in a normal model, the proposed methods still outperform the competitors by a visible margin.

{supplement}

\stitle

6 Preliminaries

Assume we observe independent data $\{(y_{i},\bm{x}_{i})\}_{i=1}^{n}$ from the linear model $y_{i}=\bm{x}_{i}^{{\rm T}}\bm{\beta}^{*}+\varepsilon_{i}$. Let $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{d})^{{\rm T}}$ be a $d$-vector of regularization parameters with $\lambda_{j}\geq 0$. Consider the optimization problem

[TABLE]

where $\widehat{\mathcal{L}}_{\tau}(\bm{\beta})=(1/n)\sum_{i=1}^{n}\ell_{\tau}(y_{i}-\bm{x}_{i}^{{\rm T}}\bm{\beta})$ and $\bm{\lambda}\circ\bm{\beta}=(\lambda_{1}\beta_{1},\ldots,\lambda_{d}\beta_{d})^{{\rm T}}$. Moreover, define the population loss $\mathcal{L}_{\tau}(\bm{\beta})=\mathbb{E}\widehat{\mathcal{L}}_{\tau}(\bm{\beta})$.

The following result provides conditions under which an $\epsilon$-optimal solution to the convex program (41) falls in an $\ell_{1}$-cone. Recall that ${\mathcal{S}}={\rm supp}(\bm{\beta}^{*})$ and ${\mathcal{S}}^{{\rm c}}=[d]\setminus{\mathcal{S}}$. Moreover, define

[TABLE]

which are, respectively, the centered score function and the approximation bias. First, we characterize the magnitude of the bias $b_{\tau}^{*}:=b(\bm{\beta}^{*})$, as a function of $\tau$.

Lemma 6.1**.**

Assume $\mu_{1}=\sup_{\bm{u}\in\mathbb{S}^{d-1}}\mathbb{E}|\bm{u}^{{\rm T}}\bm{x}|<\infty$, $\mathbb{E}(\varepsilon|\bm{x})=0$ and $\mathbb{E}(\varepsilon^{2}|\bm{x})\leq\sigma_{2}^{2}$ almost surely. Then $|b_{\tau}^{*}|\leq\mu_{1}\sigma^{2}_{2}\tau^{-1}$ for any $\tau>0$, and $\tau|b^{*}_{\tau}|\to 0$ as $\tau\to\infty$.

Proof of Lemma 6.1.

Note that $\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})=\mathbb{E}\{\ell_{\tau}^{\prime}(\varepsilon)\bm{x}\}$, where $\ell_{\tau}^{\prime}(u)=uI(|u|\leq\tau)+\tau\mathop{\mathrm{sign}}(u)I(|u|>\tau)$. Recalling $\mathbb{E}(\varepsilon|\bm{x})=0$, it follows that

[TABLE]

By the variational representation of the $\ell_{2}$-norm, we have

[TABLE]

as claimed. The second claim follows from the fact that $\mathbb{E}\{\varepsilon^{2}I(|\varepsilon|>\tau)|\bm{x}\}\to 0$ as $\tau\to\infty$. ∎

Lemma 6.2**.**

Let $\mathcal{E}$ be a subset of $[d]$ that contains ${\mathcal{S}}$. For any $\bm{\beta}\in\mathbb{R}^{d}$ satisfying $\bm{\beta}_{\mathcal{E}^{{\rm c}}}=\textbf{0}$ and $\epsilon>0$, provided $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{d})^{{\rm T}}$ satisfies $\|\bm{\lambda}_{\mathcal{E}^{{\rm c}}}\|_{\min}>\|\bm{w}(\bm{\beta})\|_{\infty}+\epsilon$, any $\epsilon$-optimal solution $\widetilde{}\bm{\beta}$ to (41) satisfies

[TABLE]

Proof of Lemma 6.2.

For any $\bm{\xi}\in\partial\|\widetilde{}\bm{\beta}\|_{1}$, define $\bm{u}_{\bm{\xi}}=\nabla\widehat{\mathcal{L}}_{\tau}(\widetilde{}\bm{\beta})+\bm{\lambda}\circ\bm{\xi}\in\mathbb{R}^{d}$. Note that

[TABLE]

Moreover, we have

[TABLE]

Together, the last two displays imply

[TABLE]

Since the right-hand side of this inequality does not depend on $\bm{\xi}$, taking the infimum with respect to $\bm{\xi}\in\partial\|\widetilde{}\bm{\beta}\|_{1}$ on both sides to reach

[TABLE]

By definition, $\widetilde{}\bm{\beta}$ is an $\epsilon$-optimal solution so that $\omega_{\bm{\lambda}}(\widetilde{}\bm{\beta})\leq\epsilon$. Putting together the pieces, we obtain

[TABLE]

Decompose $\|\widetilde{}\bm{\beta}-\bm{\beta}\|_{1}$ as $\|(\widetilde{}\bm{\beta}-\bm{\beta})_{\mathcal{E}}\|_{1}+\|(\widetilde{}\bm{\beta}-\bm{\beta})_{\mathcal{E}^{{\rm c}}}\|_{1}$, the stated result follows immediately. ∎

Lemma 6.3**.**

Consider some $\bm{\beta}\in\mathbb{R}^{d}$ satisfying $\bm{\beta}_{{\mathcal{S}}^{{\rm c}}}=\textbf{0}$, and let $\mathcal{E}\subseteq[d]$ be a subset that contains ${\mathcal{S}}$ and has cardinality $|\mathcal{E}|=k$. Assume that $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{d})^{{\rm T}}$ satisfies $\|\bm{\lambda}\|_{\infty}\leq\lambda$ and $\|\bm{\lambda}_{\mathcal{E}^{{\rm c}}}\|_{\min}\geq\rho\lambda>0$ for some $\rho\in(0,1]$ and $\lambda\geq s^{-1/2}b(\bm{\beta})$. Conditioned on event $\{\|\bm{w}(\bm{\beta})\|_{\infty}+\epsilon\leq 0.5\rho\lambda\}$, any $\epsilon$-optimal solution $\widetilde{}\bm{\beta}$ to (41) satisfies $\widetilde{}\bm{\beta}\in\bm{\beta}+\mathbb{C}(l)$, where $l=(2+\frac{2}{\rho})k^{1/2}+\frac{2}{\rho}s^{1/2}$. Moreover, let $r,\kappa>0$ satisfy

[TABLE]

Then, conditioned on the event $\mathcal{E}_{1}(r,l,\kappa)\cap\{\|\bm{w}(\bm{\beta})\|_{\infty}+\epsilon\leq 0.5\rho\lambda\}$,

[TABLE]

Proof of Lemma 6.3.

For some $r>0$ to be specified, define $\eta=\sup\{u\in[0,1]:(1-u)\bm{\beta}+u\widetilde{}\bm{\beta}\in\mathbb{B}(r)\}$, where $\mathbb{B}(r)=\{\bm{\delta}\in\mathbb{R}^{d}:\|\bm{\delta}\|_{2}\leq r\}$. Note that $\eta=1$ if $\widetilde{}\bm{\beta}\in\bm{\beta}+\mathbb{B}(r)$ and $\eta\in(0,1)$ otherwise. Then, the intermediate estimate $\widetilde{}\bm{\beta}_{\eta}:=\eta\widetilde{}\bm{\beta}+(1-\eta)\bm{\beta}$ satisfies (i) $\widetilde{}\bm{\beta}_{\eta}\in\bm{\beta}+\mathbb{B}(r)$, (ii) $\widetilde{}\bm{\beta}_{\eta}$ lies on the boundary of $\bm{\beta}+\mathbb{B}(r)$ with $0<\eta<1$ if $\widetilde{}\bm{\beta}\notin\bm{\beta}+\mathbb{B}(r)$, and (iii) $\widetilde{}\bm{\beta}_{\eta}=\widetilde{}\bm{\beta}$ with $\eta=1$ if $\widetilde{}\bm{\beta}\in\bm{\beta}+\mathbb{B}(r)$.

By the convexity of Huber loss and Lemma F.2 in Fan et al. (2018), we have

[TABLE]

First we bound the left-hand side of (44) from below. Conditioned on the stated event, Lemma 6.2 indicates

[TABLE]

from which it follows that $\|\widetilde{}\bm{\beta}-\bm{\beta}\|_{1}\leq(2+\frac{2}{\rho})k^{1/2}\|\widetilde{}\bm{\beta}-\bm{\beta}\|_{2}+\frac{2}{\rho}\lambda^{-1}b(\bm{\beta})\|\widetilde{}\bm{\beta}-\bm{\beta}\|_{2}$. Provided that $\lambda\geq s^{-1/2}b(\bm{\beta})$, this implies $\widetilde{}\bm{\beta}\in\bm{\beta}+\mathbb{C}(l)$ with $l=(2+\frac{2}{\rho})k^{1/2}+\frac{2}{\rho}s^{1/2}$. Since $\widetilde{}\bm{\beta}_{\eta}-\bm{\beta}=\eta(\widetilde{}\bm{\beta}-\bm{\beta})$, we have $\widetilde{}\bm{\beta}_{\eta}\in\bm{\beta}+\mathbb{B}(r)\cap\mathbb{C}(l)$ and conditioned on event $\mathcal{E}_{1}(r,l,\kappa)$,

[TABLE]

Next we upper bound the right-hand side of (44). For any $\bm{\xi}\in\partial\|\widetilde{}\bm{\beta}\|_{1}$, write

[TABLE]

where $\bm{u}=\nabla\widehat{\mathcal{L}}_{\tau}(\widetilde{}\bm{\beta})+\bm{\lambda}\circ\bm{\xi}$. For $\Pi_{3}=\langle\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta}),\widetilde{}\bm{\beta}-\bm{\beta}\rangle$, by the decomposition $\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta})=\bm{w}(\bm{\beta})+\nabla\mathcal{L}_{\tau}(\bm{\beta})$ we have

[TABLE]

Turning to $\Pi_{2}$, decompose $\bm{\lambda}\circ\bm{\xi}$ and $\widetilde{}\bm{\beta}-\bm{\beta}$ according to ${\mathcal{S}}\cup(\mathcal{E}\setminus{\mathcal{S}})\cup\mathcal{E}^{{\rm c}}$ to reach

[TABLE]

Since $\bm{\beta}_{\mathcal{E}^{{\rm c}}}=\textbf{0}$ and $\bm{\xi}\in\partial\|\widetilde{}\bm{\beta}\|_{1}$, we have $\langle(\bm{\lambda}\circ\bm{\xi})_{\mathcal{E}^{{\rm c}}},(\widetilde{}\bm{\beta}-\bm{\beta})_{\mathcal{E}^{{\rm c}}}\rangle=\langle\bm{\lambda}_{\mathcal{E}^{{\rm c}}},|\widetilde{}\bm{\beta}_{\mathcal{E}^{{\rm c}}}|\rangle=\langle\bm{\lambda}_{\mathcal{E}^{{\rm c}}},|(\widetilde{}\bm{\beta}-\bm{\beta})_{\mathcal{E}^{{\rm c}}}|\rangle$. Also, $\langle(\bm{\lambda}\circ\bm{\xi})_{\mathcal{E}\setminus{\mathcal{S}}},(\widetilde{}\bm{\beta}-\bm{\beta})_{\mathcal{E}\setminus{\mathcal{S}}}\rangle=\langle(\bm{\lambda}\circ\bm{\xi})_{\mathcal{E}\setminus{\mathcal{S}}},\widetilde{}\bm{\beta}_{\mathcal{E}\setminus{\mathcal{S}}}\rangle\geq 0$. Therefore,

[TABLE]

Similarly, $\Pi_{1}$ satisfies the bound

[TABLE]

Together, (46)–(49) yield

[TABLE]

Taking the infimum over $\bm{\xi}\in\partial\|\widetilde{}\bm{\beta}\|_{1}$ on both sides, it follows that

[TABLE]

It follows from (44), (45) and (50) that conditioned on $\mathcal{E}_{1}(r,l,\kappa)\cap\{\|\bm{w}(\bm{\beta})\|_{\infty}+\epsilon\leq 0.5\rho\lambda\}$,

[TABLE]

On the same event, note that

[TABLE]

Moreover, recall that $\|\bm{\lambda}_{{\mathcal{S}}}\|_{2}\leq s^{1/2}\lambda$ and $b(\bm{\beta})\leq s^{1/2}\lambda$. Plugging these bounds into (51) yields

[TABLE]

Hence, $\widetilde{}\bm{\beta}_{\eta}$ falls in the interior of $\bm{\beta}+\mathbb{B}(r)$. Via proof by contradiction, we must have $\eta=1$ and $\widetilde{}\bm{\beta}_{\eta}=\widetilde{}\bm{\beta}$. Consequently, (42) and (43) follow, respectively, from (51) and the last display. ∎

In Lemma 6.3, we need $\lambda$ to be sufficiently large in the sense that for some $s$-sparse vector $\bm{\beta}$,

[TABLE]

where $\bm{w}(\bm{\beta})=\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta})-\nabla\mathcal{L}_{\tau}(\bm{\beta})$ is the centered gradient at $\bm{\beta}$ and $b(\bm{\beta})=\|\nabla\mathcal{L}_{\tau}(\bm{\beta})\|_{2}$ quantifies the bias. To prove the weak oracle property (Proposition 2.2), we will take $\bm{\beta}=\bm{\beta}^{*}$ and control the stochastic term $\|\bm{w}^{*}\|_{\infty}=\|\bm{w}(\bm{\beta}^{*})\|_{\infty}$ and bias term $b(\bm{\beta}^{*})$ separately. For some $\bm{\beta}$ which has vanishing or negligible bias, we will only focus on the stochastic term $\|\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta})\|_{\infty}$, as described in the next lemma. Recall the event $\mathcal{E}_{2}(r,l,\kappa)$ defined in (21), and $\mathbb{B}_{\Sigma}(r)=\{\bm{\beta}\in\mathbb{R}^{d}:\|\bm{\beta}\|_{\Sigma}\leq r\}$.

Lemma 6.4**.**

Consider some $\bm{\beta}\in\bm{\beta}^{*}+\mathbb{B}_{\Sigma}(r)$ ($r>0$) satisfying $\bm{\beta}_{{\mathcal{S}}^{{\rm c}}}=\textbf{0}$, and let $\mathcal{E}\subseteq[d]$ be a subset that contains ${\mathcal{S}}$ and has cardinality $|\mathcal{E}|=k$. Assume that $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{d})^{{\rm T}}$ satisfies $\|\bm{\lambda}\|_{\infty}\leq\lambda$ and $\|\bm{\lambda}_{\mathcal{E}^{{\rm c}}}\|_{\min}\geq\rho\lambda>0$ for some $\rho\in(0,1]$. Conditioned on event $\{\|\widehat{\mathcal{L}}_{\tau}(\bm{\beta})\|_{\infty}\leq 0.5\rho\lambda\}$, any optimal solution $\widehat{}\bm{\beta}$ to (41) satisfies $\widehat{}\bm{\beta}\in\bm{\beta}+\mathbb{C}(l)$, where $l=(2+\frac{2}{\rho})k^{1/2}+\frac{2}{\rho}s^{1/2}$. Moreover, let $r,\kappa>0$ satisfy $r>\kappa^{-1}(0.5\rho k^{1/2}+2s^{1/2})\lambda$. Then, conditioned on the event $\mathcal{E}_{2}(r,l,\kappa)\cap\{\|\widehat{\mathcal{L}}_{\tau}(\bm{\beta})\|_{\infty}\leq 0.5\rho\lambda\}$,

[TABLE]

The proof of Lemma 6.4 is based on the same arguments from the proof of Lemmas 6.2 and 6.3, and thus is omitted. To prove the strong oracle property (Proposition 2.3), we will apply Lemma 6.4 with $\bm{\beta}=\widehat{}\bm{\beta}^{{{\rm ora}}}$, the oracle estimator defined in (20).

7 Proofs of Propositions

7.1 Proof of Proposition 2.1

With the initial estimate $\widetilde{}\bm{\beta}^{(0)}=\textbf{0}_{d}$, we have $\bm{\lambda}^{(0)}=p_{\lambda}^{\prime}(\textbf{0}_{d})=(\lambda,\ldots,\lambda)^{{\rm T}}\in\mathbb{R}^{d}$. Then (14) follows immediately from Lemma 6.3 with $\bm{\beta}=\bm{\beta}^{*}$, $\mathcal{E}={\mathcal{S}}$ and $\rho=1$. ∎

7.2 Proof of Proposition 2.2

In order to improve the statistical rate at step $\ell\geq 1$, we need to control the magnitude of the spurious discoveries from the last step, that is, $\max_{j\in{\mathcal{S}}^{{\rm c}}}|\widetilde{\beta}_{j}^{(\ell-1)}|$. Recall that $\bm{\lambda}^{(\ell-1)}=(\lambda_{1}^{(\ell-1)},\ldots,\lambda_{d}^{(\ell-1)})^{{\rm T}}=(p_{\lambda}^{\prime}(|\widetilde{\beta}_{1}^{(\ell-1)}|),\ldots,p_{\lambda}^{\prime}(|\widetilde{\beta}_{d}^{(\ell-1)}|))^{{\rm T}}$ and $p_{\lambda}(t)=\lambda^{2}p(t/\lambda)$ for $t\in\mathbb{R}$. Intuitively, the larger $|\widetilde{\beta}_{j}^{(\ell-1)}|$ is, the smaller $\lambda^{(\ell-1)}_{j}$ is. Motivated by this observation, we construct an augmented set $\mathcal{E}_{\ell}$ of ${\mathcal{S}}$ in each step and control the magnitude of $\|\bm{\lambda}^{(\ell-1)}_{\mathcal{E}_{\ell}^{{{\rm c}}}}\|_{\min}$.

Starting from $\widetilde{}\bm{\beta}^{(0)}=\textbf{0}$, we have $\bm{\lambda}^{(0)}=(\lambda,\ldots,\lambda)^{{\rm T}}\in\mathbb{R}^{d}$. Recall from (15) that $\lambda\geq s^{-1/2}b^{*}_{\tau}$, or equivalently, $\lambda^{-1}b^{*}_{\tau}\leq s^{1/2}$. Then, applying Lemma 6.3 with $\mathcal{E}={\mathcal{S}}$ and $l_{0}=6s^{1/2}$ we obtain that, conditioning on $\mathcal{E}_{1}(r,l_{0},\kappa)\cap\{\lambda\geq\frac{2}{p^{\prime}(\gamma)}(\|\bm{w}^{*}\|_{\infty}+\epsilon_{1})\}$,

[TABLE]

where the last inequality is due to (15). For $\ell\geq 1$, define the augmented set

[TABLE]

which depends on the solution $\widetilde{}\bm{\beta}^{(\ell-1)}$ from the previous step. We claim that the above constructed sets satisfy

[TABLE]

where $c$ is the constant determined by (16). If these were true, it follows from Lemma 6.3 with $\rho=p^{\prime}(\gamma)$, $k=(c^{2}+1)s$ and $l=(2+\frac{2}{\rho})(c^{2}+1)^{1/2}s^{1/2}+\frac{2}{\rho}s^{1/2}$ that, conditioned on $\mathcal{E}_{1}(r,l,\kappa)\cap\{\lambda\geq\frac{2}{p^{\prime}(\gamma)}(\|\bm{w}^{*}\|_{\infty}+\epsilon_{1})\}$,

[TABLE]

where the last two steps follow from (16) and (17).

We prove the earlier claim (54) by the method of induction. For $\ell=1$, we have $\bm{\lambda}^{(0)}=(\lambda,\ldots,\lambda)^{{\rm T}}\in\mathbb{R}^{d}$. Thus, (54) holds with $\mathcal{E}_{1}={\mathcal{S}}$. Next, assume (54) holds for some $\ell\geq 1$, from which (56) follows. To bound the cardinality of $\mathcal{E}_{\ell+1}$, note that for any $j\in\mathcal{E}_{\ell+1}\setminus{\mathcal{S}}$, $p^{\prime}_{\lambda}(|\widetilde{\beta}_{j}^{(\ell)}|)=\lambda_{j}^{(\ell)}<p^{\prime}(\gamma)\lambda=p_{\lambda}^{\prime}(\gamma\lambda)$. This, together with the monotonicity of $p_{\lambda}^{\prime}$ on $\mathbb{R}_{+}$, implies $|\widetilde{\beta}_{j}^{(\ell)}|>\gamma\lambda$. Recalling that $\beta^{*}_{j}=0$ for $j\in\mathcal{E}_{\ell+1}\setminus{\mathcal{S}}$, we obtain

[TABLE]

where inequality (i) applies the bound (56). Hence, $|\mathcal{E}_{\ell+1}|\leq|{\mathcal{S}}|+|\mathcal{E}_{\ell+1}\setminus{\mathcal{S}}|<(c^{2}+1)s$. By (53) and the property $p_{\lambda}^{\prime}(t)=\lambda p^{\prime}(t/\lambda)$, we are guaranteed that

[TABLE]

The two hypotheses in (54) then hold for $\ell+1$, which completes the induction step. Consequently, the bounds (55) and (56) hold for any $\ell\geq 1$.

We have shown that under proper conditions, all the estimates $\widetilde{}\bm{\beta}^{(\ell)}$ fall in a local neighborhood of $\bm{\beta}^{*}$, i.e.,$\|\widetilde{}\bm{\beta}^{(\ell)}-\bm{\beta}^{*}\|_{\Sigma}\leq r^{{{\rm crude}}}=c\gamma s^{1/2}\lambda$. To further refine this bound as signal strengthens, on the right-hand side of (55), we need to establish sharper bounds on

[TABLE]

and maintain the bias term $b^{*}_{\tau}$, instead of replacing it with an upper bound $s^{1/2}\lambda$. For each $j\in[d]$, $\lambda_{j}^{(\ell-1)}=p^{\prime}_{\lambda}(|\widetilde{\beta}_{j}^{(\ell-1)}|)$. If $|\widetilde{\beta}_{j}^{(\ell-1)}-\beta_{j}^{*}|\geq\gamma\lambda$, then $\lambda_{j}^{(\ell-1)}\leq\lambda\leq\gamma^{-1}|\widetilde{\beta}_{j}^{(\ell-1)}-\beta_{j}^{*}|$; otherwise if $|\widetilde{\beta}_{j}^{(\ell-1)}-\beta_{j}^{*}|\leq\gamma\lambda$, $\lambda_{j}^{(\ell-1)}\leq p_{\lambda}^{\prime}(|\beta_{j}^{*}|-\gamma\lambda)$ due to monotonicity of $p_{\lambda}^{\prime}$. Putting together the pieces, we conclude that

[TABLE]

For the remaining terms that involve $\mathcal{E}_{\ell}$, by the triangle inequality and (57) we obtain that

[TABLE]

Plugging the above refined bounds into (55) yields

[TABLE]

Taking $\delta=p^{\prime}(\gamma)/(2\gamma\kappa)\in(0,1)$, the contraction inequality (2.2) follows immediately. Finally, (19) is a direct consequence of (2.2) and (52). ∎

7.3 Proof of Proposition 2.3

By construction, the oracle estimator $\widehat{}\bm{\beta}^{{\rm ora}}$ is such that $\widehat{}\bm{\beta}^{\rm ora}_{{\mathcal{S}}^{\rm c}}=\textbf{0}\in\mathbb{R}^{d-s}$ and $\nabla\widehat{\mathcal{L}}_{\tau}(\widehat{}\bm{\beta}^{{\rm ora}})_{{\mathcal{S}}}=\textbf{0}\in\mathbb{R}^{s}$. With $\bm{w}^{{\rm ora}}=\nabla\widehat{\mathcal{L}}_{\tau}(\widehat{}\bm{\beta}^{\rm ora})$, the proof strategy is similar to that in the proof of Proposition 2.2 with $\epsilon_{\ell}=0$, because $\widehat{}\bm{\beta}^{(\ell)}$ are optimal solutions to $({\rm P}_{\ell})$.

Recall that $l=\{2+\frac{2}{p^{\prime}(\gamma_{0})}\}(c_{0}^{2}+1)^{1/2}s^{1/2}$ with $c_{0}>0$ determined by (23). Conditioned on the event

[TABLE]

following the proof of Proposition 2.2 and applying Lemma 6.4 with $\bm{\beta}=\widehat{}\bm{\beta}^{{\rm ora}}$ and $\mathcal{E}=\mathcal{E}_{\ell}$, it can be similarly shown that

[TABLE]

where similarly to (53) and (54), $\mathcal{E}_{\ell}={\mathcal{S}}\cup\{1\leq j\leq d:\lambda_{j}^{(\ell-1)}<p^{\prime}(\gamma_{0})\lambda\}$ is such that $|\mathcal{E}_{\ell}\setminus{\mathcal{S}}|<c_{0}^{2}s$ and thus $|\mathcal{E}_{\ell}|<(c_{0}^{2}+1)s$. In this case, the approximation bias is hidden in $\|\bm{w}^{{\rm ora}}\|_{\infty}$. Moreover, define a sequence of subsets

[TABLE]

Starting with $\widehat{}\bm{\beta}^{(0)}=\textbf{0}$, it holds under the minimum signal strength condition that ${\mathcal{S}}_{0}={\mathcal{S}}$.

To obtain a refined upper bound on $\|\bm{\lambda}^{(\ell-1)}_{{\mathcal{S}}}\|_{2}$, note that if $j\in{\mathcal{S}}\cap{\mathcal{S}}_{\ell-1}^{{\rm c}}$, $\lambda_{j}^{(\ell-1)}=p_{\lambda}^{\prime}(|\widehat{\beta}^{(\ell-1)}_{j}|)\leq p_{\lambda}^{\prime}(|\beta_{j}^{*}|-\gamma_{0}\lambda)$ due to monotonicity; otherwise if $j\in{\mathcal{S}}\cap{\mathcal{S}}_{\ell-1}$, $\lambda_{j}^{(\ell-1)}\leq\lambda$. Therefore,

[TABLE]

Since $\|\bm{\beta}_{{\mathcal{S}}}^{*}\|_{\min}\geq(\gamma_{0}+\gamma_{1})\lambda$ and $p_{\lambda}^{\prime}(t)=0$ for all $t\geq\gamma_{1}\lambda$, $\|p_{\lambda}^{\prime}(|\bm{\beta}^{*}_{{\mathcal{S}}}|-\gamma_{0}\lambda)\|_{2}$ vanishes. Turning to $\|\bm{w}^{\rm ora}_{\mathcal{E}_{\ell}}\|_{2}$, by the first-order condition of minimizing $\bm{\beta}_{{\mathcal{S}}}\mapsto(1/n)\sum_{i=1}^{n}\widehat{\mathcal{L}}_{\tau}(y_{i}-\bm{x}_{i,{\mathcal{S}}}^{{\rm T}}\bm{\beta}_{{\mathcal{S}}})$, we have $\bm{w}^{\rm ora}_{{\mathcal{S}}}=\textbf{0}$ and hence

[TABLE]

For each $j\in\mathcal{E}_{\ell}\setminus{\mathcal{S}}$, $\beta^{*}_{j}=0$ and $\lambda_{j}^{(\ell-1)}=p_{\lambda}^{\prime}(|\widehat{\beta}_{j}^{(\ell-1)}|)<p^{\prime}(\gamma_{0})\lambda=p_{\lambda}^{\prime}(\gamma_{0}\lambda)$. Hence, $|\widehat{\beta}_{j}^{(\ell-1)}-\beta^{*}_{j}|=|\widehat{\beta}_{j}^{(\ell-1)}|>\gamma_{0}\lambda$ so that $j\in{\mathcal{S}}_{\ell-1}\setminus{\mathcal{S}}$. Therefore, $\mathcal{E}_{\ell}\setminus{\mathcal{S}}\subseteq{\mathcal{S}}_{\ell-1}\setminus{\mathcal{S}}$. Combined with the earlier bound, we arrive at

[TABLE]

Since $p^{\prime}(\gamma_{0})\leq 1/2$, substituting the above estimates into (58) yields

[TABLE]

Next we bound $|{\mathcal{S}}_{\ell}|$ $(\ell\geq 1)$, the cardinality of ${\mathcal{S}}_{\ell}$. By (24), it holds for any $j\in{\mathcal{S}}_{\ell}$ that

[TABLE]

In conjunction with (59), this implies

[TABLE]

Recall that ${\mathcal{S}}_{0}={\mathcal{S}}$, we have $|{\mathcal{S}}_{\ell}|^{1/2}<\delta^{\ell}s^{1/2}$ for any $\ell\geq 1$. As long as $\ell\geq T:=\lceil\log(s^{1/2})/\log(1/\delta)\rceil$, we are guaranteed that $|{\mathcal{S}}_{\ell}|<1$, i.e. ${\mathcal{S}}_{\ell}=\emptyset$. Consequently, it follows from (59) that $\widehat{}\bm{\beta}^{(\ell)}=\widehat{}\bm{\beta}^{{\rm ora}}$ for all $\ell\geq T+1$. This completes the proof. ∎

7.4 Proof of Proposition 2.4

The proof is based on similar arguments that were used in the proof of Lemmas C.3 and C.4 in Sun, Zhou and Fan (2020). We only present the necessary steps in order to slightly relax the sub-Gaussian condition on $\bm{x}_{i}=(x_{i1},\ldots,x_{id})^{{\rm T}}$.

For any $\bm{\beta}\in\bm{\beta}^{*}+\mathbb{B}(r)\cap\mathbb{C}(l)$, write $\bm{\delta}=\bm{\beta}-\bm{\beta}^{*}$. Following the proof of Lemma C.3 in Sun, Zhou and Fan (2020), it can be shown under Condition 2.1 that

[TABLE]

where the second inequality uses the bound

[TABLE]

and $\Delta(r,l)=\sup_{\bm{\delta}\in\mathbb{B}(r)\cap\mathbb{C}(l)}(1/n)\sum_{i=1}^{n}|f_{\bm{\delta}}(\bm{x}_{i},\varepsilon_{i})-\mathbb{E}f_{\bm{\delta}}(\bm{x}_{i},\varepsilon_{i})|$ with

[TABLE]

and $\varphi_{c}(u):=u^{2}I(|u|\leq c/2)+(|u|-c)^{2}I(c/2<|u|\leq c)$ for $u\in\mathbb{R}$ and $c\geq 0$. We thus let $\tau\geq\max\{4\sqrt{2}\sigma_{2},8\log(8\rho_{\bm{x}}^{2})\sigma_{\bm{x}}r\}$ so that $(2\sigma_{2}/\tau)^{2}\leq 1/8$ and $\rho_{\bm{x}}^{2}e^{-\tau/(8\sigma_{\bm{x}}r)}\leq 1/8$. It then follows from (61) that

[TABLE]

holds uniformly over $\bm{\beta}\in\bm{\beta}^{*}+\mathbb{B}(r)\cap\mathbb{C}(l)$.

It remains to control the supremum $\Delta(r,l)$. Note that $\varphi_{c}(\cdot)$ is $c$-Lipschitz continuous, and satisfies $\varphi_{bc}(bu)=b^{2}\varphi_{c}(u)$ for $b,c>0$ and $u\in\mathbb{R}$. Thus, the above $f_{\bm{\delta}}(\bm{x}_{i},\varepsilon_{i})$ can be simplified as

[TABLE]

Moreover, since $0\leq\varphi_{c}(u)\leq\min\{(c/2)^{2},u^{2}\}$, we have $0\leq f_{\bm{\delta}}(\bm{x}_{i},\varepsilon_{i})\leq\tau^{2}/(4r)^{2}$ and by (26), $\mathbb{E}f^{2}_{\bm{\delta}}(\bm{x}_{i},\varepsilon_{i})\leq\mathbb{E}(\bm{x}_{i}^{{\rm T}}\bm{\delta}/\|\bm{\delta}\|_{2})^{4}\leq\rho_{\bm{x}}^{4}\{\mathbb{E}(\bm{x}_{i}^{{\rm T}}\bm{\delta}/\|\bm{\delta}\|_{2})^{2}\}^{2}\leq(\rho_{u}\rho_{\bm{x}}^{2})^{2}$ for all $\bm{\delta}\in\mathbb{R}^{d}$. Then, applying Bousquet’s version of Talagrand’s inequality (see, e.g. Theorem 7.3 in Bousquet (2003)), we obtain that for every $t\geq 0$,

[TABLE]

with probability at least $1-e^{-t}$. For $\mathbb{E}\Delta(r,l)$, using Rademacher symmetrization gives

[TABLE]

where $e_{1},\ldots,e_{n}$ are independent Rademacher random variables. Since $\varphi_{c}(\cdot)$ is $c$-Lipshitz, $f_{\bm{\delta}}(\bm{x}_{i},\varepsilon_{i})$ is a $(\tau/2r)$-Lipschitz function in $\bm{x}_{i}^{{\rm T}}\bm{\delta}/\|\bm{\delta}\|_{2}$, i.e., for any sample $(\bm{x}_{i},\varepsilon_{i})$ and parameters $\bm{\delta},\bm{\delta}^{\prime}\in\mathbb{R}^{d}$,

[TABLE]

Moreover, observe that $f_{\bm{\delta}}(\bm{x}_{i},\varepsilon_{i})=0$ for any $\bm{\delta}$ such that $\bm{x}_{i}^{{\rm T}}\bm{\delta}/\|\bm{\delta}\|_{2}=0$, and $I(|\varepsilon_{i}|\leq\tau/2)\in\{0,1\}$. Then, applying Talagrand’s contraction principle (see, e.g. Theorem 4.4, Theorem 4.12 and (4.20) in Ledoux and Talagrand (1991)) yields

[TABLE]

where the last inequality uses the cone constraint that $\|\bm{\delta}\|_{1}\leq l\|\bm{\delta}\|_{2}$. Next, we apply a maximal inequality for sub-exponential random variables to bound the last term on the right-hand side of (64). For $j=1,\ldots,d$, define partial sums $S_{j}=\sum_{i=1}^{n}e_{i}x_{ij}$, of which each summand satisfies $\mathbb{E}(e_{i}x_{ij})=0$ and $\mathbb{E}(e_{i}x_{ij})^{2}=\sigma_{jj}$. More over, for $k=2,4,\ldots$, $\mathbb{E}|e_{i}|^{k}=1$ and

[TABLE]

By the symmetry of Rademacher random variables, and Bernstein’s inequality (see, e.g. Theorems 2.10 in Boucheron, Lugosi and Massart (2013)), we obtain that

[TABLE]

where $\nu=2n\sigma_{\bm{x}}^{2}$ and $c=\sigma_{\bm{x}}$. Following the proof of Theorems 2.5 in Boucheron, Lugosi and Massart (2013), it can be shown that

[TABLE]

Re-arranging terms and using (64), we find that

[TABLE]

Combining this with (62) and (63) yields the bound (2.4). ∎

7.5 Proof of Proposition 2.5

Write $S_{j}=(1/n)\sum_{i=1}^{n}(\xi_{i}x_{ij}-\mathbb{E}\xi_{i}x_{ij})$ with $\xi_{i}:=\ell_{\tau}^{\prime}(\varepsilon_{i})$ for $j=1,\ldots,d$, so that $\|\bm{w}^{*}\|_{\infty}=\max_{1\leq j\leq d}|S_{j}|$. Note that $\ell_{\tau}^{\prime}(u)=uI(|u|\leq\tau)+\tau\mathop{\mathrm{sign}}(u)I(|u|>\tau)$, we have $\mathbb{E}(\xi_{i}^{2}|\bm{x}_{i})\leq\sigma_{2}^{2}$ and $|\xi_{i}|\leq\tau$. It follows that $\mathbb{E}(\xi_{i}x_{ij})^{2}\leq\sigma_{2}^{2}\sigma_{jj}$ and under Condition 2.1,

[TABLE]

for $k=3,4,\ldots.$ Bernstein’s inequality, in conjunction with the union bound, implies that for any $x\geq 0$,

[TABLE]

with probability at least $1-2de^{-x}$. Taking $x=\log(2d)+t$ proves (28). Next we use a standard covering argument to prove (29). For any $\epsilon\in(0,1)$, there exists an $\epsilon$-net $\mathcal{N}_{\epsilon}$ of the unit sphere in $\mathbb{R}^{s}$ with cardinality $|\mathcal{N}_{\epsilon}|\leq(1+2/\epsilon)^{s}$ such that

[TABLE]

For every $\bm{u}\in\mathcal{N}_{\epsilon}$, Bernstein’s condition holds: $\mathbb{E}(\xi_{i}\bm{u}^{{\rm T}}\bm{x}_{i,{\mathcal{S}}})^{2}\leq 2\sigma_{2}^{2}\sigma_{\bm{x}}^{2}$ and for $k=3,4,\ldots$,

[TABLE]

Again, applying (one-sided) Bernstein’s inequality yields, for any $x>0$,

[TABLE]

with probability at least $1-e^{-x}$. Consequently, from the union bound and (65), we have

[TABLE]

with probability at least $1-e^{\log(1+2/\epsilon)s-x}$. Taking $\epsilon=1/3$ and $x=2s+t$ proves (29). ∎

8 Proofs of Theorems

8.1 Proof of Theorem 2.1

We will apply Propositions 2.1, 2.4 and 2.5 to prove Theorem 2.1. To begin with, let $\lambda,r,\kappa>0$ satisfy (13), that is, $\lambda\geq s^{-1/2}b^{*}_{\tau}$ and $r>2.5\kappa^{-1}s^{1/2}\lambda$. Applying Proposition 2.1 with $\epsilon_{1}=0$ yields that, conditioned on the event $\mathcal{E}_{1}(r,6s^{1/2},\kappa)\cap\{\lambda\geq 2\|\bm{w}^{*}\|_{\infty}\}$, the Huber-Lasso estimator defined in (6) satisfies

[TABLE]

It remains to control the above event of interest. Taking $l=6s^{1/2}$ and $\kappa=\rho_{l}/2$, it follows from Proposition 2.4 that $\mathbb{P}\{\mathcal{E}_{1}(r,l,\kappa)\}\geq 1-e^{-t}$ as long as $\tau\gtrsim\max(\sigma_{2},r)$ and $n\gtrsim(\tau/r)^{2}(s\log d+t)$. Furthermore, Proposition 2.5 ensures that event $\{\lambda\geq 2\|\bm{w}^{*}\|_{\infty}\}$ occurs with probability at least $1-e^{-t}$ if the regularization parameter satisfies

[TABLE]

Finally, we take $\tau\asymp\sigma_{2}\sqrt{n/(\log d+t)}$ and $r\asymp\tau$. By Lemma 6.1 and Condition 2.1, $b^{*}_{\tau}\leq\sigma_{\bm{x}}\sigma_{2}^{2}\tau^{-1}\asymp\sigma_{2}\sqrt{(\log d+t)/n}$. Putting together the pieces, we conclude that under the scaling $n\gtrsim s\log d+t$, and if $\lambda$ has magnitude of the order within the range

[TABLE]

the event $\mathcal{E}_{1}(r,6s^{1/2},\rho_{l}/2)\cap\{\lambda\geq 2\|\bm{w}^{*}\|_{\infty}\}$ occurs with probability at least $1-2e^{-t}$. Combined with (66), this proves the claimed bound. ∎

8.2 Proof of Theorem 2.2

We will apply Propositions 2.2, 2.4 and 2.5 to prove (32). The key is to control the random events from Proposition 2.2, which relies on a delicate combination of all the parameters. Similarly to the proof of Theorem 2.1, we take $\kappa=\rho_{l}/2$, and let $\lambda\geq s^{-1/2}b^{*}_{\tau}$, where $b^{*}_{\tau}=\|\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})\|_{2}\leq\sigma_{\bm{x}}\sigma_{2}^{2}\tau^{-1}$ due to Condition 2.1. Given $\gamma_{1}>\gamma_{0}>0$ satisfying (30), define $\delta=p^{\prime}(\gamma_{0})/(\rho_{l}\gamma_{0})\in(0,1)$, and let $c>0$ be the constant determined by $0.5p^{\prime}(\gamma_{0})(c^{2}+1)^{1/2}+4=c\rho_{l}\gamma_{0}$, which verifies (16) with $\kappa=\rho_{l}/2$. Moreover, set $l=\{2+\frac{2}{p^{\prime}(\gamma_{0})}\}(c^{2}+1)^{1/2}s^{1/2}+\frac{2}{p^{\prime}(\gamma_{0})}s^{1/2}$, and let

[TABLE]

Consequently, we apply Proposition 2.2 to conclude that, conditioned on event $\mathcal{E}_{1}(r,l,\rho_{l}/2)\cap\{\|\bm{w}^{*}\|_{\infty}+n^{-1/2}\leq 0.5p^{\prime}(\gamma_{0})\lambda\}$, the $\epsilon_{\ell}$-optimal solutions $\widetilde{}\bm{\beta}^{(\ell)}$ satisfy

[TABLE]

for all $\ell\geq 2$. Under the minimum signal strength condition $\|\bm{\beta}^{*}_{\mathcal{S}}\|_{\min}\geq(\gamma_{0}+\gamma_{1})\lambda$, and due to the fact that $p^{\prime}_{\lambda}(t)=0$ for all $t\geq\gamma_{1}\lambda$, the deterministic term $\|p_{\lambda}^{\prime}(|\bm{\beta}^{*}_{\mathcal{S}}|-\gamma_{0}\lambda)\|_{2}$ vanishes, thus implying

[TABLE]

for all $\ell\geq 2$.

Next, we control the event $\mathcal{E}_{1}(r,l,\rho_{l}/2)\cap\{\|\bm{w}^{*}\|_{\infty}+n^{-1/2}\leq 0.5p^{\prime}(\gamma_{0})\lambda\}$ and the oracle error term $\|\bm{w}^{*}_{{\mathcal{S}}}\|_{2}$. Given $t\geq 0$, it follows from Proposition 2.4 with the above $l\asymp s^{1/2}$ and $\kappa=\rho_{l}/2$ that event $\mathcal{E}_{1}(r,l,\kappa)$ occurs with probability at least $1-e^{-t}$ as long as $\tau\gtrsim\max(\sigma_{2},r)$ and $n\gtrsim(\tau/r)^{2}(s\log d+t)$. We therefore take $r\asymp\tau$ throughout the proof. Turning to the gradient vector $\bm{w}^{*}\in\mathbb{R}^{d}$, applying Proposition 2.5 yields that with probability at least $1-2e^{-t}$,

[TABLE]

Based on the above analysis, we choose the regularization parameter $\lambda=C\sigma_{2}\sqrt{(\log d+t)/n}$ for a sufficiently large $C$. Under the scaling $n\gtrsim s\log d+t$, and if $\tau$ has magnitude of the order within the range of $\sigma_{2}$ to $\sigma_{2}\sqrt{n/(\log d+t)}$, it follows from (67) that with probability at least $1-3e^{-t}$,

[TABLE]

This leads to the claimed bound by letting $\ell\gtrsim\lceil\log(\log d+t)/\log(1/\delta)\rceil$. ∎

8.3 Proof of Theorem 2.3

The proof is based primarily on Proposition 2.3, combined with complementary probabilistic analysis. For $t\geq 0$ and a prescribed $q\geq\max(s,\log d)$, set $\tau=\sigma_{2}\sqrt{n/(q+t)}$. In order to apply the high-level result in Proposition 2.3, we need the following two technical lemmas to control the events in (24). The former controls the event $\mathcal{E}_{2}(r,l,\kappa)$ defined in (21) under proper sample size requirement, and the latter provides upper bounds on the $\ell_{\infty}$-error terms $\|\widehat{}\bm{\beta}^{{\rm ora}}-\bm{\beta}^{*}\|_{\infty}$ and $\|\bm{w}^{{\rm ora}}\|_{\infty}=\|\nabla\widehat{\mathcal{L}}_{\tau}(\widehat{}\bm{\beta}^{{\rm ora}})\|_{\infty}$.

Lemma 8.1**.**

Under the conditions of the theorem, let $\tau,r,l>0$ satisfy

[TABLE]

where $C_{0}$ is an absolute constant and $C_{1}$ depends only on $\sigma_{\bm{x}}$. Then, with probability at least $1-e^{-t}$,

[TABLE]

holds uniformly over $(\bm{\beta}_{1},\bm{\beta}_{2})\in\mathcal{C}(r,l)$, where $\mathcal{C}(r,l)$ is defined in (21).

The next lemma provides statistical properties of the oracle estimator $\widehat{}\bm{\beta}^{{\rm ora}}$ defined in (20) with $\tau=\sigma_{2}\sqrt{n/(q+t)}$. Since the oracle $\widehat{}\bm{\beta}^{{\rm ora}}$ has access to the true active set ${\mathcal{S}}$, it is essentially an unpenalized Huber estimator based on $\{(y_{i},\bm{x}_{i,{\mathcal{S}}})\}_{i=1}^{n}$.

Lemma 8.2**.**

Under the sample size scaling $n\gtrsim q+t$, the following bounds

[TABLE]

and

[TABLE]

hold with probability at least $1-7e^{-t}$.

Compared to Propositions 2.4 and 2.5, the proofs of Lemmas 8.1 and 8.2, which are placed in the following two subsections, require a more delicate analysis of the local behavior of the gradient process $\{\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta}),\bm{\beta}\in\mathbb{R}^{d}\}$ around both the underlying vector $\bm{\beta}^{*}$ and the oracle estimator $\widehat{}\bm{\beta}^{{\rm ora}}$, with the latter being random itself.

With the above preparations, we are ready to prove the result. In Proposition 2.3, we set $\kappa=\rho_{l}/2$, $\delta=2.5/(\rho_{l}\gamma_{0})\in(0,1)$ and $l=\{2+\frac{2}{p^{\prime}(\gamma_{0})}\}(c_{0}^{2}+1)^{1/2}s^{1/2}$, where $c_{0}$ is determined by (23). Taking $r=\tau/C_{1}$ in Lemma 8.1, we obtain that event $\mathcal{E}_{2}(r,l,\kappa)$ happens with probability at least $1-e^{-t}$ as long as $n\gtrsim\sigma_{2}$ and $n\gtrsim s\log d+t$. Next, let $\lambda=C\sigma_{2}\sqrt{(\log d+t)/n}$ for a sufficient large constant $C$. Then, it follows from Lemma 8.2 that the event $\{\|\bm{w}^{{\rm ora}}\|_{\infty}\leq 0.5p^{\prime}(\gamma_{0})\lambda\}\cap\{\|\widehat{}\bm{\beta}^{{\rm ora}}-\bm{\beta}^{*}\|_{\infty}\leq\lambda/(5\delta\kappa)=\gamma_{0}\lambda/6.25\}\cap\{\|\widehat{}\bm{\beta}^{{\rm ora}}-\bm{\beta}^{*}\|_{\Sigma}\leq r\}$ occurs with probability at least $1-7e^{-t}$ as long as $n\gtrsim(q+t)^{1+1/\eta}(\log d)^{-1/\eta}$. Finally, the strong oracle property is a direct consequence of Proposition 2.3. ∎

8.3.1 Proof of Lemma 8.1

By the convexity of the loss function, for any $\bm{\beta}_{1},\bm{\beta}_{2}\in\mathbb{R}^{d}$,

[TABLE]

where $I_{\mathcal{E}_{i}}$ is the indicator function of the event $\mathcal{E}_{i}:=\{|\varepsilon_{i}|\leq\tau/4\}\cap\{|\bm{x}_{i}^{{\rm T}}(\bm{\beta}_{2}-\bm{\beta}^{*})|\leq\tau/4\}\cap\{|\bm{x}_{i}^{{\rm T}}(\bm{\beta}_{1}-\bm{\beta}_{2})|\leq\tau\|\bm{\beta}_{1}-\bm{\beta}_{2}\|_{2}/(2r)\}$ on which $|y_{i}-\bm{x}_{i}^{{\rm T}}\bm{\beta}_{2}|\leq|\varepsilon_{i}|+|\bm{x}_{i}^{{\rm T}}(\bm{\beta}_{2}-\bm{\beta}^{*})|\leq\tau/2$ and $|y_{i}-\bm{x}_{i}^{{\rm T}}\bm{\beta}_{1}|\leq|\bm{x}_{i}^{{\rm T}}(\bm{\beta}_{1}-\bm{\beta}_{2})|+|\bm{x}_{i}^{{\rm T}}(\bm{\beta}_{2}-\bm{\beta}^{*})|+|\varepsilon_{i}|\leq\tau$ for all $\bm{\beta}_{1}\in\bm{\beta}_{2}+\mathbb{B}(r)$. Similarly to the proof of Proposition 2 in Loh (2017), for any $R>0$, define Lipschitz continuous functions

[TABLE]

which are smoothed versions of $u\mapsto u^{2}I(|u|\leq R)$ and $u\mapsto I(|u|\leq R)$, respectively. Moreover, $\varphi_{R}(u)\leq u^{2}I(|u|\leq R)$ and $\phi_{R}(u)\leq I(|u|\leq R)$. By (72) and the fact that $\ell_{\tau}^{\prime\prime}(u)=1$ for $|u|\leq\tau$,

[TABLE]

In what follows, we deal with $\mathbb{E}D_{0}(\bm{\beta}_{1},\bm{\beta}_{2})$ and $D_{0}(\bm{\beta}_{1},\bm{\beta}_{2})-\mathbb{E}D_{0}(\bm{\beta}_{1},\bm{\beta}_{2})$, separately.

Noting that $\varphi_{R}(u)\geq u^{2}I(|u|\leq R/2)$ and $\phi_{R}(u)\geq I(|u|\leq R/2)$, we have

[TABLE]

Write $\bm{\delta}=\bm{\beta}_{1}-\bm{\beta}_{2}$ for $\bm{\beta}_{1}\in\bm{\beta}_{2}+\mathbb{B}(r)$ and $\bm{\beta}_{2}\in\bm{\beta}^{*}+\mathbb{B}_{\Sigma}(r)$. By Hölder’s inequality and (26),

[TABLE]

[TABLE]

and $\mathbb{E}(\bm{x}_{i}^{{\rm T}}\bm{\delta})^{2}I_{|\varepsilon_{i}|>\frac{\tau}{4}}\leq(4\sigma_{2}/\tau)^{2}\|\bm{\delta}\|_{\Sigma}^{2}$. Substituting these into (74) yields

[TABLE]

Let $\tau\geq 8\max(2\sigma_{2},\{\log(16\rho_{\bm{x}}^{2})\}^{1/2}\sigma_{\bm{x}}r,4\rho_{\bm{x}}^{2}r)$, so that

[TABLE]

It thus follows that

[TABLE]

uniformly over $\bm{\beta}_{1}\in\bm{\beta}_{2}+\mathbb{B}(r)$ and $\bm{\beta}_{2}\in\bm{\beta}^{*}+\mathbb{B}_{\Sigma}(r)$.

Next, we will establish a high probability bound for the supremum

[TABLE]

Note that $\varphi_{cR}(cu)=c^{2}\varphi_{R}(u)$ for any $c>0$ and $u\in\mathbb{R}$. For each pair $(\bm{\beta}_{1},\bm{\beta}_{2})$, we write $\bm{\delta}=\bm{\beta}_{1}-\bm{\beta}_{2}$ and define

[TABLE]

so that $\Delta(r,l)=\sup_{(\bm{\beta}_{1},\bm{\beta}_{2})\in\mathcal{C}(r,l)}|(1/n)\sum_{i=1}^{n}f_{\bm{\beta}_{1},\bm{\beta}_{2}}(\bm{x}_{i},\varepsilon_{i})-\mathbb{E}f_{\bm{\beta}_{1},\bm{\beta}_{2}}(\bm{x}_{i},\varepsilon_{i})|$. Since $0\leq\varphi_{R}(u)\leq\min\{(R/2)^{2},u^{2}\}$ and $0\leq\phi_{R}(u)\leq 1$ for all $u\in\mathbb{R}$, we have

[TABLE]

By Bousquet’s version of Talagrand’s inequality (Bousquet, 2003) and (63), for any $t\geq 0$,

[TABLE]

holds with probability at least $1-e^{-t}$.

It suffices to bound the expected value $\mathbb{E}\Delta(r,l)$. Applying the symmetrization inequality for empirical processes and the connection between Gaussian complexity and Rademacher complexity (see, e.g. Lemma 4.5 in Ledoux and Talagrand (1991)), we obtain that

[TABLE]

where

[TABLE]

with $\bm{\delta}=\bm{\beta}_{1}-\bm{\beta}_{2}$ and $g_{i}$’s are independent standard normal random variables that are independent of the observations. In particular, $\mathbb{G}_{\bm{\beta}^{*},\bm{\beta}^{*}}$ is defined as zero. Let $\mathbb{E}^{*}$ be the conditional expectation given $\{(y_{i},\bm{x}_{i})\}_{i=1}^{n}$. By symmetry,

[TABLE]

Next, we apply the Gaussian comparison theorem to bound $\mathbb{E}^{*}\{\sup_{(\bm{\beta}_{1},\bm{\beta}_{2})\in\mathcal{C}(r,l)}\mathbb{G}_{\bm{\beta}_{1},\bm{\beta}_{2}}\}$, from which an upper bound on $\mathbb{E}\{\sup_{(\bm{\beta}_{1},\bm{\beta}_{2})\in\mathcal{C}(r,l)}\mathbb{G}_{\bm{\beta}_{1},\bm{\beta}_{2}}\}$ follows immediately. For another pair $(\bm{\beta}_{1}^{\prime},\bm{\beta}_{2}^{\prime})\in\mathcal{C}(r,l)$, write $\bm{\delta}^{\prime}=\bm{\beta}_{1}^{\prime}-\bm{\beta}_{2}^{\prime}$, and note that

[TABLE]

By the Lipschitz properties of $\phi_{R}$ and $\varphi_{R}$, i.e., $|\phi_{R}(u)-\phi_{R}(v)|\leq\frac{2}{R}|u-v|$ and $|\varphi_{R}(u)-\varphi_{R}(v)|\leq R|u-v|$, and recall that $\varphi_{R}(u)\leq(R/2)^{2}$, we have

[TABLE]

and

[TABLE]

Motivated by (79), (80) and the inequality that

[TABLE]

we define another (conditional) Gaussian process $\{\mathbb{Z}_{\bm{\beta}_{1},\bm{\beta}_{2}},(\bm{\beta}_{1},\bm{\beta}_{2})\in\mathcal{C}(r,l)\}$ as

[TABLE]

where $g_{1}^{\prime},g_{1}^{\prime\prime},\ldots,g_{n}^{\prime},g_{n}^{\prime\prime}$ are independent standard normal random variables that are independent of all the other variables. We have established that $\mathbb{E}^{*}(\mathbb{G}_{\bm{\beta}_{1},\bm{\beta}_{2}}-\mathbb{G}_{\bm{\beta}^{\prime}_{1},\bm{\beta}^{\prime}_{2}})^{2}\leq\mathbb{E}^{*}(\mathbb{Z}_{\bm{\beta}_{1},\bm{\beta}_{2}}-\mathbb{Z}_{\bm{\beta}^{\prime}_{1},\bm{\beta}^{\prime}_{2}})^{2}$. Then, applying Sudakov-Fernique’s Gaussian comparison inequality (see, e.g. Theorem 7.2.11 in Vershynin (2018)) yields

[TABLE]

which remains valid if $\mathbb{E}^{*}$ is replaced by $\mathbb{E}$. For the supremum of $\mathbb{Z}_{\bm{\beta}_{1},\bm{\beta}_{2}}$, it is easy to see that

[TABLE]

Together, (77), (78), (81) and (82) deliver the bound

[TABLE]

Finally we bound the maximum under expectation on the right-hand side of (83). Write $S_{j}=\sum_{i=1}^{n}g_{i}x_{ij}$ for $j=1,\ldots,d$. Under Condition 2.3, for each $1\leq j\leq d$ and $m\geq 2$ we have

[TABLE]

Let $g\sim{\sf N}(0,1)$ be independent of $\bm{x}$. Using the Legendre duplication formula, i.e., $\Gamma(s)\Gamma(s+1/2)=2^{1-2s}\sqrt{\pi}\,\Gamma(2s)$, and some algebra, we get

[TABLE]

Hence, using Bernstein’s inequality and the symmetry of normal distribution yields

[TABLE]

for all $\lambda\in(0,1/\sigma_{\bm{x}})$. Combined with Theorem 2.5 in Boucheron, Lugosi and Massart (2013), this implies

[TABLE]

Combining (83), (84) with the concentration inequality (76), we determine that with probability at least $1-e^{-t}$, $\Delta(r,l)\leq\rho_{l}/4$ as long as $n\gtrsim(\tau/r)^{2}(s+l^{2}\log d+t)$. This, together with (73) and (75), proves the claim (69). ∎

8.3.2 Proof of Lemma 8.2

To begin with, consider the decomposition

[TABLE]

where $\bm{w}(\bm{\beta})=\nabla\widehat{\mathcal{L}}_{\tau}(\bm{\beta})-\nabla\mathcal{L}_{\tau}(\bm{\beta})$, $\mathcal{L}_{\tau}(\bm{\beta})=\mathbb{E}\widehat{\mathcal{L}}_{\tau}(\bm{\beta})$ and $\bm{w}^{*}=\bm{w}(\bm{\beta}^{*})$. In the following, we control the $\ell_{\infty}$-norms of the three terms, $\bm{w}(\widehat{}\bm{\beta}^{{\rm ora}})-\bm{w}(\bm{\beta}^{*})$, $\nabla\mathcal{L}_{\tau}(\widehat{}\bm{\beta}^{{\rm ora}})$ and $\bm{w}^{*}$, separately. Throughout the proof, we take $\tau=\sigma_{2}\sqrt{n/(q+t)}$ for some $q\geq\max(s,\log d)$ and $t\geq 0$.

Applying Proposition 2.5 to the centered gradient $\bm{w}^{*}=((\bm{w}^{*}_{{\mathcal{S}}})^{{\rm T}},(\bm{w}^{*}_{{\mathcal{S}}^{{\rm c}}}))^{{{\rm T}}})^{{\rm T}}\in\mathbb{R}^{d}$ with slight modifications, we obtain that with probability at least $1-2e^{-t}$,

[TABLE]

thus implying $\|\bm{w}^{*}\|_{\infty}\lesssim\sigma_{2}\sqrt{(\log d+t)/n}$ with the same probability.

Recall that $\widehat{}\bm{\beta}^{{\rm ora}}$ and $\bm{\beta}^{*}$ have the same support ${\mathcal{S}}\subseteq[d]$. Define the oracle local neighborhood $\Theta^{*}(r)=\{\bm{\beta}\in\bm{\beta}^{*}+\mathbb{B}_{\Sigma}(r):\bm{\beta}_{{\mathcal{S}}^{\rm c}}=\textbf{0}\}$. Then, conditioned on the event $\{\|\widehat{}\bm{\beta}^{{\rm ora}}-\bm{\beta}^{*}\|_{\Sigma}\leq r\}$,

[TABLE]

We thus focus on the supremum on the right-hand side of (86). For every $s$-sparse vector $\bm{\beta}\in\Theta^{*}(r)$, we write $\bm{\delta}=(\bm{\beta}-\bm{\beta}^{*})_{{\mathcal{S}}}\in\mathbb{R}^{s}$. For $j=1,\ldots,d$, let $\bm{e}_{j}\in\mathbb{R}^{d}$ be the coordinate vector that has 1 on its $j$-th coordinate and 0 elsewhere, and define $\Delta^{0}_{j}(\bm{\delta})=\langle\bm{w}(\bm{\beta})-\bm{w}(\bm{\beta}^{*}),\bm{e}_{j}\rangle=(1/n)\sum_{i=1}^{n}(\eta_{ij}-\mathbb{E}\eta_{ij})$ for $\bm{\delta}\in\mathbb{R}^{s}$, where $\eta_{ij}=x_{ij}\{\ell_{\tau}^{\prime}(\varepsilon_{i}-\bm{x}_{i,{\mathcal{S}}}^{{\rm T}}\bm{\delta})-\ell_{\tau}^{\prime}(\varepsilon_{i})\}$. Consequently, we have

[TABLE]

where ${\rm S}=\Sigma_{{\mathcal{S}}{\mathcal{S}}}\in\mathbb{R}^{s\times s}$. In order to bound the local fluctuation $\sup_{\bm{\delta}:\|\bm{\delta}\|_{{\rm S}}\leq r}\Delta^{0}_{j}(\bm{\delta})$, we need to control the moment generating function of $\Delta^{0}_{j}(\bm{\delta})$ for each $\bm{\delta}\in\mathbb{R}^{s}$. By the Lipschitz continuity of $\ell_{\tau}^{\prime}(\cdot)$, $|\mathbb{E}(\eta_{ij})|\leq\mathbb{E}|x_{ij}\bm{x}_{i,{\mathcal{S}}}^{{\rm T}}\bm{\delta}|\leq\sigma_{jj}^{1/2}\|\bm{\delta}\|_{{\rm S}}$, $\mathbb{E}(\eta_{ij}^{2}|\bm{x}_{i})\leq x_{ij}^{2}(\bm{x}_{i,{\mathcal{S}}}^{{\rm T}}\bm{\delta})^{2}$ and

[TABLE]

The above moment inequalities, combined with the elementary inequality $|e^{u}-1-u|\leq(u^{2}/2)e^{u\vee 0}$, imply that for any $\lambda\in\mathbb{R}$ and $\lambda^{*}=\lambda/(\sigma_{jj}^{1/2}\|\bm{\delta}\|_{{\rm S}})$,

[TABLE]

Applying Hölder’s inequality to the exponential moments on the right-hand side of (88), we have

[TABLE]

and

[TABLE]

Substituting these bounds into the earlier inequality (88), we find that for any $|\lambda|\leq\sqrt{n}/C_{1}$,

[TABLE]

where $C_{1},C_{2}>0$ depend only on $\upsilon_{1}$ in Condition 2.3. A similar argument can be used to establish the same bound for each pair $(\bm{\delta},\bm{\delta}^{\prime})$, that is,

[TABLE]

The above inequality certifies condition ($\mathcal{E}d$) in Spokoiny (2012) (see Section 2 in the supplement), so that Corollary 2.2 therein applies to the process $\{\Delta^{0}_{j}(\bm{\delta})\}_{\bm{\delta}\in\mathbb{R}^{s}:\|\bm{\delta}\|_{{\rm S}}\leq r}$: with probability at least $1-e^{-x}$,

[TABLE]

as long as $n\gtrsim s+x$. Combined with (87) and the union bound, we find that

[TABLE]

with probability at least $1-2de^{-x}$ provided $n\gtrsim s+x$. Taking $x=\log(2d)+t$, it follows from (86) that conditioned on $\{\|\widehat{}\bm{\beta}^{{\rm ora}}-\bm{\beta}^{*}\|_{\Sigma}\leq r\}$,

[TABLE]

holds with probability at least $1-e^{-t}$ as long as $n\gtrsim s+\log d+t$.

Tuning to $\|\nabla\mathcal{L}_{\tau}(\widehat{}\bm{\beta}^{{\rm ora}})\|_{\infty}$, again, we control this term conditioned on the same event above. Following the proof of Lemma 6.1, it can be similarly shown that

[TABLE]

For any $\bm{\beta}\in\Theta^{*}(r)$, write $\bm{\delta}=(\bm{\beta}-\bm{\beta}^{*})_{{\mathcal{S}}}\in\mathbb{R}^{s}$, and note that

[TABLE]

Let $\mathbb{E}_{\bm{x}}$ and $\mathbb{P}_{\bm{x}}$ be the conditional expectation and probability given $\bm{x}$, respectively. By the anti-concentration property (34) of the distribution of $\varepsilon$ given $\bm{x}$, we see that for any $u\in\mathbb{R}$,

[TABLE]

Together, the last two displays imply

[TABLE]

where ${\rm H}_{\cdot{\mathcal{S}}}:=\mathbb{E}\{\ell^{\prime\prime}_{\tau}(\varepsilon)\bm{x}\bm{x}_{{\mathcal{S}}}^{{\rm T}}\}\in\mathbb{R}^{d\times s}$ is the submatrix of ${\rm H}=\nabla^{2}\mathcal{L}_{\tau}(\bm{\beta}^{*})=\mathbb{E}\{\ell^{\prime\prime}_{\tau}(\varepsilon)\bm{x}\bm{x}^{{\rm T}}\}$. For the linear term ${\rm H}_{\cdot{\mathcal{S}}}\bm{\delta}$, write $\Sigma_{\cdot{\mathcal{S}}}=\mathbb{E}(\bm{x}\bm{x}^{{\rm T}}_{{\mathcal{S}}})$ and note that

[TABLE]

Together, (90), (91) and (92) imply that conditioned on $\{\|\widehat{}\bm{\beta}^{{\rm ora}}-\bm{\beta}^{*}\|_{\Sigma}\leq r\}$,

[TABLE]

Next we consider the oracle estimator $\widehat{}\bm{\beta}^{{\rm ora}}$ with $\tau=\sigma_{2}\sqrt{n/(q+t)}$. Following an argument similar to that used to prove Theorem 2.1 in Chen and Zhou (2020), it can be shown that with probability at least $1-3e^{-t}$,

[TABLE]

and

[TABLE]

where ${\rm S}=\mathbb{E}(\bm{x}_{{\mathcal{S}}}\bm{x}_{{\mathcal{S}}}^{{\rm T}})=\Sigma_{{\mathcal{S}}{\mathcal{S}}}$. Note that the linear term $(1/n)\sum_{i=1}^{n}\ell^{\prime}_{\tau}(\varepsilon_{i})\bm{x}_{i,{\mathcal{S}}}$ in the Bahadur representation bound (95) can be written as $(1/n)\sum_{i=1}^{n}\ell^{\prime}_{\tau}(\varepsilon_{i})\bm{x}_{i,{\mathcal{S}}}=-\bm{w}_{{\mathcal{S}}}^{*}-\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})_{{\mathcal{S}}}$. It follows that

[TABLE]

Similarly to Lemma 6.1, we obtain that $\|{\rm S}^{-1/2}\nabla\mathcal{L}_{\tau}(\bm{\beta}^{*})_{{\mathcal{S}}}\|_{2}\leq(a_{\eta}\sigma_{2})^{2+\eta}\tau^{-1-\eta}$. For $\|{\rm S}^{-1}\bm{w}_{{\mathcal{S}}}^{*}\|_{\infty}$, following the proof of Proposition 2.5, it can be similarly shown that with probability at least $1-e^{-t}$,

[TABLE]

Putting together the pieces, we conclude that the $\ell_{2}$-error bound (94) and the $\ell_{\infty}$-error bound

[TABLE]

hold with probability $1-4e^{-t}$ as long as $n\gtrsim q+t$. Combined with (94), this proves (70).

Finally, it remains to deal with $\|\Sigma_{\cdot{\mathcal{S}}}(\widehat{}\bm{\beta}^{{\rm ora}}-\bm{\beta}^{*})_{{\mathcal{S}}}\|_{\infty}$. Under condition (33),

[TABLE]

Using the previous bounds (85), (90) and (95), we obtain that

[TABLE]

Combining this bound with (89), (93) and (94) yields the claim (71). ∎

8.4 Proof of Theorem 3.1

For simplicity, we write $\bm{\beta}^{(k)}=\bm{\beta}^{(1,k)}$, $\phi^{(k)}=\phi^{(1,k)}$ and $\bm{\lambda}=\bm{\lambda}^{(0)}$ throughout this section.

8.4.1 Technical lemmas

We first present three technical lemmas, which are the key ingredients to the proof. The first lemma provides an alternative to the stopping rule.

Lemma 8.3**.**

$\omega_{\bm{\lambda}}(\bm{\beta}^{(k)})\leq L(1+\gamma_{u})\|\bm{\beta}^{(k)}-\bm{\beta}^{(k-1)}\|_{2}.$

Proof of Lemma 8.3.

For simplicity, we write $\mathcal{L}(\cdot)=\widehat{\mathcal{L}}_{\tau}(\cdot)$ as the loss function of interest. Since $\bm{\beta}^{(k)}$ is the exact solution at the $k$-th iteration when $\ell=1$, the first-order optimality condition holds: there exists some $\bm{\xi}^{(k)}\in\partial\|\bm{\beta}^{(k)}\|_{1}$ such that

[TABLE]

For any $\bm{u}\in\mathbb{R}^{d}$ such that $\|\bm{u}\|_{1}=1$, we have

[TABLE]

where the last inequality is due to the Lipschitz continuity of $\nabla\mathcal{L}(\cdot)$. Taking the supremum over all $\bm{u}$ satisfying $\|\bm{u}\|_{1}\leq 1$, we obtain

[TABLE]

It remains to show that $\phi^{(k)}\leq L\gamma_{u}$ for any $k$. This is guaranteed by the iterative LAMM algorithm. Otherwise, if $\phi^{(k)}>L\gamma_{u}$, then $\phi^{\prime}\equiv\phi^{(k)}/\gamma_{u}>L$ is the quadratic parameter in the previous iteration for searching $\phi$ such that

[TABLE]

where $\widetilde{}\bm{\beta}^{(k)}$ is the new updated parameter vector under the quadratic coefficient $\phi^{\prime}$. On the other hand, it follows from the definition of $F$ and the Lipschitz continuity of $\nabla\mathcal{L}$ that

[TABLE]

This leads to a contradiction, indicating that $\phi^{(k)}\leq L\gamma_{u}$. ∎

The second lemma is a modified version of Lemma E.4 in Fan et al. (2018). We reproduce its proof here for completeness. Let $\Psi(\bm{\beta},\bm{\lambda})=\mathcal{L}(\bm{\beta})+\|\bm{\lambda}\circ\bm{\beta}\|_{1}$ with $\bm{\lambda}=\bm{\lambda}^{(0)}$.

Lemma 8.4**.**

For any $\bm{\beta}\in\mathbb{R}^{d}$, we have

[TABLE]

Proof of Lemma 8.4.

Since $F(\bm{\beta};\phi^{(k)},\bm{\beta}^{(k-1)})$ majorizes $\mathcal{L}(\bm{\beta})$ at $\bm{\beta}^{(k)}$, we have

[TABLE]

By the convexity of $\mathcal{L}(\cdot)$ and $\bm{\beta}\mapsto\|\bm{\lambda}\circ\bm{\beta}\|_{1}$,

[TABLE]

for any $\bm{\xi}^{(k)}\in\partial\|\bm{\beta}^{(k)}\|_{1}$. This further implies

[TABLE]

Plugging the expression of $F(\bm{\beta}^{(k)};\phi^{(k)},\bm{\beta}^{(k-1)})$ in (38) and (97) into (96), we obtain

[TABLE]

By the first-order optimality condition, there exists some $\bm{\xi}\in\partial\|\bm{\beta}^{(k)}\|_{1}$ such that

[TABLE]

Substituting this into (8.4.1) proves the claimed bound. ∎

Recall that $\Psi(\bm{\beta},\bm{\lambda})=\mathcal{L}(\bm{\beta})+\|\bm{\lambda}\circ\bm{\beta}\|_{1}$ and $\widehat{}\bm{\beta}^{(1)}\in\min_{\bm{\beta}\in\mathbb{R}^{d}}\Psi(\bm{\beta},\bm{\lambda})$ denotes the optimal solution in the contraction stage.

Lemma 8.5**.**

For any $k\geq 1$, we have

[TABLE]

Proof of Lemma 8.5.

For simplicity, we write $\widehat{}\bm{\beta}=\widehat{}\bm{\beta}^{(1)}$, and define $\phi_{\max}=\max_{1\leq j\leq k}\phi^{(j)}$ and $\phi_{\min}=\min_{1\leq j\leq k}\phi^{(j)}>0$. Taking $\bm{\beta}=\widehat{\bm{\beta}}$ in Lemma 8.4 gives

[TABLE]

for all $j\geq 1$. Summing over $j$ from 1 to $k$ yields

[TABLE]

which further implies

[TABLE]

Again, by Lemma 8.4 with $\bm{\beta}=\bm{\beta}^{(j-1)}$ and $k=j$,

[TABLE]

Multiplying both sides of the above inequality by $j-1$ and summing over $j$, we obtain

[TABLE]

or equivalently,

[TABLE]

Together, (99) and (100) imply

[TABLE]

from which it follows immediately that

[TABLE]

This completes the proof. ∎

8.4.2 Proof of the theorem

Recall that $\bm{\beta}^{(k)}=\bm{\beta}^{(1,k)}$ and $\phi^{(k)}=\phi^{(1,k)}$. By Lemma 8.3 and its proof,

[TABLE]

Next, taking $\bm{\beta}=\bm{\beta}^{(k-1)}$ in Lemma 8.4 yields

[TABLE]

Together, the last two displays lead to a bound for the suboptimality measure

[TABLE]

Recall that $\{\Psi(\bm{\beta}^{(k)},\bm{\lambda})\}_{k=0}^{\infty}$ is a non-increasing sequence, i.e.,

[TABLE]

Then, it follows from (101) and Lemma 8.5 that

[TABLE]

where we used the fact that $\bm{\beta}^{(0)}={\bf 0}$. By the triangle inequality,

[TABLE]

Therefore, in the contraction stage, we need $k\gtrsim\{L(1+\gamma_{u})(\|\bm{\beta}^{*}\|_{2}+\|\widehat{}\bm{\beta}-\bm{\beta}^{*}\|_{2})/\epsilon_{{\rm c}}\}^{2}$ to ensure $\omega_{\bm{\lambda}^{(0)}}(\bm{\beta}^{(k)})\leq\epsilon_{{\rm c}}$. This proves the stated result. ∎

8.5 Proof of Theorem 3.2

For convenience, we omit the index $\ell$, and use $\widehat{}\bm{\beta}$, $\bm{\beta}^{(k)}$, $\bm{\lambda}$ and $\mathcal{E}$ to denote $\widehat{}\bm{\beta}^{(\ell)}$, $\bm{\beta}^{(\ell,k)},\bm{\lambda}^{(\ell-1)}$ and $\mathcal{E}_{\ell}$, respectively, where $\mathcal{E}_{\ell}$ is the subset defined in (53) satisfying ${\mathcal{S}}\subseteq\mathcal{E}_{\ell}$ and $|\mathcal{E}_{\ell}|\leq C_{0}s$ for some constant $C_{0}>1$. Moreover, write $\mathcal{L}(\cdot)=\widehat{\mathcal{L}}_{\tau}(\cdot)$, and define $\Psi(\bm{\beta},\bm{\lambda})=\mathcal{L}(\bm{\beta})+\|\bm{\lambda}\circ\bm{\beta}\|_{1}=\mathcal{L}(\bm{\beta})+\|\bm{\lambda}^{(\ell-1)}\circ\bm{\beta}\|_{1}$, so that $\widehat{}\bm{\beta}\in\min_{\bm{\beta}}\Psi(\bm{\beta},\bm{\lambda})$.

8.5.1 Technical lemmas

We first provide several technical lemmas along with the proofs.

Lemma 8.6**.**

For any $k$-sparse ($k\geq 1$) vectors $\bm{\beta}_{1},\bm{\beta}_{2}\in\bm{\beta}^{*}+\mathbb{B}(r)$, we have

[TABLE]

where $D_{\mathcal{L}}(\bm{\beta}_{1},\bm{\beta}_{2}):=\mathcal{L}(\bm{\beta}_{1})-\mathcal{L}(\bm{\beta}_{2})-\langle\nabla\mathcal{L}(\bm{\beta}_{2}),\bm{\beta}_{1}-\bm{\beta}_{2}\rangle$ is the Bregman divergence.

Proof of Lemma 8.6.

By a second-order Taylor series expansion, there exists some $\gamma\in[0,1]$ such that $\widetilde{}\bm{\beta}=\gamma\bm{\beta}_{1}+(1-\gamma)\bm{\beta}_{2}\in\bm{\beta}^{*}+\mathbb{B}(r)$ and $D_{\mathcal{L}}(\bm{\beta}_{1},\bm{\beta}_{2})=(1/2)(\bm{\beta}_{1}-\bm{\beta}_{2})^{{\rm T}}\nabla^{2}\mathcal{L}(\widetilde{}\bm{\beta})(\bm{\beta}_{1}-\bm{\beta}_{2})$. The stated bounds then follow directly from Definition 3.1. ∎

The next lemma converts the bound on $\Psi(\bm{\beta},\bm{\lambda})-\Psi(\bm{\beta}^{*},\bm{\lambda})$ to that on $\|\bm{\beta}-\bm{\beta}^{*}\|_{2}$. Recall that for any subset $\mathcal{E}\subseteq[d]$, we write $\bm{\beta}_{\mathcal{E}}$ as a subvector of $\bm{\beta}$ indexed by $\mathcal{E}$.

Lemma 8.7**.**

Assume LSE$(1)$ condition holds. Let $\mathcal{E}\subseteq[d]$ be a subset satisfying ${\mathcal{S}}\subseteq\mathcal{E}$ and $|\mathcal{E}|\leq C_{0}s$ for some $C_{0}\geq 1$. Assume further that $\lambda\geq\max\{4\|\nabla\mathcal{L}(\bm{\beta}^{*})\|_{\infty},\|\bm{\lambda}\|_{\infty}\}$ and $\|\bm{\lambda}_{\mathcal{E}^{{\rm c}}}\|_{\min}\geq\lambda/2$. Then, for any $\bm{\beta}\in\bm{\beta}^{*}+\mathbb{B}(r)$ satisfying $\|\bm{\beta}_{{\mathcal{S}}^{{\rm c}}}\|_{0}\leq s^{\prime}$ and $\Psi(\bm{\beta},\bm{\lambda})-\Psi(\bm{\beta}^{*},\bm{\lambda})\leq Cs\lambda^{2},$ we have

[TABLE]

where $C_{1},C_{2}>0$ depend only on $C_{0},C$ and localized sparse eigenvalues.

Proof of Lemma 8.7.

We omit the arguments in $\kappa_{-}(s+s^{\prime},r,\tau)$ and $\kappa_{+}(s+s^{\prime},r,\tau)$ whenever there is no ambiguity. For any $\bm{\beta}\in\bm{\beta}^{*}+\mathbb{B}(r)$ satisfying $\|\bm{\beta}_{{\mathcal{S}}^{{\rm c}}}\|_{0}\leq s^{\prime}$, note that $\|\bm{\beta}\|_{0}\leq s+s^{\prime}$ and $\|\bm{\beta}-\bm{\beta}^{*}\|_{0}\leq s+s^{\prime}$. Using Lemma 8.6 yields

[TABLE]

Since $\Psi(\bm{\beta})-\Psi(\bm{\beta}^{*})\leq Cs\lambda^{2},$ or equivalently,

[TABLE]

it follows

[TABLE]

After some simple algebra, it can be derived that

[TABLE]

Combining the above bounds gives

[TABLE]

which further implies

[TABLE]

To bound the right-hand side of the above inequality, we discuss two cases regarding the magnitude of $\|(\bm{\beta}-\bm{\beta}^{*})_{\mathcal{E}}\|_{1}$ as compared to $s\lambda$:

• •

If $5\lambda\|(\bm{\beta}-\bm{\beta}^{*})_{\mathcal{E}}\|_{1}/4\leq Cs\lambda^{2}$, we have

[TABLE]

• •

If $5\lambda\|(\bm{\beta}-\bm{\beta}^{*})_{\mathcal{E}}\|_{1}/4>Cs\lambda^{2}$, we have

[TABLE]

thus implying

[TABLE]

Combining (103) and (104), we obtain

[TABLE]

Since $\bm{\beta}-\bm{\beta}^{*}$ is at most $(s+s^{\prime})$-sparse, $\|\bm{\beta}-\bm{\beta}^{*}\|_{1}\leq(s+s^{\prime})^{1/2}\|\bm{\beta}-\bm{\beta}^{*}\|_{2}$. The stated results then follow immediately. ∎

Recall that $\mathcal{E}_{\ell}$ is the subset defined in (53) satisfying ${\mathcal{S}}\subseteq\mathcal{E}_{\ell}$ and $|\mathcal{E}_{\ell}|\leq C_{0}s$ for some $C_{0}>1$.

Lemma 8.8**.**

Assume LSE$(C_{0})$ condition holds and $4\{\|\nabla\mathcal{L}(\bm{\beta}^{*})\|_{\infty}+\epsilon_{{\rm c}}\vee\epsilon_{{\rm t}}\}\leq\lambda\lesssim s^{-1/2}r$. For any $\ell\geq 2$, the solution sequence $\{\bm{\beta}^{(\ell,k)}\}_{k\geq 0}$ satisfies

[TABLE]

where $C_{1},C_{2}>0$ are constants depending only on the localized sparse eigenvalues.

Proof of Lemma 8.8.

We prove the theorem by the method of induction on $(\ell,k)$. Throughout, $C$ denotes a constant independent of $(n,d,s)$ and may take different values at each appearance. For the first subproblem, directly applying Proposition 4.1 and Lemma 5.4 in Fan et al. (2018) we obtain that $\|\widetilde{}\bm{\beta}^{(1)}-\bm{\beta}\|_{2}\leq C\kappa_{*}^{-1}s^{1/2}\lambda<r$, $\|\widetilde{}\bm{\beta}^{(1)}-\bm{\beta}\|_{1}\leq C\kappa_{*}^{-1}s\lambda$ and $\widetilde{}\bm{\beta}^{(1)}$ is $(s+s^{\prime})$-sparse, where $s^{\prime}\leq Cs$. It follows that $\bm{\beta}^{(2,0)}=\widetilde{}\bm{\beta}^{(1)}$ falls in a localized sparse set.

To apply the method of induction, first we assume that for any $k$, $\bm{\beta}^{(2,k)}$ falls in a localized sparse set such that (105) holds. We then use Lemma E.13 in Fan et al. (2018) to show that $\bm{\beta}^{(2,k+1)}$ also falls in a localized sparse set. To this end, we need to verify two conditions. The first one, $\|\bm{\lambda}^{(\ell)}_{\mathcal{E}_{\ell}^{{\rm c}}}\|_{\min}\geq\lambda/2$ is guaranteed by Claim (54) in the proof of Proposition 2.2, when $\gamma$ is such that $p^{\prime}(\gamma)=1/2$ and $|\mathcal{E}_{l}|\leq C_{0}s$ for some $C_{0}>1$. For the second condition, it suffices to show

[TABLE]

where $\zeta=\kappa^{*}/\kappa_{*}$. Using the mean value theorem, there exists some convex combination of $\bm{\beta}^{(2,k)}$ and $\bm{\beta}^{*}$, say $\widetilde{}\bm{\beta}$, such that

[TABLE]

With above preparations, it follows from Lemma E.13 in Fan et al. (2018) with slight modification that $\|\bm{\beta}_{\mathcal{E}_{\ell}^{{\rm c}}}^{(2,k+1)}\|_{0}\leq s^{\prime}$ and hence $\|\bm{\beta}^{(2,k+1)}\|_{0}\leq C_{0}s+s^{\prime}$.

Next, we show that $\|\bm{\beta}^{(2,k+1)}-\bm{\beta}^{*}\|_{2}\lesssim\kappa_{*}^{-1}s^{1/2}\lambda$. Again, by Lemma 8.4,

[TABLE]

This implies that $\{\Psi(\bm{\beta}^{(2,k)},\bm{\lambda}^{(1)})-\Psi(\bm{\beta}^{*},\bm{\lambda}^{(1)})\}_{k\geq 1}$ is a non-increasing sequence. By induction, it follows that

[TABLE]

Combining this with Lemma 8.7 gives the desired bounds on $\|\bm{\beta}^{(2,k+1)}-\bm{\beta}^{*}\|_{2}$ and $\|\bm{\beta}^{(2,k+1)}-\bm{\beta}^{*}\|_{1}$.

Finally, by an argument similar to that in the proof of Lemma 5.4 in Fan et al. (2018), we can derive the stated results for all $\ell\geq 3$. ∎

For $\epsilon>0$, let $\widetilde{}\bm{\beta}$ be an $\epsilon$-optimal solution to the program $\min_{\bm{\beta}}\{\mathcal{L}_{\tau}(\bm{\beta})+\|\bm{\lambda}\circ\bm{\beta}\|_{1}\}.$ The following lemma provides conditions under which $\widetilde{}\bm{\beta}$ falls in an $\ell_{1}$-cone.

Lemma 8.9**.**

Let $\mathcal{E}\subseteq[d]$ be a subset satisfying ${\mathcal{S}}\subseteq\mathcal{E}$, and assume $\lambda\geq\|\bm{\lambda}\|_{\infty}\vee 4\{\|\nabla\mathcal{L}(\bm{\beta}^{*})\|_{\infty}+\epsilon\}$ and $\|\bm{\lambda}_{\mathcal{E}^{{\rm c}}}\|_{\min}\geq\lambda/2$. Then, any $\epsilon$-optimal solution $\widetilde{}\bm{\beta}$ satisfies the cone constraint

[TABLE]

Proof of Lemma 8.9.

For any $\bm{\xi}\in\partial\|\widetilde{}\bm{\beta}\|_{1}$, let $\bm{u}=\nabla\mathcal{L}(\widetilde{\bm{\beta}})+\bm{\lambda}\circ\bm{\xi}$. By the convexity of $\mathcal{L}(\cdot)$, $\langle\nabla\mathcal{L}(\widetilde{}\bm{\beta})-\nabla\mathcal{L}(\bm{\beta}^{*}),\widetilde{}\bm{\beta}-\bm{\beta}^{*}\rangle\geq 0$. This, together with the inequality $\langle\nabla\mathcal{L}(\widetilde{}\bm{\beta})+\bm{\lambda}\circ\bm{\xi},\widetilde{}\bm{\beta}-\bm{\beta}^{*}\rangle\leq\|\bm{u}\|_{\infty}\|\widetilde{}\bm{\beta}-\bm{\beta}^{*}\|_{1}$, implies

[TABLE]

For I and II, note that $\text{I}\geq-\|\nabla\mathcal{L}(\bm{\beta}^{*})\|_{\infty}\|\widetilde{\bm{\beta}}-\bm{\beta}\|_{1},$ and

[TABLE]

Substituting the above bounds into (106) and taking the infimum over $\bm{\xi}\in\partial\|\widetilde{}\bm{\beta}\|_{1}$ yields

[TABLE]

or equivalently,

[TABLE]

This proves the stated result. ∎

8.5.2 Proof of the theorem

Restricting our attention to the $\ell$-th subproblem, we write $\phi^{(k)}=\phi^{(\ell,k)}$ for simplicity. Define the subset $\mathbb{L}=\{\alpha\widehat{\bm{\beta}}+(1-\alpha)\bm{\beta}^{(k-1)}:0\leq\alpha\leq 1\}$. Due to local majorization, we have

[TABLE]

where we used the convexity of $\mathcal{L}(\cdot)$ in the second inequality. Since $\Psi(\bm{\beta},\bm{\lambda})=\mathcal{L}(\bm{\beta})+\|\bm{\lambda}\circ\bm{\beta}\|_{1}$ is minimized at $\widehat{\bm{\beta}}$, by convexity we have

[TABLE]

Next, we bound the right-hand side of (107). By Lemma 8.8,

[TABLE]

Similarly, it can be shown the the optimum $\widehat{}\bm{\beta}$ satisfies the same properties. Hence,

[TABLE]

By the first-order optimality condition, there exists some $\widehat{\bm{\xi}}\in\partial\|\widehat{}\bm{\beta}\|_{1}$ such that $\nabla\mathcal{L}(\widehat{\bm{\beta}})+\bm{\lambda}\circ\widehat{\bm{\xi}}={\bf 0}$. Moreover, define $D_{\mathcal{L}}(\bm{\beta}_{1},\bm{\beta}_{2})=\mathcal{L}(\bm{\beta}_{1})-\mathcal{L}(\bm{\beta}_{2})-\langle\nabla\mathcal{L}(\bm{\beta}_{2}),\bm{\beta}_{1}-\bm{\beta}_{2}\rangle$. Using Definition 3.1, Lemma 8.6, and the convexity of $\mathcal{L}(\cdot)$ and $\ell_{1}$-norm, we obtain that

[TABLE]

where $\kappa_{-}=\kappa_{-}(C_{0}s+2s^{\prime},r,\tau)$. Plugging this bound into (107) yields

[TABLE]

Following the proof of Lemma 8.3, it can be similarly shown that $\phi^{(k)}\leq\gamma_{u}\kappa^{*}$ under Condition 3.2. Consequently,

[TABLE]

where $\zeta={\kappa^{*}}/{\kappa_{*}}$.

By an argument similar to that in the proof of Lemma 8.3, we can show that, for $\ell\geq 2$,

[TABLE]

Further, using Lemma 8.4 to bound $\|\bm{\beta}^{(\ell,k)}-\bm{\beta}^{(\ell,k-1)}\|_{2}$ from above and noting that $\phi^{(k)}\geq\kappa_{*}$, we obtain

[TABLE]

where the last step applies Lemmas 8.4 and 8.8.

To make the right-hand side of the above inequality smaller than $\epsilon_{{\rm t}}$, we need $k$ to be sufficiently large that $k\geq C_{1}\log(C_{2}s^{1/2}\lambda/\epsilon_{{\rm t}})$, where $C_{1},C_{2}>0$ are constants depending only on localized sparse eigenvalues and $\gamma_{u}$. This completes the proof. ∎