Inference robust to outliers with l1-norm penalization

TL;DR

This paper introduces a robust linear regression method using l1-norm penalization that effectively handles outliers, achieves standard asymptotic properties, and simplifies inference procedures.

Contribution

It develops a square-root lasso-based estimator for outlier-robust inference with theoretical guarantees and practical computational advantages.

Findings

Estimator achieves asymptotic normality with standard variance.

Method handles growing number of outliers with sample size.

Provides a practical alternative to OLS for robust inference.

Abstract

This paper considers the problem of inference in a linear regression model with outliers where the number of outliers can grow with sample size but their proportion goes to 0. We apply the square-root lasso estimator penalizing the l1-norm of a random vector which is non-zero for outliers. We derive rates of convergence and asymptotic normality. Our estimator has the same asymptotic variance as the OLS estimator in the standard linear model. This enables to build tests and confidence sets in the usual and simple manner. The proposed procedure is also computationally advantageous as it amounts to solving a convex optimization program. Overall, the suggested approach constitutes a practical robust alternative to the ordinary least squares estimator.

Tables

Table 1. Table 1: . bias, variance, mean squared error (MSE) and coverage of 95 % percent 95 95\% confidence intervals for λ = 2.01 ​ 2 ​ log ⁡ ( n ) 𝜆 2.01 2 𝑛 \lambda=2.01\sqrt{2\log(n)} .

value	p	n	${\hat{\beta }}_1$	${\tilde{\beta }}_1^{OLS}$	${\hat{\beta }}_2$	${\tilde{\beta }}_2^{OLS}$
bias	0.025	100	0.127	0.301	0.278	0.671
variance	0.025	100	0.060	0.130	0.097	0.221
MSE	0.025	100	0.076	0.221	0.174	0.671
coverage	0.025	100	0.82	0.47	0.75	0.20
bias	0.01	1000	0.045	0.120	0.133	0.361
variance	0.01	1000	0.002	0.004	0.007	0.018
MSE	0.001	1000	0.004	0.018	0.025	0.148
coverage	0.001	1000	0.74	0.16	0.24	0.00
bias	0.001	10000	0.005	0.015	0.017	$0.057$
variance	0.001	10000	1.08 $\times {10}^{-4}$	1.28 $\times {10}^{-4}$	2.21 $\times {10}^{-4}$	5.23 $\times {10}^{-4}$
MSE	0.001	10000	1.33 $\times {10}^{-4}$	3.53 $\times {10}^{-4}$	5.10 $\times {10}^{-4}$	3.772 $\times {10}^{-3}$
coverage	0.001	10000	0.93	0.66	0.68	0.03

Equations

$$\begin{array}[]{cc}y_{i}=x_{i}^{\top}\beta+\alpha_{i}+\epsilon_{i}&\forall i=1,\dots,n,\end{array}$$

$${\alpha }_i=y_i-{\tilde{y}}_i+({\tilde{x}}_i-x_i)\beta .$$

$$y_i=\beta +{\alpha }_i+{ϵ}_i,$$

$$\hat{\beta }\in \underset{\beta \in {\mathbb{R}}^K}{\mathrm{argmin}}\sum _{i=1}^n(y_i-\beta -{\alpha }_i{)}^2+\lambda \sum _{i=1}^n|{\alpha }_i|,\lambda >0,$$

$$(\hat{\beta },\hat{\alpha })\in \underset{\beta \in {\mathbb{R}}^K,\alpha \in {\mathbb{R}}^n}{\mathrm{argmin}}\frac{1}{\sqrt{n}}||y-X\beta -\alpha |{|}_2+\frac{\lambda }{n}||\alpha |{|}_1,$$

$$\frac{1}{\sqrt{n}}||y-X\beta -\hat{\alpha }|{|}_2+\frac{\lambda }{n}||\hat{\alpha }|{|}_1\le \frac{1}{\sqrt{n}}||y-Xb-\hat{\alpha }|{|}_2+\frac{\lambda }{n}||\hat{\alpha }|{|}_1,$$

$$\hat{\beta }={\left(X^{\top }X\right)}^{-1}{X}^{\top }(y-\hat{\alpha }).$$

$$\frac{1}{\sqrt{n}}||M_X(y-\alpha )|{|}_2+\frac{\lambda }{n}||\alpha |{|}_1\le \frac{1}{\sqrt{n}}||y-Xb-\alpha |{|}_2+\frac{\lambda }{n}||\alpha |{|}_1.$$

$$\hat{\alpha }\in \underset{\alpha \in {\mathbb{R}}^N}{\mathrm{argmin}}\frac{1}{\sqrt{n}}||M_X(y-\alpha )|{|}_2+\frac{\lambda }{n}||\alpha |{|}_1.$$

$$\frac{1}{n}||\hat{\alpha }-\alpha |{|}_1=O_P\left(p\lambda \right).$$

$$\hat{\beta }-\beta ={\left(X^{\top }X\right)}^{-1}{X}^{\top }ϵ+(X^{\top }X{)}^{-1}{X}^{\top }(\alpha -\hat{\alpha }).$$

$${\left|\left|{\left(X^{\top }X\right)}^{-1}{X}^{\top }(\alpha -\hat{\alpha })\right|\right|}_2$$

$$=O_P\left(\frac{1}{n}||X|{|}_{\infty }||\alpha -\hat{\alpha }|{|}_1\right)\text{(By Ho¨lder’s inequality)}.$$

$$||(X^{\top }X{)}^{-1}{X}^{\top }(\alpha -\hat{\alpha })|{|}_2=O_P\left(p\lambda ||X|{|}_{\infty }\right).$$

$$\frac{\hat{\beta }-\beta }{\max \left(\frac{1}{\sqrt{n}},p\lambda ||X|{|}_{\infty }\right)}=O_P\left(1\right).$$

$$\sqrt{n}(\hat{\beta }-\beta )\stackrel{d}{\to }\mathcal{N}(0,{\sigma }^2{\Sigma }^{-1}).$$

$$(\hat{\beta },\hat{\alpha },\hat{\sigma })\in \arg \underset{\beta \in {\mathbb{R}}^K,\alpha \in {\mathbb{R}}^n,\sigma \in {\mathbb{R}}^{+}}{\min }\frac{\sigma }{2}+\frac{1}{2\sigma }{\left|\left|y-X\beta -\alpha \right|\right|}_2^2+\frac{\lambda }{2\sqrt{n}}{\left|\left|\alpha \right|\right|}_1.$$

$$\alpha_{i}=\left\{\begin{array}[]{cc}0&\text{if $x_{2i}<q_{1-p}$}\\ 5x_{2i}&\text{if $x_{2i}\geq q_{1-p}$},\end{array}\right.$$

$${\tilde{\beta }}^{OLS}\in \underset{\beta \in {\mathbb{R}}^K}{\mathrm{argmin}}||y-X\beta -\alpha |{|}_2^2$$

$$||M_Xϵ|{|}_2^2$$

$$=||ϵ|{|}_2^2-{ϵ}^{\top }X(X^{\top }X{)}^{-1}{X}^{\top }ϵ.$$

$$2\sqrt{n}\frac{||M_X(ϵ)|{|}_{\infty }}{||M_Xϵ|{|}_2}$$

$$\le \frac{1}{\sigma }||ϵ|{|}_{\infty }+\left|\frac{\sqrt{n}}{||M_Xϵ|{|}_2}-\frac{1}{\sigma }\right|||ϵ|{|}_{\infty }+\frac{1}{\sigma }||P_Xϵ|{|}_{\infty }+\left|\frac{\sqrt{n}}{||M_Xϵ|{|}_2}-\frac{1}{\sigma }\right|||P_Xϵ|{|}_{\infty }$$

$$2\sqrt{n}\frac{||M_Xϵ|{|}_{\infty }}{||M_Xϵ|{|}_2}=2\frac{||ϵ|{|}_{\infty }}{\sigma }+o_P(||ϵ|{|}_{\infty })+O_P(1)$$

$$\underset{n\to \infty }{\lim }\mathbb{P}\left(2c\sqrt{2\log (n)}\ge 2\frac{||ϵ|{|}_{\infty }}{\sigma }\right)=1.$$

$$\mathbb{P}\left(\frac{\left|{ϵ}_i\right|}{\sigma }\ge t\right)\le 2e^{-\frac{t^2}{2}}.$$

$$\mathbb{P}\left(2c\sqrt{2\log (n)}\ge 2\frac{||ϵ|{|}_{\infty }}{\sigma }\right)$$

$$\le n{e}^{-c\log (n)}\to 0.$$

$$\underset{n\to \infty }{\lim }\mathbb{P}\left(c\sqrt{\log (n)}\ge 2\frac{||ϵ|{|}_{\infty }}{\sigma }\right)=1.$$

$$\mathbb{P}\left(\frac{\left|{ϵ}_i\right|}{\sigma }\ge t\right)\le 2e^{-\frac{t^2}{2b}}.$$

$$\mathbb{P}\left(4\sqrt{b}\sqrt{\log (n)}\ge 2\frac{||ϵ|{|}_{\infty }}{\sigma }\right)$$

$$\le 2n{e}^{-2\log (n)}\to 0.$$

$$\underset{n\to \infty }{\lim }\mathbb{P}\left(c\log (n)\ge 2\frac{||ϵ|{|}_{\infty }}{\sigma }\right)=1.$$

Full text

Inference robust to outliers with $\bm{\ell_{1}}$-norm penalization††thanks: I thank my PhD supervisor Professor Eric Gautier for his availability and great help. I am also grateful to Anne Ruiz-Gazen, Jean-Pierre Florens, Thierry Magnac and Nour Meddahi for useful comments. I acknowledge financial support from the ERC POEMH 337665 grant.

Jad Beyhum 111jad.beyhum@gmail.com

Toulouse School of Economics, Université Toulouse Capitole

This paper considers the problem of inference in a linear regression model with outliers where the number of outliers can grow with sample size but their proportion goes to [math]. We apply the square-root lasso estimator penalizing the $\ell_{1}$-norm of a random vector which is non-zero for outliers. We derive rates of convergence and asymptotic normality. Our estimator has the same asymptotic variance as the OLS estimator in the standard linear model. This enables to build tests and confidence sets in the usual and simple manner. The proposed procedure is also computationally advantageous as it amounts to solving a convex optimization program. Overall, the suggested approach constitutes a practical robust alternative to the ordinary least squares estimator.

KEYWORDS: Machine learning, high-dimensional statistics, square-root lasso, outliers, robust inference.

MSC 2010 Subject Classification: Primary 62F35; secondary 62J05, 62J07.

1 Introduction

This paper considers a linear regression model with outliers. The statistican observes a dataset of $n$ i.i.d. realizations of an outcome scalar random variables $y_{i}$ and a random vector of covariates $x_{i}$ with support in $\mathbb{R}^{K}$, such that $\Sigma=\mathbb{E}[x_{i}x_{i}^{\top}]$ is positive definite. We place ourselves in the Huber’s contamination framework, that is the distribution of $(y_{i},x_{i})$ is a mixture between two distributions. With probability $1/2<1-p\leq 1$, it corresponds to a linear regression model with conditionally homoscedastic errors, that is there exists $\beta\in\mathbb{R}^{K}$ and scalar i.i.d. random variables $\epsilon_{i}$ such that $y_{i}=x_{i}^{\prime}\beta+\epsilon_{i}$, $\mathbb{E}[x_{i}\epsilon_{i}]=\mathbb{E}[\epsilon_{i}]=0$ and $0<\text{var}[\epsilon_{i}^{2}|x_{i}]=\sigma^{2}<\infty$. With probability $p$, the distribution is unspecified. An observation $(y_{i},x_{i})$ is called an outlier when it was generated according to this unspecified distribution $G$. The goal of the statistician is to estimate the parameter $\beta$. This model can be rewritten as

[TABLE]

where $\alpha_{i}$ is scalar random variable which is equal to [math] when an observation is not an outlier and which dependence with $x_{i}$ and $\epsilon_{i}$ is left unrestricted. The probability that $\alpha_{i}$ is different from [math] is hence, $p=\mathbb{P}\left(\alpha_{i}\neq 0\right)=\mathbb{E}\left[\left|\left|\alpha\right|\right|_{0}/n\right]$. We derive estimation results in an asymptotic where $p$ goes to [math] as a function of the sample size $n$.

This model can represent various situations of practical interest. First, the statistician could be interested in $\beta$ because it corresponds to the slope of the best linear predictor of $y_{i}$ given $x_{i}$ for the observations for which $\alpha_{i}=0$. These coefficients are of interest because, in the presence of outliers, the slope of the best linear predictor of $y_{i}$ given $x_{i}$ for the whole population may differ greatly from $\beta$ and hence a statistical analysis based on the whole population may lead to a poor prediction accuracy for the large part of the population that are not outliers.

Second, if $\beta$ is given a causal interpretation, then it represents the causal effect of the regressors for the population of "standard" individuals. That is, for instance, if the aim is evaluate a program, it could be that the treatment effect is negative for most of the population but strongly positive for a small fraction of the individuals, the outliers. The policy maker may not be willing to implement a policy that has a negative effect on most of the population, giving interest to a statistical procedure that estimates the treatment effect of the majority of the population robustly.

Finally, $\beta$ could represent the true coefficient of the best linear predictor of $y_{i}$ given $x_{i}$ in a measurement errors model. Indeed, assume that our population follows the model $\tilde{y}_{i}=\tilde{x}_{i}\beta+\tilde{\epsilon}_{i}$ with $\mathbb{E}[\tilde{x}_{i}\tilde{\epsilon}_{i}]=0$ but that we do not observe $(\tilde{y}_{i},\tilde{x}_{i})$ but $(y_{i},x_{i})$, this fits in our framework with $\epsilon_{i}=\tilde{\epsilon}_{i}$ and

[TABLE]

Hence, $\alpha_{i}$ allows for both measurement errors in $x_{i}$ - called outliers in the $x$-direction - and in $y_{i}$, the outliers in the $y$-direction, for a small fraction of the population, see Rousseeuw and Leroy (2005) for a precise discussion.

This paper develop results on the estimation of $\beta$ when the vector $\alpha=\left(\alpha_{1},\dots,\alpha_{n}\right)^{\top}$ is sparse in the sense that $p$ goes to [math] with $n$. We rely on a variant of the square-root lasso estimator of Belloni et al. (2011a) which penalizes the $\ell_{1}$-norm of the vector $\alpha$. The advantages of our estimator are that the penalty parameter does not depend on the variance of the error term and is computationally tractable. If the vector $\alpha$ is sparse enough, we show that our estimator is $\sqrt{n}$-consistent and asymptotically normal. It has the same asymptotic variance as the OLS estimator in the standard linear model without outliers.

Related literature. This paper is connected to at least two different research fields. First, it draws on the literature on inference in the high-dimensional linear regression model and closely related variants of this model. A series of papers from Belloni et al. (2011b, 2012, 2014a, 2014b, 2016, 2017) study a variety of models ranging from panel data models to quantile regression in an high-dimensional setting. Gautier et al. (2011) proposes inference procedures in an high-dimensional IV model with a large number of both regressors and instrumental variables. Javanmard and Montanari (2014); Van de Geer et al. (2014); Zhang and Zhang (2014) suggest debiasing strategies of the lasso estimator to obtain confidence intervals in a high-dimensional linear regression model. We borrow from this literature by using an $\ell_{1}$-penalized estimator and complete existing research by deriving inference results for the linear regression model with outliers.

Next, our work is related to the literature on robust regression. For detailed accounts of this field, see Rousseeuw and Leroy (2005); Hampel et al. (2011); Maronna et al. (2018). The literature identifies a trade-off between efficiency and robustness, as explicited below. Indeed, $M$-estimators (such as the Ordinary Least-Squares (OLS) estimator) are efficient when data is generated by the standard linear model without outliers and Gaussian errors but this comes at the price of a breakdown point - the maximum proportion of the data that can be contaminated without the estimator performing arbitrarily poorly - of [math]. By contrast, $S$-estimators such as the Least Median of Squares (LMS) and the Least Trimmed Squares (LTS) have a strictly positive and fixed breakdown point. They are also asymptotically normal in the model without outliers but are not efficient and have computational issues because of the non-convexity of their objective functions (see Rousseeuw and Leroy (2005)). Our estimator is efficient under certain conditions, because it attains the same asymptotic variance as the OLS estimator in the standard linear model. Unlike this literature, our procedure relies on a convex program and is computationally tractable, see Belloni et al. (2011a) for a detailed analysis. The proposed approach therefore provides a simple efficient alternative to the rest of the literature.

Within the robust regression literature some authors have considered the application of $\ell_{1}$-norm penalization to robust estimation. In particular, our model nests the Huber’s contamination model for location estimation introduced in Huber et al. (1964). Indeed, if there is a single constant regressor, our model nests the following framework:

[TABLE]

where $\epsilon_{i}\sim\mathcal{N}(0,1)$ i.i.d., $\beta\in\mathbb{R}$ is the mean of $y_{i}$ for non-outlying coefficients while $\mathbb{E}[y_{i}|\alpha_{i}\neq 0]$ is left unrestricted. Chen et al. (2018) show that the minimax lower bound for the squared $\ell_{2}$-norm estimation error is of order greater than $\max(1/n,p^{2})$ under gaussian errors, where $\lvert\lvert\alpha\rvert\rvert_{0}$ is the number of outliers in the sample. When $p\sqrt{\log(n)}\to 0$, we attain this lower bound up to a factor $\log(n)^{2}$. Several strategies have been proposed to tackle this location estimation problem. The one which is the closest to ours is soft-thresholding using a lasso estimator, that is use

[TABLE]

see for instance Collier and Dalalyan (2017). We substitute this estimator with a square-root lasso that has the advantage to provide guidance on the choice of the penalty level that is independent from the variance of the noise (see Belloni et al. (2011a)). We extend the analysis of this type of estimators to the linear regression model and add inference results to the literature. Other $\ell_{1}$-norm penalized estimators for robust linear regression have been studied in the literature such as in Lambert-Lacroix et al. (2011); Dalalyan (2012); Li (2012); Alfons et al. (2013), but the authors do not provide inference results. Fan et al. (2017) considers robust estimation in the case where $\beta$ is a high-dimensional parameter. Its estimator penalizes the Huber loss function by a term proportional to the $\ell_{1}$-norm of $\beta$.

Notations. We use the following notations. For a matrix $M$, $M^{\top}$ is its transpose, $\lvert\lvert M\rvert\rvert_{2}$ is its $\ell_{2}$-norm, $\lvert\lvert M\rvert\rvert_{1}$ is the $\ell_{1}$-norm, $\lvert\lvert M\rvert\rvert_{\infty}$ is its sup-norm, $\lvert\lvert M\rvert\rvert_{\text{op}}$ is its operator norm and $\lvert\lvert M\rvert\rvert_{0}$ is the number of non-zero coefficients in $M$, that is its $\ell_{0}$-norm. For a probabilistic event $\mathcal{E}$, the fact that it happens w.p.a. $1$ (with probability approaching $1$) signifies that $\mathbb{P}\left(\mathcal{E}\right)\xrightarrow[n\to\infty]{}1$. Then, for $k=1,\dots,K$, $x_{k}$ is the vector $((x_{1})_{k},\dots,(x_{n})_{k})^{\top}$ and $X$ is the matrix $(x_{1},\dots,x_{n})^{\top}$. $P_{X}$ is the projector on the vector space spanned by the columns of the matrix $X$ and $M_{X}=I_{n}-P_{X}$, where $I_{n}$ is the identity matrix of size $n$. We introduce $y=(y_{1},\dots,y_{n})^{\top}$ and $\epsilon=(\epsilon_{1},\dots,\epsilon_{n})^{\top}$.

2 Low-dimensional linear regression

2.1 Framework

The probabilistic framework consists of a sequence of data generating processes (henceforth, DGPs) that depend on the sample size $n$. The joint distribution of $(x_{i},\epsilon_{i})$ is independent from the sample size. We consider an asymptotic where $n$ goes to $\infty$ and where $p$, the contamination level, depends on $n$ while the number of regressors remains fixed.

Our estimation strategy is able to handle models where $\alpha$ is sparse, that is $\left|\left|\alpha\right|\right|_{0}/n=o_{P}(1)$ or, in other words, $p\to 0$. Potentially, every individual’s $y_{i}$ can be generated by a distribution that does not follow a linear model but the difference between the distribution of $y_{i}$ and the one yielded by a linear model can only be important for a negligible proportion of individuals. Our subsequent theorems will help to quantify these previous statements.

2.2 Estimation procedure

We consider an estimation procedure that estimates both the coefficients $\alpha_{i}$ and the effects of the regressors $\beta$ by a square-root lasso that penalizes only the coefficients $\alpha_{i}$, that is

[TABLE]

where $\lambda$ is a penalty level whose choice is discussed later. The advantage of the square-root lasso over the lasso estimator is that the penalty level does not depend on an estimate of the variance of $\epsilon_{i}$. Hence, our procedure is simple in that it does not make use of any tuning parameter unlike the least median of quares and least trimmed squares estimators. An important remark is that if $\beta$ is such that $X\beta=P_{X}(y-\widehat{\alpha})$, then

[TABLE]

for any $b\in\mathbb{R}^{K}$. Therefore, if $X^{\top}X$ is positive definite, $\widehat{\beta}$ is the OLS estimator of the regression of $y-\widehat{\alpha}$ on $X$, that is

[TABLE]

Then, notice also that for all $\alpha\in\mathbb{R}^{n}$ and $b\in\mathbb{R}^{K}$, we have

[TABLE]

Hence, because $\frac{1}{\sqrt{n}}\lvert\lvert M_{X}(y-\alpha)\rvert\rvert_{2}+\frac{\lambda}{n}\lvert\lvert\alpha\rvert\rvert_{1}$ is feasible, it holds that

[TABLE]

Under assumptions developed below, this procedure yields consistent estimation and asymptotic normality for $\widehat{\beta}$. Remark that model (1) can be seen as a standard linear model with the coefficient $\alpha_{i}$ corresponding to the slope parameter of a dummy variable which value is $1$ for the individual $i$ and [math] otherwise. Hence, our analysis of the square-root lasso fits in the framework of Belloni et al. (2011a). However, our approach is met with additional technical difficulties because we penalize only a subset of the variables and there is no hope to estimate $\alpha$ consistently as each of its entries is indirectly observed only once. As a result, we develop new assumptions and theorems that are better suited for the purposes of this paper.

2.3 Assumptions and results

The main assumption concerns the choice of the penalty level:

Assumption 2.1

We have $\lim\limits_{n\to\infty}\mathbb{P}\left(\lambda\geq 2\sqrt{n}\frac{\lvert\lvert M_{X}\epsilon\rvert\rvert_{\infty}}{\lvert\lvert M_{X}\epsilon\rvert\rvert_{2}}\right)=1.$

The tuning of $\lambda$ prescribed by this penalty level depends on the distributional assumptions made on $\epsilon$, in particular on the tails. The next lemma provides guidance on how to choose the regularization parameter according to assumptions on $\epsilon$:

Lemma 2.1

It holds that $2\sqrt{n}\frac{\lvert\lvert M_{X}\epsilon\rvert\rvert_{\infty}}{\lvert\lvert M_{X}\epsilon\rvert\rvert_{2}}=2\frac{\lvert\lvert\epsilon\rvert\rvert_{\infty}}{\sigma}+o_{P}(\lvert\lvert\epsilon\rvert\rvert_{\infty})+O_{P}(1)$. Additionally, if $\psi$ is such that $\lim\limits_{n\to\infty}\mathbb{P}\left(\psi\geq 2\frac{\lvert\lvert\epsilon\rvert\rvert_{\infty}}{\sigma}\right)=1$ and $\varphi\to\infty$, then for any $c>1$, $\lambda=c\psi+\varphi$ satisfies Assumption 2.1.

The proof is given in Appendix A. This lemma suppresses the role of the matrix $X$ in the choice of the penalty and simplifies the decision procedure. It leads to the subsequent corollary:

Corollary 2.1

The following hold:

• (i)

If $\epsilon_{i}$ are gaussian random variables, then $\lambda=2c\sqrt{2\log(n)}$ satisfies Assumption 2.1 for any $c>1$;

• (ii)

If $\epsilon_{i}$ are sub-gaussian random variables, then there exists a constant $c>0$ such that $\lambda=c\sqrt{\log(n)}$ satisfies Assumption 2.1;

• (iii)

If $\epsilon_{i}$ are sub-exponential random variables, then then there exists a constant $c>0$ such that $\lambda=c\log(n)$ satisfies Assumption 2.1.

The proof is given in Appendix A. The statistician can use Corollary 2.1 to decide on the penalization parameter given how heavy she expects the tails of the error term to be in her data. In practice, it is advised to choose the smallest penalty verifying Assumption 2.1. This can be done by Monte-Carlo simulations. Notice that our approach allows for heavy-tailed distributions such as sub-exponential random variables.

To derive the convergence rate of our estimator, we first bound the estimation error on $\alpha$ and obtain the following result:

Lemma 2.2

Under Assumption 2.1 and if $p\max\left(\lambda,\sqrt{\left|X\right|_{\infty}}\right)=o_{P}(1)$ (and , it holds that

[TABLE]

The proof is given in Appendix B. The rate of convergence of $\lvert\lvert\widehat{\alpha}-\alpha\rvert\rvert_{1}/n$ therefore is lower than $p\sqrt{\log(n)}$ if the errors are gaussian or sub-gaussian and we choose the penalty level as in Lemma 2.1. Note that, as standard in works related to the lasso estimator (see Bühlmann and Van De Geer (2011)), in our proof we make use of a compatibility condition that states that a compatibility constant is bounded from below with probability approaching one. The condition that $p\lvert\lvert X\rvert\rvert_{\infty}=o_{P}(1)$ is enough to show that this property holds as shown in Lemma B.1 in Appendix B. It is possible to find other sufficient conditions but it is outside the scope of this paper. Remark that if $\{x_{i}\}_{i}$ are i.i.d. sub-Gaussian random variables then $\lvert\lvert X\rvert\rvert_{\infty}=O_{P}\left(\sqrt{\log(n)}\right)$ allowing for the sparsity level $p=o_{P}(1/\sqrt{\log(n)})$.

Here, we show how to derive the rate of convergence of $\widehat{\beta}$ from Lemma 2.2. Assume that $p\max\left(\lambda,\left|X\right|_{\infty}\right)=o_{P}(1)$. Substituting $y$ by $X\beta+\alpha+\epsilon$ in (2), we obtain

[TABLE]

Now, notice that $\left(X^{\top}X\right)^{-1}X^{\top}(\alpha-\widehat{\alpha})=(X^{\top}X/n)^{-1}X^{\top}(\alpha-\widehat{\alpha})/n$. By the law of large numbers, $\left(X^{\top}X/n\right)^{-1}=O_{P}(1)$, which implies that

[TABLE]

By Lemma 2.2, this implies that

[TABLE]

Finally, by the central limit theorem and Slutsky’s lemma, we have that $\sqrt{n}(X^{\top}X)^{-1}X^{\top}\epsilon\xrightarrow[]{d}\mathcal{N}(0,\sigma\Sigma^{-1})$. This leads to Theorem 2.2.

Theorem 2.2

Under Assumption 2.1 and if $p\max\left(\lambda,\sqrt{\left|X\right|_{\infty}}\right)=o_{P}(1)$, it holds that

[TABLE]

This result allows to derive the rates of convergence under different assumptions on the tails of the distributions of the regressors and the error term. For instance, if $\{x_{i}\}_{i}$ and $\{\epsilon_{i}\}_{i}$ are i.i.d. sub-Gaussian random variables, then $\widehat{\beta}$ is consistent as long as $p\log(n)\to 0$ for the choice of $\lambda$ proposed in Lemma 2.1. In this case, this implies that our estimator reaches (up to a logarithmic factor) the minimax lower bound for the Huber’s contamination location model under gaussian errors, which is $\max(1/n,p^{2})$ in $\ell_{2}$-norm according to Chen et al. (2018). We attain the rate $\max(1/n,p^{2}\log(n))$. Remark also that equation (5) explains the role of $\lvert\lvert X\rvert\rvert_{\infty}$ in the convergence rate of $\widehat{\beta}$. For an individual $i$, if $x_{i}$ is large then an error in the estimation of $\alpha_{i}$ can contribute to an error in the estimation of $\beta$ via the term $(X^{\top}X)^{-1}X^{\top}(\alpha-\widehat{\alpha})$ in (4). $\lvert\lvert X\rvert\rvert_{\infty}$ measures the maximum influence that an observation can have.

To show that our estimator is asymptotically normal, it suffices to assume that the bias term $(X^{\top}X)^{-1}X^{\top}(\alpha-\widehat{\alpha})$ in (4) vanishes asymptotically:

Theorem 2.3

Under Assumption 2.1, assuming that $p\lambda\lvert\lvert X\rvert\rvert_{\infty}\sqrt{n}=o_{P}(1)$ (and $p\max\left(\lambda,\sqrt{\left|X\right|_{\infty}}\right)=o_{P}(1)$), we have

[TABLE]

Moreover, $\widehat{\sigma}^{2}=\frac{1}{n}\sum_{i=1}^{n}(y_{i}-x_{i}^{\top}\hat{\beta}-\widehat{\alpha})^{2}$ and $\widehat{\Sigma}=\frac{1}{n}\sum_{i=1}^{n}x_{i}x_{i}^{\top}$ are consistent estimators of, respectively, $\sigma^{2}$ and $\Sigma$.

The proof that $\widehat{\sigma}^{2}\xrightarrow{\mathbb{P}}\sigma^{2}$ is given in Appendix C.When the entries of $X$ and $\epsilon$ are sub-Gaussian, for the choice of the penalty prescribed in Lemma 2.1, the contamination level needs to satisfy $p\log(n)\sqrt{n}\to 0$ to be able to use 2.3 to prove asymptotic normality. Notice that the asymptotic variance of our estimator corresponds to the one of the OLS estimator in the standard linear model under homoscedasticity. Hence, confidence sets and tests can be built in the same manner as in the theory of the OLS estimator.

An important last remark concerns the meaning of confidence intervals developed using Theorem 2.3. Note that they are obtained under an asymptotic with triangular array data under which the number of outliers is allowed to go to infinity while the proportion of outliers goes to [math]. The interpretation of a 95% confidence interval $I$ built with Theorem 2.3 is as follows: if the number of outliers in our data is low enough and the sample size is large enough, then there is a probability of approximatively $0.95$ that $\beta$ belongs to $I$.

3 Computation and simulations

3.1 Iterative algorithm

We propose the following algorithm to compute our estimator. Because $u=\min_{\sigma>0}\left\{\frac{\sigma}{2}+\frac{1}{2\sigma}u^{2}\right\}$, as long as $\left|\left|y-X\widehat{\beta}-\widehat{\alpha}\right|\right|_{2}^{2}>0$, we have that

[TABLE]

This is a convex objective and we propose to iteratively minimize over $\beta$, $\alpha$, and $\sigma$. Let us start from $\left(\beta^{(0)},\alpha^{(0)},\sigma^{(0)}\right)$ and compute the following sequence for $t\in\mathbb{N}^{*}$ until convergence:

$\beta^{(t+1)}\in\arg\min_{\beta\in\mathbb{R}^{K}}\left|\left|y-X\beta-\alpha^{(t)}\right|\right|_{2}^{2};$ 2. 2.

$\alpha^{(t+1)}\in\arg\min_{\alpha\in\mathbb{R}^{n}}\left|\left|y-X\beta^{(t+1)}-\alpha\right|\right|_{2}^{2}+\frac{\lambda\sigma^{(t)}}{\sqrt{n}}\left|\left|\alpha\right|\right|_{1};$ 3. 3.

$\sigma^{(t+1)}=\left|\left|y-X\beta^{(t+1)}-\alpha^{(t+1)}\right|\right|_{2}$.

The following lemma explains how to perform step 2:

Lemma 3.1

For $i=1,\dots,n$, if $\left|y_{i}-(X\beta^{(t+1)})_{i}\right|\leq\frac{\lambda\sigma^{(t)}}{\sqrt{n}}$ then $\alpha^{(t+1)}_{i}=0$. If $\left|y_{i}-(X\beta^{(t+1)})_{i}\right|>\frac{\lambda\sigma^{(t)}}{\sqrt{n}}$ then $\alpha^{(t+1)}_{i}=y_{i}-\left(X\beta^{(t+1)}\right)_{i}-\text{sign}\left(y_{i}-\left(X\beta^{(t+1)}\right)_{i}\right)\frac{\lambda\sigma^{(t)}}{\sqrt{n}}$.

The proof is given in Appendix D.

3.2 Simulations

We apply this computation approach in a small simulation exercise. The data generating process is as follows: there are two regressors $x_{1i}$ and $x_{2i}$, with $x_{1i}=1$ for all $i$ and $x_{2i}$ are i.i.d. $\mathcal{N}(0,1)$ random variables. $\epsilon_{i}$ are i.i.d. $\mathcal{N}(0,1)$ random variables. Then, we set

[TABLE]

where $q_{1-p}$ is such that $\mathbb{P}(x_{2i}\geq q_{1-p})=p$. In table 1, we present the bias, the variance, the mean squared error (MSE) and the coverage of $95\%$ confidence intervals for our estimator $\widehat{\beta}$ computed using the algorithm of Section 3.1, where we use $100$ iterations and with $\lambda=2.01\sqrt{2\log(n)}$. This choice corresponds to the one outlined in Corollary 2.1. The bias, the variance and the coverage of $95\%$ confidence intervals for the naive OLS estimator:

[TABLE]

are also reported. For the OLS estimator, the confidence intervals correspond to the ones of the standard linear model. The presented results are averages among $8000$ replications. We observe that our estimator brings a substantial improvement in estimation precision with respect to the OLS estimator.

Appendix A Choice of the penalization parameter

A.1 Proof of Lemma 2.1

We start by proving the next two technical lemmas:

Lemma A.1

It holds that $\lvert\lvert P_{X}\epsilon\rvert\rvert_{\infty}=O_{P}(1)$.

Proof. Because of the assumptions on the joint distribution of $(x_{i},\epsilon_{i})$, we have that $\sqrt{n}(X^{\top}X)^{-1}X^{\top}\epsilon\xrightarrow[]{d}\mathcal{N}(0,\sigma\Sigma^{-1})$, therefore $\sqrt{n}\lvert\lvert(X^{\top}X)^{-1}X^{\top}\epsilon\rvert\rvert_{2}=O_{P}(1)$. Because $X(X^{\top}X)^{-1}X^{\top}\epsilon=\frac{X}{\sqrt{n}}\sqrt{n}(X^{\top}X)^{-1}X^{\top}\epsilon$, we obtain that $\lvert\lvert P_{X}\epsilon\rvert\rvert_{2}\leq\frac{\lvert\lvert X\rvert\rvert_{2}}{\sqrt{n}}\sqrt{n}\lvert\lvert(X^{\top}X)^{-1}X^{\top}\epsilon\rvert\rvert_{2}=O_{P}\left(\frac{\lvert\lvert X\rvert\rvert_{2}}{\sqrt{n}}\right)=O_{P}(1)$, by the law of large numbers. $\Box$

Lemma A.2

It holds that $\frac{\sqrt{n}}{\lvert\lvert M_{X}\epsilon\rvert\rvert_{2}}-\frac{1}{\sigma}=o_{P}(1)$.

Proof. First, remark that, by the theorem of Pythagore,

[TABLE]

Now, this leads to $\frac{1}{n}\lvert\lvert M_{X}\epsilon\rvert\rvert_{2}^{2}=\frac{1}{n}\lvert\lvert\epsilon\rvert\rvert_{2}^{2}-\frac{1}{n}\epsilon^{\top}X(X^{\top}X)^{-1}X^{\top}\epsilon$. Because $\{x_{i}\}_{i}$ and $\{\epsilon_{i}\}_{i}$ are i.i.d. and $\mathbb{E}[x_{i}\epsilon_{i}]=0$, we have that $\sqrt{n}(X^{\top}X)^{-1}X^{\top}\epsilon\xrightarrow[]{d}\mathcal{N}(0,\sigma\Sigma^{-1})$ and $\frac{1}{\sqrt{n}}X^{\top}\epsilon\xrightarrow[]{d}\mathcal{N}(0,\sigma\Sigma)$. This implies that $\epsilon^{\top}X(X^{\top}X)^{-1}X^{\top}\epsilon=O_{P}(1/n)$. We also have that $\frac{1}{n}\lvert\lvert\epsilon\rvert\rvert_{2}^{2}\xrightarrow[]{\mathbb{P}}\sigma^{2}$, which leads to $\frac{1}{n}\lvert\lvert M_{X}\epsilon\rvert\rvert_{2}^{2}\xrightarrow[]{\mathbb{P}}\sigma^{2}$. We conclude by the continuous mapping theorem. $\Box$

Now, we proceed with the proof of Lemma 2.1. Notice that

[TABLE]

Using lemmas A.1 and A.2, we obtain

[TABLE]

The rest of the lemma is a direct consequence of (7) and the pigeonhole principle.

A.2 Proof of Corollary 2.1

Proof of (i) By Lemma 2.1 it is sufficient to show that for $c>1$,

[TABLE]

Let us remember the gaussian bound (see Lemma B.1 in Giraud (2014)): for $t\geq 0$, we have

[TABLE]

Then, we have

[TABLE]

Proof of (ii) By Lemma 2.1 it is sufficient to show that there exists $c>0$ such that

[TABLE]

Let us remember the sub-gaussian bound (see Proposition 2.5.2 in Vershynin (2018)): for $t\geq 0$, there exists $b>0$ such that

[TABLE]

Then, we have

[TABLE]

Proof of (iii) By Lemma 2.1 it is sufficient to show that there exists $c>0$ such that

[TABLE]

Let us remember the sub-exponential bound (see Proposition 2.7.1 in Vershynin (2018)): for $t\geq 0$, there exists $b>0$ such that

[TABLE]

Then, we have, for $n$ large enough,

[TABLE]

Appendix B Proof of lemma 2.2

B.1 Compatibility constant

For $\delta\in\mathbb{R}^{n}$, we denote by $\delta_{J}\in\mathbb{R}^{n}$ the vector for which $(\delta_{J})_{i}=\delta_{i}$ if $\alpha_{i}\neq 0$ and $(\delta_{J})_{i}=0$ otherwise. Let us also define $\delta_{J^{c}}=\delta-\delta_{J}$. We introduce the following cone:

[TABLE]

We work with the following compatibility constant (see Bühlmann and Van De Geer (2011) for a discussion of the role of compatibility conditions in the lasso literature) corresponding to

[TABLE]

We use the following lemma:

Lemma B.1

If $p^{2}\lvert\lvert X\rvert\rvert_{\infty}=o_{P}(1)$, there exists $\kappa_{*}>0$ such that $\kappa>\kappa_{*}$ w.p.a. $1$.

Proof. Take $\delta\in C$, to show this result, notice that

[TABLE]

Therefore, we have

[TABLE]

where $X_{k}=(x_{1k},...,x_{nK})^{\top}$. Next, we have that

[TABLE]

Now, because we have $\frac{1}{n}\sum_{i=1}^{n}x_{i}x_{i}^{\top}\xrightarrow[]{\mathbb{P}}\Sigma$ by the law of large numbers, we obtain that $\left|\left|\left(X^{\top}X/n\right)^{-1}\right|\right|_{\text{op}}=O_{P}(1)$ and that $\sum_{k=1}^{K}\lvert\lvert X_{k}\rvert\rvert_{2}/\sqrt{n}=\sum_{k=1}^{K}\sqrt{\left(X^{\top}X/n\right)_{kk}}=O_{P}(1)$, both implying that$\frac{1}{\sqrt{n}}\sum_{k=1}^{K}\lvert\lvert X_{k}\rvert\rvert_{2}\left|\left|\left(X^{\top}X/n\right)^{-1}\right|\right|_{\text{op}}=O_{P}(1)$. We conclude the proof using that $p^{2}\lvert\lvert X\rvert\rvert_{\infty}=o_{P}(1)$. $\Box$

B.2 End of the proof of Lemma 2.2

Throughout this proof, we work on the event

[TABLE]

which has probability approaching $1$ according to Assumption 2.1, Lemma B.1, and the condition that $p\lambda\to 0$. Let us define $\Delta=\widehat{\alpha}-\alpha$. Now, remark that

[TABLE]

Next, we use the fact that $\left|\left|\alpha+\Delta_{J^{c}}\right|\right|_{1}=\lvert\lvert\alpha\rvert\rvert_{1}+\left|\left|\Delta_{J^{c}}\right|\right|_{1}$. Combining this and (9), we get

[TABLE]

By definition of $\widehat{\alpha}$ and concentrating our objective function in $\beta$, we have

[TABLE]

By convexity, if $M_{X}\epsilon\neq 0$, it holds that

[TABLE]

where (12) comes from $\lambda\geq 2\sqrt{n}\lvert\lvert M_{X}\epsilon\rvert\rvert_{2}/\lvert\lvert M_{X}\epsilon\rvert\rvert_{\infty}$. This last inequality is also straightforwardly true when $M_{X}\epsilon=0$. This and (11) imply

[TABLE]

Using (10), we get

[TABLE]

Then, as $\lvert\lvert\Delta\rvert\rvert_{1}=\left|\left|\Delta_{J^{c}}\right|\right|_{1}+\left|\left|\Delta_{J}\right|\right|_{1}$, we obtain

[TABLE]

which implies that $\Delta\in C$. Using $y=X\beta+\alpha+\epsilon$, we get

[TABLE]

By Hölder’s inequality, this results in

[TABLE]

Because $\lambda\geq 2\sqrt{n}\frac{\lvert\lvert M_{X}\epsilon\rvert\rvert_{\infty}}{{\lvert\lvert M_{X}\epsilon\rvert\rvert_{2}}}$ , we obtain

[TABLE]

This implies that

[TABLE]

By equations (10) and (11), we have $\frac{1}{\sqrt{n}}\lvert\lvert M_{X}(y-\widehat{\alpha})\rvert\rvert_{2}-\frac{1}{\sqrt{n}}\lvert\lvert M_{X}(y-\alpha)\rvert\rvert_{2}\leq\frac{\lambda}{n}\left(\left|\left|\Delta_{J}\right|\right|_{1}-\left|\left|\Delta_{J^{c}}\right|\right|_{1}\right)$. Using the fact that $\Delta\in C$ and (12), this yields

[TABLE]

Next, notice that

[TABLE]

This implies

[TABLE]

Combining (15) and (16) and remarking that $\lvert\lvert M_{X}\epsilon\rvert\rvert_{2}=\lvert\lvert M_{X}(y-\alpha)\rvert\rvert_{2}$, we obtain

[TABLE]

Now, as $\Delta\in C$, this implies that

[TABLE]

From now on assume that $\lvert\lvert M_{X}\Delta\rvert\rvert_{2}\neq 0$, we get

[TABLE]

which implies again that

[TABLE]

Finally, as $\Delta\in C$, we have

[TABLE]

The last inequality also holds if $M_{X}\Delta=0$ because, as $\kappa>\kappa_{*}$, this implies that $\Delta_{J}=0$ and hence $\Delta=0$ using the fact that $\Delta$ belongs to $C$. To conclude the proof, use (17), the fact that $\lvert\lvert M_{X}\epsilon\rvert\rvert_{2}/\sqrt{n}\leq\lvert\lvert\epsilon\rvert\rvert_{2}/\sqrt{n}=o_{P}(1)$ by the law of large numbers and the condition $p\max\left(\lambda,\lvert\lvert X\rvert\rvert_{\infty}\right)=o_{P}(1)$.

Appendix C Proof that $\mathbf{\widehat{\sigma}^{2}\xrightarrow{\mathbb{P}}\sigma^{2}}$ in Theorem 2.3

We have

[TABLE]

Next, remark that there exists $c,\eta>0$ such that $\mathbb{P}\left(\left|X\right|_{\infty}\geq c\right)\geq\delta$ for $n$ large enough. Indeed otherwise, this would imply that the law of large numbers cannot hold for $x_{i}$. This yields that $\sqrt{n}p\lambda=o(1)$. Then, because of Lemma 2.2, Theorem 2.3 and $p\lambda\left|X\right|_{\infty}=o_{P}(1)$ it holds that

[TABLE]

Next, we have

[TABLE]

by the law of large numbers. Then, by Hölder’s inequality, we obtain that

[TABLE]

By the inequality of Cauchy-Schwartz, this also leads to $\frac{\left<X\left(\beta-\widehat{\beta}\right),\alpha-\widehat{\alpha}\right>}{n}=o_{P}(1)$. Then, the law of large numbers implies that $\frac{\left<X\left(\beta-\widehat{\beta}\right),\epsilon\right>}{n}=o_{P}(1)$ and $\frac{\left<\alpha-\widehat{\alpha},\epsilon\right>}{n}=o_{P}(1)$, which concludes the proof.

Appendix D Proof of Lemma 3.1

By Lemma D.5 in Giraud (2014), there exists $\widehat{z}$, a random vector in $\mathbb{R}^{n}$, such that the first-order conditions of step 2 are

[TABLE]

where, for $i=1,\dots,n$, $\widehat{z}_{i}\in[-1,1]$ if $\alpha^{(t+1)}_{i}=0$ and $\widehat{z}_{i}=\text{sign}\left(\alpha^{(t+1)}_{i}\right)$ if $\alpha^{(t+1)}_{i}\neq 0$. This yields that, if $\alpha^{(t+1)}_{i}\neq 0$,

[TABLE]

Hence, if $\alpha^{(t+1)}_{i}>0$, we obtain

[TABLE]

and, therefore, $y_{i}-\left(X\beta^{(t+1)}\right)_{i}>\frac{\lambda\sigma^{(t)}}{\sqrt{n}}\geq 0$. Similarly, if $\alpha^{(t+1)}_{i}<0$, we have

[TABLE]

and, therefore, $y_{i}-\left(X\beta^{(t+1)}\right)_{i}<-\frac{\lambda\sigma^{(t)}}{\sqrt{n}}\leq 0$. This shows that, if $\alpha^{(t+1)}_{i}\neq 0$, we have

[TABLE]

and

[TABLE]

Next, if $\alpha^{(t+1)}_{i}=0$, (18) implies that

[TABLE]

This concludes the proof.