Robust Regression via Heuristic Hard Thresholding

1 Introduction

The presence of noises and corruptions in real-world data can be inevitably caused by the experimental errors, accidental outliers, or even adversarial data attacks. As the traditional least squares regression methods are vulnerable to outlier observations [Maronna et al., 2006], we study Robust Least Square Regression (RLSR) to handle the problem of learning a reliable set of regression coefficients given the presence of several adversarial corruptions in its response vector. A commonly adopted model from existing methods assumes that the observed response is obtained from the generative model y = XT β∗ +u, where β∗ is the true regression coefficients that we wish to recover and u ∈ Rn is the corruption vector with arbitrary values. Due to the ubiquitousness of data corruptions and popularity of regression methods, RLSR has become a critical component of several important realworld applications in various domains such as signal processing [Zoubir et al., 2012; Studer et al., 2012], economics [Rousseeuw and Leroy, 2005], bioinformatics [Lourenco et al., 2011] and image processing [Naseem et al., 2012; Wright et al., 2009]. A large body of literature on robust regression problem has been built up over the last few decades, but most of them [Smolic and Ohm, 2000; Huber and Ronchetti, 2009; She and Owen, 2011; Jung et al., 2016] lack the theoretical guarantee of regression coefficients recovery. To theoretically guarantee the recovery performance, [Chen et al., 2013] proposed a trimmed inner product based algorithm, but the recovery boundary of their method is not tight in a massive dataset. Also, [Mcwilliams et al., 2014] proposed a sub-sampling algorithm for large scale corrupted linear regression, the recovery result they provide is not close to a exact recovery [Bhatia et al., 2015]. To pursuit the exact recovery results for RLSR problem, some work focused on L1 penalty based convex formulations [Wright and Ma, 2010; Nguyen and Tran, 2013]. However, these methods imposed severe restrictions on the data distribution such as rowsampling from an incoherent orthogonal matrix [Nguyen and Tran, 2013].

Several studies requires the corruption ratio parameter which is difficult to be determined manually. For instance, instead of the exact corruption ratio, [Chen et al., 2013] requires the upper bound of outliers number which is also difficult to estimate. [She and Owen, 2011] relies on a regularization parameter to control the size of uncorrupted set based on soft-thresholding. Recently, [Bhatia et al., 2015] proposed a hard-thresholding algorithm for the RLSR problem: argmin ∑i∈S (yi - xTi β)2 with the constraint: |S| ≥ (1 - γ)n, where |S| is the size of the uncorrupted set and γ is the corruption ratio. Although the method guarantees an exact recovery of β under a mild assumption for covariate matrix, their results are highly dependent on the corruption ratio parameter γ inputed by users. Specifically, the parameter is required to be larger than the exact corruption ratio γ ∗ to ensure its convergence. Unfortunately, it is seldom practical to estimate the corruption ratio under the assumption that the response vector is arbitrarily corrupted. Furthermore, empirical results show that if the parameter γ is more than 50% off the true value, the recovery error can be more than double in size.

To address the lack of rigorous analysis on corruption ratio and theoretical guarantee on exact recovery, we proposed a new model, Robust Least squares regression algorithm via Heuristic Hard thresholding (RLHH). The main contributions of our study are summarized as follows: 1) The design of an efficient algorithm to address the RLSR problem without parameterizing its corruption. The algorithm RLHH is proposed to recover the regression coefficients and uncorrupted set efficiently. Unlike with a fixed corruption ratio, our method alternately estimates the optimal corruption ratio based on residual errors using optimized regression coefficients in each iteration. 2) An exact recovery guarantee under a mild assumption regarding input variables. We prove that our RLHH algorithm converges at a geometric rate and recovers β∗ exactly under the assumption that the least squares function satisfies both the Subset Strong Convexity (SSC) and Subset Strong Smoothness (SSS) properties. Specifically, we prove that our heuristic hard thresholding function ensures that the residual of the estimated uncorrupted set in each iteration has a tight upper error bound for the true uncorrupted set. 3) Empirical effectiveness and efficiency. Our proposed algorithm was evaluated with 6 competing methods in synthetic data. The results demonstrate that our approach consistently outperforms existing methods in both regression coefficients and uncorrupted set recovery, delivering a competitive running time.

The reminder of this paper is organized as follows. Section 2 gives a formal problem formulation. The proposed RLHH algorithm is presented in Section 3. Section 4 presents the proof for the recovery guarantees. In Section 5, the experimental results are analyzed and the paper concludes with a summary of our work in Section 6.

2 ProblemFormulation

In this study, we consider the problem of RLSR with adversarially corrupted data. Given a covariate matrix X = [x1, ..., xn], where each column xi ∈ Rp×1 and β∗ represents the ground truth coefficients of the regression model, we assume the corresponding response vector y ∈ Rn×1 is generated using the following model:

where y∗ = XT β∗ and u is the unbounded corruption vector introduced by an adversary. ε represents the additive dense noise, where εi ∼ N(0, σ2 ). The goal of our problem is to learn a new problem, which is to recover the regression coefficients β∗ and simultaneously determine the uncorrupted point set ˆ S. The problem is formally defined as follows:

Given a subset S ⊂ [n], yS restricts the row of y to indices in S and XS signifies that the columns of X are restricted to indices in S. Therefore, we have yS ∈ R|S|×1 and XS ∈ Rp×|S| . We use the notation S∗ = supp(u) to denote the set of uncorrupted points. Also, for any vector v ∈ Rn , the notation vS represents the |S|-dimensional vector containing the components in S. The function G(·) is to determine the size of set S according to the regression coefficients β, which is explained in Section 3.

To prove the theoretical recovery of regression coefficients, we require that the least squares function satisfies the Subset Strong Convexity (SSC) and Subset Strong Smoothness (SSS), which are defined as follows:

Definition 1. SSC and SSS Property. The least squares function f (β) = ‖yS -XT S β‖22 satisfies 2ζα-Subset Strong Convexity Property and 2κα-Subset Strong Smoothness if the following holds:

Equation (3) is equivalent to:

where λmin and λmax are denoted as the smallest and largest eigenvalues of matrix X, respectively.

The optimization problem in Equation (2) is non-convex (jointly in β and S) in general and existing methods cannot guarantee the exact recovery and efficient convergence rate.

3 The ProposedMethodology

In order to solve the problem in Equation (2) efficiently with the guarantee on the exact recovery of regression coefficients, we propose a novel heuristic hard thresholding based robust regression algorithm, RLHH. The algorithm iteratively optimizes the regression coefficients β and uncorrupted set S until convergence. The optimization of S is very challenging because it mounts to a non-convex discrete optimization problem. To handle it, we propose a heuristic corruption estimator to determine the size of set S and then apply the estimated uncorrupted size into heuristic hard thresholding method for the optimization of S elements.

Denoting residual vector r = y-XT β and rδ(k) be the kth elements of r in ascending order of magnitude, the heuristic estimator G(·) determines the size of optimal uncorrupted set τe by optimizing the following problem:

where the function L is defined as

The variable τo in the constraint is defined as follows:

where τ ′ = τ - ⌈n/2⌉ and Sτ ′ is the position set containing the smallest τ ′ elements in residual r. The constraint is imposed to avoid the case when τ is close to n, where the residual becomes so arbitrary that the denominator can become very large, making L much smaller than the value of the estimated threshold τe.

Applying the optimal uncorrupted set size generated by G(·), the heuristic hard thresholding is defined as follows:

Definition 2. Heuristic Hard Thresholding. Denoting δ-1 r (i) as the position of the ith element in residual vector r's ascending order of magnitude. The heuristic hard thresholding of r is defined as

Algorithm 1: RLHH ALGORITHM

Input: Corrupted training data {xi, yi}, i = 1...n, tolerance ϵ Output: solution ˆ β 1 S0 = [n], t ←0 2 repeat 3 βt+1 ← (XStXT St)-1XStySt 4 for i = 1..n do 5 rt+1 i ← |ySti -xTStiβ| 6 St+1 ←Hτ∗ (rt+1 ), where τ∗ is solved by Equation (4). 7 t ← t + 1 8 until ‖rt+1 St+1 - rtSt ‖2 < ϵn 9 return βt+1

The optimization of S is formulated as solving the Equation (7). The set returned by HG(β) will be used as the estimated uncorrupted set.

We will first present the reasoning behind our choice of function in this section and then show that our heuristic function can indeed ensure a rigorous recovery of regression coefficients β in Section 4. Basically, the function follows a natural intuition that data points with unbounded corruption always have a higher residual ri = yi -Xiβ in magnitude compared to uncorrupted data. Moreover, as this corruption is arbitrary and unbounded, when the residual vector r is sorted in ascending order, the slope of overall corrupted data is always much larger than the slope of the uncorrupted data. As Figure 1 shows, point p∗ has the minimum L value in the feasible domain. Therefore, we can estimate the corresponding threshold τ∗ of point p∗ as the optimal threshold. To avoid a zero value for the numerator of L(τ), we add 1 to all the values in the residual vector.

In Algorithm 1, an efficient robust regression algorithm, RLHH, based on heuristic hard thresholding is proposed to solve the RLSR problem. It follows an intuitive strategy of updating β to provide a better fit for the current estimated set S in Line 3, and updating the residual vector r in Line

Point

Figure 2: Residual r in ascending order for the 1 st (left) and 5 th (right) iterations.

4-5. It then estimates an active set S of uncorrupted points via heuristic hard thresholding in Line 6 based on the residual vector r = y-Xβ in the current iteration. The active set is initialized using the entire data samples in Line 1. The algorithm continues until the change in the residual vector falls within a small range. Figure 2 shows that the residual of the uncorrupted set in the 1 st and 5 th iteration, respectively. It intuitively explains the convergence progress of our algorithm: the optimization steps of β based on St makes the rSt smaller than its previous iteration, and it leads to smaller L values for items in St. Then these items in St have much higher possibility to kept in St+1 than items in [n] \ St. This progress continues until the active set is fixed.

4 Theoretical Recovery Analysis

For convenience, the convergence analyses for the case without dense noise will be presented, i.e. y = XT β + u. The convergence proof relies on the optimality of two steps carried out by the algorithm, the β optimization step that selects the best coefficients based on the uncorrupted set, and the heuristic hard threshold step that automatically discovers the best active set based on the current regression coefficients.

Lemma 1. For a given residual vector r ∈ Rn , let δ(k) be the k-th position of the ascending order in vector r, i.e. rδ(1) ≤ rδ(2) ≤ ... ≤ rδ(n). For any 1 ≤ τ1 < τ2 ≤ n, let S1 = {δ(i)|1 ≤ i ≤ τ1} and S2 = {δ(i)|1 ≤ i ≤ τ2}. We then have ‖rS1 ‖22 ≤ τ1 τ2 ‖rS2 ‖22 ≤ ‖rS2 ‖22 .

Proof. Let S3 = {δ(i) : τ1 + 1 ≤ i ≤ τ2}. Clearly , we have ‖rS2 ‖22 = ‖rS1 ‖22 + ‖rS3 ‖22 . Moreover, since each element in S3 is larger than any of the element in S1, we have ‖rS1 ‖22 ≤ ‖rS2 ‖22 + |S3| |S1| ‖rS1 ‖22 ≤ |S1| |S1|+|S3| ‖rS2 ‖22 = τ1 τ2 ‖rS2 ‖22 ≤ ‖rS2 ‖22 .

Lemma 2. Let St be the estimated uncorrupted set at the tth iteration. If τt ≥ τ∗ = γn, then |S∗ ∩ St| ≥ τt - n2 .

Lemma 3. Let τt be the estimated uncorrupted threshold at the t-th iteration. If τt > τ∗ = γn, then ‖rtSt ‖22 ≤ [1 + 128(1-γ) 2γ-1 ]‖rtS∗ ‖22 .

Proof. To simplify the notation, we will omit all the subscripts t that signify the t-th iteration in the explanation below and assumes the residual vector r is sorted in ascending order of magnitude. According to the optimization step in equation (4), we have the following properties:

According to the optimization step in equation (4), we have the following properties:

The inequality (a) follows τo ≥ τ/4 and inequality (b) follows the definition of τo in equation (6) and the fact that |S∗ ∩ St| ≥ τ ′ in Lemma 2. The inequality (c) follows |St \ S∗| ≤ (1-γ) · n and ‖rSt\S∗ ‖22 ≤ |St \S∗|r2 τ . Then we have

(11) The inequality (d) follows ‖rS∗ ‖22 = ‖rS∗\St ‖22 + ‖rS∗∩St ‖22 . The inequality (e) follows τ ′ = τt - n2 .

Theorem 4. Let X = [x1, ..., xn] ∈ Rp be the given data matrix and y = XT β∗ + u be the corrupted output with ‖u‖0 = γn. Let Σ0 be an invertible matrix such that ˜X = Σ-1/2 0 X, f(β) = ‖yS - ˜XSβ‖22 satisfies the SSC and SSS properties at level α, γ with 2ζα,γ and 2κα,γ If the data satisfies κγ ζ1-α < 1 √λ (√2-1), then after t = O (log 1η µ‖u‖2 ϵ ) iterations, Algorithm 1 yields an ϵ-accurate solution βt .

Proof. Let Gt = (XSt XT St )-1XSt , the t-th iteration of Algorithm 1 satisfies

Thus, the residual in t+1-th iteration for any set S ⊂ [n], yields

For each iteration, we have two conditions when choosing different values of τ t+1 . For condition 1, τ t+1 ≤ τ∗, we have ‖rt+1 St+1 ‖22 ≤ ‖rt+1 S∗ ‖22 (see Lemma 1).

where α1 = maxt{1 - τt n }. The inequality (a) follows ‖rt+1 St+1 ‖22 ≤ ‖rt+1 S∗ ‖22 , and inequality (b) follows from the setting ˜X = Σ-1/2 0 X, SSC/SSS properties, |St| ≤ (1-γ) · n and |S∗ \ St+1| ≤ α1 · n. Solving the quadratic equation for the corruption vector gives us

For condition 2, τ t+1 > τ∗. According to Lemma 3, ‖rtSt ‖22 ≤ λ‖rtS∗ ‖22 where λ = 1 + 128(1-γ) 2γ-1 , we have

where α2 = maxt{1 - τt n }. The inequality (c) follows Lemma 3, inequality (d) follows from the definition of SSC/SSS properties, |St| ≤ (1-α2)·n and |S∗\St+1| ≤ γ·n. Inequality (e) notices the fact that √λ ≥ 1. Solving the quadratic equation in Equation (14) gives us

Combine these two conditions and let t1 be the iterations for the case of condition 1. We get

5 Experimental Results

In this section, we report the extensive experimental evaluation carried out to verify the robustness and efficiency of the proposed method. All the experiments were conducted on a 64-bit machine with Intel(R) core(TM) quad-core processor (i7CPU@3.6GHz) and 32.0GB memory. Details of both the source code and sample data used in the experiment can be downloaded here 1 .

1 https://github.com/xuczhang/RLHH

[Wright and Ma, 2010] [Nguyen and Tran, 2013]. For extensive L1 minimization solvers, [Yang et al., 2010] showed that the Homotopy and DALM solvers outperform other proposed methods both in terms of recovery properties and running time. Both of the L1 solver methods are parameter free. Another recently proposed hard thresholding method, Torrent (Abbr. Torr), developed for robust regression [Bhatia et al., 2015] was also compared to our method. As the method requires a parameter for corruption ratio, which is difficult to estimate in practice, we chose 4 versions with different parameter settings: TORR*, TORR25, TORR50, and TORR80. TORR* uses the true corruption ratio as its parameter, and the others apply parameters that are uniformly distributed across the range of ±25%, ±50%, and ±80% off the true value, respectively. All the results will be averaged over 10 runs.

5.3 Recovery of regression coefficients

We selected 6 competing methods with which to evaluate the recovery performance of β: OLS, DALM, Homotopy, TORR*, TORR25, TORR50. As the recovery error for the OLS method is almost 10 times larger than those of the other methods, its result is not shown in Figure 3 in order to present the other results properly. Figures 3(a), 3(b), and 3(c) show the recovery performance for different data sizes when the feature number is fixed. Looking at the results, we can conclude that: 1) the RLHH method outperforms all the competing methods except for TORR*, whose parameter is hardly given in practice. 2) The results of the TORR methods are significantly affected by their corruption ratio parameters; TORR50 performs almost twice as badly as TORR* and yields worse results than one of the L1-Solver methods, DALM. However, RLHH performs consistently throughout, with no impact of the parame-

5.1 Datasets andMetrics

To demonstrate the performance of our proposed method, we carried out comprehensive experiments in synthetic datasets. Specifically, the simulation samples were randomly generated according to the model in Equation (1) for RLSR problem, sampling the regression coefficients β∗ ∈ Rp as a random unit norm vector. The covariance data X was drawn independently and identically distributed from xi ∼ N(0, Ip) and the uncorrupted response variables were generated as y∗ i = xTi β∗ . The set of corrupted points S was selected as a uniformly random (n-τ∗)-sized subset of [n], where τ∗ is the size of the uncorrupted set. The corrupted response vector was generated as y = y∗ + u + ε, where the corruption vector u was sampled from the uniform distribution [-5‖y∗ ‖∞, 5‖y∗ ‖∞] and the additive dense noise was εi ∼ N(0, σ2 ).

Following the setting in [Bhatia et al., 2015], we measured the performance of the regression coefficients recovery using the standard L2 error e = ‖ ˆ β - β∗ ‖2 , where ˆ β represents the recovered coefficients for each method and β∗ is the true regression coefficients. To validate the performance for corrupted set discovery, the F1-score is measured by comparing the discovered corrupted sets with the actual ones. To compare the scalability of each method, the CPU running time for each of the competing methods was also measured.

5.2 ComparisonMethods

The following methods are included in the performance comparison presented here: Ordinary least squares (OLS). The OLS method ignores the corruption of data and trains the model based on the whole dataset. We also compared our method to the regularized L1 algorithm for robust regression ter. 3) The L1-Solver methods generally exhibit worse performance than the hard thresholding based algorithms. Specifically, compared to DALM, Homotopy is more sensitive to the number of corrupted instances in the data. Figure 3(d) shows its similar performance when the feature number increases. Figures 3(e) and 3(f) show that RLHH performs equally as well as TORR* without dense noise, with both achieving an exact recovery of regression coefficients β.

5.4 Recovery of Uncorrupted Sets

As most competing methods do not explicitly estimate uncorrupted sets, we compared our proposed method with the TORR algorithm using a number of different parameter settings ranging from the true corrupted ratio up to a deviation of 30%. As the results shown in Table 1, we found that: 1) The F1 score of RLHH is 1.1% less than that of TORR* on average, although it is important to note that the latter uses the true corruption ratio, which cannot be estimated exactly in practice. This indicates that our method provides very close to an optimal estimation result for an uncorrupted set. 2) The RLHH method significantly outperforms the other methods, especially when the true corruption ratio is high. 3) The results of the TORR methods are highly dependent on the corruption ratio parameter: the results for a 25% corruption estimation error are much better than those for a 50% error. However, RLHH is a parameter free method that is capable of obtaining a good result consistently. 4) When increasing the feature number and corruption ratio, the F1 scores slightly

(a) p=400, n=4000, no dense

noise

Figure 4: Running time for different corruption ratios and data sizes

(b) p=100, cr=0.1, dense noise increase for all the methods. 5) In a no dense noise setting, RLHH performs a near optimal recovery result, while TORR* exactly recovers the result only because it is using the true corruption ratio.

5.5 Efficiency

To evaluate the efficiency of our proposed method, we compared the performances of all the competing methods for two different data settings: different corruption ratios and data sizes. In general, as Figure 4 shows, the hard thresholding based methods significantly outperformed the L1-Solver based methods. Also, the running time for RLHH increases slowly as either the corruption ratio or the data size increases, just as in the TORR methods. In addition, even though RLHH performs the additional step of estimating the uncorrupted set in each optimization iteration, the efficiency of RLHH still outperforms TORR, which indicates that 1) the heuristic hard thresholding step in RLHH always performs efficiently, even under in different circumstances; and 2) The RLHH algorithm converges more quickly than TORR.

6 Conclusion

In this paper, a novel robust regression algorithm, RLHH, is proposed to recover the regression coefficients and the uncorrupted set in the presence of adversarial corruption in the response vector. To achieve this, we designed a heuristic hard thresholding method with which to estimate the optimal uncorrupted set that is alternately updated with the optimized regression coefficients. We demonstrate that our algorithm can recover true regression coefficients exactly, with a geometric convergence rate. Extensive experiments on massive simulation data demonstrated that the proposed algorithm outperforms other comparable methods in both effectiveness and efficiency.