Robust Estimation of High-Dimensional Mean Regression

Introduction

Our era has witnessed the massive explosion of data and a dramatic improvement of technology in collecting and processing large data sets. We often encounter huge data sets that the number of features greatly surpasses the number of observations. It makes many traditional statistical analysis tools infeasible and poses great challenge on developing new tools. Regularization methods have been widely used for the analysis of high-dimensional data. These methods penalize the least squares or the likelihood function with the $L_{1}$-penalty on the unknown parameters (Lasso, Tibshirani (1996)), or a folded concave penalty function such as the SCAD (Fan and Li, 2001) and the MCP(Zhang, 2010). However, these penalized least-squares methods are sensitive to the tails of the error distributions, particularly for ultrahigh dimensional covariates, as the maximum spurious correlation between the covariates and the realized noises can be large in those cases. As a result, theoretical properties are often obtained under light-tailed error distributions (Bickel, Ritov and Tsybakov, 2009; Fan and Lv, 2011).

To tackle the problem of heavy-tailed errors, robust regularization methods have been extensively studied. Li and Zhu (2008), Wu and Liu (2009) and Zou and Yuan (2008) developed robust regularized estimators based on quantile regression for the case of fixed dimensionality. Belloni and Chernozhukov (2011) studied the $L_{1}$-penalized quantile regression in high dimensional sparse models. Fan, Fan, and Barut (2014) further considered an adaptively weighted $L_{1}$ penalty to alleviate the bias problem and showed the oracle property and asymptotic normality of the corresponding estimator. Other robust estimators were developed based on Least Absolute Deviation (LAD) regression. Wang (2013) studied the $L_{1}$-penalized LAD regression and showed that the estimator achieves near oracle risk performance under the high dimensional setting.

The above methods essentially estimate the conditional median (or quantile) regression, instead of the conditional mean regression function. In the applications where the mean regression is of interest, these methods are not feasible unless a strong assumption is made that the distribution of errors is symmetric around zero. A simple example is the heteroscedastic linear model with asymmetric noise distribution. Another example is to estimate the conditional variance function such as ARCH model (Engle, 1982). In these cases, the conditional mean and conditional median are very different. Another important example is to estimate large covariance matrix without assuming light-tails. We will explain this more in details in Section 5. In addition, LAD-based methods tend to penalize strongly on small errors. If only a small proportion of samples are outliers, they are expected to be less efficient than the least squares based method.

A natural question is then how to conduct ultrahigh dimensional mean regression when the tails of errors are not light? How to estimate the sample mean with very fast concentration when the distribution has only bounded second moment? These simple questions have not been carefully studied. LAD-based methods do not intend to answer these questions as they alter the problems of the study. This leads us to consider Huber loss as another way of robustification. The Huber loss (Huber, 1964) is a hybrid of squared loss for relatively small errors and absolute loss for relatively large errors, where the degree of hybridization is controlled by one tuning parameter. Unlike the traditional Huber loss, we allow the regularization parameter to diverge (or converge if its reciprocal is used) in order to reduce the bias induced by the Huber loss for estimating conditional mean regression function. In this paper, we consider the regularized approximate quadratic (RA-Lasso) estimator with an $L_{1}$ penalty and show that it admits the same $L_{2}$ error rate as the optimal error rate in the light-tail situation. In particular, if the distribution of errors is indeed symmetric around 0 (where the median and mean agree), this rate is the same as the regularized LAD estimator obtained in Wang (2013). Therefore, the RA-Lasso estimator does not lose efficiency in this special case. In practice, since the distribution of errors is unknown, RA-Lasso is more flexible than the existing methods in terms of estimating the conditional mean regression function.

A by-product of our method is that the RA-Lasso estimator of the population mean has the exponential type of concentration even in presence of the finite second moment. Catoni (2012) studied this type of problem and proposed a class of losses to result in a robust $M$-estimator of mean with exponential type of concentration. We further extend his idea to the sparse linear regression setting.

As done in many other papers, estimators with nice sampling properties are typically defined through the optimization of a target function such as the penalized least-squares. The properties that are established are not the same as the ones that are computed. Following the framework of Agarwal, Negahban, and Wainwright (2012), we propose the composite gradient descent algorithm for solving the RA-Lasso estimator and develop the sampling properties by taking computational error into consideration. We show that the algorithm indeed produces a solution that admits the same optimal $L_{2}$ error rate as the theoretical estimator after sufficient number of iterations.

This paper is organized as follows. First, in Section 2, we introduce the RA-Lasso estimator and show that it has the same $L_{2}$ error rate as the optimal rate under light-tails. In Section 3, we study the property of the composite gradient descent algorithm for solving our problem and show that the algorithm produces a solution that performs as well as the solution in theory. In Section 4, we show the connection between Huber loss and Catoni loss and establish an concentration inequality for robust estimation of mean. The estimation of the error’s variance is investigated in Section 5. Numerical studies are given in Section 6 and 7 to compare our method with two competitors. All technical proofs are presented in Section 8.

RA-Lasso estimator

where $\{\mathbf{x}_{i}\}_{i=1}^{n}$ are independent and identically distributed (i.i.d) $p$-dimensional covariate vectors, $\{\epsilon_{i}\}_{i=1}^{n}$ are i.i.d errors, and $\boldsymbol{\beta}^{\ast}$ is a $p$-dimensional regression coefficient vector. We consider the high-dimensional setting, where $\log(p)=o(n^{b})$ for some constant $0<b<1$. We assume the distributions of $\mathbf{x}$ and $\epsilon$ are independent and both have mean 0. Under this assumption, $\boldsymbol{\beta}^{\ast}$ is the mean effect of $y$ conditioning on $\mathbf{x}$, which is assumed to be of interest.

To adapt for different magnitude of errors and robustify the estimation, we propose to use the Huber loss (Huber, 1964):

To estimate $\boldsymbol{\beta}^{\ast}$, we propose to solve the following convex optimization problem:

In the following, we give the rate of the approximation and estimation error, respectively. We show that $\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\ast}\rVert_{2}$ admits the same rate as the optimal rate under light tails, as long as the two tuning parameters $\alpha$ and $\lambda_{n}$ are properly chosen. We first give the rate of the approximation error under some moment conditions on $\mathbf{x}$ and $\epsilon$. We assume both $\boldsymbol{\beta}^{\ast}$ and $\boldsymbol{\beta}_{\alpha}^{\ast}$ are interior points of an $L_{2}$ ball with sufficiently large radius.

It holds that $\lVert\boldsymbol{\beta}_{\alpha}^{\ast}-\boldsymbol{\beta}^{\ast}\rVert_{2}=O(\alpha^{k-1})$, under the following conditions:

$\operatorname{E}|\epsilon|^{k}\leq M_{k}<\infty$, for some $k\geq 2$.

$0<\kappa_{l}\leq\lambda_{\min}(\operatorname{E}[\mathbf{x}\mathbf{x}^{T}])\leq\lambda_{\max}(\operatorname{E}[\mathbf{x}\mathbf{x}^{T}])\leq\kappa_{u}<\infty$,

Theorem 1 reveals that the approximation error vanishes faster if higher moments of error distribution exists. We next give the rate of the estimation error $\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{\alpha}^{\ast}\rVert_{2}$. This part differs from the existing work regarding the estimation error of high dimensional regularized $M$-estimator (Negahban, et al., 2012; Agarwal, Negahban, and Wainwright, 2012) as the population minimizer $\boldsymbol{\beta}_{\alpha}^{\ast}$ now varies with $\alpha$. However, we will show that the estimation error rate does not depend on $\alpha$, given a uniform sparsity condition.

In order to be solvable in the high-dimensional setting, $\boldsymbol{\beta}^{\ast}$ is usually assumed to be sparse or weakly sparse, i.e. many elements of $\boldsymbol{\beta}^{\ast}$ are zero or small. By Theorem 1, $\boldsymbol{\beta}_{\alpha}^{\ast}$ converges to $\boldsymbol{\beta}^{\ast}$ as $\alpha$ goes to 0. In view of this fact, we assume that $\boldsymbol{\beta}_{\alpha}^{\ast}$ is uniformly weakly sparse when $\alpha$ is sufficiently small. In particular, we assume that there exists a small constant $r>0$, such that $\boldsymbol{\beta}_{\alpha}^{\ast}$ belongs to an $L_{q}$-ball with a uniform radius $R_{q}$ that

for all $\alpha\in(0,r]$, and some $q\in$. When the conditional distribution of $\varepsilon_{i}$ is symmetric, $\beta_{\alpha,j}^{\ast}=\beta_{j}^{*}$ for all $\alpha$ and $j$. Therefore the condition reduces to that $\boldsymbol{\beta}^{\ast}$ is in the $L_{q}$ ball. In a special case where $q=1$, it follows from Theorem 1 that if $\boldsymbol{\beta}^{\ast}$ belongs to the $L_{1}$-ball with radius $R_{1}/2$ and $r\leq[R_{1}/(2c_{0}\sqrt{p})]^{\frac{1}{k-1}}$, where $c_{0}$ is a generic constant, then $\boldsymbol{\beta}_{\alpha}^{\ast}$ belongs to the $L_{1}$-ball with radius $R_{1}$ for all $\alpha\in(0,r]$. For a general $q\in[0,1)$, we assume a uniform upper bound $R_{q}$ as in (2.4), which is allowed to diverge to infinity.

Since the RA-quadratic loss is convex, we show that with high probability the estimation error $\widehat{\boldsymbol{\Delta}}=\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{\alpha}^{\ast}$ belongs to a star-shaped set, which depends on $\alpha$ and the threshold level $\eta$ of signals.

Under Conditions (C1) and (C3), by choosing $\lambda_{n}=\kappa_{\lambda}\sqrt{\frac{\log p}{n}}$ and $\alpha\geq\frac{L\lambda_{n}}{4v}$, where $\kappa_{\lambda}$, $v$ and $L$ are some constants, with probability greater than $1-2p^{-c_{0}}$,

where $c_{0}=\kappa_{\lambda}^{2}/(32v)-1$, $\eta$ is a positive constant, $S_{\alpha\eta}=\{j:|\beta_{\alpha,j}^{\ast}|>\eta\}$ and $\boldsymbol{\Delta}_{S_{\alpha\eta}}$ denotes the subvector of $\boldsymbol{\Delta}$ with indices in set $S_{\alpha\eta}$.

We further verify a restricted strong convexity (RSC) condition, which has been shown to be critical in the study of high dimensional regularized $M$-estimator (Negahban, et al., 2012; Agarwal, Negahban, and Wainwright, 2012). Let

The loss function ${\cal L}_{n}$ satisfies RSC condition on a set $S$ with curvature $\kappa_{{\cal L}}>0$ and tolerance $\tau_{{\cal L}}$ if

Under conditions (C1)-(C3), for all $\lVert\boldsymbol{\Delta}\rVert_{2}\leq 1$, there exist uniform positive constants $\kappa_{1}$ and $\kappa_{2}$, such that

with probability at least $1-c_{1}\exp(-c_{2}n)$ for some positive constants $c_{1}$ and $c_{2}$.

Suppose conditions (C1)-(C3) hold and assume that

Under conditions of Lemma 1 and (2.7), with probability at least $1-2p^{-c_{0}}-c_{1}\exp(-c_{2}n)$,

Finally, Theorem 1 and 2 together lead to the following main result, which gives the rate of the statistical error $\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\ast}\rVert_{2}$.

Under conditions of Lemma 1 and (2.7), with probability at least $1-2p^{-c_{0}}-c_{1}\exp(-c_{2}n)$,

Next, we compare our result with the existing results regarding the robust estimation of high dimensional linear regression model.

When the distribution of $\epsilon$ is symmetric around 0, then $\boldsymbol{\beta}_{\alpha}^{\ast}=\boldsymbol{\beta}^{\ast}$ for any $\alpha$, which has no approximation error. If $\epsilon$ has heavy tails in addition to being symmetric, we would like to choose $\alpha$ sufficiently large to robustify the estimation. It then follows from Theorem 2 that $\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\ast}\rVert_{2}=O_{P}(\sqrt{R_{q}}[(\log p)/n]^{1/2-q/4})$, where $R_{q}=\sum_{j=1}^{p}|\beta_{j}^{\ast}|^{q}$. The rate is the same as the minimax rate (Raskutti, Wainwright, and Yu, 2011) for weakly sparse model under the light tails. In a special case that $q=0$, it gives $\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\ast}\rVert_{2}=O_{P}(\sqrt{s(\log p)/n})$, where $s$ is the number of nonzero elements in $\boldsymbol{\beta}^{\ast}$. This is the same rate as the regularized LAD estimator in Wang (2013) and the regularized quantile estimator in Belloni and Chernozhukov (2011). It suggests that our method does not lose efficiency for symmetric heavy-tailed errors.

If the distribution of $\epsilon$ is asymmetric around 0, the quantile and LAD based methods are inconsistent, since they estimate the median instead of the mean. Theorem 3 shows that our estimator still achieves the optimal rate given that $\alpha=o(\{R_{q}[(\log p)/n]^{1-\frac{q}{2}}\}^{\frac{1}{2(k-1)}})$ even though the $k$-th moment of $\varepsilon$ is assumed. Recall from conditions in Lemma 1 that we also need to choose $\alpha>c\sqrt{(\log p)/n}$ for some constant $c$. Given the sparsity condition (2.7), $\alpha$ can be chosen to meet the above two requirements. In terms of estimating the conditional mean effect, errors with heavy but asymmetric tails give the case where the RA-Lasso has the biggest advantage over the other estimators.

In practice, the distribution of errors is unknown. However, we proved that our method is no worse than the existing methods for any type of errors, as long as the tuning parameters are chosen properly. Hence, our method is more flexible.

Geometric convergence of computational error

The gradient descent algorithm (Nesterov, 2007; Agarwal, Negahban, and Wainwright, 2012) is usually applied to solve the convex problem (2.3). For example, we can replace the RA-quadratic loss with its local quadratic approximation (LQA) and iteratively solve the following optimization problem:

where $\gamma_{u}$ is a fixed constant at each iteration, and the side constraint “$\lVert\boldsymbol{\beta}\rVert_{1}\leq\rho$” is introduced to guarantee good performance in the first few iterations and $\rho$ is allowed to be sufficiently large. To solve (3.1), the update can be computed by a two-step procedure. We first solve (3.1) without the norm constraint by soft-thresholding the vector $\widehat{\boldsymbol{\beta}}^{t}-\frac{1}{\gamma_{u}}\nabla{\cal L}_{n}(\widehat{\boldsymbol{\beta}}^{t})$ at level $\lambda_{n}$ and call the solution $\check{\boldsymbol{\beta}}$. If $\lVert\check{\boldsymbol{\beta}}\rVert_{1}\leq\rho$, set $\widehat{\boldsymbol{\beta}}^{t+1}=\check{\boldsymbol{\beta}}$. Otherwise, $\widehat{\boldsymbol{\beta}}^{t+1}$ is obtained by further project $\check{\boldsymbol{\beta}}$ onto the $L_{1}$-ball $\{\boldsymbol{\beta}:\lVert\boldsymbol{\beta}\rVert_{1}\leq\rho\}$. The projection can be done (Duchi, et al., 2008) by soft-thresholding $\check{\boldsymbol{\beta}}$ at level $\pi_{n}$, where $\pi_{n}$ is given by the following procedure: (1) sort $\{|\check{\beta}_{j}|\}_{j=1}^{p}$ into $b_{1}\geq b_{2}\geq\ldots\geq b_{p}$; (2) find $J=\max\{1\leq j\leq p:b_{j}-(\sum_{r=1}^{j}b_{r}-\rho)/j>0\}$ and let $\pi_{n}=(\sum_{r=1}^{J}b_{j}-\rho)/J$.

The key is that the RA-quadratic loss function ${\cal L}_{n}$ satisfies the restricted strong convexity (RSC) condition and the restricted smoothness condition (RSM) with some uniform constants, namely $\delta{\cal L}_{n}(\boldsymbol{\Delta},\boldsymbol{\beta})$ as defined in (2.5) satisfies the following conditions:

for all $\boldsymbol{\beta}$ and $\boldsymbol{\Delta}$ in some set of interest, with parameters $\gamma_{l}$, $\tau_{l}$, $\gamma_{u}$ and $\tau_{u}$ that do not depend on $\alpha$. We show that such conditions hold with high probability.

We further show in Theorem 4 that, whenever $R_{q}(\frac{\log p}{n})^{1-(q/2)}=o(1)$, which is required for consistency of any method over the weak sparse $L_{q}$ ball by the known minimax results (Raskutti, Wainwright, and Yu, 2011), it holds that $\lVert\widehat{\boldsymbol{\beta}}^{t}-\widehat{\boldsymbol{\beta}}\rVert_{2}=o(\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{\alpha}^{\ast}\rVert_{2})$ for sufficiently many iterations with lower bound specified in Theorem 4. Hence,

which has the same rate as $\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\ast}\rVert_{2}$. Hence, from a statistical point of view, there is no need to iterate beyond $t$ steps.

Under conditions of Theorem 3, suppose we choose $\lambda_{n}$ as in Lemma 1 and also satisfying

where $|S_{\alpha\eta}|$ denotes the cardinality of set $S_{\alpha\eta}$ and $\bar{\gamma}_{l}=\gamma_{l}-64\tau_{l}|S_{\alpha\eta}|$, then with probability at least $1-p^{-c_{0}}-c_{1}\exp(-c_{2}n)$, we have

where $\phi_{n}(\boldsymbol{\beta})={\cal L}_{n}(\boldsymbol{\beta})+\lambda_{n}\lVert\boldsymbol{\beta}\rVert_{1}$ and $\widehat{\boldsymbol{\beta}}^{0}$ is the initial value, $\delta=\varepsilon^{2}/(1-\kappa)$ is the tolerance level, $\kappa$ and $\varepsilon$ are some constants as will be defined in (8.21) and (8.22), respectively.

Connection with Catoni loss

Catoni (2012) considered the estimation of the mean of heavy-tailed distribution with fast concentration. He proposed an $M$-estimator by solving

where the influence function $\psi_{c}(x)$ is chosen such that $-\log(1-x+x^{2}/2)\leq\psi_{c}(x)\leq\log(1+x+x^{2}/2)$. He showed that this $M$-estimator has the exponential type of concentration by only requiring the existence of the variance. It performed as well as the sample mean under the light-tail case. In Section 2, we essentially showed the same type of concentration for the RA-quadratic loss under the linear regression setting.

The estimation of mean can be regarded as a univariate linear regression where the covariate equals to 1. In that special case, we have a more explicit concentration result for the RA-mean estimator, which is the estimator that minimizes the RA-quadratic loss. Let $\{y_{i}\}_{i=1}^{n}$ be an i.i.d sample from some unknown distribution with $\operatorname{E}(y_{i})=\mu$ and $\mbox{var}(y_{i})=\sigma^{2}$. The RA-mean estimator $\widehat{\mu}_{\alpha}$ of $\mu$ is the solution of

for parameter $\alpha\to 0$, where the influence function $\psi(x)=x$ if $|x|\leq 1$, $\psi(x)=1$, if $x>1$ and $\psi(x)=-1$ if $x<-1$. The following theorem gives the exponential type of concentration of $\widehat{\mu}_{\alpha}$ around $\mu$.

Assume $\frac{\log(1/\delta)}{n}\leq 1/8$ and let $\alpha=\sqrt{\frac{\log(1/\delta)}{nv^{2}}}$ where $v\geq\sigma$. Then,

The above result provides fast concentration of the mean estimation with only two moments assumption. This is very useful for last scale hypothesis testing (Efron, 2010; Fan, Han, and Gu, 2012) and covariance matrix estimation (Bickel and Levina, 2008; Fan, Liao and Mincheva, 2013), where uniform convergence is required. Taking the estimation of large covariance matrix as an example, in order for the elements of the sample covariance matrix to converge uniformly, the aforementioned authors require the underlying multivariate distribution be sub-Gaussian. This restrictive assumptions can be removed if we apply the robust estimation with concentration bound. Regarding $\sigma_{ij}=\operatorname{E}X_{i}X_{j}$ as the expected value of the random variable $X_{i}X_{j}$ (it is typically not the same as the median of $X_{i}X_{j}$), it can be estimated with accuracy

where $v\geq\max_{i,j\leq p}\sqrt{\mbox{var}(X_{i}X_{j})}$ and $\widehat{\sigma}_{ij}$ is RA-mean estimator using data $\{X_{ik}X_{jk}\}_{k=1}^{n}$. Since there are only $O(p^{2})$ elements, by taking $\delta=p^{-3}$ and the union bound, we have

when $\operatorname{E}X_{i}^{4}<\infty$. This robustified covariance estimator requires much weaker condition than the sample covariance and has far wide applicability than the sample covariance. It can be regularized further in the same way as the sample covariance matrix.

where the influence function $\psi_{c}(t)$ is given by

Let $\widehat{\boldsymbol{\beta}}^{c}$ be the corresponding solution. Then, $\widehat{\boldsymbol{\beta}}^{c}$ has the same convergence rate as the RA-Lasso, when the second or the third moment of errors exists.

Suppose condition (C1) holds for $k=2$ or 3, (C2), (C3) and (2.7) hold, then with probability at least $1-2p^{-c_{0}}-c_{1}\exp(-c_{2}n)$,

Unlike the RA-lasso, the order of bias of $\widehat{\boldsymbol{\beta}}^{c}$ cannot be further improved, even when higher moments of errors exist beyond the third order. The reason is that the Catoni loss is not exactly the quadratic loss over any finite intervals. Similar results regarding the computational error of $\widehat{\boldsymbol{\beta}}^{c}$ could also be established as in Theorem 4, since the RSC/RSM conditions also hold for Catoni loss with uniform constants.

Variance Estimation

We estimate $\sigma^{2}=E\epsilon^{2}$ based on the RA-Lasso estimator and a cross-validation scheme. To ease the presentation, we assume the data set can be evenly divided into $K$ folds with $m$ observations in each fold. Then, we estimate $\sigma^{2}$ by

where $\widehat{\boldsymbol{\beta}}^{(-k)}$ is the RA-Lasso estimator obtained by using data points outside the $k$-th fold. We show that $\widehat{\sigma}^{2}$ is asymptotically efficient.

Under conditions of Theorem 3, if $R_{q}(\log p)^{1-q/2}/n^{(1-q)/2}\to 0$ for $q\in[0,1)$, and $\alpha=o\left(\{R_{q}[(\log p)/n]^{1-\frac{q}{2}}\}^{\frac{1}{2(k-1)}}\right)$, then

Simulation Studies

In this section, we assess the finite sample performance of the RA-Lasso and compare it with other methods through various models. We simulated data from the following high dimensional model

where we generated $n=100$ observations and the number of parameters was chosen to be $p=400$. We chose the true regression coefficient vector as

where the first 20 elements are all equal to 3 and the rest are all equal to 0. To involve various shapes of error distributions, we considered the following five scenarios:

Normal with mean 0 and variance 4 (N(0,4));

Two times the t-distribution with degrees of freedom 3 ($2t_{3}$);

Mixture of Normal distribution(MixN): $0.5N(-1,4)+0.5N(8,1)$;

Log-normal distribution (LogNormal): $\epsilon=e^{1+1.2Z}$, where $Z$ is standard normal.

Weibull distribution with shape parameter = 0.3 and scale parameter = 0.5 (Weibull).

In order to meet the model assumption, the errors were standardized to have mean 0. Table 1 categorizes the five scenarios according to the shapes and tails of the error distributions.

To obtain our estimator, we iteratively applied the gradient descent algorithm. We compared RA-Lasso with another two methods in high-dimensional setting: (a) Lasso: the penalized least-squares estimator with $L_{1}$-penalty as in Tibshirani (1996); and (b) R-Lasso: the R-Lasso estimator in Fan, Fan, and Barut (2014), which is the same as the regularized LAD estimator with $L_{1}$-penalty as in Wang (2013). Their performance under the five scenarios was evaluated by the following four measurements:

$L_{2}$ error, which is defined as $\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\ast}\rVert_{2}$.

$L_{1}$ error, which is defined as $\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{\ast}\rVert_{1}$.

Number of false positives (FP), which is number of noise covariates that are selected.

Number of false negatives (FN), which is number of signal covariates that are not selected.

We also measured the relative gain of RA-Lasso with respect to R-Lasso and Lasso, in terms of the difference to the oracle estimator. The oracle estimator $\widehat{\boldsymbol{\beta}}_{\text{oracle}}$ is defined to be the least square estimator by using the first 20 covariates only. Then, the relative gain of RA-Lasso with respect to Lasso ($\text{RG}_{\text{A,L}}$) in $L_{2}$ and $L_{1}$ norm are defined as

The relative gain of RA-Lasso with respect to R-Lasso ($\text{RG}_{\text{A,R}}$) is defined similarly.

For RA-Lasso, the tuning parameters $\lambda_{n}$ and $\alpha$ were chosen optimally based on 100 independent validation datasets. We ran a 2-dimensional grid search to find the best $(\lambda_{n},\alpha)$ pair that minimizes the mean $L_{2}$-loss of the 100 validation datasets. Such an optimal pair was then used in the simulations. Similar method was applied in choosing the tuning parameters in Lasso and R-Lasso.

The above simulation model is based on the additive model (6.1), in which error distribution is independent of covariates. However, this homoscedastic model makes the conditional mean and the conditional median differ only by a constant. To further examine the deviations between the mean regression and median regression, we also simulated the data from the heteroscedastic model

where the constant $c=\sqrt{3}\|\boldsymbol{\beta}^{\ast}\|^{2}$ makes $\operatorname{E}[c^{-1}(\mathbf{x}_{i}^{T}\boldsymbol{\beta}^{\ast})^{2}]^{2}=1$. Note that $\mathbf{x}_{i}^{T}\boldsymbol{\beta}^{\ast}\sim N(0,\|\boldsymbol{\beta}^{\ast}\|^{2})$ and therefore $c$ is chosen so that the average noise levels is the same as that of $\epsilon_{i}$. For both the homoscedastic and the heteroscedastic models, we ran 100 simulations for each scenario. The mean of each performance measurement is reported in Table 2 and 3, respectively.

Tables 2 and 3 indicate that our method had the biggest advantage when the errors were asymmetric and heavy-tailed (LogNormal and Weibull). In this case, R-Lasso had larger $L_{1}$ and $L_{2}$ errors due to the bias for estimating the conditional median instead of the mean. Even though Lasso did not have bias in the loss component, it did not perform well due to its sensitivity to outliers. The advantage of our method is more pronounced in the heteroscedastic model than in the homoscedastic model. Both of them clearly indicate that if the errors come from asymmetric and heavy-tailed distributions, our method is better than both Lasso and R-Lasso. When the errors were symmetric and heavy-tailed ($2t_{3}$), our estimator performed closely as R-Lasso, both of which outperformed Lasso. The above two cases evidently showed that RA-Lasso was robust to the outliers and did not lose efficiency when the errors were indeed symmetric. Under the light-tailed scenario, if the errors were asymmetric (MixN), our method performed similarly as Lasso. R-Lasso performed worse, since it had bias. For the regular setting (N(0, 4)), where the errors were light-tailed and symmetric, the three methods were comparable with each other.

In conclusion, RA-Lasso is more flexible than Lasso and R-Lasso. The tuning parameter $\alpha$ automatically adapts to errors with different shapes and tails. It enables RA-Lasso to render consistently satisfactory results under all scenarios.

Real Data Example

In this section, we use a microarray data to illustrate the performance of Lasso, R-Lasso and RA-Lasso. Huang, et al. (2011) studied the role of innate immune system on the development of atherosclerosis by examining gene profiles from peripheral blood of 119 patients. The data were collected using Illumina HumanRef8 V2.0 Bead Chip and are available on Gene Expression Omnibus. The original study showed that the toll-like receptors (TLR) signaling pathway plays an important role on triggering the innate immune system in face of atherosclerosis. Under this pathway, the “TLR8” gene was found to be a key atherosclerosis-associated gene. To further study the relationship between this key gene and the other genes, we regressed it on another 464 genes from 12 different pathways (TLR, CCC, CIR, IFNG, MAPK, RAPO, EXAPO, INAPO, DRS, NOD, EPO, CTR) that are related to the TLR pathway. We applied Lasso, R-Lasso and RA-Lasso to this data. The tuning parameters for all methods were chosen by using five-fold cross validation. Figure 1 shows our choice of the penalization parameter based on the cross validation results. For RA-Lasso, the choice of $\alpha$ was insensitive to the results and was fixed at 5. We then applied the three methods with the above choice of tuning parameters to select significant genes. The QQ-plots of the residuals from the three methods are shown in Figure 2. The selected genes by the three methods are reported in Table 4. After the selection, we regressed the expression of TLR8 gene on the selected genes, the $t$-values from the refittings are also reported in Table 4.

Table 4 shows that Lasso only selected one gene. R-Lasso selected 17 genes. Our proposed RA-Lasso selected 34 genes. Eight genes (CSF3, IL10, AKT1, TOLLIP, TLR1, SHC1, EPOR, and TJP1) found by R-Lasso were also selected by RA-Lasso. Compared with Lasso and R-Lasso, our method selected more genes, which could be useful for a second-stage confirmatory study. It is clearly seen from Figure 2 that the residuals from the fitted regressions had heavy right tail and skewed distribution. We know from the simulation studies in Section 6 that RA-Lasso tends to perform better than Lasso and R-Lasso in this situation. For further investigation, we randomly chose 24 patients as the test set; applied three methods to the rest patients to obtain the estimated coefficients, which in return were used to predict the responses of 24 patients. We repeated the random splitting 100 times, the boxplots of the Mean Absolute/Squared Error of predictions are shown in Figure 3. RA-Lasso has better predictions than Lasso and R-Lasso.

Proofs

This together with (8.1) completes the proof. ∎

where $\psi(x)=x$, for $|x|\leq 1$; $\psi(x)=1$, for $x>1$; and $\psi(x)=-1$, for $x<-1$. Using $\alpha^{-1}|\psi(\alpha x)|\leq|x|$ and assumption (C3), we have

where $v$ is a constant depending on $\kappa_{0}$ and $M_{2}$ and the last inequality follows from a similar argument as in the proof of Theorem 1. By (C3) and that $|\psi(x)|\leq 1$, $\psi[\alpha(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}_{\alpha}^{\ast})]x_{ij}$ is also sub-Gaussian. For any $k\geq 3$, using the relation between the $k$th moment and the second moment of sub-Gaussian random variables (Rivasplata, 2012),

where $L$ is a constant depending on $\kappa_{0}$ only. Hence,

By Bernstein inequality (Proposition 2.9 of Massart and Picard (2007)) and note that $\operatorname{E}(\frac{2}{\alpha}\psi[\alpha(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}_{\alpha}^{\ast})]\mathbf{x}_{i})=\boldsymbol{0}$, we have

Let $t=n\lambda_{n}^{2}/(32v)$ and observe that $\frac{2Lt}{\alpha n}\leq\sqrt{\frac{2vt}{n}}$ by the choice of $\lambda_{n}$ and $\alpha$. We have

It then follows from union inequality that

where $c_{0}=\kappa_{\lambda}^{2}/(32v)-1$. This completes the proof. ∎

where $\psi^{\prime}(x)=1$ for $|x|\leq 1$, and $\psi^{\prime}(x)=0$ otherwise. Note that each term in (8.5) is nonnegative. However, the quadratic component in (8.5) is not Lipschitz continuous with a bounded Lipschitz coefficient. In order to apply the contraction theorem of Ledoux and Talagrand (1991), we introduce a truncation function that is Lipschitz and bound (8.5) from below. Let

where $I(\cdot)$ is the indicator function. Clearly, $\varphi_{t}(u)\leq u^{2}$ and satisfies the Lipschitz condition with Lipschitz coefficient bounded by $2t$. We first show

for $0<\alpha\leq 1/(T+\tau)$, where the thresholds $T$ and $\tau$ will be chosen as in (8.10).

Let $A_{i}=\psi^{\prime}[\alpha(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}_{\alpha}^{\ast}+v\mathbf{x}_{i}^{T}\boldsymbol{\Delta})](\mathbf{x}_{i}^{T}\boldsymbol{\Delta})^{2}$. We need only to show that

When $|y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}_{\alpha}^{\ast}|>T$ or $|\mathbf{x}_{i}^{T}\boldsymbol{\Delta}|>\tau$, the right hand side is zero and the inequality holds trivially. Thus, we need only to consider the case $|y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}_{\alpha}^{\ast}|\leq T$ and $|\mathbf{x}_{i}^{T}\boldsymbol{\Delta}|\leq\tau$. In this case,

and hence $\psi^{\prime}[\alpha(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}_{\alpha}^{\ast}+v\mathbf{x}_{i}^{T}\boldsymbol{\Delta})]=1$. Using this,

Using (8.7), to prove the Lemma, we need to show that, for any fixed $\delta\in(0,1]$, with high probability

where constants $\kappa_{1}$ and $\kappa_{2}$ do not depend on $\alpha$ and

which equals to the event (8.8) with $\tau=\delta\tau_{1}$.

To establish (8.9), let us consider its complementary event. Define

Our goal is to show that the probability of this event is very small, which is demonstrated through the following three steps.

First, we show that with the following choice of truncation

Finally, we use a standard peeling argument (Alexander, 1987; Van de Geer, 2000) to establish

The result (c) together with (8.11) show that the probability of the complementary event of (8.9) with $\kappa_{1}=\kappa_{l}/4$ and $\kappa_{2}=40\tau^{2}\kappa_{0}\kappa_{1}^{-1}$ is bounded by $\exp(-c_{1}n-c_{2}\log p)$, which completes the proof.

Note that, $g_{\boldsymbol{\Delta}}(\mathbf{x})=(\mathbf{x}^{T}\boldsymbol{\Delta})^{2}$ for all $\mathbf{x}$ such that $|y-\mathbf{x}^{T}\boldsymbol{\beta}_{\alpha}^{\ast}|\leq T$ and $|\mathbf{x}^{T}\boldsymbol{\Delta}|\leq\tau/2$. Therefore, we have

To bound the first term on the right hand side of (8.13), it follows from the Cauchy-Schwartz inequality that

Since $\mathbf{x}^{T}\boldsymbol{\Delta}$ is sub-Gaussian with parameter at most $\kappa_{0}^{2}$ by assumption (C3), we have $\operatorname{E}(\mathbf{x}^{T}\boldsymbol{\Delta})^{4}\leq 16\kappa_{0}^{4}$. Meanwhile, for any $0<\alpha\leq 1/(T+\tau)$, it follows from the Chebyshev inequality and Theorem 1 that

To bound the second term on the right hand side of (8.13), by the concentration inequality of sub-Gaussian variables, we have

Then, by the choice of $T$ and $\tau$ in (8.10),

Next, we bound $\operatorname{E}Z(t)$. Let $\{\omega_{i}\}_{i=1}^{n}$ be an i.i.d. sequence of Rademacher variables. A symmetrization theorem (Theorem 14.3 of Bühlmann and Van De Geer (2011)) yields

By definition, the function $\varphi_{\tau}$ is Lipschitz with parameter at most $2\tau\leq 2\tau^{2}$ and $\varphi_{\tau}(0)=0$. Therefore, by the Ledoux-Talagrand contraction theorem (Ledoux and Talagrand (1991), p.112), we have

Since the variables $\{x_{ij}\}_{i=1}^{n}$ are zero-mean i.i.d. sub-Gaussian with parameter at most $\kappa_{0}^{2}$, so are $\{\omega_{i}x_{ij}I(|y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}_{\alpha}^{\ast}|\leq T)\}_{i=1}^{n}$. Since $\operatorname{E}\left\lVert\frac{1}{n}\sum_{i=1}^{n}\omega_{i}\mathbf{x}_{i}I(|y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}_{\alpha}^{\ast}|\leq T)\right\rVert_{\infty}$ is the maxima of $p$ such terms, known bounds on the expectation of sub-Gaussian maxima (e.g. see Ledoux and Talagrand (1991), p.79) yield

where constants $c_{1}^{\prime}$ and $c_{2}^{\prime}$ depends on $\kappa_{l}$ and $\kappa_{0}$ only. This result holds for each given $t$.

since $Z(\lVert\boldsymbol{\Delta}\rVert_{1})\geq 2^{m-3}\kappa_{l}$ for $\boldsymbol{\Delta}\in B_{m}$. By letting $2^{m-3}\kappa_{l}={\kappa_{l}}/{4}+40\tau^{2}\kappa_{0}t\sqrt{(\log p)/n}$ as in (8.12) and solving for $t$, by (8.12), we obtain

where the last inequality follows from sum of geometric series. ∎

Therefore, $|S_{\alpha\eta}|\leq\eta^{-q}R_{q}$. Let $S_{\alpha\eta}^{c}=\{1,2,\ldots,p\}\backslash S_{\alpha\eta}$, we have

By the Cauchy-Schwartz inequality and (8.17), we can bound further that

With $\lambda_{n}=\kappa_{\lambda}\sqrt{(\log p)/n}$ and $\eta=\lambda_{n}$, it holds that

which is no larger than $\kappa_{1}/2$ under assumption (2.7). On the other hand,

Therefore, RSC holds with $\kappa_{{\cal L}}=\frac{\kappa_{1}}{2}$ and $\tau^{2}_{{\cal L}}=4R_{q}\kappa_{2}{\kappa_{\lambda}}^{1-q}(\frac{\log p}{n})^{1-(q/2)}$. ∎

Let $A_{1}$ and $A_{2}$ denote the events that Lemma 1 and Lemma 3 hold, respectively. By Theorem 1 of Negahban, et al. (2012), within $A_{1}\cap A_{2}$, it holds that

where (i) follows from the choice of $\eta=\lambda_{n}$. On the other hand, by Lemma 1 and 3, $P(A_{1}\cap A_{2})\geq 1-2p^{-c_{0}}-c_{1}\exp(-c_{2}n)$. ∎

From the proof of Lemma 2, we can see that (2.6) indeed holds for all $\boldsymbol{\beta}$ and $\boldsymbol{\Delta}\in\{\boldsymbol{\Delta}:\lVert\boldsymbol{\Delta}\rVert_{2}\leq 1\}$ that

Using the fact that $ab\leq(a^{2}+b^{2})/2$, we conclude that

Therefore, (3.2) holds with $\gamma_{l}=\kappa_{1}$ and $\tau_{l}=\kappa_{1}\kappa_{2}^{2}(\log p)/(2n)$. Meanwhile, since $|\psi^{\prime}(\cdot)|\leq 1$, it follows from (8.5) that

Under the sub-Gaussianity assumption (C3), it follows from some existing work (e.g. page 18 of Loh and Wainwright (2013)) that, with probability great than $1-c_{1}\exp(-c_{2}n)$, it holds that

where $c_{1}$ and $c_{2}$ are some generic constants. Hence, (3.3) holds with $\gamma_{u}={3\kappa_{u}}$ and $\tau_{u}=\kappa_{u}(\log p)/n$. ∎

We prove the theorem by the following two steps:

(a) We first show that, for any $\delta^{2}\geq\varepsilon^{2}/(1-\kappa)$, $\phi(\widehat{\boldsymbol{\beta}}^{t})-\phi(\widehat{\boldsymbol{\beta}})\leq\delta^{2}$, for all $t$ greater than the right hand side of (8.20), where $\kappa\in[0,1)$ is a contraction constant and $\varepsilon$ is a tolerance parameter, which will be given in (8.21) and (8.22), respectively.

(b) We use RSC condition (3.2) to transform the upper bound of $\phi(\widehat{\boldsymbol{\beta}}^{t})-\phi(\widehat{\boldsymbol{\beta}})$ into the upper bound of $\lVert\widehat{\boldsymbol{\beta}}^{t}-\widehat{\boldsymbol{\beta}}\rVert_{2}$.

For step (a), since our loss function is convex, we apply Theorem 2 of Agarwal, Negahban, and Wainwright (2012). In order for our proof to be self-contained, we cite their theorem as the follows:

[Theorem 2 of Agarwal, Negahban, and Wainwright (2012)] Suppose for any data set $Z_{1}^{n}$, the loss function ${\cal L}_{n}(.;Z_{1}^{n})$ is convex and differentiable and the regularizer ${\cal R}$ is a norm. Consider the optimization problem of $\widehat{\theta}=\operatorname{argmin}_{{\cal R}(\theta)\leq\rho}\{{\cal L}_{n}(\theta;Z_{1}^{n})+\lambda_{n}{\cal R}(\theta)\}$ for a radius $\rho$ such that $\theta^{\ast}$ is feasible, where $\theta^{\ast}=\operatorname{argmin}\operatorname{E}{\cal L}_{n}(\theta;Z_{1}^{n})$, and a regularization parameter $\lambda_{n}$ satisfying bound

where ${\cal R}^{\ast}$ is the dual norm of the regularizer. In addition, suppose that the loss function ${\cal L}_{n}$ satisfies the RSC/RSM condition with parameters ($\gamma_{l}$, $\tau_{l}$) and ($\gamma_{u}$, $\tau_{u}$), respectively. Let $(\mathcal{M},\bar{\mathcal{M}}^{\perp})$ be any ${\cal R}$-decomposable pair of subspaces such that

where $\Psi(\bar{\mathcal{M}})=\sup_{\theta\in\bar{\mathcal{M}}\backslash\{0\}}{{\cal R}(\theta)}/{\lVert\theta\rVert_{2}}$, $\bar{\gamma}_{l}=\gamma_{l}-64\tau_{l}\Psi^{2}(\bar{\mathcal{M}})$, $\xi=(1-{64\tau_{u}\bar{\gamma}_{l}^{-1}\Psi^{2}(\bar{\mathcal{M}})})^{-1}$, and $\chi=2\left({\bar{\gamma}_{l}}/(4\gamma_{u})+128\tau_{u}\bar{\gamma}_{l}^{-1}\Psi^{2}(\bar{\mathcal{M}})\right)\tau_{l}+8\tau_{u}+2\tau_{l}$. Denote $\varepsilon^{2}=8\xi\chi\left(6\Psi(\bar{\mathcal{M}})\lVert\widehat{\theta}-\theta^{\ast}\rVert_{2}+8{\cal R}(\Pi_{\mathcal{M}^{\perp}}(\theta^{\ast}))\right)^{2}$, where $\Pi_{\mathcal{M}^{\perp}}(\theta^{\ast})$ is the projection of $\theta^{\ast}$ onto $\mathcal{M}^{\perp}$. Then for any $\delta^{2}\geq\varepsilon^{2}/(1-\kappa)$, we have $\phi_{n}(\widehat{\theta}^{t})-\phi_{n}(\widehat{\theta})\leq\delta^{2}$ for all

where $\phi_{n}(\theta)={\cal L}_{n}(\theta;Z_{1}^{n})+\lambda_{n}{\cal R}(\theta)$, $\widehat{\theta}^{t}$ is the solution by the gradient descent algorithm after $t^{\text{th}}$ iteration, and $\theta^{0}$ is the initial value of $\theta$.

In fact, Theorem 2 of Agarwal, Negahban, and Wainwright (2012) is a deterministic statement for all choices of pairs $(\mathcal{M},\bar{\mathcal{M}}^{\perp})$. From Lemma 1 and Lemma 4, we have shown that with our choice of $\lambda_{n}$, the RA-quadratic loss function satisfy (8.18) and RSC/RSM with probability at least $1-2p^{-c_{0}}$ and $1-c_{1}\exp(-c_{2}n)$, respectively. Hence, Theorem 2 of Agarwal, Negahban, and Wainwright (2012) applies to our problem with high probability. We further choose the pair $(\mathcal{M},\bar{\mathcal{M}}^{\perp})=(S_{\alpha\eta},S_{\alpha\eta}^{c})$ and give the explicit expression of constants for our problem as the follows:

where $\bar{\gamma}_{l}=\kappa_{1}-32\kappa_{1}\kappa_{2}^{2}|S_{\alpha\eta}|{(\log p)}/{n}$, $\xi=\{1-64\kappa_{u}|S_{\alpha\eta}|(\log p)/(n\bar{\gamma}_{l})\}^{-1}$, and $\chi=2\{{\bar{\gamma}_{l}}/(4\gamma_{u})+{128\tau_{u}|S_{\alpha\eta}|}/{\bar{\gamma}_{l}}+1\}\tau_{l}+8\tau_{u}$. It remains to check (8.19). By (8.21), $\kappa\in[0,1)$ is equivalent to requiring

With $\eta=\lambda_{n}$, it follows from (8.16) that

Hence, (8.23) holds when $n$ is sufficiently large. Moreover, from (8.19) we need

which is satisfied under the stated assumption. It then follows from Theorem 2 of Agarwal, Negahban, and Wainwright (2012) that, for any $\delta^{2}\geq{\varepsilon^{2}}/(1-\kappa)$, $\phi(\widehat{\boldsymbol{\beta}}^{t})-\phi(\widehat{\boldsymbol{\beta}})\leq\delta^{2}$, for all iterations $t$ greater than the right hand side of (8.20).

For step (b), it follows from RSC condition that

Since $\widehat{\boldsymbol{\beta}}$ is the minimizer of $\phi(\boldsymbol{\beta})$, by the first-order condition, $[\nabla{\cal L}_{n}(\widehat{\boldsymbol{\beta}})+\lambda_{n}\nabla\lVert\widehat{\boldsymbol{\beta}}\rVert_{1}]^{T}(\widehat{\boldsymbol{\beta}}^{t}-\widehat{\boldsymbol{\beta}})\geq 0$. Therefore,

By the convexity of the $L_{1}$-norm, $\lVert\widehat{\boldsymbol{\beta}}^{t}\rVert_{1}-\lVert\widehat{\boldsymbol{\beta}}\rVert_{1}-[\nabla\lVert\widehat{\boldsymbol{\beta}}\rVert_{1}]^{T}(\widehat{\boldsymbol{\beta}}^{t}-\widehat{\boldsymbol{\beta}})\geq 0$. Hence,

Next, we bound $\lVert\widehat{\boldsymbol{\beta}}^{t}-\widehat{\boldsymbol{\beta}}\rVert_{1}$. It follows from Lemma 3 of Agarwal, Negahban, and Wainwright (2012) that

where $\delta$ is defined as in (a). Then, by the Cauchy-Schwartz inequality,

Equations (8.24) and (8.25) together with results in (a) imply that,

We now bound the second term in (8.26). By (8.16) and (8.17), we have

Since $\bar{\gamma}_{l}\asymp 1$, $\kappa\asymp 1$, $\xi\asymp 1$, $\chi\asymp\frac{\log p}{n}$, and $\tau_{l}\asymp\frac{\log p}{n}$, it follows from (8.26), (8.27) and (8.28) that

The proof follows the same spirit of the proof of Proposition 2.4 of Catoni (2012). The influence function $\psi(x)$ satisfies

Using this and independence, with $r(\theta)=\frac{1}{\alpha n}\sum_{i=1}^{n}\psi[\alpha(Y_{i}-\theta)]$, we have

Similarly, $\operatorname{E}\left\{\exp[-\alpha nr(\theta)]\right\}\leq\exp\left\{-n\alpha(\mu-\theta)+n\alpha^{2}[v^{2}+(\mu-\theta)^{2}]\right\}$. Define

Similarly, $P(r(\theta)<B_{-}(\theta))\leq\delta$.

Let $\theta_{+}$ be the smallest solution of the quadratic equation $B_{+}(\theta_{+})=0$ and $\theta_{-}$ be the largest solution of the equation $B_{-}(\theta_{-})=0$. Under the assumption that $\frac{\log(1/\delta)}{n}\leq 1/8$ and the choice of $\alpha=\sqrt{\frac{\log(1/\delta)}{nv^{2}}}$, we have $\alpha^{2}v^{2}+\frac{\log(1/\delta)}{n}\leq 1/4$. Therefore,

With $\alpha=\sqrt{\frac{\log(1/\delta)}{nv^{2}}}$, $\theta_{+}\leq\mu+4v\sqrt{\frac{\log(1/\delta)}{n}}$, $\theta_{-}\geq\mu-4v\sqrt{\frac{\log(1/\delta)}{n}}$. Since the map $\theta\mapsto r(\theta)$ is non-increasing, under event $\{B_{-}(\theta)\leq r(\theta)\leq B_{+}(\theta)\}$

i.e. $|\widehat{\mu}_{\alpha}-\mu|\leq 4v\sqrt{\frac{\log(1/\delta)}{n}}$. Meanwhile, $P(B_{-}(\theta)\leq r(\theta)\leq B_{+}(\theta))>1-2\delta$. ∎

Hence, by the Cauchy-Schwartz inequality,

Hence, $\lVert\boldsymbol{\beta}_{\alpha}^{c\ast}-\boldsymbol{\beta}^{\ast}\rVert_{2}=O(\alpha^{2})$. If $\epsilon$ only has the second moment exist, by a first-order Taylor expansion of $g_{\alpha}^{\prime}(x)$ similarly as in (8.29), we have $\lVert\boldsymbol{\beta}_{\alpha}^{c\ast}-\boldsymbol{\beta}^{\ast}\rVert_{2}=O(\alpha)$. Next, since $(\psi_{c})^{\prime}(0)=1$, by the same argument as in the proof of Lemma 3, RSC holds for Catoni’s loss with probability no less than $1-c_{1}\exp(-c_{2}n)$. Hence, similarly as in Theorem 2, $\lVert\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{\alpha}^{c\ast}\rVert_{2}=O(\sqrt{R_{q}}[(\log p)/n]^{1/2-q/4})$. This together with $\lVert\boldsymbol{\beta}_{\alpha}^{c\ast}-\boldsymbol{\beta}^{\ast}\rVert_{2}=O(\alpha^{k-1})$ completes the proof. ∎

Given that $\operatorname{E}\epsilon^{4}$ exists, by Central Limit Theorem, $\sqrt{n}(\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}^{2}-\sigma^{2})\stackrel{{\scriptstyle\mathcal{D}}}{{\to}}\mathcal{N}(0,\operatorname{E}\epsilon^{4}-\sigma^{4})$. Let $z_{i}=\mathbf{x}_{i}^{T}(\widehat{\boldsymbol{\beta}}^{(-k)}-\boldsymbol{\beta}^{\ast})$. We now need to prove that the last two terms are negligible. Conditioning on data outside the $k$th fold,

Hence, $m^{-1/2}\sum_{i\in\text{fold }k}\epsilon_{i}\mathbf{x}_{i}^{T}(\widehat{\boldsymbol{\beta}}^{(-k)}-\boldsymbol{\beta}^{\ast})=O_{P}\left(\lVert\widehat{\boldsymbol{\beta}}^{(-k)}-\boldsymbol{\beta}^{\ast}\rVert_{2}\right)=o_{P}(1)$, where the last equality follows from Theorem 3. By an analogous argument, we have

This completes the proof of the Theorem. ∎

References