On Robust Estimation of High Dimensional Generalized Linear Models

1 Introduction

Statistical models in machine learning allows us to make strong predictions even from limited data, by leveraging specific assumptions imposed on the model space. However, on the flip side, when the specific model assumptions do not exactly hold, these standard methods could deteriorate severely. Constructing estimators that are robust to such departures from model assumptions is thus an important problem, and forms the main focus of Robust Statistics [Huber, 1981; Hampel et al., 1986; Maronna et al., 2006]. In this paper, we focus on the robust estimation of high-dimensional generalized linear models (GLMs). GLMs are a very general class of models for predicting a response given a covariate vector, and include many classical conditional distributions such as Gaussian, logistic, etc. In classical GLMs, the data points are typically low dimensional and are all assumed to be actually drawn from the assumed model. In our setting of high dimensional robust GLMs, there are two caveats: (a) the true parameter vector can be very high dimensional and furthermore, (b) certain observations are outliers, and could have arbitrary values with no quantitative relationship to the assumed generalized linear model.

Existing Research: Robust Statistics. There has been a long line of work [Huber, 1981; Hampel et al., 1986; Maronna et al., 2006] on robust statistical estimators. These are based

Pradeep Ravikumar

Department of Computer Science University of Texas, Austin pradeepr@cs.utexas.edu

on the insight that the typical log-likelihood losses, such as the squared loss for the ordinary least squares estimator for linear regression, are very sensitive to outliers, and that one could devise surrogate losses instead that are more resistant to such outliers. Rousseeuw [1984] for instance proposed the least median estimator as a robust variant of the ordinary least squares estimator. Another class of approaches fit trimmed estimators to the data after an initial removal of candidate outliers Rousseeuw and Leroy [1987]. There have also been estimators that model the outliers explicitly. [Gelman et al., 2003] for instance model the responses using a mixture of two Gaussian distributions: one for the regular noise, and the other for the outlier noise, typically modeled as a Gaussian with high variance. Another instance is where the outliers are modeled as being drawn from heavy-tailed distributions such as the t distribution [Lange et al., 1989].

Existing Research: Robustness and Sparsity. The past few years have actually led to an understanding that outlier robust estimation is intimately connected to sparse signal recovery [Candes and Tao, 2005; Antoniadis, 2007; Jin and Rao, 2010; Mitra et al., 2010; She and Owen, 2011]. The main insight here is that if the number of outliers is small, it could be cast as a sparse error vector that is added to the standard noise. The problem of sparse signal recovery itself has seen a surge of recent research, and where a large body of work has shown that convex and tractable methods employing the likes of ℓ1 regularization enjoy strong statistical guarantees [Donoho and Elad, 2003; Ng, 2004; Candes and Tao, 2006; Meinshausen and B¨uhlmann, 2006; Tropp, 2006; Zhao and Yu, 2007; Wainwright, 2009; Yang et al., 2012]. Intriguingly, Antoniadis [2007]; She and Owen [2011] show that even classical robust statistics methods could be cast as sparsity encouraging M-estimators that specifically use non-convex regularization. Jin and Rao [2010]; Mitra et al. [2010] have also suggested the use of non-convex penalization based methods such as SCAD [Fan and Li, 2001] for robust statistical estimation. Convex regularization based estimators however have been enormously successful in highdimensional statistical estimation, and in particular provide tractable methods that scale to very high-dimensional problems, and yet come with rigorous guarantees. To complete the story on the connection between robustness and sparsity, it is thus vital to obtain bounds on the performance of the convex regularization based estimators for general highdimensional robust estimation. For the task of high dimensional robust linear regression, there has been some interesting recent work [Nguyen and Tran, 2011] that have provided precisely such bounds. In this paper, we provide such an analysis for GLMs beyond the standard Gaussian linear model.

It turns out that the story for robust GLMs beyond the standard Gaussian linear model is more complicated. In particular, outlier modeling in GLMs could be done in two ways: (a) in the parameter space of the GLM, or (b) in the output space. For the linear model these two approaches are equivalent, but significant differences emerge in the general case. We show that the former approach always leads to convex optimization problems but only works under rather stringent conditions. On the other hand, the latter approach can lead to a non-convex M-estimator, but enjoys better guarantees. However, we show that all global minimizers of the M-estimation problem arising in the second approach are close to each other, so that the non-convexity in the problem is rather benign. Leveraging recent results [Agarwal et al., 2010; Loh and Wainwright, 2011], we can then show that projected gradient descent will approach one of the global minimizers up to an additive error that scales with the statistical precision of the problem. Our main contributions are thus as follows:

1. For robust estimation of GLMs, we show that there are two distinct ways to use the connection between robustness and sparsity.
2. For each of these two distinct approaches, we provide M-estimators, that use ℓ1 regularization, and in addition, appropriate constraints. For the first approach, the M-estimation problem is convex and tractable. For the second approach, the M-estimation problem is typically non-convex. But we provide a projected gradient descent algorithm that is guaranteed to converge to a global minimum of the corresponding M-estimation problem, up to an additive error that scales with the statistical precision of the problem.
3. One of the main contributions of the paper is to provide ℓ2 error bounds for each of the two M-estimators, for general values of the tuple (n, p, s, k). The analysis of corrupted general non-linear models in highdimensional regimes is highly non-trivial: it combines the twin difficulties of high-dimensional analysis of nonlinear models, and analysis given corrupted observations. The presence of both these elements, specifically the interactions therein, required a subtler analysis, as well as slightly modified M-estimators.

2 Problem Statement and Setup

We consider generalized linear models (GLMs) where the response variable has an exponential family distribution, conditioned on the covariate vector,

Examples. The standard linear model with Gaussian noise, the logistic regression and the Poisson model are typical examples of this model. In case of standard linear model, the domain of variable y, Y, is the set of real numbers, R, and with known scale parameter σ, the probability of y in (1) can be rewritten as

where the normalization function A(a) in (1) in this case becomes a2 /2. Another very popular example in GLM models is logistic regression given a categorical output variable:

where Y is {0, 1}, and the normalization function A(a) = log(1 + exp(a)). We can also derive the Poisson regression model from (1) as follows:

where Y is {0, 1, 2, ...}, and the normalization function A(a) = exp(a). Our final example is the case where the variable y follows an exponential distribution:

where Y is the set of non-negative real numbers, and the normalization function A(a) = -log(-a). Any distribution in the exponential family can be written as the GLM form (1) where the canonical parameter of exponential family is 〈θ⋆ , x〉. Note however that some distributions such as Poisson or exponential place restrictions on 〈θ⋆ , x〉 to be valid parameter, so that the density is normalizable, or equivalently the normalization function A(〈θ⋆ , x〉) < +∞.

In the GLM setting, suppose that we are given n covariate vectors, xi ∈ Rp , drawn i.i.d. from some distribution, and corresponding response variables, yi ∈ Y, drawn from the distribution P(y|xi, θ⋆ ) in (1). A key goal in statistical estimation is to estimate the parameters θ∗ , given just the samples Z n 1 := {(xi, yi)}ni=1 . Such estimation becomes particularly challenging in a high-dimensional regime, where the dimension p is potentially even larger than the number of samples n. In this paper, we are interested in such high dimensional parameter estimation of a GLM under the additional caveat that some of the observations yi are arbitrarily corrupted. We can model such corruptions by adding an "outlier error" parameter e⋆i in two ways: (i) we consider e⋆i in the "parameter space" to the uncorrupted parameter 〈θ⋆ , xi〉, or (ii) introduce e⋆i in the output space, so that the output yi is actually the sum of e⋆i and the uncorrupted output ¯ yi. For the specific case of the linear model, both these approaches are exactly the same. We assume that only some of the examples are corrupted, which translates to the error vector e⋆ ∈ Rn being sparse. We further assume that the parameter θ⋆ is also sparse. We thus assume:

with support sets S and T, respectively. We detail the two approaches (modeling outlier errors in the parameter space and output space respectively) in the next two sections.

3 Modeling gross errors in the parameter space

In this section, we discuss a robust estimation approach by modeling gross outlier errors in the parameter space. Specifically, we assume that the i-th response yi is drawn from the conditional distribution in

(1) but with a "corrupted" parameter 〈θ⋆ , xi〉 + √ne⋆i , so that the samples are distributed as P(yi|xi, θ⋆ , e⋆i ) = exp { h(yi)+yi(〈θ⋆ ,xi〉+√ne⋆i )-A(〈θ⋆ ,xi〉+√ne⋆i ) c(σ) }. We can

then write down the resulting negative log-likelihood as,

We thus arrive at the following ℓ1 regularized maximum likelihood estimator:

In the sequel, we will consider the following constrained version of the MLE (a0, b0 are constants independent of n, p):

The additional regularization provided by the constraints allow us to obtain tighter bounds for the resulting M-estimator.

We note that the M-estimation problem in (6) is convex: adding the outlier variables e does not destroy the convexity of the original problem with no gross errors. On the other hand, as we detail below, extremely stringent conditions are required for consistent statistical estimation. (The strictness of these conditions is what motivated us to also consider output space gross errors in the next section, where the conditions required are more benign).

3.1 ℓ2 Error Bound

We require the following stringent condition:

We assume the covariates are multivariate Gaussian distributed as described in the following condition:

Assumption 2. Let X be the n × p design matrix, with the n samples {xi} along the n rows. We assume that each sample xi is independently drawn from N(0,Σ). Let λmax and λmin > 0 be the maximum and minimum eigenvalues of the covariance matrix Σ, respectively, and let ξ be the maximum diagonal entry of Σ. We assume that ξλmax = Θ(1).

Additionally, we place a mild restriction on the normalization function A(·) that all examples of GLMs in Section 2 satisfy:

Assumption 3. The double-derivative A′′ (·) of the normalization function has at most exponential decay: A′′ (η) ≥ exp(-cη) for some c > 0.

Theorem 1. Consider the optimal solution (̂ θ, ̂ e) of (6) with the regularization parameters:

where c1 and c2 are some known constants. Then, there exist positive constants K, c3 and c4 such that with probability at least 1- Kn , the error (̂∆, ̂Γ) := (̂ θ - θ⋆ , ̂ e-e⋆ ) is bounded by

Note that the theorem requires the assumption that the outlier errors are bounded as ‖e⋆‖2 ≤ b0 √n . Since the "corrupted" parameters are given by 〈θ⋆ , xi〉+ √ne⋆i , the gross outlier error scales as √n‖e⋆‖2, which the assumption thus entails be bounded by a constant (independent of n). Our search to find a method that can tolerate larger gross errors led us to introduce the error in the output space in the next section.

4 Modeling gross errors in the output space

In this section, we investigate the consequences of modeling the gross outlier errors directly in the response space. Specifically, we assume that a perturbation of the i-th response, yi - √ne⋆i is drawn from the conditional distribution in (1) with parameter 〈θ⋆ , xi〉, so that the samples are distributed as P(yi|xi, θ⋆ , e⋆i ) =

then write down the resulting likelihood as, Lo(θ, e;Z n 1 ) := 1n ∑ni=1 [B (yi - √nei) - (yi - √nei)〈θ, xi〉 + A(〈θ, xi〉)], where B(y) = -h(y), and the resulting ℓ1 regularized maximum likelihood estimator as:

Note that when B(y) is set to -h(y) as above, the estimator has the natural interpretation of maximizing regularized loglikelihood, but in the sequel we allow it to be an arbitrary function taking the response variable as an input argument. As we will see, setting this to a function other than -h(y) will allow us to obtain stronger statistical guarantees.

Just as in the previous section, we consider a constrained version of the MLE in the sequel:

A key reason we introduce these constraints will be seen in the next section: these constraints help in designing an efficient iterative optimization algorithm to solve the above optimization problem (by providing bounds on the iterates right from the first iteration). One unfortunate facet of the Mestimation problem in (8), and the allied problem in (7), is that it is not convex in general. We will nonetheless show that the computationally tractable algorithm we provide in the next section is guaranteed to converge to a global optimum (up to an additive error that scales at most with the statistical error of the global optimum).

We require the following bounds on the ℓ2 norms of θ⋆ , e⋆ .

Assumption 4. ‖θ⋆‖2 ≤ a0 and ‖e⋆‖2 ≤ b0 for some constants a0, b0.

When compared with Assumption 1 in the previous section, the above assumption imposes a much weaker restriction on the magnitude of the gross errors. Specifically, with the √n scaling included, the above Assumption 4 allows the ℓ2 norm of the gross error to scale as √n, whereas Assumption 1 in the previous section required the norm ‖θ⋆‖2 to be bounded above by a constant.

4.1 ℓ2 Error bound

It turned out, given our analysis, that a natural selection of the function B(·) is to use the quadratic function (we defer discussion due to lack of space). Thus, in the spirit of classical robust statistics, we considered the modified log-likelihood objective in (7) with the above setting of B(·): Lo(θ, e; Z n 1 ) := 1n ∑ni=1 [ 12 (yi - √nei)2 - (yi - √nei)〈θ, xi〉 + A(〈θ, xi〉)].

Similarly here, we assume the random design matrix has rows sampled from a sub-Gaussian distribution:

Assumption 5. Let X be the n× p design matrix which has each sample xi in its ith row. Let λmax and λmin > 0 be the maximum and minimum eigenvalues of the covariance matrix of x, respectively. For any v ∈ Rp , the variable 〈θ, xi〉 is sub-Gaussian with parameter at most κ2u ‖v‖22 .

Theorem 2. Consider the optimal solution (̂ θ, ̂ e) of (8) with the regularization parameters:

where c′ , c′′ , c1, c2, c3 are some known constants. Then, there exist positive constants K, L and c4 such that if n ≥ Lmax(s, k) log p, then with probability at least 1 - nK , the error (̂∆, ̂Γ) := (̂ θ - θ⋆ , ̂ e - e⋆ ) is bounded by

Remarks. Nguyen and Tran [2011] analyze the specific case of the standard linear regression model (which nonetheless is a member of the GLM family), and provide the bound:

which is asymptotically equivalent to the bound in Theorem 2. As we noted earlier, for the linear regression model, both approaches of modeling outlier errors in the parameter space or the output space are equivalent, so that we could also compare the linear regression bound to our bound in Theorem 1.There too, the bounds can be seen to be asymptotically equivalent. We thus see that the generality of the GLM family does not adversely affect ℓ2 norm convergence rates even when restricted to the simple linear regression model.

5 A Tractable OptimizationMethod for the Output Space Modeling Approach

In this section we focus on the M-estimation problem (8) that arises in the second approach where we model errors in the output space. Unfortunately, this is not a tractable optimization problem: in particular, the presence of the bilinear term ei〈θ, xi〉 makes the objective function Lo non-convex. A tractable seemingly-approximate method would be to solve for a local minimum of the objective, by using a gradient descent based method. In particular, projected gradient descent (PGD) applied to the M-estimation problem (8) produces the iterates:

where η > 0 is a step-size parameter. Note that even though Lo is non-convex, the problem above is convex, and decouples in θ, e. Moreover, minimizing a composite objective over the ℓ1 ball can be solved very efficiently by performing two projections onto the ℓ1-ball (see Agarwal et al. [2010] for instance for details).

While the projected gradient descent algorithm with the iterates above might be tractable, one concern might be that these iterates would atmost converge to a local minimum, which might not satisfy the consistency and ℓ2 convergence rates as outlined in Theorem 2. However, the following theorem shows that the concern is unwarranted: the iterates converge to a global minimum of the optimization problem in (8), up to an additive error that scales at most as the statistical error, (‖̂ θ - θ⋆‖22 + ‖̂ e - e⋆ ‖22).

Theorem 3. Suppose all conditions of Theorem 2 hold and that n > c20 (k + s)2 log(p). Let F(θ, e) denote the objective function in (8) and let (̂ θ, ̂ e) be a global optimum of the problem. Then, when we apply the PGD steps above with appropriate step-size η, there exist universal constants C1,C2 > 0 and a contraction coefficient γ < 1, independent of (n, p, s, k), such that ‖θt - ̂ θ‖22 + ‖et - ̂ e‖22 ≤ C1 (‖̂ θ - θ⋆‖22 + ‖̂ e - e⋆‖22) for all iterates t ≥ T where

6 Experimental Results

In this section, we provide experimental validation, over both simulated as well as real data, of the performance of our Mestimators.

6.1 Simulation Studies

In this section, we provide simulations corroborating Theorems 1 and 2. The theorems are applicable to any distribution in the GLM family (1), and as canonical instances, we consider the cases of logistic regression (3), Poisson regression (4), and exponential regression (5). (The case of the standard linear regression model under gross errors has been previously considered in Nguyen and Tran [2011].)

We instantiated our models as follows. We first randomly selected a subset S of {1, . . . , p} of size √p as the support set (indexing non-zero values) of the true parameter θ∗ . We then set the nonzero elements, θ∗S , to be equal to ω, which we vary as noted in the plots. We then randomly generated n i.i.d samples, {x1, ..., xn}, from the normal distribution N(0, σ2 Ip×p). Given each feature vector xi, we drew the corresponding true class label ¯ yi from the corresponding GLM distribution. To simulate the worst instance of gross errors, we selected the k samples with the highest value of 〈θ∗ , xi〉 and corrupted them as follows. For logistic regression, we just flipped their class labels, to yi = (1 - ¯ yi). For the Poisson and exponential regression models, the corrupted response yi is obtained by adding a gross error term δi to ¯ yi. The learning algorithms were then given the corrupted dataset {xi, yi}ni=1 . We scaled the number of corrupted samples k with the total number of samples n in three different ways: logarithmic scaling with k = Θ(logn), square root scaling with k = Θ(√n), and linear scaling with k = Θ(n). For each tuple of (n, p, s, k), we drew 50 batches of n samples, and plot the average.

Figure 1 plots the ℓ2 norm error ‖ˆ θ - θ∗‖2 of the parameter estimates, against the number of samples n. We compare three methods: (a) the standard ℓ1 penalized GLMMLE (e.g. ℓ1 logistic reg.), which directly learns a GLM regression model over the corrupted data; (b) our first M-estimator (6), which models error in the parameter space (Gross error in param); and (c) our second M-estimator, which models error in the output space (8) (Gross error in output). As a gold standard, we also include the performance of the standard ℓ1 penalized GLM regression over the uncorrupted version of the dataset, {xi, ¯ yi}ni=1 (w/o gross error). Note that the ℓ2 norm error is just on the parameter estimates, and we exclude the error in estimating the outliers e themselves, so that we could compare against the gold-standard GLM regression on the uncorrupted data.

While the M-estimation problem with gross errors in the output space is not convex, it can be seen that the proximal gradient descent (PGD) iterates converge to the true θ∗ , corroborating Theorem 3. In the figure, the three rows correspond to the three different GLMs, and the three columns correspond to different outlier scalings, with logarithmic (first column), square-root (second column), and linear (third column) scalings of the number of outliers k as a function of the number of samples n. As the figure shows, the approaches modeling the outliers in the output and parameter spaces perform overwhelmingly better than the baseline ℓ1 penalized GLM regression estimator, and their error even approaches the estimator that is trained from uncorrupted data, even under settings where the number of outliers is a linear fraction of the number of samples. The approach modeling outliers in the output space seems preferable in some cases (logistic, exponential), while the approach modeling outliers in the parameter space seems preferable in some cases (Poisson).

6.2 Real Data Examples

In this section, we evaluate the performance of our estimators on some real binary classification datasets, obtained from LIBSVM (http://www.csie.ntu.edu.tw/∼cjlin/libsvmtools/ datasets/). We focused on the logistic regression case, and compared our two proposed approaches against the standard logistic regression. Note that the datasets we consider have p < n so that we can set λn,θ = 0, thus not adding further sparsity encouraging ℓ1 regularization to the parameters. We created variants of the datasets by adding varying number of outliers (by randomly flipped the values of the responses y as in the experiments on the simulated datasets). Given each dataset, we split it into three groups; 0.2 of training dataset, 0.4 of validation dataset and 0.4 of test dataset. The validation dataset is used to choose the best performance of λn,e and the constraint constant ρ where we solved the optimization problem under the constraint ‖e‖ ≤ ρ.

Figure 2 plots performance comparisons on 3 datasets, one row for each dataset. We varied the fraction of training dataset provided to each algorithm, and columns correspond to these varying fractions: 10% (Left column), 50% (Center) and 100% (Right column). Each graph has a group of four bar-plots corresponding to the four different types of outliers: original dataset without adding artificial outliers (w/o), and where the number of outliers scales as log(n) (Log), √n (Sqrt) or 0.1(n) (Linear), given n training examples. Our proposed robust methods perform as well or better, with particularly strong performance, with more outliers, and/or where less samples are used for the training. We found the latter phenomenon interesting, and worthy of further research: that robustness might help the performance of regression models even in the absence of outliers by preventing overfitting.

7 Conclusion

We have provided a comprehensive analysis of statistical estimation of high dimensional GLMs with grossly corrupted observations. We detail two distinct approaches for modeling sparse outlier errors in GLMs: incidentally these are equivalent in the linear case, though distinct for general GLMs. For both approaches, we provide tractable M-estimators, and analyze their consistency by providing ℓ2 error bounds. The parameter space approach is nominally more intuitive and computationally tractable, but requires stronger conditions for the error bounds to hold (and in turn for ℓ2 consistency). In contrast, the second output space based approach leads to a nonconvex problem, which makes statistical and computational analyses harder. Nonetheless, we show that this second approach is better than the first on the statistical front, since we obtain better bounds that require weaker conditions to hold, and on the computational front it is comparable, as we show a simple projected gradient descent algorithm converges to one of the global optima up to statistical precision.

Acknowledgments. We acknowledge the support of NSF via IIS-1149803, and DoD via W911NF-12-1-0390.