Robust Estimation via Robust Gradient Estimation

Introduction

In this paper, we present a general class of estimators that are computationally tractable, and have strong robustness guarantees. The estimators we propose are obtained by robustifying iterative updates of risk minimization, and are broadly applicable to a wide-range of parametric statistical models. In the risk minimization framework, the target parameter $\theta^{*}$ is defined as the solution to an optimization problem:

where $\overline{\mathcal{L}}$ is an appropriate loss-function, $R$ is the population risk and $\Theta$ is the set of feasible parameters. The statistical inference problem within the risk minimization framework is then to compute an approximate minimizer to the above program when given access to samples $\mathcal{D}_{n}=\{z_{1},\ldots,z_{n}\}$. A classical approach to do so is via empirical risk minimization (ERM), where we substitute the empirical expectation given the samples for the population expectation in the specification of the risk objective. While most modern statistical estimators use the above empirical risk minimization framework, a standard assumption that is imposed on $\mathcal{D}_{n}$ is that the data has no outliers, and has no arbitrary deviations from model assumptions; i.e., it is typically assumed that each of the $z_{i}$’s are independent and identically distributed according to the distribution $P$. Moreover, many analyses of risk minimization further assume that $P$ follows a sub-gaussian distribution, or has otherwise well-controlled tails in order to appropriately control the deviation between the population risk and its empirical counterpart. Due in part to these caveats with ERM, the seminal work of $M$-estimation replaces the risk minimization objective with a robust counterpart, so that the minimizer of the empirical expectation of the robust counterpart is more robust than the ERM minimizer. As noted above, for strong statistical guarantees, these in turn require solving computationally intractable non-convex optimization programs.

In contrast to this classical work, we propose a class of estimators that have a shift in perspective: rather than specify a robust objective, we consider an algorithm, namely projected gradient descent, that directly optimizes the population risk objective in Eq. (1), and focus on making this algorithm robust. Thus, in contrast to specifying the robust parameter estimate as the solution to an optimization program as in $M$-estimation, which in turn could be computationally intractable, we specify the robust parameter estimate as the limit of a sequence of iterative updates that are individually robust as well as computationally tractable. We find that this shift in perspective leads to estimators that are both computationally tractable as well with strong robustness guarantees, that are as broadly applicable as ERM or $M$-estimators, and moreover with a unified statistical treatment for varied statistical models.

In addition to being applicable to a variety of statistical models, our general results are also applicable to a variety of notions of robustness. In this paper, we derive corollaries in particular for two canonical robustness settings:

Robustness to arbitrary outliers: In this setting, we focus on Huber’s $\epsilon$-contamination model, where rather than observe samples directly from $P$ in (1) we instead observe samples drawn from $P_{\epsilon}$ which for an arbitrary distribution $Q$ is defined as:

The distribution $Q$ allows for arbitrary outliers, which may correspond to gross corruptions or more subtle deviations from the assumed model. This model can be equivalently viewed as model mis-specfication in the Total Variation (TV) metric.

Robustness to heavy-tails: In this setting, we are interested in developing estimators under weak moment assumptions. We assume that the distribution $P$ from which we obtain samples only has finite low-order moments (see Section 5.3 for a precise characterization). Such heavy tailed distributions arise frequently in the analysis of financial data and large-scale biological datasets (see for instance examples in ). In contrast to classical analyses of empirical risk minimization , in this setting the empirical risk is not uniformly close to the population risk, and methods that directly minimize the empirical risk perform poorly (see Section 4).

While we provide corollaries demonstrating robustness with respect to the above deviations, we emphasize that our framework is more general. Below, we provide an outline of our results and contributions.

Estimators. Our first contribution is to introduce a new class of robust estimators for risk minimization (1). These estimators are based on robustly estimating gradients of the population risk to then plug in to a projected gradient descent algorithm, and are computationally tractable by design. A crucial ingredient of our framework is the design of robust gradient estimators for the population risk in (1). Our main insight is that in this general risk minimization setting, the gradient of the population risk is simply a multivariate mean vector, and we can leverage prior work on mean estimation to design robust gradient estimators. Thus, for our two canonical robustness cases, we develop such robust gradient estimators building on prior work for robust mean estimation in the Huber model , and in the heavy-tailed model . Another perspective of our framework is that it significantly generalizes the applicability of mean estimation methods to general parametric models.

Empirical Investigations. Our estimators are computationally practical, and accordingly, our second contribution is to conduct extensive numerical experiments on real and simulated data with our proposed class of estimators. We provide guidelines for tuning parameter selection, and compare the proposed estimators with several competitive baselines . We find that our estimators consistently perform well across different settings, and across various metrics.

Statistical Guarantees. Finally, we provide rigorous robustness guarantees for the estimators we propose for a variety of classical statistical models: linear regression, logistic regression, and exponential family models. Our contributions in this direction are two-fold: building on prior work we provide a general result on the stability of gradient descent for risk minimization, and show that under certain conditions, gradient descent can be quite tolerant to inaccurate gradient estimates. Subsequently, in concrete settings, we provide a careful analysis of the quality of gradient estimation afforded by our proposed gradient estimators, and combine these results to obtain guarantees on our final parameter estimates.

Broadly, as we discuss in the sequel, our work suggests that our class of estimators based on robust gradient estimation offer a variety of practical, conceptual, statistical and computational advantages for robust estimation. They provide the general applicability of classical $M$-estimators, together with computational practicality even for large-scale models, as well as strong robustness guarantees.

There has been extensive work in the broad area of robust statistics (see for instance and references therein); we focus this section on some lines of work that are most related to this paper. For the robustness setting of $\epsilon$-contaminated models, several classical estimators have been developed that are optimally robust for a variety of inferential tasks, including hypothesis testing , mean estimation , general parametric estimation , and non-parametric estimation . However, a major drawback with this classical line of work has been that most of the estimators with strong robustness guarantees are computationally intractable , while the remaining ones use heuristics and are consequently not optimal . A complementary line of recent research has focused on providing minimax upper and lower bounds on the performance of estimators under $\epsilon$-contamination model, without the constraint of computational tractability. Recently, there has been a flurry of research in theoretical computer science designing provably robust estimators which are computationally tractable while achieving near-optimal contamination dependence, for special classes of problems such as computing means and covariances. Some of the proposed algorithms are however not computationally practical as they rely on the ellipsoid algorithm or require solving semi-definite programs which can be slow for modern problem sizes.

While in the general $\epsilon$-contamination setting, the contamination distribution could be arbitrary, there has been a lot of work in settings where the contamination distribution is restricted in various ways. For example, recent work in high-dimensional statistics (for instance ) have studied problems like principal component analysis and linear regression under the assumption that the corruptions are evenly spread throughout the dataset.

For the robustness setting of heavy tailed distributions, robust estimators aim to relax the sub-gaussian or sub-exponential distributional assumptions that are typically imposed on the target distribution, and allow it to be a heavy tailed distribution. Most approaches in this category substitute the empirical mean of the risk objective in risk minimization with robust mean estimators such as that exhibit sub-gaussian type concentration around the true mean for distributions satisfying mild moment assumptions. The median-of-means estimator and Catoni’s mean estimator are two popular examples of such robust mean estimators. In particular, Hsu and Sabato use the median-of-means estimator to solve the corresponding robust variant of ERM. Although this estimator has strong theoretical guarantees, and is computationally tractable, as noted by the authors in it performs poorly in practice. In recent work Brownlees et al. use the Catoni’s mean estimator to solve the corresponding robust variant of ERM. The authors provide risk bounds similar to bounds one can achieve under sub-gaussian distributional assumptions. However, their estimator is not easily computable and the authors do not provide a practical algorithm to compute the estimator. Other recent work by Lerasle and Oliveira , Lugosi and Mendelson use similar ideas to derive estimators that perform well theoretically, in heavy-tailed situations. However, these approaches involve optimization of complex objectives for which no computationally tractable algorithms exist. We emphasize that in contrast to our work, these works focus on robustly estimating the population risk which does not directly lead to a computable estimator. In contrast, we consider robustly estimating the gradient of the population risk, and embedding these estimates within the iterative algorithm of projected gradient descent, which leads naturally to a computionally practical estimator.

2 Outline

We conclude this section with a brief outline of the remainder of the paper. In Section 2, we provide some background on risk minimization and the running robustness settings of Huber contamination, and heavy-tailed noise models. In Section 3, we introduce our class of robust estimators, and provide concrete algorithms for the $\epsilon$-contaminated and heavy-tailed settings. In Section 4 we study the empirical performance of our estimator on a variety of tasks and datasets. We complement our empirical results with theoretical guarantees in Sections 5, 6 and 7. We defer technical details to the Appendix. Finally, we conclude in Section 8 with a discussion of some open problems.

Background and Problem Setup

In this section we provide the necessary background on risk minimization, gradient descent, and introduce two notions of robustness that we consider in this work.

The goal of risk minimization is to minimize the population risk $\mathcal{R}(\theta)$, given only $n$ samples $\mathcal{D}_{n}=\{z_{i}\}_{i=1}^{n}$, in order to estimate the unknown parameter $\theta^{*}$.

In this work we assume that the population risk is convex to ensure tractable minimization. Moreover, in order to ensure identifiability of the parameter $\theta^{*}$, we impose two standard regularity conditions on the population risk. These properties are defined in terms of the error of the first-order Taylor approximation of the population risk, i.e. defining, $\tau(\theta_{1},\theta_{2}):=\mathcal{R}(\theta_{1})-\mathcal{R}(\theta_{2})-\langle{\nabla\mathcal{R}(\theta_{2})},\;{\theta_{1}-\theta_{2}}\rangle$, we assume that

2 Illustrative Examples of Risk Minimization

The framework of risk minimization is a central paradigm of statistical estimation and is widely applicable. In this section, we provide illustrative examples that fall under this framework.

For this setting we use the squared loss as our loss function, which induces the following population risk:

2.2 Generalized Linear Models

Once again, the true parameter $\theta^{*}$ is the minimizer of the resulting population risk $\mathcal{R}(\theta)$. It is easy to see that Linear Regression with Gaussian Noise lies in the family of generalized linear models. A popular instance of such GLMs is a logistic regression model.

In this case the $(x,y)$ pairs are linked as:

This corresponds to setting $\Phi(t)=\log(1+\exp(t))$ and $c(\gamma)=1$ in (5). The hessian of the population risk is given by

Note that as $\theta$ diverges, the minimum eigenvalue of the hessian approaches and the loss is no longer strongly convex. To prevent this, in this case we take the parameter space $\Theta$ to be bounded.

2.3 Exponential Families and Canonical Parameters.

In more details, we can write the true distribution $P$ in this case as

where $h(z)$ is some base measure. The negative log-likelihood gives us the following loss function:

3 Empirical Risk Minimization

Given data $\mathcal{D}_{n}=\{z_{i}\}_{i=1}^{n}$, empirical risk minimization (ERM) substitutes the empirical expectation of the risk for the population risk in the risk minimization objective:

Most modern statistical estimators follow this ERM recipe above. When the loss is the log-likelihood of the statistical model, this reduces to the classical Maximum Likelihood Estimation (MLE) principle. The empirical risk minimizer is however a poor estimator of $\theta^{*}$ in the presence of outliers in the data: since ERM depends on the sample mean, outliers in the data can effect the sample mean and lead ERM to sub-optimal estimates. This observation has led to a large body of research that focuses on developing robust M-estimators, where we substitute in the empirical expectation of a robust counterpart of the loss function $\overline{\mathcal{L}}$; the resulting estimators have favorable statistical properties, but are often computationally intractable.

4 Projected Gradient Descent

A popular approach for solving the empirical risk minimization problem is projected gradient descent. Projected gradient descent generates a sequence of iterates $\{\theta^{t}\}_{t=0}^{\infty}$, by refining an initial parameter $\theta_{0}\in\Theta$ via the update:

5 Robust Estimation

One of the goals of this work is to develop general statistical estimation methods that are robust under varied robustness settings. We derive corollaries in particular for two robustness models: Huber’s $\epsilon$-contamination model, and the heavy-tailed model. We now briefly review these two notions of robustness.

Huber’s $\epsilon$-contamination model: Huber proposed the $\epsilon$-contamination model where we observe samples that are obtained from a mixture of the form

Heavy-tailed model: In the heavy-tailed model it is assumed that the data follows a heavy-tailed distribution (i.e, $P$ is heavy-tailed). While heavy-tailed distributions have various possible characterizations: in this paper we consider a characterization via gradients. For a fixed $\theta\in\Theta$ we let $P_{g}^{\theta}$ denote the multivariate distribution of the gradient of population loss, i.e. $\nabla\overline{\mathcal{L}}(\theta;z).$ We refer to a potentially heavy-tailed distribution as one for which our only assumption on $P_{g}^{\theta}$ is that it has finite second moments for any $\theta\in\Theta$. As we illustrate in Section 7, in various concrete examples this translates to relatively weak low-order moment assumptions on the data distribution $P$.

Given $n$ i.i.d observations from $P$, our objective is to estimate the minimizer of the population risk. From a conceptual standpoint, the classical analysis of risk-minimization which relies on uniform concentration of the empirical risk around the true risk, fails in the heavy-tailed setting necessitating new estimators and analyses .

Robust Gradient Descent via Gradient Estimation

Gradient descent and its variants are at the heart of modern optimization and are well-studied in the literature. Suppose we have access to the true distribution $P_{\theta^{*}}$. Then to minimize the population risk $\mathcal{R}(\theta)$, we can use projected gradient descent, where starting at some initial $\theta^{0}$ and for an appropriately chosen step-size $\eta$, we update our estimate according to:

However, we only have access to $n$ samples $\mathcal{D}_{n}=\{z_{i}\}_{i=1}^{n}$. The key technical challenges are then to estimate the gradient of $\mathcal{R}(\theta)$ from samples $\mathcal{D}_{n}$, and to ensure that an appropriate modification of gradient descent is stable to the resulting estimation error.

To address the first challenge we observe that the gradient of the population risk at any point $\theta$ is the mean of a multivariate distribution, i.e. $\nabla\mathcal{R}(\theta)=E_{z\sim P}\left[\nabla\overline{\mathcal{L}}(\theta;z)\right]$. Accordingly, the problem of gradient estimation can be reduced to a multivariate mean estimation problem, where our goal is to robustly estimate the true mean $\nabla\mathcal{R}(\theta)$ from $n$ samples $\{\nabla\overline{\mathcal{L}}(\theta;z_{i})\}_{i=1}^{n}$. For a given sample-size $n$ and confidence parameter $\delta\in(0,1)$ we define a gradient estimator:

A function $g(\theta;\mathcal{D}_{n},\delta)$ is a gradient estimator, if for functions $\alpha$ and $\beta$, with probability at least $1-\delta$, at any fixed $\theta\in\Theta$, the estimator satisfies the following inequality:

In subsequent sections, we will develop conditions under which we can obtain gradient estimators with strong control on the functions $\alpha(n,\delta)$ and $\beta(n,\delta)$ in the Huber and heavy-tailed models. Furthermore, by investigating the stability of gradient descent we will develop sufficient conditions on these functions such that gradient descent with an inaccurate gradient estimator still returns an accurate estimate.

To minimize $\mathcal{R}(\theta)$, we replace $\nabla\mathcal{R}(\theta)$ in equation (10) with the gradient estimator $g(\theta;\mathcal{D}_{n},\delta)$ and perform projected gradient descent. In order to avoid complex statistical dependency issues that can arise in the analysis of gradient descent, for our theoretical results we consider a sample-splitting variant of the algorithm where each iteration is performed on a fresh batch of samples. We summarize the overall robust gradient descent algorithm via gradient estimation in Algorithm 1. In contrast to $M$-estimation where we use robust estimates of the overall loss function, here we use robust estimates of the gradient, a small shift in perspective, but which has strong statistical and computational consequences: we obtain a computationally practical algorithm, and moreover with strong robustness guarantees via careful statistical analyses of the stability of the resulting biased and inexact gradient descent iterates.

We further assume that the number of gradient iterations $T$ is specified a-priori, and accordingly we define:

We discuss methods for selecting $T$, and the impact of sample-splitting in later sections. As confirmed in our experiments (see Section 4), sample-splitting should be viewed as a device introduced for theoretical convenience which can likely be eliminated via more complex uniform arguments (see for instance the work ).

It can be seen that the key ingredient in the robust gradient descent Algorithm 1 is a robust estimator of the gradients. Next, we consider the two notions of robustness described in Section 2, and derive specific gradient estimators for each of the models using the framework described above. Although we derive corollaries of our general results for these two settings of Huber contamination and heavy-tailed models, we emphasize that our class of estimators are more general and are not restricted to these two notions of robustness.

There has been a flurry of recent interest in designing mean estimators which, under the Huber contamination model, can robustly estimate the mean of a random vector. While some of these results are focused on the case where the uncorrupted distribution is Gaussian, or isotropic, we are more interested in robust mean oracles for more general distributions. Lai et al. proposed a robust mean estimator for general distributions, satisfying weak moment assumptions, and we leverage the existence of such an estimator to design a Huber gradient estimator $g(\theta;\mathcal{D}_{n},\delta)$ which works in the Huber contamination model.

The estimator builds upon the fact that with a single dimension, it is relatively easy to estimate the gradient robustly. In higher dimensions, the crucial insight of Lai et al. is that the effect of the contamination distribution $Q$ on the mean of uncontaminated distribution $P$ is effectively one-dimensional provided we can accurately estimate the direction along which the mean is shifted. In our context, if we can compute the gradient shift direction, i.e. the direction of the difference between the sample (corrupted) mean gradient and the true (population) gradient, then the true gradient can be estimated by using a robust 1D-mean algorithm along the gradient-shift direction and a non-robust sample-gradient in the orthogonal direction since the contamination has no effect on the gradient in this orthogonal direction. In order to identify this gradient shift direction, we follow Lai et al. and use a recursive Singular Value Decomposition (SVD) based algorithm. In each stage of the recursion, we first remove gross-outliers via a truncation algorithm (described in more detail in the Appendix, and termed HuberOutlierGradientTruncation in Algorithm 2). We subsequently identify two subspaces using an SVD – a clean subspace where the contamination has a small effect on the mean and another subspace where the contamination has a potentially larger effect. We use a simple sample-mean estimator in the clean subspace and recurse our computation on the other subspace. Building on the work of Lai et al. , in Lemma 1 and Appendix K we provide a careful non-asymptotic analysis of this gradient estimator.

Algorithm 2 presents the overall Huber gradient estimator $g(\theta;\mathcal{D}_{n},\delta)$.

2 Gradient Estimation in the Heavy-Tailed model

To design gradient estimators for the heavy-tailed model, we leverage recent work on designing robust mean estimators in this setting. These robust mean estimators build on the classical work of Alon et al. , Nemirovski and Yudin and Jerrum et al. on the so-called median-of-means estimator. For the problem of one-dimensional mean estimation, Lerasle and Oliveira , Catoni et al. propose robust mean estimators that achieve exponential concentration around the true mean for any distribution with bounded second moment. In this work we require mean estimators for multivariate distributions. Several recent works () extend the median-of-means estimator of to general metric spaces. In this paper we use the geometric median-of-means estimator (Gmom), which was originally proposed and analyzed by Minsker , to design the gradient estimator $g(\theta;\mathcal{D}_{n},\delta)$.

Algorithm 3 presents the gradient estimator $g(\theta;\mathcal{D}_{n},\delta)$ obtained using Gmom as the mean estimator.

3 Choice of Hyper-Parameters

In this section, we discuss how to tune the hyperparameters for our algorithms. In particular, note that the gradient estimators described in Algorithms 2, 3 depend on corruption level $\epsilon$, and on confidence $\delta$, which are not known in advance.

Since the standard hyper-parameter selection techniques such as cross validation, hold-out validation, pick hyper-parameters that minimize the empirical mean of the loss on validation data, they can’t be used in the presence of outliers in the data. One criteria we could use in such cases is to choose hyper-parameters that minimize a robust estimate of the population risk on validation data. However, we can’t use any of the existing robust mean estimators to estimate the population risk because they themselves depend on hyper-parameters such as corruption level $\epsilon$.

We now consider the Huber contamination model and propose a heuristic based on Scheffe’s tournament estimator for hyper-parameter selection. In particular, we consider the gradient descent procedure described in Algorithm 2 and explain our technique for choosing $\epsilon,\delta$ using hold out cross validation. Note that our goal is to pick hyper-parameters that minimize the population risk $\mathcal{R}(\theta)$. Under the assumption of strong convexity of $\mathcal{R}(\theta)$, this is equivalent to picking hyper-parameters that minimize the parameter error $\|\theta-\theta^{*}\|_{2}$.

We begin with the problem of density estimation, where we are given $n$ i.i.d samples $\{z_{i}\}_{i=1}^{n}$ from $(1-\epsilon)P_{\theta^{*}}+\epsilon Q$, where $P_{\theta^{*}}$ belongs to the class of distributions $\{P_{\theta}:\theta\in\Theta\}$, and $Q$ is an arbitrary distribution. Our goal is to estimate $\theta^{*}\in\Theta$ from the samples. Suppose $\left\{P_{\widehat{\theta}_{1}},P_{\widehat{\theta}_{2}},\ldots,P_{\widehat{\theta}_{m}}\right\}$ are the candidate solutions returned by Algorithm 2 for different settings of $\epsilon,\delta$. Consider the following pairwise test function:

Following , it can be shown that the above procedure picks a $j^{*}$ such that $P_{\widehat{\theta}_{j^{*}}}$ is close to $P_{\theta^{*}}$ in TV metric. For distributions $\{P_{\theta}:\theta\in\Theta\}$ whose TV metric is roughly equivalent to the parameter error, the above procedure results in hyper-parameters which minimize the parameter error. This procedure can be extended to supervised learning problems such as regression and classification.

For the Heavy-Tailed setting we experimented with (a) empirical mean of the risk on validation data: $\frac{1}{n_{val}}\sum_{i=1}^{n_{val}}\overline{\mathcal{L}}(\theta;z^{\prime}_{i})$ where $\{z^{\prime}_{i}\}_{i=1}^{n_{val}}$ is the validation data, which does not require any tuning parameters, as well as (b) median of means based mean of the risk on validation data, for various confidence levels $\delta$. However, both the techniques in the context of hold-out validation resulted in models with similar performance. So, in our experiments with heavy tailed distributions, we present results obtained using the empirical risk as in (a) on hold-out validation data.

Experiments

In this section we demonstrate our proposed methods for the Huber contamination and heavy-tailed models, on a variety of simulated and real data examples.

We first consider the Huber contamination model and demonstrate the practical utility of gradient-descent based robust estimator described in Algorithms 1 and 2.

Recall the linear regression model described in (4) where we observe paired samples $\{(x_{i},y_{i})\}_{i=1}^{n}$. We assume that the $(x,y)$ pairs sampled from the true distribution $P$ are linked via a linear model:

We now describe the experiment setup, the data model and present the results.

We fix the contamination level $\epsilon=0.1$ and $\sigma^{2}=0.1$. Next, we generate $(1-\epsilon)n$ clean covariates from a multivariate Gaussian $x\sim\mathcal{N}(0,\mathcal{I}_{p})$, and we generate the corresponding clean responses using $y=\left\langle x,\theta^{*}\right\rangle+w$ where $\theta^{*}=[1,\ldots,1]^{T}$ and $w\sim\mathcal{N}(0,\sigma^{2})$. We simulate an outlier distribution by drawing the covariates from $\mathcal{N}(0,p^{2}\mathcal{I}_{p})$, and setting the responses to . The total number of samples is set to be $(10\frac{p}{\epsilon^{2}})$. We note that the sample size we choose increases with the dimension. This scaling is used to ensure that the statistical (minimax) error, in the absence of any contamination, is roughly 0.001. An optimally robust method should have error close to 0.1 (roughly equal to corruption level), which ours does (see Figure 1).

As our baselines, we use OLS, TORRENT , the Huber-loss M-estimator, RANSAC and a plugin estimator (detailed further in Section sec:linregtheory, and which in a nutshell robustly estimates the sufficient statistics required for the OLS estimator). TORRENT is an iterative hard-thresholding based alternating minimization algorithm, where in one step, it calculates an active set of examples by keeping only $(1-\epsilon)n$ samples which have the smallest absolute values of residual $r=y-\left\langle x,\theta^{t}\right\rangle$, and in the other step it updates the current estimates by solving OLS on the active set. Bhatia et al. have shown the superiority of TORRENT over a variety of other convex-penalty based outlier techniques, hence, we do not compare against those methods. The plugin estimator is implemented using Algorithm 2 to estimate both the mean vector $\frac{1}{n}\sum_{i=1}^{n}y_{i}x_{i}$ and the covariance matrix $\frac{1}{n}\sum_{i=1}^{n}x_{i}x_{i}^{T}$, which are the required sufficient statistics for the OLS estimator.

Error vs dimension $p$: All estimators except our proposed algorithm perform poorly with increasing dimension, as shown in Figure 1(a). Note that the TORRENT algorithm has strong guarantees when only the response $y$ is corrupted but performs poorly in the Huber contamination model where both $x$ and $y$ may be contaminated. We find that the error for the robust plugin estimator increases with dimension. We investigate this theoretically in Section 6.1, where we find that the error of the plugin estimator grows with the norm of $\theta^{*}$. In our experiments, we choose $\|\theta^{*}\|_{2}=\sqrt{p}$, and thus Figure 1(a) corroborates Corollary 3 in Section 6.1.

Error vs $\epsilon$: In Figure 1(b) we find that the parameter error $\|\widehat{\theta}-\theta^{*}\|_{2}$ increases linearly with the contamination rate $\epsilon$ and we study this further in Section 6.1.

Error vs iteration $t$: Finally, Figure 1(c) shows that the convergence rate decreases with increasing contamination $\epsilon$ and after $\epsilon$ is high enough, the algorithm remains stuck at $\theta_{0}$, corroborating Lemma 8 (in the Appendix).

Hyper-parameter Tuning: In Figures 1(d) and 1(e), we find the final solution chosen by our tournament based heuristic for hyper-parameter selection (TournamentGD) has roughly the same performance as the algorithm which knows the true value of $\epsilon$ (OracleGD). In particular, our final error does not scale with $p$.

Next, we study the performance of our proposed method in the context of classification.

1.2 Synthetic Experiments: Logistic Regression

We simulate a linearly separable classification problem, where the clean covariates are sampled from $\mathcal{N}(0,\mathcal{I}_{p})$, the corresponding clean responses are computed as $y=\text{sign}(\left\langle x,\theta^{*}\right\rangle)$ where $\theta^{*}=[1/\sqrt{p},\ldots,1/\sqrt{p}]^{T}$. We simulate the outlier distribution by adding asymmetric noise, i.e. we flip the labels of one class, and increase the variance of the corresponding covariates by multiplying them by $p^{2}$. The total number of samples are set to be $(10\frac{p}{\epsilon^{2}})$.

We measure the 0-1 classification error on a held-out (clean) test set. We study how the 0-1 error changes with $p$ and $\epsilon$ and the parameter estimation error of our proposed method for different contamination levels $\epsilon$.

We use the logistic regression MLE and the linear Support Vector Machine (SVM) as our baselines.

0/1 Error vs dimension $p$: In Figure 2(a) we observe that both the SVM and logistic regression MLE perform poorly with increasing dimension. The logistic regression MLE completely flips the labels and has a 0-1 error close to 1, whereas the linear SVM outputs a random hyperplane classifier that flips the label for roughly half of the dataset.

0/1 Error vs $\epsilon$ and $t$: Figures 2(b) and 2(c) show qualitatively similar results to the linear regression setting, i.e. that the error of our proposed estimator degrades gracefully (and grows linearly) with the contamination level $\epsilon$ and that the gradient descent iterates converge linearly.

1.3 Robust Face Reconstruction

In this experiment, we show the efficacy of our algorithm by attempting to reconstruct face images that have been corrupted with heavy occlusion, where the occluding pixels play the role the outliers. We use the data from the Cropped Yale Dataset . The dataset contains 38 subjects, and each image has $192\times 168$ pixels. Following the methodology of Wang et al. , we choose 8 face images per subject, taken under mild illumination conditions and computed an eigenface set with 20 eigenfaces. Then given a new corrupted face image of a subject, the goal is to get the best reconstruction/approximation of the true face. To remove scaling effects, we normalized all images to $$range. One image per person was used to test reconstruction. Occlusions were simulated by randomly placing$10$blocks of size$30\times 30$. We repeated this 10 times for each test image. In this setting, each image of a subject corresponds to a different task; i.e.$X$is a common fixed eigenface basis,${y}$ is an observed(occluded) image, and the goal is to reconstruct(de-noise) the given image using the given basis. Note that in this example, we use a linear regression model as the uncontaminated statistical model, which is almost certainly not an exact match for the unknown ground truth distribution. Despite this model misspecification, as our results show, that robust mean based gradient algorithms do well.

We use Root Mean Square Error (RMSE) between the original and reconstructed image to evaluate the performance of the algorithms. We also compute the best possible reconstruction of the original face image by using the 20 eigenfaces.

Wang et al. implemented popular robust estimators such as RANSAC, Huber Loss etc. and showed their poor performance. Wang et al. then proposed an alternate robust regression algorithm called Self Scaled Regularized Robust Regression(SCRRR). Hence, use TORRENT, SCRRR and OLS as baselines. We also compare against the best possible RMSE obtained by reconstructing the un-occluded image using the eigenfaces.

Table 1 shows that the mean RMSE is best for our proposed gradient descent based method and that the recovered images are in most cases closer to the un-occluded original image. (Figure 3). Figure 3(c) shows a case when none of the methods succeed in reconstruction.

2 Heavy-tailed Estimation

We now consider the heavy-tailed model and present experimental results on synthetic datasets comparing the gradient descent based robust estimator described in Algorithms 1 and 3 (which we call RobustGD) with ERM and several other recent proposals. In these experiments we focus on the problem of linear regression which is described in Section 4.1 and work with heavy-tailed noise distributions.

We compare RobustGD with several baselines. Since we are always in the low-dimensional ($n\geq p$) setting, the solution to ERM has a closed form expression and is simply the OLS solution. We also study OLS-GD, which performs a gradient descent on ERM and is equivalent to using empirical mean as the gradient oracle in our framework. We also compare against the robust estimation techniques of Hsu and Sabato and Duchi and Namkoong , which we refer to as RobustHS, RobustDN and two classical techniques namely the LASSO and ridge regression. In our experiments, all the iterative techniques are run until convergence.

We use two metrics to compare the performance of various approaches: a) parameter error which is defined as $\|\theta-\theta^{*}\|_{2}$ and b) to compare the performance of two estimators $\widehat{\theta}_{1}$, $\widehat{\theta}_{2}$, we define the notion of relative efficiency:

Roughly, this corresponds to the percentage improvement in the parameter error obtained using $\widehat{\theta}_{1}$ over $\widehat{\theta}_{2}$. Whenever $\text{\sc RelEff}(\widehat{\theta}_{1},\widehat{\theta}_{2})>0$, $\widehat{\theta}_{1}$ has a lower parameter error, and higher the value, the more the fractional improvement.

To reduce the variance in the plots presented here, we averaged results over $20$ repetitions. Figure 4 shows the benefits of using RobustGD over other baselines.

Error vs number of iterations: In Figures 4(a), 4(b) we plot the excess risk of various approaches against the number of iterations (for OLS, LASSO, ridge regression and the method of Hsu and Sabato we only plot the excess risk of the final iterate). We see that upon convergence RobustGD has a much lower parameter error. As expected, OLS-GD converges to OLS.

Error vs number of samples: Next, in Figures 4(c), 4(d) we plot the parameter error as $n/p$ increases. We see that RobustGD is always better than other baselines, even when the number of samples is 12 times the dimension $p$.

Relative Efficiency vs $\beta$, and $\sigma$: In Figure 4(e), we plot the relative efficiency against $\beta$, the moment bound of Pareto distribution. This shows that the percentage improvement in the excess risk by RobustGD decreases as the moment bound $\beta$ increases. This behavior is expected because as we increase the moment bound the tails of the noise distribution become lighter. This shows that there is more benefit in using RobustGD in the heavy tailed setting. We do a similar study to see the relative efficiency against the variance of the noise distribution. Figure 4(f) plots relative efficiency against standard deviation of the noise distribution.

Theoretical Preliminaries

In this section we develop some theoretical preliminaries. We first develop a general theory on convergence of projected gradient descent in Section 5.1. Next we analyze the gradient estimators defined in Algorithms 2 and 3 in Sections 5.2 and 5.3 respectively. Finally in Sections 6 and 7 we present consequences of our general theory for the canonical examples of risk minimization described in Section 2.2, under Huber contamination and heavy-tailed models.

In this section we develop a general theory for the convergence of the projected gradient descent described in Algorithm 1. Note that our gradient estimators could be biased and are not guaranteed to be consistent estimators of the true gradient $\nabla\mathcal{R}(\theta)$. This is especially true in the Huber contamination model where it is impossible to obtain consistent estimators of the gradient of the risk because of the non-vanishing bias caused by the contaminated samples. Hence, we turn our attention to understanding the behavior of projected gradient descent with a biased, inexact, gradient estimator of the form in (11). Before we present our main result, we define the notion of stability of a gradient estimator, which plays a key role in the convergence of gradient descent.

We denote by $\kappa$ the following contraction parameter:

and note that $\kappa<1$. With these definitions in place we state our main result on the stability of gradient descent:

We defer a proof of this result to the Appendix. For the bound (15), the first term is decreasing in $T$, while the second term is increasing in $T$. This suggests that for a given $n$ and $\delta$, we need to run just enough iterations for the first term to be bounded by the second. Hence, we can fix the number of iterations $T^{*}$ as the smallest positive integer such that:

Since we obtain linear convergence, i.e. $\kappa<1$, typically a logarithmic number of iterations suffice to obtain an accurate estimate.

Theorem 1 provides a general result for risk minimization and parameter estimation, and requires bounds on $\alpha(\widetilde{n},\widetilde{\delta}),\beta(\widetilde{n},\widetilde{\delta})$ which capture the the error suffered by the gradient estimator for a given risk minimization problem. In any concrete instantiation for a given gradient estimator, risk pair, we first estimate these gradient estimator error bounds by studying the distribution of the gradient of the risk, and then apply Theorem 1. In the next two sections, we provide some general analyses of the gradient estimator in Algorithm 2 for the Huber contamination model, and the gradient estimator in Algorithm 3 for the heavy-tailed model, and which apply to any risk minimization problem. In Sections 6,7 we then instantiate these gradient estimator error results for various illustrative statistical models such as linear regression, logistic regression, and general exponential families. Plugging these into Theorem 1, we then get consequences of our robustness guarantees for various statistical model, robustness setting pairs.

2 General Analysis of Huber Contamination Gradient Estimator in Algorithm 2

We now analyze the gradient estimator described in Algorithm 2 for Huber contamination model and study the error suffered by it. As stated before, Algorithm 2 uses the robust mean estimator of Lai et al. . Hence, while our proof strategy mimics that of Lai et al. , we present a different result which is obtained by a more careful non-asymptotic analysis of the algorithm.

and with this definition in place we have the following result:

We note in particular, if $n\rightarrow\infty$ (with other parameters held fixed) then $\gamma(n,p,\delta,\epsilon)\rightarrow 0$ and the error of our gradient estimator satisfies

and has only a weak dependence on the dimension $p$.

3 General Analysis of Heavy-tailed Model Gradient Estimator in Algorithm 3

In this section we analyze the gradient estimator for heavy-tailed setting, described in Algorithm 3. The following result shows that the gradient estimate has exponential concentration around the true gradient, under the mild assumption that the gradient distribution has bounded second moment. Its proof follows directly from the analysis of geometric median-of-means estimator of Minsker . We use $\text{tr}\left(\Sigma_{\theta}\right)$ to denote the trace of the matrix $\Sigma_{\theta}$.

The results of the Lemmas 1 and 2 effectively ensure that under relatively mild moment assumptions we can robustly estimate multivariate mean vectors and in subsequent sections we show how to leverage these strong guarantees for robust parametric estimation.

Consequences for Estimation under ϵitalic-ϵ\epsilon-Contaminated Model

We now turn our attention to the examples introduced earlier, and present specific applications of Theorem 1, for parametric estimation under Huber contamination model. As shown in Lemma 1, we need the added assumption that the true gradient distribution has bounded fourth moments, which suggests the need for additional assumptions. We make our assumptions explicit and defer the technical details to the Appendix.

In the asymptotic setting when the number of samples $n\rightarrow\infty$ (and other parameters are held fixed), we see that for the Huber Gradient Estimator, the corresponding maximum allowed contamination level is

i.e. the better conditioned the covariance matrix $\Sigma$, the higher the contamination level we can tolerate.

Comparing bounds (17) and (18), we see that even when it does not have to estimate the covariance matrix, the error of the plugin estimator depends on $\|\theta^{*}\|_{2}$, which would make the estimator vacuous if $\|\theta^{*}\|_{2}$ scales with the dimension $p$. On the other hand, the asymptotic rate of our robust gradient estimator is independent of $\|\theta^{*}\|_{2}$. This disadvantage of plugin estimation is inescapable, and is seen for instance in known minimax results for robust mean estimation that show that the dependence on $\|\theta^{*}\|_{2}$ is unavoidable for any oracle which estimates the mean of $xy$ in the $\epsilon$-contaminated setting. Next, we apply our estimator to generalized linear models.

2 Generalized Linear Models

Here we assume that the covariates have bounded 8${}^{\text{th}}$ moments. Additionally, we assume smoothness of $\Phi^{\prime}(\cdot)$ around $\theta^{*}$. In particular, we assume that there exist universal constants $L_{\Phi,2k}$, $B_{2k}$ such that

Consider the statistical model in equation (5), and suppose that the number of samples $n$ is large enough such that

and the contamination level is such that,

for some contraction parameter $\kappa<1$.

Note that for the case of linear regression with gaussian noise, it is relatively straightforward to see that $L_{\Phi,2k}=C_{2k}\|\Sigma\|_{2}^{k}$, $B_{\Phi,2k}=0$, $M_{\Phi,t,k}=1\ \ \forall(t\geq 2,k\in\mathcal{N})$ and $M_{\Phi,t,k}=0\ \ \forall(t\geq 3,k\in\mathcal{N})$ under the assumption of bounded $8^{th}$ moments of the covariates; which essentially leads to an equivalence between Theorem 2 and Theorem 4 for this setting. In the following section, we instantiate the above Theorem for logistic regression and compare and contrast our results to other existing methods.

By observing that $\Phi^{(t)}(\cdot)$ is bounded for logistic regression for all $t\geq 1$, we can see that $L_{\Phi,2k}=0$, and that there exists a universal constant $C>0$ such that $B_{\Phi,2k}<C$ and $M_{\Phi,t,k}<C\ \ \forall(t\geq 1,k\in\mathcal{N})$.

for some contraction parameter $\kappa<1$.

Under the restrictive assumption that $x\sim\mathcal{N}(0,\mathcal{I}_{p})$, Du et al. exploited Stein’s trick to derive a plugin estimator for logistic regression. However, similar to the linear regression, the error of the plugin estimator scales with $\|\theta^{*}\|_{2}$, which is avoided in our robust gradient descent algorithm. We also note that our algorithm extends to general covariate distributions.

3 Exponential Family

Here we assume that the random vector $\phi(z),z\sim P$ has bounded 4${}^{\text{th}}$ moments.

for some contraction parameter $\kappa<1$.

4 Discussion and Limitations

Consequences for Heavy-Tailed Estimation

In this section we present specific applications of Theorem 1 for parametric estimation, under heavy tailed setting. The proofs of the results can be found in the Appendix.

Consider the statistical model in equation (4). There are universal constants $C_{1},C_{2}>0$ such that if

for some contraction parameter $\kappa<1$.

2 Generalized Linear Models

In this section we consider generalized linear models described in Equation (5), where the covariate $x$ is allowed to have a heavy tailed distribution. Here we assume that the covariates have bounded 4${}^{\text{th}}$ moment. Additionally, we assume smoothness of $\Phi^{\prime}(\cdot)$ around $\theta^{*}$. Specifically, we assume that there exist universal constants $L_{\Phi,2k}$, $B_{2k}$ such that

Consider the statistical model in equation (5). There are universal constants $C_{1},C_{2}>0$ such that if

for some contraction parameter $\kappa<1$.

We now instantiate the above Theorem for logistic regression model.

Consider the model in equation(7). There are universal constants $C_{1},C_{2}>0$ such that if

for some contraction parameter $\kappa<1$.

3 Exponential Family

We now instantiate Theorem 1 for parameter estimation in heavy-tailed exponential family distributions. Here we assume that the random vector $\phi(z),z\sim P$ has bounded 2${}^{\text{nd}}$ moments, and we obtain the following result:

for some contraction parameter $\kappa<1$ and universal constant $C$.

Discussion

In this paper we introduced a broad class of robust estimators, that leverage the inherent robustness of gradient descent, together with the observation that for risk minimization in most statistical models, the gradient of the risk takes the form of a simple multivariate mean, which can be robustly estimated using recent work on robust mean estimation. In contrast to classical $M$-estimators that use robust estimates of the risk, our class of estimators employ a shift in perspective, and use robust estimates of gradients of the risk instead, which can then be embedded into a simple projected gradient descent iterative algorithm. Our class of robust gradient descent estimators work well in practice and in many cases outperform other robust (and non-robust) estimators. We also show that these estimators have strong robustness guarantees under varied robustness settings, including Huber’s $\epsilon$-contamination model and for heavy-tailed distributions.

There are several avenues for future work, including a better understanding of robust mean estimation, any improvement in which would immediately translate to improved guarantees for our robust gradient descent estimators. Finally, it would also be of interest to understand the extent to which we could replace gradient descent with other optimization methods such as accelerated gradient descent or Newton’s method. We note however, that although these methods may have faster rates of convergence in the classical risk minimization settings, in our setup their stability to using inexact gradients is far more crucial and warrants further investigation.

Acknowledgements

A.P., A.S., P.R. acknowledge the support of PNC, and NSF via IIS-1149803, IIS-1664720, DMS-1264033. S.B. acknowledges the support of NSF via DMS-1713003. We thank Larry Wasserman and Ankit Pensia for helpful comments on the paper.

References

Appendix A Proof of Theorem 1

In this section, we present the proof of our main result on projected gradient descent with an inexact gradient estimator. To ease the notation we will often omit $\{D_{n},\delta\}$ from $g(\theta;D_{n},\delta)$.

At any iteration step $t\in\{1,2,\ldots,T\}$, by assumption we have that with probability at least $1-\frac{\delta}{T}$,

Taking union bound, (27) holds over all iteration steps $t\in\{1\dots T\}$, with probability at least $1-\delta$. For the remainder of the analysis, we assume this event to be true.

Let $g(\theta^{k})=\nabla\mathcal{R}(\theta^{k})+e_{k}$ be the noisy gradient. Let $\alpha=\alpha(n/T,\delta/T)$ and $\beta=\beta(n/T,\delta/T)$ for brevity.

We have the following Lemma from Bubeck .

where Equation (28) follows from contraction property of projections. Now, we can write $\|\theta^{k}-\theta^{*}-\eta(\nabla\mathcal{R}(\theta^{k})-\nabla R(\theta^{*}))\|_{2}$ as

for some $\kappa<1$. Solving the induction,we get:

Appendix B Proof of Theorem 2

The proof of Theorem 2 follows from Theorem 4, where we study Generalized Linear Models, which include linear regression as a special case. For the case of linear regression with gaussian noise, it is relatively straightforward to see that the smoothness parameters satisfy $L_{\Phi,2k}=C_{2k}\|\Sigma\|_{2}^{k}$, $B_{\Phi,2k}=0$, $M_{\Phi,t,k}=1\ \ \forall(t\geq 2,k\in\mathcal{N})$ and $M_{\Phi,t,k}=0\ \ \forall(t\geq 3,k\in\mathcal{N})$ under the assumption of bounded $8^{th}$ moments of the covariates. Substituting these values in Theorem 4 gives us the required result.

Appendix C Proof of Theorem 4

To prove our result on Robust Generalized Linear Models, we first study the distribution of gradients of the corresponding risk function.

Consider the model in Equation (5), then there exist universal constants $C_{1},C_{2}>0$ such that

The gradient $\nabla\overline{\mathcal{L}}(\theta)$ and it’s expectation can be written as:

where the last line follows from our assumption of smoothness.

Using the $C_{r}$ inequality, we have that

where the last line follows from our assumption that $P_{\theta^{*}}(y|x)$ is in the exponential family, hence, the cumulants are higher order derivatives of the log-normalization function.

Bounded Fourth Moment. To show that the fourth moment of the gradient distribution is bounded, we have

where the last step follows from the fact that the 8th central moment can be written as a polynomial involving the lower cumulants, which in turn are the derivatives of the log-normalization function.

By assumption $L_{\Phi,k},B_{\Phi,k},M_{\Phi,t,k}$ are all bounded for $k,t\leq 8$, which implies that there exist constants $c_{1},c_{2}>0$ such that

Having studied the distribution of the gradients, we use Lemma 1 to characterize the stability of Huber Gradient estimator. Using Lemma 1, we know that at any point $\theta$, the Huber Gradient Estimator $g(\theta,\delta/T)$ satisfies that with probability $1-\delta/T$,

Appendix D Proof of Corollary 3

We begin by studying the distribution of the random variable $xy=xx^{T}\theta^{*}+x.w$.

Consider the model in Equation (4), with $x\sim\mathcal{N}(0,\mathcal{I}_{p})$ and $w\sim\mathcal{N}(0,1)$ then there exist universal constants $C_{1},C_{2}$ such that

Now, $Z=(xx^{T}-I)\theta^{*}$ can be written as:

Hence the covariance matrix can be written as:

Control of $A$. Using Cauchy Schwartz, and normality of 1D projections of normal distribution

Control of $B$, $B\leq\|\theta^{*}\|_{2}^{4}$.

Control of $C$, $C=O(1)$, using independence of $w$ and normality of 1D projections of normal distribution.

Appendix E Proof of Theorem 6

To prove our result on Robust Exponential Family, we first study the distribution of gradients of the corresponding risk function.

Consider the model in Equation (8), then there exists a universal constant $C_{1}$ such that

By Fisher Consistency of the negative log-likelihood, we know that

Bounded moments follows from our assumption that the sufficient statistics have bounded 4th moments. ∎

Having studied the distribution of the gradients, we use Lemma 1 to characterize the stability of Huber Gradient estimator. Using Lemma 1, we know that at any point $\theta$, the Huber Gradient Estimator $g(\theta,\delta/T)$ satisfies that with probability $1-\delta/T$,

Appendix F Proof of Corollary 7

Using the contraction property of projections, we know that

By Fisher Consistency of the negative log-likelihood, we know that

The true parameter $\theta^{*}$ can be obtained by inverting the $\nabla A$ operator whenever possible.

where $A^{*}$ is the convex conjugate of $A$. We can use the following result to control the Lipschitz smoothness $A^{*}$.

(Strong/Smooth Duality) Assume $f(\cdot)$ is closed and convex. Then $f(\cdot)$ is smooth with parameter $M$ if and only if its convex conjugate $f(\cdot)$ is strongly convex with parameter $m=\frac{1}{M}$.

A proof of the above theorem can be found in . Hence, we have that:

Combining the above with Equation (63) recovers the result of Corollary 7.

Appendix G Proof of Theorem 8

Before we present the proof of Theorem 8, we first study the distribution of gradients of the loss function. This will help us bound the error in the gradient estimator.

where $\Delta=\theta-\theta^{*}$ and $E[xx^{T}]=\Sigma$.

Next, we bound the operator norm of the covariance of the gradients $\nabla\overline{\mathcal{L}}(\theta)$ at any point $\theta$.

where the second last step follows from Cauchy-Schwartz and the last step follows from our assumption of bounded 4th moments (see Equation (14)). ∎We now proceed to the proof of Theorem 8. From Lemma 2, we know that at any point $\theta$, the gradient estimator described in Algorithm 3, $g(\theta;D_{\widetilde{n}},\widetilde{\delta})$, satisfies the following with probability at least $1-\delta$,

To complete the proof of this theorem, we use the results from Theorem 1. Note that the gradient estimator satisfies the stability condition if $\alpha(\widetilde{n},\widetilde{\delta})<\tau_{l}$. This holds when

Now suppose $\widetilde{n}$ satisfies the above condition, then plugging $\beta(\widetilde{n},\widetilde{\delta})$ into Theorem 1 gives us the required result.

Appendix H Proof of Theorem 9

To prove the Theorem we use the result from Lemma 4, where we derived the following expression for covariance of $\nabla\overline{\mathcal{L}}(\theta)$

From Lemma 2, we know that at any point $\theta$, the gradient estimator described in Algorithm 3, $g(\theta;D_{\widetilde{n}},\widetilde{\delta})$, satisfies the following with probability at least $1-\delta$,

We now use the results from Theorem 1. The gradient estimator satisfies the stability condition if $\alpha(\widetilde{n},\widetilde{\delta})<\tau_{l}$. This holds when

Now suppose $\widetilde{n}$ satisfies the above condition, then plugging $\beta(\widetilde{n},\widetilde{\delta})$ into Theorem 1 gives us the required result.

Appendix I Proof of Theorem 11

Since $\alpha(\widetilde{n},\widetilde{\delta})=0$, the stability condition is always satisfied, as long as $\tau_{l}>0$. Substituting $\beta(\widetilde{n},\widetilde{\delta})$ into Theorem 1 gives us the required result.

Appendix J Upper bound on Contamination Level

We provide a complementary result, which gives an upper bound for the contamination level $\epsilon$ based on the initialization point $\theta^{0}$, above which, Algorithm 1 would not work. The key idea is that the error incurred by any mean estimation oracle is lower bounded by the variance of the distribution, and that if the zero vector lies within that error ball, then any mean oracle can be forced to output $\mathbf{0}$ as the mean. For Algorithm 1, this implies that, in estimating the mean of the gradient, if the error is high, then one can force the mean to be $\mathbf{0}$ which forces the algorithm to converge. For the remainder of the section we consider the case of linear regression with $x\sim\mathcal{N}(0,\mathcal{I}_{p})$ in the asymptotic regime of $n\rightarrow\infty$.

Consider the model in equation(4) with $x\sim\mathcal{N}(0,\mathcal{I}_{p})$ and $w\sim\mathcal{N}(0,1)$, then there exists a universal constant $C_{1}$ such that if $\epsilon>C_{1}\frac{\|\theta^{0}-\theta^{*}\|_{2}}{\sqrt{1+2\|\theta^{0}-\theta^{*}\|_{2}^{2}}}$, then for every gradient oracle, there exists a contamination distribution $Q$ such that, Algorithm 1 will converge to $\theta^{0}$ even when the number of samples $n\rightarrow\infty$.

Using Lemma 5, we know that for any point $\theta$,

where $\Delta=\theta-\theta^{*}$. Let $P_{\nabla\overline{\mathcal{L}}(\theta)}$ represent the distribution $\nabla\overline{\mathcal{L}}\theta$. Similarly, let $P_{\epsilon,\nabla\overline{\mathcal{L}}(\theta),Q}$ represent the corresponding $\epsilon$-contaminated distribution. Then, using Theorem 2.1 , we know that the minimax rate for estimating the mean of the distribution of gradients is given by:

Suppose that the contamination level $\epsilon$ is such that,

then for every oracle there exists a corresponding $Q$ such that Algorithm 1 will remain stuck at $\theta^{0}$.

Chen et al. provide a general minimax lower bound of $\Omega(\epsilon)$ for $\epsilon$-contamination models in this setting. In contrast, using Algorithm 1 with as oracle, we can only $O(\sqrt{\epsilon\log p})$ close to the true parameter even when the contamination is small, which implies that our procedure is not minimax optimal. Our approach is nonetheless the only practical algorithm for robust estimation of general statistical models.

Appendix K Proof of Lemma 1

In this section we present a refined, non-asymptotic analysis of the robust mean estimator of , described in Algorithm 2. We begin by introducing some preliminaries. We subsequently analyze the algorithm in 1-dimension and finally turn our attention to the general algorithm.

Unless otherwise stated, we assume throughout that the random variable $X$ has bounded fourth moments, i.e. for every unit vector $v$,

We summarize some useful results from , which bound the deviation of the conditional mean/covariance from the true mean/covariance.

For random variables with bounded fourth moments we can use Chebyshev’s inequality to obtain tail bounds.

[Lemma 3.14 ] Let $X$ have bounded fourth moments, then for every unit vector $v$ we have that,

With these preliminaries in place we use the following result from .

Equivalently, with probability at least $1-\delta$,

We now turn our attention to an analysis of Algorithm 2 for the 1-dimensional case.

K.2 The case when p=1𝑝1p=1

Firstly, we analyze Algorithm 2 when $p=1$.

where $t=C_{8}\sqrt{\frac{1}{n}\log\left(\frac{n}{\delta}\right)}$. which can be further simplified to,

By an application of Hoeffding’s inequality we obtain that with probability at least $1-\delta/3$, the fraction of corrupted samples (i.e. samples from the distribution $Q$) is less than $\epsilon+\sqrt{\frac{\log(3/\delta)}{2n}}$. We condition on this event through the remainder of this proof. We let $\eta$ denote the fraction of corrupted samples. Further, we let $S_{P}$ be the samples from the true distribution. Let $n_{P}$ be the cardinality of this set, i.e. $n_{P}:=|S_{P}|$.

Let $I_{1-\eta}$ be the interval around $\mu$ containing $1-\eta$ mass of $P_{\theta}^{*}$. Then, using Lemma 11, we have that:

Using Lemma 13 we obtain that with probability at least $1-\delta/3$ the number of samples from the distribution $P$ that fall in the interval $I_{1-\eta}$ is at least $1-\eta-t$ where $t$ is upper bounded as:

Now we let $\widetilde{S}$ be the set of points in the smallest interval containing $(1-\eta-t)(1-\eta)$ fraction of all the points.

This can be re-written as, that with probability at least $(1-\delta/3)$, there exists a universal constant $C_{0}$ such that,

Using Equation (78), we know that $(1-\eta-t)$ fraction of $S_{D}$ lie in $\mathcal{I}_{1-\eta}$.

Let $\widetilde{S}$ be the set of points in the smallest interval containing $(1-\eta-t)(1-\eta)$ fraction of the points.

We know that the length of minimum interval containing $(1-\eta-t)(1-\eta)$ fraction of the points of $S$ is less than length of smallest interval containing $(1-\eta-t)$ fraction of points of $S_{D}$, which in turn is less than length of $I_{1-\eta}$.

Now, $I_{1-\eta}$ and minimum interval containing $(1-\eta-t)$ fraction of points of $S_{D}$ need to overlap. This is because, $n$ is large enough such that $t<\frac{1}{2}-\eta$ hence, the extreme points for such an interval can be atmost $2\text{length}\left(I_{1-\eta}\right)$ away.

Hence, the distance of all chosen noise-points from $\mu$ will be within the $\text{length}\left(I_{1-\eta}\right)$.

Moreover, the interval of minimum length with $(1-\eta-t)(1-\eta)$ fraction of $S$ will contain at least $1-3\eta-t$ fraction of $S_{D}$.

Hence, we can bound the error of $mean(\widetilde{S})$ by controlling the sources of error.

All chosen noise points are within $\text{length}\left(I_{1-\eta}\right)$, and there are atmost $\eta$ of them, hence the maximum error can be $\eta\text{length}\left(I_{1-\eta}\right)$.

Next, the mean of chosen good points will converge to the mean of the conditional distribution. i.e. points sampled from $D$ but conditioned to lie in the minimum length interval. The variance of these random variables is upper bounded using Lemma 10.

To control the distance between the mean($E(X)$ and the conditional mean($E(X|A)$), where $A$ is the event that a sample $x$ is in the chosen interval. We know that $P(A)\geq 1-3\eta-t$, hence, using Lemma 3.11, we get that there exists a constant $C_{13}$ such that,

Hence, with probability at least $1-\delta/3$, the mean of $\widetilde{S}$ will be within

Taking union-bound over all conditioning statements, and upper bounding, $\eta$ with $\epsilon+\sqrt{\frac{\log(3/\delta)}{2n}}$, we recover the statement of the lemma.

K.3 The case when p>1𝑝1p>1

To prove the case for $p>1$, we use a series of lemmas. Lemma 15 proves that the outlier filtering constrains the points in a ball around the true mean. Lemma 17 controls the error in the mean and covariance the true distribution after outlier filtering ($\widetilde{D}$). Lemma 18 controls the error for the mean of $\widetilde{S}$ when projected onto the bottom span of the covariance matrix $\Sigma_{\widetilde{S}}$.

Pick $n$ orthogonal directions $v_{1},v_{2},\ldots,v_{n}$, and use method for one-dimensions, and using union bound, we can recover the result. ∎

Next, we prove the case when $p>1$. Firstly, we prove that after the outlier step,

After the outlier removal step, there exists universal constants $C_{11}>0$ such that with probability at least $1-\delta$, every remaining point $x$ satisfies,

where $r_{1}^{*}=C_{10}\frac{C_{4}^{\frac{1}{4}}\sqrt{p\|\Sigma\|_{2}}}{\eta^{\frac{1}{4}}}$ and $r_{2}^{*}=C_{1}C_{4}^{\frac{1}{4}}\sqrt{p\|\Sigma\|_{2}}(\eta+t)^{\frac{3}{4}}+C_{2}\sqrt{p\|\Sigma\|_{2}}(\eta+t)^{\frac{1}{2}}\sqrt{\frac{\log(1/\delta)}{n}}$ and $t=C_{8}\sqrt{\frac{1}{n}\log\left(\frac{np}{\delta}\right)}$. Here $\eta\leq\epsilon+\sqrt{\frac{\log(1/\delta)}{2n}}$ is the fraction of samples corrupted.

Let $\widetilde{S}$ be the set of points chosen after the outlier filtering. Let $\widetilde{S}_{D}$ be set of good points chosen after the outlier filtering. Let $\widetilde{S}_{N}$ be the set of bad points chosen after the outlier filtering.

Using VC theory we know that for every closed ball $\mathcal{B}(\mu,r)=\{x|\|x-\mu\|_{2}\leq r\}$, there exists a constant $C_{9}$ such that with probability at least $1-\delta$

Let $B^{*}=B(\mu,r^{*})$ for $r_{1}^{*}=C_{10}\frac{C_{4}^{\frac{1}{4}}}{(\eta)^{\frac{1}{4}}}\sqrt{p\|\Sigma\|_{2}}$. Then, we claim that

To see this, suppose we have some $x\in D$. Let $z=x-\mu$. Let $z_{i}=z^{T}v_{i}$ for some orthogonal directions $v_{1},v_{2},\ldots,v_{p}$. Let $Z^{2}=\sum z_{i}^{2}=\|z\|_{2}^{2}$.

Now, $E(Z^{4})\leq p^{2}\max_{i}E(z_{i}^{4})\leq C_{4}p^{2}\|\Sigma\|_{2}^{2}$. Plugging this in the above, we have that $P(x\in B^{*}|x\sim D)\geq 1-\eta$.

Hence, we have that $P(x\in B^{*}|x\in_{u}S_{D})\geq 1-\eta-t_{2}$.

Using Lemma 15, we have that at least $(1-\eta-t_{2})$ fraction of good points are $r_{1}^{*}+r_{2}^{*}$ away from $a$. Hence, we have that the minimum radius of the ball containing all the $(1-\eta-t_{2})(1-\eta)$ has a radius of atmost $r_{1}^{*}+r_{2}^{*}$, which when combined with the triangle inequality recovers the statement of lemma.

As before, let $\widetilde{S}$ be the set of points after outlier filtering. Let $\mu_{\widetilde{S}}=mean(\widetilde{S})$, $\mu_{\widetilde{S}_{D}}=mean(\widetilde{S}_{D})$, $\mu_{\widetilde{S}_{N}}=mean(\widetilde{S}_{N})$.

Let $\widetilde{S}_{D}$ be the set of clean points remaining after the outlier filtering. Then, with probability at least $1-\delta$, we have that

We first prove the bounds on the mean shift.

Control of B. We use Lemma 9 on $X=x^{T}\frac{\mu_{\widetilde{D}}-\mu}{\|\mu_{\widetilde{D}}-\mu\|_{2}}$ for $x\sim D$, and $A$ be the event that $x$ is not removed by the outlier filtering.

Control of A. Using Lemma 10,we have that $\|\Sigma_{\widetilde{D}}\|_{2}\leq(1+2(\eta+t_{2}))\|\Sigma\|_{2}$. Now, we use Bernstein’s inequality . Lemma 12 with $R=C(r_{1}^{*}+2r_{2}^{*}+B)$, we get that, with probability at least $1-\delta$,

Next, we prove the bound for covariance matrix.

To control $\|\Sigma_{\widetilde{S}_{D}}-\Sigma_{\widetilde{D}}\|_{2}$, we use Bernstein’s inequality, with $Z_{k}=\frac{(x_{k}-\mu_{\widetilde{D}})(x_{k}-\mu_{\widetilde{D}})^{T}-\Sigma_{\widetilde{D}}}{n}$. From, Lemma 16, we know that the points are constrained in a ball. Plugging this into Lemma 12,

where $R^{2}=C\left(r_{1}^{*^{2}}+r_{2}^{*^{2}}+B^{2}\right)$. ∎

Let $W$ be the bottom $p/2$ principal components of the covariance matrix after filtering $\Sigma_{\widetilde{S}}$. Then there exists a universal constant $C>0$ such that with probability at least $1-\delta$, we have that

where $\delta_{\mu}=\mu_{\widetilde{S}_{N}}-\mu_{\widetilde{S}_{D}}$, $P_{W}$ is the projection matrix on the bottom $p/2$-span of $\Sigma_{\widetilde{S}}$, $\beta(n,\delta)$ is as defined in Lemma 17 and $\gamma(n,\delta)=\left(\eta^{\frac{1}{2}}+(\eta+t)^{5/2}+\eta(\eta+t)\frac{\log(1/\delta)}{n}\right)$

where $t=C_{8}\sqrt{\frac{1}{n}\log\left(\frac{np}{\delta}\right)}$.

Using that $W$ is the space spanned by the bottom $p/2$ eigenvectors of $\Sigma_{\widetilde{S}}$ and $P_{W}$ is corresponding projection operator, we have that:

Following some algebraic manipulation in , we get that,

Having established all required results, we are now ready to prove Lemma 1. We first present a result for general mean estimation. The proof of Lemma 1 then follows directly from this result.

where $\eta=\epsilon+\sqrt{\frac{\log(p)\log(\log p/\delta)}{2n}}$ and $t_{2}=\sqrt{\frac{p\log(p)\log(n/(p\delta))}{n}}$.

We divide $n$ samples into $\left\lfloor\log(p)\right\rfloor$ different sets. We choose the first set and keep that as our active set of samples. We run our outlier filtering on this set, and let the remaining samples after the outlier filtering be $\widetilde{S}_{D}$. By orthogonality of subspaces spanned by eigenvectors, coupled with triangle inequality and contraction of projection operators, we have that

where $V$ is the span of the top $p/2$ principal components of $\Sigma_{\widetilde{S}}$ and where $\widehat{\mu}_{V}$ is the mean vector of returned by the running the algorithm on the reduced dimensions $dim(V)=p/2$. From Lemma 18, both $\beta(n,\delta)$ and $\gamma(n,\delta)$ are monotonically increasing in the dimension; moreover the upper bound in Lemma 17 is also monotonically increasing in the dimension $p$, hence, the error at each step of the algorithm can be upper bounded by error incurred when running on dimension $p$, with $n/\log(p)$ samples, and probability of $\delta/\log p$. Hence, the overall error for the recursive algorithm can be upper bounded as,

Combining Lemma 17 and Lemma 18 which are instantiated for $n/\log p$ samples and probability $\delta/\log p$, we get,