Outlier Detection Using Nonconvex Penalized Regression

Introduction

Outliers are a pervasive problem in statistical data analysis. Nonrigorously, outliers refer to one or more observations that are different from the bulk of the data. ? estimate that a routine data set may contain about 1-10% (or more) outliers. Unfortunately, outliers often go unnoticed [Roussbook], although they may have serious effects in estimation, inference, and model selection [weisb]. Perhaps the most popular statistical modeling method is ordinary least squares (OLS) regression. OLS is very sensitive to outliers — a single unusual observation may break it down completely. Our goal in this work is outlier identification for regression models, together with robust coefficient estimation.

Suspected outliers are most commonly found by looking at residuals $r_{i}=y_{i}-{\boldsymbol{x}}_{i}^{\mathsf{T}}\hat{\boldsymbol{\beta}}$ where $\hat{\boldsymbol{\beta}}$ is the OLS estimate of ${\boldsymbol{\beta}}$. It is well-known that such raw residuals can fail to detect outliers at leverage points. A better way to detect an outlier is the leave-one-out approach [weisb]. If the $i$’th case is suspected to be an outlier, then we compute the externally studentized residual

where ${\boldsymbol{X}}_{(i)}$, $\hat{\boldsymbol{\beta}}_{(i)}$ and $\hat{\sigma}_{(i)}$ are the predictor matrix, coefficient estimate and scale estimate respectively, based on $n-1$ observations, leaving out the $i$’th. Large values $|t_{i}|>\eta$ are then taken to suggest that observation $i$ is an outlier. The threshold $\eta=2.5$ [Roussbook] is a reasonable choice. If ${\boldsymbol{\epsilon}}\sim{\mathcal{N}}(0,\sigma^{2}I)$ then $t_{i}\sim t_{(n-p-1)}$ and we can even attach a significance level to $|t_{i}|$. After removing an apparent outlier from the data, one then looks for others.

Studentized residuals, and other leave-one-out methods such as Cook’s distance and DFFITS, are simple and effective when there is only one outlier. When there are multiple outliers, these simple methods can fail. Two phenomena have been remarked on. In masking, when an outlying $i$’th case has been left out, the remaining outliers cause either a large value of $\hat{\sigma}_{(i)}$ or a small value of $|y_{i}-{\boldsymbol{x}}_{i}^{\mathsf{T}}\hat{\boldsymbol{\beta}}_{(i)}|$, or both, and as a result observation $i$ does not look like an outlier. Therefore, multiple outliers may mask each other and go undetected. In swamping, the effect of outliers is to make $|y_{i}-{\boldsymbol{x}}_{i}^{\mathsf{T}}\hat{\boldsymbol{\beta}}_{(i)}|$ large for a non-outlying case $i$. Swamping could lead one to delete good observations and becomes more serious in the presence of multiple outliers.

In this paper we take the studentized residual as our starting point. The $t$-test for whether observation $i^{\prime}$ is an outlier is the same as testing whether the parameter $\gamma$ is zero in the regression ${\boldsymbol{y}}={\boldsymbol{X}}{\boldsymbol{\beta}}+\gamma 1_{i=i^{\prime}}+{\boldsymbol{\epsilon}}$. Because we don’t know which observations might be outliers, we use a model

in which the parameter $\gamma_{i}$ is nonzero when observation $i$ is an outlier. This formulation was earlier used by ? and ?. This mean-shift model allows any combination of observations to be outliers. It has $n+p$ regression parameters and only $n$ data points. Our approach is to fit (1.2) imposing sparsity on ${\boldsymbol{\gamma}}$ in order to avoid the trivial estimate $\hat{\boldsymbol{\gamma}}={\boldsymbol{y}}$ and to get a meaningful estimate of ${\boldsymbol{\beta}}$. The resulting algorithm is called thresholding (denoted by $\Theta$) based iterative procedure for outlier detection, or $\Theta$-IPOD for short.

All of our proposals (apart from one exception noted where it arises) require a preliminary robust regression to be run. This practice is in line with the best current robust regression methods. The preliminary regression supplies a robust estimate of ${\boldsymbol{\beta}}$, and usually a robust estimate of $\sigma$ as well. The robust regression methods that we compare to are known to outperform the preliminary regressions that they use as starting points. We will compare $\Theta$-IPOD to those best performing methods.

The rest of the paper is organized as follows. Section 2 surveys the literature on robust regression. Section 3 develops the soft-IPOD algorithm which fits (1.2) using an $L_{1}$ penalty on ${\boldsymbol{\gamma}}$. This algorithm minimizes a convex criterion, but it is not robust. Section 4 develops a family of algorithms replacing soft-thresholding by a general thresholding rule $\Theta$. We find that some nonconvex criteria properly identify multiple outliers in some standard challenging test cases. Section 5 investigates the computational efficiency of $\Theta$-IPOD in comparison to iteratively reweighted least squares (IRLS). IRLS requires a QR decomposition at each iteration, while $\Theta$-IPOD needs only one. As a result, though possibly requiring more iterations, it is much faster in large problems that we investigate. In Section 6, we discuss the important problem of parameter tuning in outlier detection and carry out an empirical study to demonstrate the advantage of our penalized approach. Section 7 extends the technique to high-dimensional data with $p>n$. Our conclusions are in Section 8.

Survey of robust regression

Many methods have been developed for multiple outlier identification. Robust or resistant regressions, such as $M$-estimators [Huberbook] and Least Trimmed Squares (LTS) [Roussbook], were proposed to provide a trustworthy coefficient estimate even in the presence of multiple outliers. There are also clustering-based procedures, such as ?, and some informal methods based on graphics. Multiple outlier detection procedures usually alternate between two steps. One step scans regression output (coefficient and variance estimates) to identify seemingly clean observations. The other fits a linear regression model to those clean observations. The algorithm can be initialized with OLS, but generally it is better to initialize it with something more robust.

If a data set contains more than one outlier, masking may occur and the task of outlier detection is much more challenging. Well-known examples of masking include the Belgian telephone data and the Hertzsprung-Russell star data [Roussbook], as well as some artificial datasets like the Hawkins-Bradu-Kass (HBK) data [HBK] and the Hadi-Simonoff (HS) data [HS]. The HS and HBK data sets both have swamping effects.

The main challenge in multiple outlier detection is to counter masking and swamping effects. ? describe two broad classes of algorithms — direct methods and indirect methods. The direct procedures include forward search algorithms [HS97, PY95, atki:rian:2000] and backward selection [menj:wels:2010] among others. Indirect methods are those that use residuals from a robust regression estimate to identify the outliers. Examples of indirect methods include Least Median of Squares (LMS) [lms], Least Trimmed Squares (LTS) [Roussbook], S-estimators [sest], MM-estimators [mmest], one-step GM estimators [*]*gmest and S1S estimators [s1s]. Almost all start with an initial high breakdown point estimate (not necessarily efficient) say from LTS, S, MTS [Nguyen2010], or ? (denoted by PY) fast procedure. Most published examples in the outlier detection literature are in $10$ or fewer dimensions. Methods with high breakdown point typically have costs that grow exponentially in the dimension. In practice, when there are a large number of predictors, PY can be applied to provide an initial estimate with certain robustness, although in theory its breakdown point property is not well established [PYfast].

It is worth mentioning that outlier identification and robust regression are two closely related but not quite identical problems [*]*Roussbook,Yohaibook. Even if we perfectly identified the regression coefficient vector, there could still be some overlap between the residual distributions for good and outlying points. That is, given the true regression coefficient ${\boldsymbol{\beta}}$ there would be type one and type two errors in trying to identify outliers. On the other hand, if we could perfectly identify all gross outliers, it is a relatively easy task to obtain a robust coefficient estimate, as will be supported by our experiments in Section 6.

Soft-IPOD

We will use the mean shift model (1.2) from the introduction which predicts ${\boldsymbol{y}}$ by the usual linear model ${\boldsymbol{X}}{\boldsymbol{\beta}}$ plus an outlier term ${\boldsymbol{\gamma}}$. If $\gamma_{i}=0$ then the $i$’th case is good, and otherwise it is an outlier. Our goals are to find a robust estimate of ${\boldsymbol{\beta}}$ as well as to estimate ${\boldsymbol{\gamma}}$ thereby identifying which cases are outliers and which are not. We assume that ${\boldsymbol{\gamma}}$ is sparse because outliers should not be the norm. We suppose at first that $n>p$ and that ${\boldsymbol{X}}=[{\boldsymbol{x}}_{1},\dots,{\boldsymbol{x}}_{n}]^{\mathsf{T}}$ has full rank $p$. Section 7 considers the case when ${\boldsymbol{\beta}}$ is sparse, too, with $p$ possibly greater than $n$. Yet the majority of our paper focuses on the outlier problem only. Let ${\boldsymbol{H}}$ be the hat matrix defined by ${\boldsymbol{H}}={\boldsymbol{H}}({\boldsymbol{X}})={\boldsymbol{X}}({\boldsymbol{X}}^{\mathsf{T}}{\boldsymbol{X}})^{-1}{\boldsymbol{X}}^{\mathsf{T}}$. The $i$’th diagonal entry of ${\boldsymbol{H}}$, denoted $h_{i}$, is called the leverage of the $i$’th observation.

The assumed sparsity of ${\boldsymbol{\gamma}}$ motivates using an $L_{1}$-penalized regression to minimize

Although (3.1) is a well-formulated model and soft-IPOD is computationally efficient, in the presence of multiple outliers with moderate or high leverage values, this method fails to remove masking and swamping effects. Take the artificial HBK data as an illustration. Using a robust estimate $\hat{\sigma}$ from LTS, and $\lambda_{i}=\hat{\sigma}\sqrt{2(1-h_{i})\log n}$, we obtained $\hat{\boldsymbol{\gamma}}=[0,\cdots,0,-8.6,-9.7,-7.6,-8.4,0,\cdots,0]^{\mathsf{T}}$ which identifies cases 11-14 as serious outliers, while the true ${\boldsymbol{\gamma}}$ is $[10,\cdots,10,0,\cdots,0]^{\mathsf{T}}$ with cases 1-10 being the actual outliers. This erroneous identification is not a matter of parameter tuning. Figure 1 plots the soft-IPOD solution over the relevant range for $\lambda$. Whatever value of $\lambda$ we choose, cases 11-14 are sure to be swamped, and cases 1-10 are very likely to be masked. The $L_{1}$ approach is not able to identify the correct outliers without swamping. Our extensive experience shows that this $L_{1}$ technique hardly works for any benchmark dataset in the outlier detection literature. According to ?, a convex criterion is inherently incompatible with robustness.

Next we delineate the parallels between soft-IPOD and Huber’s M-estimate regression which is similarly non-robust. As pointed out by ? and ?, there is a connection between the $L_{1}$-penalized regression (3.1) and Huber’s $M$-estimate. Huber’s loss function is

Huber’s method with concomitant scale estimation, minimizes

over ${\boldsymbol{\beta}}$ and $\sigma$ jointly, where $c\geq 0$ and $\lambda>0$ are given constants. We can prove the following result, which is slightly more general than the version given by ? and ?.

For any $c\geq 0$ and $\lambda>0$ suppose that we minimize

over ${\boldsymbol{\beta}}$, ${\boldsymbol{\gamma}}$ and $\sigma$. Then the minimizing ${\boldsymbol{\beta}}$ and $\sigma$ match those from minimizing Huber’s criterion (3.2). For any fixed $\sigma>0$ minimizing (3.3) over ${\boldsymbol{\beta}}$ and ${\boldsymbol{\gamma}}$ yields the minimizer of (3.2) over ${\boldsymbol{\beta}}$.

This connection is helpful to understand the inherent difficulty with $L_{1}$-penalized regression described earlier. It is well known that Huber’s method cannot even handle moderate leverage points well [Huberbook, p. 192] and is prone to masking and swamping in outlier detection. Its break-down point is $0$.

A promising way to improve (3.1) is to adopt a different penalty function, possibly nonconvex. Our launching point in this paper is, however, the operator $\Theta$ in the above iterative algorithm. Substituting an appropriate $\Theta$ for soft-thresholding in the IPOD (see Algorithm 1), we may obtain a good estimator of ${\boldsymbol{\gamma}}$ (and ${\boldsymbol{\beta}}$ as well).

Θ\Theta-IPOD

To deal with masking and swamping in the presence of multiple outliers, we consider the iterative procedure described in Section 3 using a general $\Theta$ operator, referred to as $\Theta$-IPOD. A $\Theta$-IPOD estimate is a limit point of $({\boldsymbol{\beta}}^{(j)},{\boldsymbol{\gamma}}^{(j)})$, denoted by $(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\gamma}})$. Somewhat surprisingly, simply replacing soft-thresholding by hard-thresholding (henceforth hard-IPOD) resolves the masking and swamping problem, for the challenging HBK example problem. After picking $\hat{\sigma}$ from LTS and using a zero start we obtained

and $\hat{\boldsymbol{\gamma}}_{11:75}=\boldsymbol{0}$ which perfectly detects the true outliers. The coefficient estimate $\hat{\boldsymbol{\beta}}$ from hard-IPOD directly gives the OLS estimate computed from the clean observations.

Some questions naturally arise: What optimization problem is $\Theta$-IPOD trying to solve? For an arbitrary $\Theta$, does $\Theta$-IPOD converge at all? Or, under what conditions does $\Theta$-IPOD converge? To answer these questions, we limit our discussions of $\Theta$ to thresholding rules.

We begin by defining the class of threshold functions we will study. It includes well-known thresholds such as soft and hard thresholding, SCAD, and Tukey’s bisquare.

A threshold function is a real valued function $\Theta(t;\lambda)$ defined for $-\infty<t<\infty$ with $\lambda$ as the parameter ($0\leq\lambda<\infty$) such that

$\Theta(t;\lambda)\leq\Theta(t^{\prime};\lambda)$ for $t\leq t^{\prime}$,

$\lim_{t\to\infty}\Theta(t;\lambda)=\infty$, and

$0\leq\Theta(t;\lambda)\leq t$ for $0\leq t<\infty$.

In words, $\Theta(\cdot;\lambda)$ is an odd monotone unbounded shrinkage rule for $t$, at any $\lambda$.

A vector version of $\Theta$ is defined componentwise if either $t$ or $\lambda$ are replaced by vectors. When both $t$ and $\lambda$ are vectors, we assume they have the same dimension.

Huber’s soft-thresholding rule corresponds to an absolute error, or $L_{1}$ penalty. More generally, for any thresholding rule $\Theta(\cdot;\lambda)$, a corresponding penalty function $P=P_{\Theta}$ can be defined. There may be multiple penalty functions for a given threshold as demonstrated by hard thresholding (4.7) below. The following three-step construction finds the penalty with the smallest curvature [antrev, SheTISP]:

where $u\geq 0$ holds throughout (4.1). The constructed penalty $P_{\Theta}$ is nonnegative and is continuous in $\theta$.

Let $\Theta$ be a thresholding rule as given by Definition 4.1 and let $P_{\Theta}$ be the corresponding penalty defined in (4.1). Given $\lambda_{i}\geq 0$, the objective function is defined by

where $P(\cdot;\cdot)$ is any function satisfying $P(\theta;\lambda)-P(0;\lambda)=P_{\Theta}(\theta;\lambda)+q(\theta;\lambda)$ where $q(\cdot;\lambda)$ is nonnegative and $q(\Theta(\theta;\lambda))=0$ for all $\theta$. Then the $\Theta$-IPOD iteration sequence $({\boldsymbol{\beta}}^{(j)},{\boldsymbol{\gamma}}^{(j)})$ satisfies

The function $q$ will often be zero, but we use non-zero $q$ below to demonstrate that multiple penalties (infinitely many, as a matter of fact) yield hard thresholding, including the $L_{0}$-penalty. Theorem 4.1 shows that $\Theta$-IPOD converges. Any limit point of $({\boldsymbol{\beta}}^{(j)},{\boldsymbol{\gamma}}^{(j)})$ must be a stationary point of (4.2). $\Theta$-IPOD also gives a general connection between penalized regression (4.2) and $M$-estimators. Recall that an $M$-estimator is defined to be a solution to the score equation

where $\lambda$ is a general parameter of the $\psi$ function. Although ${\boldsymbol{\beta}}$ and $\sigma$ can be simultaneously estimated by Huber’s Proposal 2 [Huberbook], a more common practice is to fix $\sigma$ at an initial robust estimate and then optimize over ${\boldsymbol{\beta}}$ [Hampelbook]. Unless otherwise specified, we consider equation (4.4) as constraining ${\boldsymbol{\beta}}$ with $\sigma$ fixed.

For any thresholding rule $\Theta(\cdot;\lambda)$, for any $\Theta$-IPOD estimate $(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\gamma}})$, $\hat{\boldsymbol{\beta}}$ is an $M$-estimate associated with $\psi$, as long as $(\Theta,\psi)$ satisfies

We have mentioned that Huber’s method or soft-IPOD behaves poorly in outlier detection. The problem is that it never rejects gross outliers that have moderate or high leverage. To reject gross outliers, redescending $\psi$-functions are advocated, corresponding to a class of thresholdings offering little shrinkage for large components, or using nonconvex penalties for solving the sparsity problem (4.2). The differences between $\Theta$-IPOD and the corresponding $M$-estimator are as follows:

$M$-estimators focus on robust estimation of ${\boldsymbol{\beta}}$. For an explicit sparse ${\boldsymbol{\gamma}}$ estimate, a cutoff value is usually needed for the residuals. Minimizing $f_{P}$ from equation (4.2) directly yields a sparse $\hat{\boldsymbol{\gamma}}$ for outlier detection and a robust $\hat{\boldsymbol{\beta}}$.

Instead of designing a robust loss function in $M$-estimators, (4.2) considers a penalty function; $\lambda$ is not a criterion (loss) parameter but a regularization parameter that we will tune in a data-dependent way.

Figure 2 illustrates some of the better known threshold functions along with their corresponding penalties, $\psi$-functions and loss functions. In this article, we make use of the following formulas

for soft and hard thresholding, respectively. The usual penalty for hard thresholding is $P(x;\lambda)=1_{|x|<\lambda}(\lambda|x|-x^{2}/2)+1_{|x|\geq\lambda}\lambda^{2}/2$. Theorem 4.1 justifies use of the $L_{0}$-penalty $P(x;\lambda)=\lambda^{2}/2\cdot I_{x\neq 0}$ using

2 Θ\Theta-IPOD and TISP

Algorithm 1 can be simplified. The $\Theta$-IPOD iteration can be carried out without recomputing ${\boldsymbol{\beta}}^{(j)}$ at each iteration. We need only update ${\boldsymbol{\gamma}}$ via

at each iteration, where $\lambda_{i}=\lambda\sqrt{1-h_{i}}$. The multiplication $({\boldsymbol{I}}-{\boldsymbol{H}}){\boldsymbol{y}}$ in (4.8) can also be precomputed. After getting the final $\hat{\boldsymbol{\gamma}}$, we can estimate ${\boldsymbol{\beta}}$ by OLS. The resulting simplified $\Theta$-IPOD algorithm is given in Algorithm 2. For any given thresholding rule $\Theta$, let $f_{P}({\boldsymbol{\gamma}})\equiv\frac{1}{2}\|({\boldsymbol{I}}-{\boldsymbol{H}})({\boldsymbol{y}}-{\boldsymbol{\gamma}})\|_{2}^{2}+\sum_{i=1}^{n}P(\gamma_{i};\lambda_{i}),$ where $P$ can be any penalty function satisfying the conditions in Theorem 4.1. Then $f_{P}({\boldsymbol{\gamma}}^{(j+1)})\leq f_{P}({\boldsymbol{\gamma}}^{(j)})$ for the $\Theta$-IPOD iterates ${\boldsymbol{\gamma}}^{(j)}$, $j\geq 0$.

Setting up for Algorithm 2 costs $O(np^{2})$ for a dense regression. The dominant cost in a given iteration comes from computing ${\boldsymbol{H}}{\boldsymbol{\gamma}}^{(j)}$. Given a QR decomposition of ${\boldsymbol{H}}$ we can compute that matrix product in $O(np)$ work as $Q(Q^{\mathsf{T}}{\boldsymbol{\gamma}}^{(j)})$. The cost of the update could be even less if ${\boldsymbol{\gamma}}^{(j)}$ has fewer than $p$ nonzero entries and one maintains the dense matrix ${\boldsymbol{H}}$.

To give (4.8) another explanation, suppose the spectral decomposition of the hat matrix ${\boldsymbol{H}}$ is given by ${\boldsymbol{H}}={\boldsymbol{U}}{\boldsymbol{D}}{\boldsymbol{U}}^{\mathsf{T}}$. Define an index set $c=\{i:D_{ii}=0\}$ and let ${\boldsymbol{U}}_{c}$ be formed by taking the corresponding columns of ${\boldsymbol{U}}$. Then a reduced model can be obtained from the mean shift outlier model (1.2)

We can build a connection between the reduced model and simplified $\Theta$-IPOD. ? proposed a class of thresholding-based iterative selection procedures (TISP) for model selection and shrinkage. $\Theta$-TISP for solving the sparsity problem (4.9) is given by

with $k_{0}$ equal to the largest singular value of the Gram matrix ${\boldsymbol{A}}^{\mathsf{T}}{\boldsymbol{A}}={\boldsymbol{H}}$. The iteration (4.10) reduces exactly to (4.8). Therefore, all TISP studies apply to the $\Theta$-IPOD algorithm. For example, the TISP convergence theorem can be used to establish a version of Theorem 4.1, and the nonasymptotic probability bounds for sparsity recovery reveal masking and swamping errors. In particular, the advocated hard-thresholding-like $\Theta$ in TISP corresponds to a redescending $\psi$ in our outlier identification problem.

The simplified procedure (4.8) is easy to implement and is computationally efficient, because the iteration does not involve complicated operations like matrix inversion. Model (4.9) is simpler than the original (1.2) because ${\boldsymbol{\beta}}$ does not appear and all observations are clean. They have non-outlying errors because we have moved the outlier variables into the regression. Using this characterization of $\hat{\boldsymbol{\gamma}}$, it is not difficult to show that the IPOD-estimate $\hat{\boldsymbol{\beta}}$ satisfies the regression, scale, and affine equivariant properties [Hampelbook] desirable for a good robust regression estimator:

$\hat{\boldsymbol{\beta}}({\boldsymbol{X}}\boldsymbol{C},{\boldsymbol{y}})=\boldsymbol{C}^{-1}\hat{\boldsymbol{\beta}}({\boldsymbol{X}},{\boldsymbol{y}})$, for any nonsingular $\boldsymbol{C}$.

We are making here a mild assumption that the initial robust estimate is equivariant, as LTS, S and PY [PYfast] are, and we’re ignoring the possible effect of convergence criteria on equivariance.

Θ\Theta-IPOD vs. IRLS

We consider some computational issues in this section. As Proposition 4.1 suggests, $\Theta$-IPOD solves an $M$-estimation problem. The standard fitting algorithm for $M$-estimates is the well-known iteratively re-weighted least squares (IRLS). Let $w(t;\lambda)=\psi(t;\lambda)/t$, taking $0/0=0$ if necessary. The $\psi$-equation (4.4) that defines an $M$-estimator can be rewritten as $\sum_{i=1}^{n}w\left({r_{i}}/\sigma;\lambda\right)r_{i}{\boldsymbol{x}}_{i}=\boldsymbol{0}$, where $r_{i}=y_{i}-{\boldsymbol{x}}_{i}^{\mathsf{T}}{\boldsymbol{\beta}}$. Accordingly, an $M$-estimate corresponds to a weighted LS estimate as is well known. These multiplicative weights can help downweight the bad observations. Iteratively updating the weights yields the IRLS algorithm, which is the most common method for computing $M$-estimates. Model (1.2) indicates that $M$-estimation can also be characterized through additive effects on all observations.

Given $n$ and $p$, we report the total cost of computing all $M$-estimates for these different combinations of the number of outliers ($O$) and leverage value ($L$), each combination simulated 10 times. The scale parameter ($\lambda$) decreased from $\|({\boldsymbol{I}}-{\boldsymbol{H}}){\boldsymbol{y}}./\sqrt{\mbox{diag}({\boldsymbol{I}}-{\boldsymbol{H}})}\|_{\infty}$ to $0.5$ with fixed step size $-0.1$, where ./ stands for elementwise division. The upper limit is the largest possible standardized residual. Empirically, the lower bound $0.5$ yields approximately half of $\gamma_{i}$ nonzero. Also for $\lambda<0.5$ IRLS often encountered a singular WLS during the iteration, or took exceptionally long time to converge, and so could not be compared to $\Theta$-IPOD. We used IRLS (with fixed $\sigma=1$, as an oracle would have) and simplified $\Theta$-IPOD. The common convergence criterion was $\|{\boldsymbol{\gamma}}^{(j+1)}-{\boldsymbol{\gamma}}^{(j)}\|_{\infty}<10^{-4}$. We studied all sample sizes $n\in\{30,50,100,200,300,400,500,600,700,800,900,1000\}$ and we took $p=n/10$. The CPU times (in seconds) are plotted against the sample size in Figure 3.

The speed advantage of $\Theta$-IPOD over IRLS in this simulation tended to increase with $n$ and $p$. At $n=1000$ the hard-IPOD algorithm required about 2 or 3 times (averaged over all lambda values) as many iterations as IRLS but was about $10$ times faster than IRLS, due to faster iterations. For small $n$, we saw little speed difference. As remarked above, IRLS was sometimes unstable, with singular WLS problems arising for redescending $\psi$, when fewer than $p$ weights were nonzero. Such cases cause only a small problem for $\Theta$-IPOD: the update ${\boldsymbol{H}}{\boldsymbol{\gamma}}^{(j)}$ cannot then take advantage of sparsity, but it is very stable.

Although we attained impressive computational gain over the popular IRLS, we consider speed to be secondary compared to robustness (and the speed advantage of $\Theta$-IPOD can be moot when the preliminary method is very expensive).

Parameter Tuning in Outlier Detection

The parameter $\lambda$ in an $M$-estimators (4.4) is often chosen to be a constant (for all $n$), based on either efficiency or breakdown of the estimator. The value $2.5\hat{\sigma}$ is popular [*]*Roussbook,Wilbook,Yohaibook. But as mentioned in Section 3, even with no outliers and ${\boldsymbol{X}}={\boldsymbol{I}}$, a constant $\lambda$ independent of $n$ is far from optimal [Donoho]. It is also hard, in robust regressions, to select the cutoff value $\eta$ at which to identify outliers, because the residual distribution is usually unknown. ? base an asymptotically efficient choice for $\eta$ on a Kolmogorov-Smirnov statistic, but they need to assume that the standardized robust residuals are IID ${\mathcal{N}}(0,1)$.

Specifically, because all candidate estimates lie along the $\Theta$-IPOD solution path, we design a local BIC to apply BIC on a proper local interval of the degrees of freedom (DF). First we generate the hard-IPOD solution path by decreasing $\lambda$ from $\|({\boldsymbol{I}}-{\boldsymbol{H}}){\boldsymbol{y}}./\sqrt{\mbox{diag}({\boldsymbol{I}}-{\boldsymbol{H}})}\|_{\infty}$ to $0$. Given $\lambda$ and the corresponding estimate $\hat{\boldsymbol{\gamma}}(\lambda)$, let $nz(\lambda)=\{i:\hat{\boldsymbol{\gamma}}_{i}(\lambda)\neq 0\}$. We rely on model (4.9) and study its variable selection to give the correct form of BIC. For hard-IPOD, $\hat{\boldsymbol{\gamma}}_{nz}$ is an OLS estimate with one parameter per detected outlier and the degrees of freedom are given by $\mbox{DF}(\lambda)=|nz(\lambda)|$. We use BIC with a slight modification:

We carried out simulation experiments to test the performance of the tuned $\Theta$-IPOD. The matrix ${\boldsymbol{X}}$ was generated the same way as in Section 5 using dimension $p\in\{15,50\}$, $n=1000$ observations of which the first $O\in\{200,150,100,50,10\}$ were outliers at the highly leveraged location given by $L\in\{15,20\}$ times a vector of $1$s. The outliers were generated by a mean shift ${\boldsymbol{\gamma}}=(\{5\}^{O},\{0\}^{n-O})^{\mathsf{T}}$ added to ${\boldsymbol{y}}$. Because $\Theta$-IPOD is affine, regression, and scale equivariant and the outliers are at a single $x$ value, we may set ${\boldsymbol{\beta}}=\boldsymbol{0}$ without loss of generality. The intercept term is always included in the modeling.

Five outlier detection methods were considered for comparison: hard-IPOD (tuned), MM-estimator, ? fully efficient one-step procedure (denoted by GY), the compound estimator S1S, and the LTS. (The direct procedures, such as ?, ?, behaved poorly and their results were not reported.) The S-PLUS Robust library provides implementations of MM, GY, and LTS; for the implementation of S1S, we refer to ?. Robust also provides a default initial estimate with high breakdown point (not necessarily efficient), which is used in the first three algorithms in our experiments: when $p=15$, it is the S-estimate computed via random resampling; when $p=50$, it is the estimate from the fast PY procedure. Since it is well known that the initial estimate is outperformed by MM in terms of estimation efficiency and robustness, both theoretically and empirically [mmest, fastS, PYfast], its results are not reported. (All of these routines are available for the R language [Rlang] as well. The package robust provides an R version of the Insightful Robust library.)

All methods apart from $\Theta$-IPOD require a cutoff value to identify which residuals are outliers. We applied the fully efficient procedure which performs at least as well as the fixed choice of $\eta=2.5$ in various situations [GYeff].

Each model was simulated $100$ times. We report outlier identification results for each algorithm using three benchmark measures:

In outlier detection, masking is more serious than swamping. The former can cause gross distortions while the latter is often just a matter of lost efficiency.

Ideally, $\mbox{M}\approx 0$, $\mbox{S}\approx 0$, and $\mbox{JD}\approx 100\%$. JD is the most important measure on easier problems while M makes the most sense for hard problems. The simulation results are summarized in Tables 1 and 2. Figures 4 and 5 present M and JD for $p=15$ and $50$ respectively. While our main purpose is to identify outliers, a robust coefficient estimate $\hat{\boldsymbol{\beta}}$ can be easily obtained from $\Theta$-IPOD. The MSE in ${\boldsymbol{\beta}}$ for $p=15$ is shown in Figure 6. All methods had small slope errors, though $\Theta$-IPOD and GY performed best. The results for $p=50$ (not shown) are similar.

MM and GY are two standard methods provided by the S-PLUS Robust library. Nevertheless, as seen from the Tables 1 and 2, the MM-estimator, though perhaps most popular in robust analysis, does not yield good identification results when the outliers have high leverage values and the number of outliers is not small (for example, $O=200$, $L=15,20$). GY improves MM a lot in this situation and gives similar results otherwise. The experiments also show that S1S behaves poorly in outlier detection and LTS works better in the presence of $\leq 5\%$ outliers. Unfortunately, all four methods have high masking probabilities and very low joint identification rates, which become worse for large $p$. The $\Theta$-IPOD method dominates them significantly for masking and joint detection.

To judge statistical significance of these MC results we constructed paired $t$ statistics based on our $100$ replicates. The numerator in each was the number of outliers missed by a competing method minus the number missed by $\Theta$-IPOD; the denominator was the standard error of the numerator. Most of the $t$-statistics were larger than $3$ and grew rapidly with $O$. LTS versus $\Theta$-IPOD had the smallest $t$-statistic ($1.2$ for $O=15$) but also the largest ($\sim 1700$ for $O=200$). While $\Theta$-IPOD does much better on masking, it is slightly worse for swamping. This is an acceptable tradeoff because masking causes far more harm.

Outlier Detection with p>np>n

Here we extend outlier detection to problems with $p>n$ including of course some high-dimensional problems. This context is more challenging but has diverse modern applications in signal processing, bioinformatics, finance, and neuroscience among others. Performing the task of outlier identification for data with $p>n$ or even $p\gg n$ goes beyond the traditional robust analysis which requires a large number of observations relative to the dimensionality.

The convergence of (7.1) is guaranteed even for $n<p$ [SheTISP], as long as $[{\boldsymbol{X}}\ {\boldsymbol{I}}]$ is scaled properly, which here amounts to taking $k_{0}\geq\sqrt{\|{\boldsymbol{X}}\|_{2}^{2}+1}$. Let $\hat{\boldsymbol{\beta}}$ and $\hat{\boldsymbol{\gamma}}$ denote the estimates from (7.1). The nonzero components of $\hat{\boldsymbol{\beta}}$ locate the relevant predictors while $\hat{\boldsymbol{\gamma}}$ identifies the outliers and measures their outlyingness. (We could also apply (7.1) when $p<n$ for simultaneous variable selection and outlier detection.) Not surprisingly, $L_{1}$ (or soft-thresholding) fails again for this challenging sparsity problem and thresholdings corresponding to redescending $\psi$’s should be used in (7.1). It is natural to consider the elastic net [ZouHas] for this problem as well. But the elastic net has a convex criterion, and it can be shown to have breakdown zero.

Since it is usually unknown how severe the outliers are, one may adopt the following hybrid of hard-thresholding and ridge-thresholding

referred to as the hybrid-thresholding in ?, or hard-ridge thresholding in this paper. The corresponding penalty from the three-step construction is

Setting $q(x;\lambda,\eta)=\frac{1+\eta}{2}(|x|-\lambda)^{2}1_{0<|x|<\lambda}$, we get an alternative penalty

From (7.3), we see that hard-ridge thresholding successfully fuses the $L_{0}$-penalty and the $L_{2}$-penalty (ridge-penalty) via $\Theta$-thresholding. The $L_{0}$ portion induces sparsity, while the $L_{2}$ portion shrinks ${\boldsymbol{\beta}}$. The difference from hard thresholding is as follows. When hard thresholding makes $\hat{\gamma}_{i}\neq 0$ we get $\hat{\gamma}_{i}=r_{i}$, removing any influence of observation $i$ on $\hat{\beta}$. With hard-ridge thresholding, the influence of the observation can be removed partially, with the extent of removal controlled by $\eta$. This is helpful for nonzero but small $|\gamma_{i}|$ corresponding to mild outlyingness. In addition, the $L_{2}$ shrinkage also plays a role in estimation and prediction.

The large $p$ setting brings a computational challenge. Because the large-$p$ sparse $\Theta$-IPOD can be very slow, we used the following proportional $\Theta$-IPOD to screen out some ‘nuisance dimensions’ with coefficients being exactly zero (due to our sparsity assumption). Concretely, at each update in (7.1), $\lambda$ was chosen to get precisely $\alpha n$ nonzero components in the new $({\boldsymbol{\beta}}^{(j+1)},{\boldsymbol{\gamma}}^{(j+1)})$. We used $\alpha=0.75$ though other choices could be made. Proportional $\Theta$-IPOD yields at most $\alpha n$ candidates predictors to have nonzero coefficients, making the problem simpler. Next, we run (7.1) with just those predictors, getting a solution path in $\lambda$ and choosing $\lambda$ by BIC∗. In proportional $\Theta$-IPOD, a variable that gets killed at one iteration may reappear in the fit at a later iteration. This is quite different from independent screenings based on marginal statistics like FDR [fdr] or SIS [fanlv].

For our example we used the proportional hard-ridge IPOD. Specifically, for $\eta$ in a small grid of values, we ran proportional $\Theta$-IPOD using the hard-ridge thresholding function $\Theta$. We chose $\eta$ by BIC∗.

We perform the joint robust variable selection and outlier detection on the sugar data of ?. The data set concerns NIR spectroscopy of compositions of three sugars in aqueous solution. We concentrate on glucose (sugar 2) as the response variable. The predictors are second derivative spectra of 700 absorbances at frequencies corresponding to wavelengths of 1100nm to 2498nm in steps of 2nm. There are 125 samples for model training. A test set with 21 test samples is available. These 21 samples were specially designed to be difficult to predict, with compositions outside the range of the other 125. See ? for details.

? used an MCMC Bayesian approach for variable selection, but only 160 equally spaced wavelengths were analyzed due to the computational cost. We used all $700$ wavelengths in fitting our sparse robust model that also estimates outliers. The problem size is $125\times(700+125)$. After the screening via proportional hard-ridge-IPOD, we ran hard-ridge-IPOD and tuned it as follows. First we set $\lambda=0$ which turns $\eta$ into an ordinary ridge parameter. We fit that ridge parameter $\eta$, obtaining $\eta^{*}$ by minimizing $\mbox{BIC}^{*}$ of Section 6. Then for each $\eta$ in the grid $\{0.5\eta^{*},0.05\eta^{*},0.005\eta^{*}\}$, we found $\lambda(\eta)$ to minimize $\mbox{BIC}^{*}$, and finally chose among the three $(\lambda(\eta),\eta)$ combinations to minimize $\mbox{BIC}^{*}$. We think there is no reason to use $\eta>\eta^{*}$ and a richer grid than the one we used would be reasonable, but we chose a small grid for computational reasons. Our experience is that even small $\eta$ values improve prediction.

Figure 7 plots the final estimates $\hat{\boldsymbol{\beta}}$ and $\hat{\boldsymbol{\gamma}}$. From the left panel, we see that the estimated model is sparse. It selects only $15$ of the wavelengths. The estimate $\hat{\boldsymbol{\gamma}}$ shown in right side of Figure 7 suggests that observation $99$ might be an outlier. We found $\hat{\gamma}_{99}=1.7706$ and $r_{99}=y_{99}-{\boldsymbol{x}}_{99}^{\mathsf{T}}\hat{\boldsymbol{\beta}}=1.7711$, which indicates that this observation is (almost) unused in the model fitting. Our model has good prediction performance; the mean-squared error is 0.219 on the test data, improving the reported MCMC results by about 39%. The robust residual plot is shown in Figure 8.

Discussion

The main contribution of this paper is to consider the class of $\Theta$-estimators [SheTISP], under the mean shift outlier model assumption, defined by the fixed point ${\boldsymbol{\gamma}}=\Theta({\boldsymbol{H}}{\boldsymbol{\gamma}}+({\boldsymbol{I}}-{\boldsymbol{H}}){\boldsymbol{y}};\boldsymbol{\lambda})$ which can be directly computed by $\Theta$-IPOD. With a good design of $\Theta$, we successfully identified the outliers as well as estimating the coefficients robustly. This technique is associated with $M$-estimators, but gives a new characterization in form of penalized regressions. Furthermore, we successfully generalized this penalized/thresholding methodology to high-dimensional problems to accommodate and identify gross outliers in variable selection and coefficient estimation.

When outliers are also leverage points, the Gram matrix of the reduced model (4.9), i.e., ${\boldsymbol{A}}^{\mathsf{T}}{\boldsymbol{A}}={\boldsymbol{I}}-{\boldsymbol{H}}$, may demonstrate high correlation between clean observations and outliers for small $h_{i}$. It is well known from recent advances of the lasso (e.g. the irrepresentable conditions [Zhao] and the sparse Riesz condition [ZhangHuang]) that the convex $L_{1}$-penalty encounters great trouble in this situation. Rather, nonconvex penalties which correspond to redescending $\psi$’s must be applied, with a high breakdown point initial estimate obtained by, say, the fast-LTS [fastLTS] or fast-S [fastS] when $p$ is small, or the fast PY procedure [PYfast] when $p$ is large.

In this framework, determining the efficiency parameter in $M$-estimators and choosing a cutoff value for outlier identification are both accomplished by tuning the choice of $\boldsymbol{\lambda}$. Our experience shows that adopting an appropriate data-dependent choice of $\boldsymbol{\lambda}$ is crucial to guarantee good detection performance.

We close with a note of caution on an important issue for the outlier detection literature. We have found that our robust regression algorithms work better than competitors on our simulated data sets and that they give the right answers on some small well studied real data sets. But our method and the others we study all rely on a preliminary robust fit. The most robust preliminary fits have a cost that grows exponentially with dimension and for them $p=15$ is already large. In high dimensional problems, the preliminary fit of choice is seemingly the PY procedure. It performed well in our examples, but it does not necessarily have high breakdown. Should the preliminary method fail, followup methods like $\Theta$-IPOD or the others, may or may not correct it. We have seen hard-IPOD work well from a non-robust start (e.g. the HBK problem) but it would not be reasonable to expect this will always hold. We expect that $\Theta$-IPOD iterations will benefit from improvements in preliminary robust fitting methods for high dimensional problems. For low dimensional problems robust preliminary methods like LTS are fast enough.

Appendix A Proofs

Since both Huber’s method and the $L_{1}$-penalized regression are convex, it is sufficient to consider the KKT equations. In this proof, we define

The KKT equations for Huber’s estimate $(\hat{\boldsymbol{\beta}},\hat{\sigma})$ are given by

Let $\hat{\boldsymbol{r}}={\boldsymbol{y}}-{\boldsymbol{X}}\hat{\boldsymbol{\beta}},\hat{G}=\{i:|\hat{r}_{i}|<\lambda\hat{\sigma}\},\mbox{ and }\hat{O}=\{i:|\hat{r}_{i}|>\lambda\hat{\sigma}\}.$ Since

(A.2) becomes $nc=\lambda^{2}|\hat{O}|+\sum_{i\in\hat{G}}\hat{r}_{i}^{2}/\hat{\sigma}^{2}$. To summarize, $(\hat{\boldsymbol{\beta}},\hat{\sigma})$ satisfies

Next, the joint KKT equations for $L_{1}$-penalized regression estimates $(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\gamma}},\hat{\sigma})$ are

from which it follows that ${\boldsymbol{X}}^{\mathsf{T}}{\psi}\bigl(({{\boldsymbol{y}}-{\boldsymbol{X}}\hat{\boldsymbol{\beta}}})/{\hat{\sigma}};\lambda\bigr)=0$, and $\hat{\sigma}^{2}={\|\hat{\boldsymbol{r}}-\hat{\boldsymbol{\gamma}}\|_{2}^{2}}/{nc}$. But $\|\hat{\boldsymbol{r}}-\hat{\boldsymbol{\gamma}}\|_{2}^{2}=\sum_{i\in\hat{G}}\hat{r}_{i}^{2}+\lambda^{2}\hat{\sigma}^{2}|\hat{O}|$, so we obtain

which are exactly the same as (A.3) and (A.4). Therefore, (3.3) leads to Huber’s method with joint scale estimation. When $\sigma$ is fixed, it is not difficult to show that the two minimizations still yield the same ${\boldsymbol{\beta}}$-estimate.

Finally, Huber notices his method behaves poorly even for moderate leverage points and considers an improvement of using $c_{1}\psi({\hat{r}_{i}}/{(c_{2}\sigma)})$ to replace $\psi(\hat{r}_{i}/\sigma)$ [Huberbook]. He claims, based on heuristic arguments, that $c_{1}=c_{2}=\sqrt{1-h_{i}}$ is a good choice. This is quite natural as seen from our new characterization. The regression matrix for ${\boldsymbol{\gamma}}$ in (3.3) is actually ${\boldsymbol{I}}-{\boldsymbol{H}}$ with the column norms given by $\sqrt{1-h_{i}}$. The weighted $L_{1}$-penalty of $\sum(\lambda\sqrt{1-h_{i}})\times|\gamma_{i}|$ then leads to Huber’s improvement, due to the fact that $K\psi(t/K;\lambda)=\psi(t;K\lambda)$ for $K>0$. (However it cannot completely avoid masking and swamping unless a redescending $\psi$ is used.) ∎

By definition, for any $\Theta$-IPOD estimate $(\hat{\boldsymbol{\beta}},\hat{\boldsymbol{\gamma}})$, $\hat{\boldsymbol{\gamma}}$ is a fixed point of ${\boldsymbol{\gamma}}=\Theta({\boldsymbol{H}}{\boldsymbol{\gamma}}+({\boldsymbol{I}}-{\boldsymbol{H}}){\boldsymbol{y}};{\lambda})$, and $\hat{\boldsymbol{\beta}}=({\boldsymbol{X}}^{\mathsf{T}}{\boldsymbol{X}})^{-1}{\boldsymbol{X}}^{\mathsf{T}}({\boldsymbol{y}}-\hat{\boldsymbol{\gamma}})$. It follows that

and so $\hat{\boldsymbol{\beta}}$ is an $M$-estimate associated with $\psi$. ∎

The second inequality in (4.3) is straightforward from the algorithm design. To show the first inequality is true, it is sufficient to prove the following lemma.

Given a thresholding rule $\Theta$, let $P$ be any function satisfying $P(\theta;\lambda)=P(0;\lambda)+P_{\Theta}(\theta;\lambda)+q(\theta;\lambda)$ where $q(\cdot;\lambda)$ is nonnegative and $q(\Theta(\theta;\lambda))=0$ for all $\theta$. Then, the minimization problem $\min_{\theta}(t-\theta)^{2}/2+P(\theta;\lambda)$ has a unique optimal solution $\hat{\theta}=\Theta(t;\lambda)$ for every $t$ at which $\Theta(\cdot;\lambda)$ is continuous.

This is a generalization of Proposition 3.2 in ?. Note that $P$ (and $P_{\Theta}$) may not be differentiable at 0 and may not be convex.

Without loss of generality, suppose $t>0$. It suffices to consider $\theta\geq 0$ since $f(\theta)\geq f(-\theta),\forall\theta\geq 0$, where $f(\theta)\equiv(t-\theta)^{2}/2+P(\theta;\lambda)$. First, $\Theta^{-1}$ given by (4.1) is well-defined. We have

References