Robust Estimators in High Dimensions without the Computational Intractability

Introduction

A central goal of machine learning is to design efficient algorithms for fitting a model to a collection of observations. In recent years, there has been considerable progress on a variety of problems in this domain, including algorithms with provable guarantees for learning mixture models [FOS08, KMV10, MV10, BS10, HK13], phylogenetic trees [CGG02, MR05], HMMs [AHK12], topic models [AGM12, AGHK13], and independent component analysis [AGMS15]. These algorithms crucially rely on the assumption that the observations were actually generated by a model in the family. However, this simplifying assumption is not meant to be exactly true, and it is an important direction to explore what happens when it holds only in an approximate sense. In this work, we study the following family of questions:

This is a natural formalization of the problem of designing robust and efficient algorithms for distribution estimation. We refer to it as (proper) agnostic distribution learning and we refer to the samples as being $\varepsilon$-corrupted. This family of problems has its roots in many fields, including statistics, machine learning, and theoretical computer science. Within computational learning theory, it is related to the agnostic learning model of Haussler [Hau92] and Kearns, Schapire, and Sellie [KSS94], where the goal is to learn a labeling function whose agreement with some underlying target function is close to the best possible, among all functions in some given class. In the even more challenging malicious noise model [Val85, KL93], an adversary is allowed to corrupt both the labels and the samples. A major difference with our setting is that these models apply to supervised learning problems, while here we will work in an unsupervised setting.

Within statistics and machine learning, inference problems like Question 1.1 are often termed “estimation under model misspecification.” The usual prescription is to use the maximum likelihood estimator [Hub67, Whi82], which is unfortunately hard to compute in general. Even ignoring computational considerations, the maximum likelihood estimator is only guaranteed to converge to the distribution $P^{\prime}$ in $\mathcal{D}$ that is closest (in Kullback-Leibler divergence) to the distribution from which the observations are generated. This is problematic because such a distribution is not necessarily close to $P$ at all.

A branch of statistics – called robust statistics [HR09, HRRS86] – aims to tackle questions like the one above. The usual formalization is in terms of breakdown point, which (informally) is the fraction of observations that an adversary would need to control to be able to completely corrupt an estimator. In low-dimensions, this leads to the prescription that one should use the empirical median instead of the empirical mean to robustly estimate the mean of a distribution, and interquartile range for robust estimates of the variance. In high-dimensions, the Tukey median [Tuk75] is a high-dimensional analogue of the median that, although provably robust, is hard to compute [JP78]. Similar hardness results have been shown [Ber06, HM13] for essentially all known estimators in robust statistics.

Is high-dimensional agnostic distribution learning even possible, algorithmically? The difficulty is that corruptions are often hard to detect in high dimensions, and could bias the natural estimator by dimension-dependent factors. In this work, we study agnostic distribution learning for a number of fundamental classes of distributions: $(1)$ a single Gaussian, $(2)$ a product distribution on the hypercube $\{0,1\}^{d}$, $(3)$ mixtures of two product distributions (under a natural balancedness condition), and $(4)$ mixtures of $k$ Gaussians with spherical covariances. Prior to our work, all known efficient algorithms (e.g., [LT15, BD15]) for these classes required the error guarantee, $f(\varepsilon,d)$, to depend polynomially in the dimension $d$. Hence, previous efficient estimators could only tolerate at most a $1/\operatorname{poly}(d)$ fraction of errors. In this work, we obtain the first efficient algorithms for the aforementioned problems, where $f(\varepsilon,d)$ is completely independent of $d$ and depends polynomially (often, nearly linearly) in the fraction $\varepsilon$ of corrupted samples. Our work is just a first step in this direction, and there are many exciting questions left to explore.

2 Our Techniques

All of our algorithms are based on a common recipe. The first question to address is the following: Even if we were given a candidate hypothesis $P^{\prime}$, how could we test if it is $\varepsilon$-close in total variation distance to $P$? The usual way to certify closeness is to exhibit a coupling between $P$ and $P^{\prime}$ that marginally samples from both distributions, where the samples produced from each agree with probability $1-\varepsilon$. However, we have no control over the process by which samples are generated from $P$, in order to produce such a coupling. And even then, the way that an adversary decides to corrupt samples can introduce complex statistical dependencies.

Unfortunately, in our agnostic setting, we cannot afford for (1) to depend on the dimension $d$ at all. Any such dependence would appear in the error guarantee of our algorithm. Instead, the starting point of our algorithms is a notion of parameter distance that satisfies

Given our notion of parameter distance satisfying (2), our main ingredient is an efficient method for robustly estimating the parameters. We provide two algorithmic approaches which are based on similar principles. Our first approach is faster, requiring only approximate eigenvalue computations. Our second approach relies on convex programming and achieves slightly better sample complexity, in some cases matching the information-theoretic limit. Notably, either approach can be used to give all of our concrete learning applications with nearly identical sample complexity and error guarantees. In what follows, we specialize to the problem of robustly learning the mean $\mu$ of a Gaussian whose covariance is promised to be the identity, which we will use to illustrate how both approaches operate. We emphasize that what is needed to learn the parameters in more general settings requires many additional ideas.

Our second algorithmic approach relies on convex programming. Here, instead of rejecting corrupted samples, we compute appropriate weights $w_{i}$ for the samples $X_{i}$, such that the weighted empirical average $\widehat{\mu}_{w}=\mathop{\textstyle\sum}_{i=1}^{N}w_{i}X_{i}$ is close to $\mu$. We work with the convex set:

We prove that any set of weights in $\mathcal{C}_{\delta}$ yields a good estimate $\widehat{\mu}_{w}=\mathop{\textstyle\sum}_{i=1}^{N}w_{i}X_{i}$ in the obvious way. The catch is that the set $\mathcal{C}_{\delta}$ is defined based on $\mu$, which is unknown. Nevertheless, it turns out that we can use the same type of spectral arguments that underlie the filtering approach to design an approximate separation oracle for $\mathcal{C}_{\delta}$. Combined with standard results in convex optimization, this yields an algorithm for robustly estimating $\mu$.

The third and final ingredient is some new concentration bounds. In both of the approaches above, at best we are hoping that we can remove all of the corrupted points and be left with only the uncorrupted ones, and then use standard estimators (e.g., the empirical average) on them. However, an adversary could have removed an $\varepsilon$-fraction of the samples in a way that biases the empirical average of the remaining uncorrupted samples. What we need are concentration bounds that show for sufficiently large $N$, for samples $X_{1},X_{2},\ldots,X_{N}$ from a Gaussian with mean $\mu$ and identity covariance, that every set of $(1-\varepsilon)N$ samples produces a good estimate for $\mu$. In some cases, we can derive such concentration bounds by appealing to known concentration inequalities and taking a union bound. However, in other cases (e.g., concentration bounds for degree-two polynomials of Gaussian random variables) the existing concentration bounds are not strong enough, and we need other arguments to prove that every set of $(1-\varepsilon)N$ samples produces a good estimate.

3 Our Results

We give the first efficient algorithms for agnostically learning several important distribution classes with dimension-independent error guarantees. Our first main result is for a single arbitrary Gaussian with mean $\mu$ and covariance $\Sigma$, which we denote by $\operatorname{\mathcal{N}}(\mu,\Sigma)$. In the previous subsection, we described our convex programming approach for learning the mean vector when the covariance is promised to be the identity. A technically more involved version of the technique can handle the case of zero mean and unknown covariance. More specifically, consider the following convex set, where $\Sigma$ is the unknown covariance matrix and $\left\|\cdot\right\|_{F}$ is the Frobenius norm:

We design an approximate separation oracle for this unknown convex set, by analyzing the spectral properties of the fourth moment tensor of a Gaussian. Combining these two intermediate results, we obtain our first main result (below). Throughout this paper, we will abuse notation and write $N\geq\widetilde{\Omega}(f(d,\varepsilon,\tau))$ when referring to our sample complexity, to signify that our algorithm works if $N\geq Cf(d,\varepsilon,\tau)\mbox{polylog}(f(d,\varepsilon,\tau))$ for a large enough universal constant $C$.

We can alternatively establish Theorem 1.2 via our filtering technique. See Section 5. In the first version of our paper, our analysis required $N\gtrsim d^{3}\log^{2}(1/\tau)/\varepsilon^{2}$ samples. In [DKK+17], we showed that a simple adaptation of our algorithm and analysis achieves the improved sample complexity above, which is information-theoretically optimal up to logarithmic factors. We have incorporated this modification (along with the analysis) into this version of the paper, for the sake of completeness.

For the sake of simplicity in the presentation, we did not make an effort to optimize the sample complexity of our robust estimators in the above setting. We note that methods similar to the analysis of the Gaussian setting can lead to near-optimal sample complexity in this setting as well. We also remark that for the case of balanced binary product distributions, our algorithm achieves an error of $O(\varepsilon\sqrt{\log(1/\varepsilon)}).$

Interestingly enough, the above two distribution classes are trivial to learn in the noiseless case, but in the agnostic setting the learning problem turns out to be surprisingly challenging. Using additional ideas, we are able to generalize our agnostic learning algorithms to mixtures of the above classes under some natural conditions. We note that even in the noiseless case, learning mixtures of the above families is non-trivial. First, we study $2$-mixtures of $c$-balanced products, which stipulates that the coordinates of the mean vector of each component are in the range $(c,1-c)$. We prove:

This generalizes the algorithm of Freund and Mansour [FM99] to the agnostic setting. An interesting open question is to improve the $\varepsilon$-dependence in the above bound to (nearly) linear, or to remove the assumption of balancedness and obtain an agnostic algorithm for mixtures of two arbitrary product distributions.

Finally, we give an agnostic learning algorithm for mixtures of spherical Gaussians.

In total, these results give new robust and computationally efficient estimators for several well-studied distribution learning problems that can tolerate a constant fraction of errors independent of the dimension. This points to an interesting new direction of making robust statistics algorithmic. The general recipe we have developed here gives us reason to be optimistic about many other problems in this domain.

4 Discussion and Related Work

Our results fit in the framework of density estimation and parameter learning which are both classical problems in statistics with a rich history (see e.g., [BBBB72, DG85, Sil86, Sco92, DL01]). While these problems have been studied for several decades by different communities, the computational complexity of learning is still not well understood, even for some surprisingly simple distribution families. Most textbook estimators are hard to compute in general, especially in high-dimensional settings. In the past few decades, a rich body of work within theoretical computer science has focused on designing computationally efficient distribution learning algorithms. In a seminal work, Kearns, Mansour, Ron, Rubinfeld, Schapire, and Sellie [KMR+94] initiated a systematic investigation of the computational complexity of distribution learning. Since then, efficient learning algorithms have been developed for a wide range of distributions in both low and high-dimensions [Das99, FM99, AK01, VW02, CGG02, MR05, BV08, KMV10, MV10, BS10, DDS12, CDSS13, DDO+13, CDSS14a, CDSS14b, HP15, ADLS17, DDS15b, DDKT16, DKS16b, DKS16a].

We will be particularly interested in efficient learning algorithms for mixtures of high-dimensional Gaussians and mixtures of product distributions, as this is the focus of our algorithmic results in the agnostic setting. In a pioneering work, Dasgupta [Das99] introduced the problem of parameter estimation of a Gaussian mixture to theoretical computer science, and gave the first provably efficient algorithms under the assumption that the components are suitably well-separated. Subsequently, a number of works improved these separation conditions [AK01, VW02, BV08] and ultimately removing them entirely [KMV10, MV10, BS10]. In another line of work, Freund and Mansour [FM99] gave the first polynomial time algorithm for properly learning mixtures of two binary product distributions. This algorithm was substantially generalized to phylogenetic trees [CGG02] and to mixtures of any constant number of discrete product distributions [FOS08]. Given the vast body of work on high-dimensional distribution learning, there is a plethora of problems where one could hope to reconcile robustness and computational efficiency. Thus far, the only setting where robust and efficient algorithms are known is on one-dimensional distribution families, where brute-force search or some form of polynomial regression often works. In contrast, essentially nothing is known about efficient agnostic distribution learning in the high-dimensional setting that we study here.

Question 1.1 also resembles learning in the presence of malicious errors [Val85, KL93]. There, an algorithm is given samples from a distribution along with their labels according to an unknown target function. The adversary is allowed to corrupt an $\varepsilon$-fraction of both the samples and their labels. A sequence of works studied the problem of learning a homogeneous halfspace with malicious noise in the setting where the underlying distribution is a Gaussian [Ser01, Ser03, KLS09], culminating in the work of Awasthi, Balcan, and Long [ABL17], who gave an efficient algorithm that finds a halfspace with agreement $O(\varepsilon)$. There is no direct connection between their problem and ours, especially since one is a supervised learning problem and the other is unsupervised. We note however that there is an interesting technical parallel in that the work [KLS09] also uses spectral methods to detect outliers. Both their work and our algorithm for agnostically learning the mean are based on the intuition that an adversary can only substantially bias the empirical mean if the corruptions are correlated along some direction. More specifically, [KLS09] produce a “hard” filter which leads to errors that scale logarithmically with the dimension, even in a weaker corruption model than ours. Our algorithms need to handle many significant conceptual and technical complications that arise when working with higher moments or other distribution families.

Another connection is to the work on robust principal component analysis (PCA). PCA is a transformation that (among other things) is often justified as being able to find the affine transformation $Y=\Sigma^{-1/2}(X-\mu)$ that would place a collection of Gaussian random variables in isotropic position. One can think of our results on agnostically learning a Gaussian as a type of robust PCA that tolerates gross corruptions, where entire samples are corrupted. This is different than other variants of the problem where random sets of coordinates of the points are corrupted [CLMW11], or where the uncorrupted points were assumed to lie in a low-dimensional subspace to begin with [ZL14, LMTZ15]. Finally, Brubaker [Bru09] studied the problem of clustering samples from a well-separated mixture of Gaussians in the presence of adversarial noise. The goal of [Bru09] was to separate the Gaussian components from each other, while the adversarial points are allowed to end up in any of clusters. Our work is orthogonal to [Bru09], since even if such a clustering is given, the problem still remains to estimate the parameters of each component.

5 Concurrent and Subsequent Work

After the initial publication of our results [DKK+16], there has been a flurry of recent work on robust high-dimensional estimation. Diakonikolas, Kane, and Stewart [DKS16c] studied the problem of learning the parameters of a graphical model in the presence of noise, when given its graph theoretic structure. Charikar, Steinhardt, and Valiant [CSV17] developed algorithms that can tolerate a fraction of corruptions greater than a half, under the weaker goal of outputting a small list of candidate hypotheses that contains a parameter set close to the true values. Balakrishnan, Du, Li, and Singh [Li17, DBS17, BDLS17] studied sparse mean and covariance estimation in the presence of noise obtaining computationally efficient robust algorithms with sample complexity sublinear in the dimension. Diakonikolas, Kane, and Stewart [DKS17] proved statistical query lower bounds providing evidence that the error guarantees of our robust mean and covariance estimation algorithms are best possible, within constant factors, for efficient algorithms. In a subsequent paper [DKK+17], we obtained improved bounds on the sample complexity of our algorithms, which are optimal up to polylogarithmic factors. For the sake of completeness, we include these improved sample bounds in the present version of this paper. In the same work [DKK+17], we showed that our algorithmic approach easily extends to obtain dimension-independent robustness guarantees under much weaker distributional assumptions, and gave a practical demonstration of the efficacy of our robust algorithms on both real and synthetic data.

Since the initial submission of the journal version of this paper, there has been a substantial amount of work on robust high-dimensional estimation in a variety of settings. Diakonikolas, Kane, and Stewart [DKS18a] studied PAC learning of geometric concept classes (including low-degree polynomial threshold functions and intersections of halfspaces) in the same corruption model as ours, obtaining the first dimension-independent error guarantees for these classes. Steinhardt, Charikar, and Valiant [SCV18] focused on deterministic conditions of a dataset which allow robust estimation to be possible. In our initial publication, we gave explicit deterministic conditions in various settings; by focusing directly on this goal, [SCV18] somewhat relaxed some of these assumptions. Meister and Valiant [MV17] studied learning in a crowdsourcing model, where the fraction of honest workers may be very small (similar to [CSV17]). Qiao and Valiant [QV18] considered robust estimation of discrete distributions in a setting where we have several sources (a fraction of which are adversarial) who each provide a batch of samples. A number of simultaneous works [KSS18, HL18, DKS18b] investigated robust mean estimation in even more general settings, and apply their techniques to learning mixtures of spherical Gaussians under minimal separation conditions. Finally, several concurrent results study robustness in supervised learning tasks [PSBR18, KKM18, DKK+18], including regression and SVM problems. Despite all of this rapid progress, there are still many interesting theoretical and practical questions left to explore.

6 Organization

The structure of this paper is as follows: In Section 2, we introduce basic notation and a number of useful facts that will be required throughout the paper, as well as the formal definition of our adversary model. In Section 3, we discuss several natural approaches to high-dimensional agnostic learning, all of which lose polynomial factors that depend on the dimension, in terms of their error guarantee.

The main body of the paper is in Sections 4–8. Sections 4 and 6 illustrate our convex programming framework, while Sections 5, 7, and 8 illustrate our filter framework. More specifically, in Sections 4 and 5, we analyze the setting of a single Gaussian with unknown mean and unknown covariance, using our convex programming and filter frameworks, respectively. In Section 6, we generalize the convex programming method to obtain an agnostic algorithm for mixtures of spherical Gaussians with unknown means. In Section 7, we apply our filter techniques to a binary product distribution, and generalize these in Section 8 to obtain an agnostic learning algorithm for a mixture of two balanced binary product distributions.

We note that for some of the more advanced applications of our frameworks, the technical details can get in the way of the fundamental ideas. For the reader who is interested in seeing the details of our most basic application of the convex programming framework, we recommend reading the case a Gaussian with unknown mean, in Section 4.3. Similarly, for the filter framework, we suggest either the Gaussian with unknown mean in Section 5.1 or the balanced product distribution in Section 7.1.

Preliminaries

Throughout this paper, if $v$ is a vector, we will let $\|v\|_{2}$ denote its Euclidean norm. If $M$ is a matrix, we will let $\|M\|_{2}$ denote its spectral norm, and $\|M\|_{F}$ denote its Frobenius norm. We will also let $\preceq$ and $\succeq$ denote the PSD ordering on matrices. For a discrete distribution $P$, we will denote by $P(x)$ the probability mass at point $x$. For a continuous distribution, let it denote the probability density function at $x$. Let $S$ be a multiset over $\{0,1\}^{d}$ . We will write $X\in_{u}S$ to denote that $X$ is drawn from the empirical distribution defined by $S$. Throughout the paper, we let $\otimes$ denote the Kronecker product of matrices.

As a measure of distance between distributions, we will use the notion of total variation distance:

2 Types of Adversaries

In this paper, we will consider a powerful model for agnostic distribution learning that generalizes many other existing models. The standard setup involves an oblivious adversary who chooses a distribution that is close in total variation distance to an unknown distribution in some class $\mathcal{D}$.

Given $\varepsilon>0$, and a class of distributions $\mathcal{D}$, the full adversary operates as follows: The algorithm specifies some number of samples $m$. The adversary generates $m$ samples $X_{1},X_{2},\ldots,X_{m}$ from some (unknown) distribution $D\in\mathcal{D}$. It then draws $m^{\prime}$ from an appropriate distribution. This distribution is allowed to depend on $X_{1},X_{2},\ldots,X_{m}$, but when marginalized over the $m$ samples satisfies $m^{\prime}\sim\mbox{Bin}(m,\varepsilon)$. The adversary is allowed to inspect the samples, removes $m^{\prime}$ of them, and replaces them with arbitrary points. The set of $m$ points is given (in any order) to the algorithm.

We rely on the following well-known fact:

Now we can describe how the full adversary can corrupt samples from $D$ to get samples distributed according to $P$.

The full adversary can simulate any oblivious adversary.

We draw $m$ samples $X_{1},X_{2},\ldots,X_{m}$ from $D$. We delete each sample $X_{i}$ independently with probability $\varepsilon_{2}$ and replace it with an independent sample from $N_{2}$. This gives a set of samples $Y_{1},Y_{2},\ldots,Y_{m}$ that are independently sampled from $(1-\varepsilon_{2})D+\varepsilon_{2}N_{2}$. Since the distributions $(1-\varepsilon_{1})P+\varepsilon_{1}N_{1}$ and $(1-\varepsilon_{2})D+\varepsilon_{2}N_{2}$ are identical, we can couple them to independent samples $Z_{1},Z_{2},\ldots,Z_{m}$ from $(1-\varepsilon_{1})P+\varepsilon_{1}N_{1}$. Now each sample $Z_{i}$ that came from $N_{1}$, we can delete and replace with an independent sample from $P$. The result is a set of samples that are independently sampled from $P$ where we have made $m^{\prime}$ edits and marginally $m^{\prime}\sim\mbox{Bin}(m,\varepsilon_{1}+\varepsilon_{2})$, although $m^{\prime}$ has and needs to have some dependence on the original samples from $D$. ∎

The challenge in working with the full adversary is that even the samples that came from $D$ can have biases. The adversary can now choose how to remove uncorrupted points in a careful way so as to compensate for certain other biases that he introduces using the corrupted points.

Throughout this paper, we will make use of the following notation and terminology:

We say a set of samples $X_{1},X_{2},\ldots,X_{m}$ is an $\varepsilon$-corrupted set of samples generated by the oblivious (resp. full) adversary if it is generated by the process described above in the definition of the oblivious (resp. full) adversary. If it was generated by the full adversary, we let $G\subseteq[m]$ denote the indices of the uncorrupted samples, and we let $E\subseteq[m]$ denote the indices of the corrupted samples.

In this paper, we will give a number of algorithms for agnostic distribution learning that work in the full adversary model. In our analysis, we will identify a set of events that ensure the algorithm succeeds and will bound the probability that any of these events does not occur when $m$ is suitably large. We will often explicitly invoke the assumption that $|E|\leq 2\varepsilon m$. We can do this even though the number of points that are corrupted is itself a random variable, because by the Chernoff bound, as long as $m\geq O\left(\frac{\log 1/\tau}{\varepsilon}\right)$, we know that $|E|\leq 2\varepsilon m$ holds with probability at least $1-O(\tau)$. Thus, making the assumption that $|E|\leq 2\varepsilon m$ costs us an additional additive $O(\tau)$ term in our union bound, when bounding the failure probability of our algorithms.

3 Distributions of Interest

One object of study in this paper is the Gaussian (or Normal) distribution.

A Gaussian distribution $\mathcal{N}(\mu,\Sigma)$ with mean $\mu$ and covariance $\Sigma$ is the distribution with probability density function

We will also be interested in binary product distributions.

A (binary) product distribution is a probability distribution over $\{0,1\}^{d}$ whose coordinate random variables are independent. Note that a binary product distribution is completely determined by its mean vector.

We will also be interested in mixtures of such distributions.

A mixture $P$ of distributions $P_{1},\dots,P_{k}$ with mixing weights $\alpha_{1},\dots,\alpha_{k}$ is the distribution defined by

where $\alpha_{j}\geq 0$ for all $j$ and $\sum_{j\in[k]}\alpha_{j}=1$.

4 Bounds on TV Distance

The Kullback-Liebler divergence (also known as relative entropy, information gain, or information divergence) is a well-known measure of distance between two distributions.

The primary interest we have in this quantity is the fact that (1) the KL divergence between two Gaussians has a closed form expression, and (2) it can be related (often with little loss) to the total variation distance between the Gaussians. The first statement is expressed in the fact below:

Let $\operatorname{\mathcal{N}}(\mu_{1},\Sigma_{1})$ and $\operatorname{\mathcal{N}}(\mu_{2},\Sigma_{2})$ be two Gaussians such that $\det(\Sigma_{1}),\det(\Sigma_{2})\neq 0$. Then

The second statement is encapsulated in the well-known Pinsker’s inequality:

With this we can show the following two useful lemmas, which allow us to relate parameter distance between two Gaussians to their total variation distance. The first bounds the total variation distance between two Gaussians with identity covariance in terms of the Euclidean distance between the means:

In the case where $\Sigma_{1}=\Sigma_{2}=I$, (3) simplifies to

Pinsker’s inequality (Theorem 2.12) then implies that

The second bounds the total variation distance between two mean zero Gaussians in terms of the Frobenius norm of the difference between their covariance matrices:

Let $\delta>0$ be sufficiently small. Let $\Sigma_{1},\Sigma_{2}$ such that $\|I-\Sigma_{2}^{-1/2}\Sigma_{1}\Sigma_{2}^{-1/2}\|_{F}=\delta$. Then,

Let $M=\Sigma_{2}^{-1/2}\Sigma_{1}\Sigma_{2}^{-1/2}$. Then (3) simplifies to

Since both terms in the last line are rotationally invariant, we may assume without loss of generality that $M$ is diagonal. Let $M=\mbox{diag}(1+\lambda_{1},\ldots,1+\lambda_{d})$. Thus, the KL divergence between the two distributions is given exactly by $\frac{1}{2}\sum_{i=1}^{d}\left(\lambda_{i}-\log(1+\lambda_{i})\right),$ where we are guaranteed that $(\sum_{i=1}^{d}\lambda_{i}^{2})^{1/2}=\delta$. By the second order Taylor approximation to $\ln(1+x)$, for $x$ small, we have that for $\delta$ sufficiently small,

Our algorithm for agnostically learning an arbitrary Gaussian will be based on solving two intermediate problems: (1) We are given samples from $\operatorname{\mathcal{N}}(\mu,I)$ and our goal is to learn $\mu$, and (2) We are given samples from $\operatorname{\mathcal{N}}(0,\Sigma)$ and our goal is to learn $\Sigma$. The above bounds on total variation distance will allow us to conclude that our estimate is close in total variation distance to the unknown Gaussian distribution in each of the two settings.

We note the following folklore sample complexity bounds for learning a Gaussian in the non-agnostic setting.

$N=\Theta\left(\frac{d+\log(1/\tau)}{\varepsilon^{2}}\right)$ samples are both necessary and sufficient to learn a $d$-dimensional Gaussian with unknown mean and known covariance to total variation distance $\varepsilon$ with probability $1-\tau$.

$N=\Theta\left(\frac{d^{2}+\log(1/\tau)}{\varepsilon^{2}}\right)$ samples are both necessary and sufficient to learn a $d$-dimensional Gaussian with unknown mean and covariance to total variation distance $\varepsilon$ with probability $1-\tau$.

We will also need the following lemma bounding the total variation distance between two product distributions:

Let $P,Q$ be binary product distributions with mean vectors $p,q\in(0,1)^{d}.$ We have that

The elementary inequality $2\sqrt{ab}=a+b-(\sqrt{a}-\sqrt{b})^{2},$ $a,b>0,$ gives that

where the last inequality follows from the union bound. ∎

5 Additional Concentration Lemmata

In this section, we list a number of standard concentration inequalities for nice random variables which we will frequently use throughout this paper. The proofs of these results are standard and omitted, see e.g., [Ver10] for a more thorough treatment of these results.

The first is a Chernoff bound for bounded random variables.

Let $Z_{1},\ldots,Z_{d}$ be independent random variables with $Z_{i}$ supported on $[a_{i},b_{i}].$ Let $Z=\sum_{i=1}^{d}Z_{i}$. Then for any $T>0$,

We will also require the following tail bounds for Gaussians and quadratic forms of Gaussians:

By standard union bound arguments (see e.g. [Ver10]), we obtain the following concentration results for the empirical mean and covariance of a set of Gaussian vectors:

Let $n$ be a positive integer. Let $D$ be a sub-gaussian distribution with mean and covariance $I$. Let $Y_{i}\sim D$ be independent, for $i=1,\ldots,n$. Then, there exist universal constants $A,B>0$ so that for all $t>0$, we have

With the same setup as in Lemma 2.21, there exist universal constants $A,B>0$ so that for all $t>0$, we have

6 Agnostic Hypothesis Selection

Several of our algorithms will return a polynomial-sized list of hypotheses at least one of which is guaranteed to be close to the target distribution. Usually (e.g., in a non-agnostic setting), one could use a polynomial number of additional samples to run a tournament to identify the candidate hypothesis that is (roughly) the closest to the target distribution. In the discussion that follows, we will refer to these additional samples as test samples. Such hypothesis selection algorithms have been extensively studied [Yat85, DL96, DL97, DL01, DK14, AJOS14, SOAJ14, DDS15a, DDS15b]. Unfortunately, against a strong adversary we run into a serious technical complication: the training samples and test samples are not necessarily independent. Moreover even if we randomly partition our samples in training and test, a priori there are an unbounded set of possible hypotheses that the training phase could output, and when we analyze the tournament we cannot condition on the list of hypotheses and assume that the test samples are sampled anew. Our approach is to require our original algorithm to return only hypotheses from some finite set of possibilities, and as we will see this mitigates the problem.

As a simple corollary of the agnostic tournament, observe that this allows us to do agnostic learning without knowing the precise error rate $\varepsilon$. Throughout the paper, we assume the algorithm knows $\varepsilon$, and guarantees that the output will have error which is at most $O(f(\varepsilon))$. However, if the algorithm is not given this information, and instead given an $\eta$ and asked to return something with error at most $O(f(\varepsilon+\eta))$, we may simply grid over $\{\eta,(1+\gamma)\eta,(1+\gamma)^{2}\eta,\ldots,1\}$ (here $\gamma$ is some arbitrary constant that governs a tradeoff between runtime and accuracy), run our algorithm with $\varepsilon$ set to each element in this set, and perform hypothesis selection via Tournament. Then it is not hard to see that we are guaranteed to output something which has error at most $O(f(\varepsilon+(1+\gamma)\eta))$.

Firstly, we randomly choose a subset of $N$ of our samples and a disjoint subset of $C(\log(|\mathcal{M}|)+\log(1/\tau))/\varepsilon^{2}$ of our samples for some sufficiently large $C$. Note that with high probability over our randomization, at most a $2\varepsilon$-fraction of samples from each subset are corrupted. Thus, we may instead consider the stronger adversary who sees a set $S_{1}$ of $N$ independent samples from $\Pi$ and another set, $S_{2}$, of $C(\log(|\mathcal{M}|)+\log(1/\tau))/\varepsilon^{2}$ samples from $\Pi$ and can arbitrary corrupt a $2\varepsilon$-fraction of each, giving sets $S^{\prime}_{1}$,$S^{\prime}_{2}$.

With probability at least $1-\tau/2$ over $S_{1}$, the original algorithm run on $S^{\prime}_{1}$ returns a set $\mathcal{L}$ satisfying the desired properties.

We claim that with probability at least $1-\tau/2$ over the choice of $S_{2}$, we have for any $P,Q\in\mathcal{M}$:

This follows by Chernoff bounds and a union bound over the $|\mathcal{M}|^{2}$ possibilities for $P$ and $Q$. Since the total variation distance between the uniform distributions over $S_{2}$ and $S^{\prime}_{2}$ is at most $2\varepsilon$, we also have for $S^{\prime}_{2}$ that

Therefore, if we throw away any $Q$ in our list for which there is a $P$ in our list such that

Some Natural Approaches, and Why They Fail

Many of the agnostic distribution learning problems that we study are so natural that one would immediately wonder why simpler approaches do not work. Here we detail some other plausible approaches, and what causes them to lose dimension-dependent factors (if they have any guarantees at all!). For the discussion that follows, we note that by Fact 2.13 in order to achieve an estimate that is $O(\varepsilon)$-close in total variation distance (for a Gaussian when $\mu$ is unknown and $\Sigma=I$) it is necessary and sufficient that $\|\hat{\mu}-\mu\|_{2}=O(\varepsilon)$.

One plausible approach for robust mean estimation in high dimensions is to agnostically learn along each coordinate separately. For instance, if our goal is to agnostically learn the mean of a Gaussian with known covariance $I$, we could try to learn each coordinate of the mean separately. But since an $\varepsilon$-fraction of the samples are corrupted, our estimate can be off by $\varepsilon$ in each coordinate and would be off by $\varepsilon\sqrt{d}$ in high dimensions.

Maximum Likelihood

Given a set of samples $X_{1},\ldots,X_{N}$, and a class of distributions $\mathcal{D}$, the maximum likelihood estimator (MLE) is the distribution $F\in\mathcal{D}$ that maximizes $\prod_{i=1}^{N}F(X_{i})$. Equivalently, $F$ minimizes the negative log likelihood (NLL), which is given by

and so the $\mu$ which minimizes $\textrm{NLL}(\operatorname{\mathcal{N}}(\mu,I),X_{1},\ldots,X_{N})$ is the mean of the samples $X_{i}$, since for any set of vectors $v_{1},\ldots,v_{N}$, the average of the $v_{i}$’s is the minimizer of the function $h(x)=\sum_{i=1}^{N}\|v_{i}-x\|_{2}^{2}$. Hence, if an adversary places an $\varepsilon$-fraction of the points at some very large distance, then the estimate for the mean would need to move considerably in that direction. By placing the corruptions further and further away, the MLE can be an arbitrarily bad estimate. That is, even though it is well known [Hub67, Whi82] that the MLE converges to the distribution $F\in\mathcal{D}$ that is closest in KL-divergence to the distribution from which our samples were generated (i.e., after the adversary has added corruptions), $F$ is not necessarily close to the uncorrupted distribution.

Geometric Median

In one dimension, it is well-known that the median provides a provably robust estimate for the mean in a number of settings. The mean of a set of points $a_{1},\ldots,a_{N}$ is the minimizer of the function $f(x)=\sum_{i=1}^{N}(a_{i}-x)^{2}$, and in contrast the median is the minimizer of the function $f(x)=\sum_{i=1}^{N}|a_{i}-x|$. In higher dimensions, there are many natural definitions for the median that generalize the one-dimensional case. The Tukey median is one such notion, but as we discussed it is hard to compute [JP78] and the best known algorithms run in time exponential in $d$. Motivated by this, the geometric median is another high-dimensional notion of a median. It often achieves better robustness than the mean, and can be computed quickly [CLM+16]. The formal definition is:

Unfortunately, this notion of median still incurs an error containing a factor of $O(\sqrt{d})$:

Given a set $S$ of $N=\Omega\left(\frac{d+\log(1/\tau)}{\varepsilon^{2}}\right)$ samples from $\operatorname{\mathcal{N}}(0,I)$, then with probability at least $1-\tau$, there exists a corruption $S^{\prime}$ of $S$, such that:

Agnostically Learning a Gaussian, via Convex Programming

In this section we give a polynomial time algorithm to agnostically learn a single Gaussian up to error $\widetilde{O}(\varepsilon)$. Our approach is based on the following ingredients: First, in Section 4.1, we define the set $S_{N,\varepsilon}$, which will be a key algorithmic object in our framework. In Section 4.2 we give key, new concentration bounds on certain statistics of Gaussians. We will make crucial use of these concentration bounds throughout this section. In Section 4.3 we give an algorithm to agnostically learn a Gaussian with unknown mean and whose covariance is promised to be the identity via convex programming. This will be an important subroutine in our overall algorithm, and it also helps to illustrate our algorithmic approach without many of the additional complications that arise in our later applications. In Section 4.4 we show how to robustly learn a Gaussian with mean zero and unknown covariance again via convex programming. Finally, in Section 4.5 we show how to combine these two intermediate results to get our overall algorithm.

An important algorithmic object for us will be the following set:

For any $\frac{1}{2}>\varepsilon>0$ and any integer $N$, let

and so we see that this set is designed to capture the notion of selecting a set of $(1-2\varepsilon)N$ samples from $N$ samples.

Given $w\in S_{N,\varepsilon}$ we will use the following notation

to denote the total weight on good and bad points respectively. The following facts are immediate from $|E|\leq 2\varepsilon N$ and the properties of $S_{N,\varepsilon}$.

If $w\in S_{N,\varepsilon}$ and $|E|\leq 2\varepsilon N$, then $w_{b}\leq\frac{2\varepsilon}{1-2\varepsilon}$. Moreover, the renormalized weights $w^{\prime}$ on good points given by $w^{\prime}_{i}=\frac{w_{i}}{w_{g}}$ for all $i\in G$, and $w^{\prime}_{i}=0$ otherwise, satisfy $w^{\prime}\in S_{N,4\varepsilon}$.

2 Concentration Inequalities

Throughout this section and in Section 6 we will make use of various concentration bounds on low moments of Gaussian random variables. Some are well-known, and others are new but follow from known bounds and appropriate union bound arguments.

We will also be interested in how well various statistics of Gaussians concentrate around their expectation, when we take the worst-case set of weights in $S_{N,\varepsilon}$. This is more subtle than standard settings such as Lemma 2.21 or Lemma 2.22 because as we take more samples, any fixed statistic (e.g. taking the uniform distribution over the samples) concentrates better but the size of $S_{N,\varepsilon}$ (e.g. the number of sets of $(1-2\varepsilon)N$ samples) grows too. We defer the proofs to Appendix A. The first concerns the behavior of the empirical covariance:

Fix $\varepsilon\leq 1/2$ and $\tau\leq 1$. There is a $\delta_{1}=O(\varepsilon\log 1/\varepsilon)$ such that if $Y_{1},\ldots,Y_{N}$ are independent samples from $\operatorname{\mathcal{N}}(0,I)$ and $N=\Omega\left(\frac{d+\log(1/\tau)}{\delta_{1}^{2}}\right)$, then

A nearly identical argument (Using Hoeffding instead of Bernstein in the proof of Theorem 5.50 in [Ver10]) yields:

Fix $\varepsilon$ and $\tau$ as above. There is a $\delta_{2}=O(\varepsilon\sqrt{\log 1/\varepsilon})$ such that if $Y_{1},\ldots,Y_{N}$ are independent samples from $\operatorname{\mathcal{N}}(0,I)$ and $N=\Omega\left(\frac{d+\log(1/\tau)}{\delta_{2}^{2}}\right)$, then

Note that by Cauchy-Schwarz, this implies:

Fix $\varepsilon$ and $\tau$ as above. There is a $\delta_{2}=O(\varepsilon\sqrt{\log 1/\varepsilon})$ such that if $Y_{1},\ldots,Y_{N}$ are independent samples from $\operatorname{\mathcal{N}}(0,I)$ and $N=\Omega\left(\frac{d+\log(1/\tau)}{\delta_{2}^{2}}\right)$, then

Fix $\tau>0$. Let $X_{1},\ldots,X_{N}\sim\operatorname{\mathcal{N}}(0,I)$. Then, with probability $1-\tau$, we have that $\|X_{i}\|_{2}\leq O\left(\sqrt{d\log(N/\tau)}\right)$ for all $i=1,\ldots,N$.

2.2 Estimation Error in the Frobenius Norm

Let $X_{1},...,X_{N}$ be $N$ i.i.d. samples from $\operatorname{\mathcal{N}}(0,I)$. In this section we demonstrate a tight bound on how many samples are necessary such that the sample covariance is close to $I$ in Frobenius norm. Let $\widehat{\Sigma}$ denote the empirical covariance, defined to be

By self-duality of the Frobenius norm, we know that

Since there is a $1/4$-net over all PSD matrices with Frobenius norm $1$ of size $9^{d^{2}}$ (see e.g. Lemma 1.18 in [RH17]), the Vershynin-type union bound argument combined with Lemma 2.20 immediately gives us the following:

There exist universal constants $A,B>0$ so that for all $t>0$, we have

By the argument as used in the proof of Lemma 4.3, we obtain:

Fix $\varepsilon,\tau>0$. There is a $\delta_{1}=O(\varepsilon\log 1/\varepsilon)$ such that if $X_{1},\ldots,X_{N}$ are independent samples from $\operatorname{\mathcal{N}}(0,I)$, with

Since the proof is essentially identical to the proof of Lemma 4.3, we omit the proof. However, we note that in fact, the proof technique there can be used to show something slightly stronger, which we will require later. The technique actually shows that if we take any set of size at most $\varepsilon N$, and take the uniform weights over that set, then the empirical covariance is not too far away from the truth. More formally:

Fix $\varepsilon,\tau>0$. There is a $\delta_{2}=O(\varepsilon\log 1/\varepsilon)$ such that if $X_{1},\ldots,X_{N}$ are independent samples from $\operatorname{\mathcal{N}}(0,I)$, with

2.3 Understanding the Fourth Moment Tensor

Our algorithms will be based on understanding the behavior of the fourth moment tensor of a Gaussian when restricted to various subspaces. Let $\otimes$ denote the Kronecker product on matrices. We will make crucial use of the following definition:

We will also require the following definitions:

and let $\Pi_{S}$ and $\Pi_{\operatorname{{\mathcal{S}}}^{\perp}}$ denote the projection operators onto $\operatorname{{\mathcal{S}}}$ and $\operatorname{{\mathcal{S}}}^{\perp}$ respectively. Finally let

In fact, the projection of $v=M^{\flat}$ onto $\operatorname{{\mathcal{S}}}$ where $M$ is symmetric can be written out explicitly. Namely, it is given by

The key result in this section is the following:

It is important to note that the two terms above are not the same; the first term is high rank, but the second term is rank one. The proof of this theorem will require Isserlis’ theorem, and is deferred to Appendix A.

2.4 Concentration of the Fourth Moment Tensor

We also need to show that the fourth moment tensor concentrates:

Fix $\varepsilon,\tau>0$. Let $Y_{i}\sim\operatorname{\mathcal{N}}(0,I)$ be independent, for $i=1,\ldots,N$, where we set

To do so will require somewhat more sophisticated techniques than the ones used so far to bound spectral deviations. At a high level, this is because fourth moments of Gaussians have a sufficiently larger variance that the union bound techniques used so far are insufficient. However, we will show that the tails of degree four polynomials of Gaussians still sufficiently concentrate such that removing points cannot change the mean by too much. The proof requires slightly fancy machinery and appears in Appendix B.

3 Finding the Mean, Using a Separation Oracle

In this section, we consider the problem of approximating $\mu$ given $N$ samples from $\operatorname{\mathcal{N}}(\mu,I)$ in the full adversary model. Our algorithm will be based on working with the following convex set:

It is not hard to show that $\mathcal{C}_{\delta}$ is non-empty for reasonable values of $\delta$ (and we will show this later). Moreover we will show that for any set of weights $w$ in $\mathcal{C}_{\delta}$, the empirical average

will be a good estimate for $\mu$. The challenge is that since $\mu$ itself is unknown, there is not an obvious way to design a separation oracle for $\mathcal{C}_{\delta}$ even though it is convex. Our algorithm will run in two basic steps. First, it will run a very naive outlier detection to remove any points which are more than $O(\sqrt{d})$ away from the good points. These points are sufficiently far away that a very basic test can detect them. Then, with the remaining points, it will use the approximate separation oracle given below to approximately optimize with respect to $C_{\delta}$. It will then take the outputted set of weights and output the empirical mean with these weights. We will explain these steps in detail below.

Our results will hold under the following deterministic conditions:

The concentration bounds we gave earlier were exactly bounds on the failure probability of either of these conditions, albeit for $S_{N,\varepsilon}$ instead of $S_{N,4\varepsilon}$.

The first step of our algorithm will be to remove points which have distance which is much larger than $O(\sqrt{d})$ from the mean. Our algorithm is very naive: it computes all pairwise distances between points, and throws away all points which have distance more than $O(\sqrt{d})$ from more than a $2\varepsilon$-fraction of the remaining points.

Suppose that (7) holds. Then NaivePrune removes no uncorrupted points, and moreover, if $X_{i}$ is not removed by NaivePrune, we have $\|X_{i}-\mu\|_{2}\leq O\left(\sqrt{d\log(N/\tau)}\right)$.

That no uncorrupted point is removed follows directly from (7) and the fact that there can be at most $2\varepsilon N$ corrupted points. Similarly, if $X_{i}$ is not removed by NaivePrune, that means there must be an uncorrupted $X_{j}$ such that $\|X_{i}-X_{j}\|_{2}\leq O(\sqrt{d\log(N/\tau)})$. Then the desired property follows from (7) and a triangle inequality. ∎

Henceforth, for simplicity we shall assume that no point was removed by NaivePrune, and that for all $i=1,\ldots,N$, we have $\|X_{i}-\mu\|_{2}<O(\sqrt{d\log(N/\tau)})$. Otherwise, we can simply work with the pruned set, and it is evident that nothing changes.

3.2 The Separation Oracle

Our main result in this section is an approximate separation oracle for $\mathcal{C}_{\delta}$. Throughout this section, let $w\in S_{N,\varepsilon}$ and set $\widehat{\mu}=\sum_{i=1}^{N}w_{i}X_{i}$. Moreover, let $\Delta=\mu-\widehat{\mu}$. Our first step is to show that any set of weights that does not yield a good estimate for $\mu$ cannot be in the set $\mathcal{C}_{\delta}$:

Suppose that (8)-(9) holds. Suppose that $\|\Delta\|_{2}=\Omega(\sqrt{\varepsilon\delta_{1}})=\Omega(\varepsilon\log 1/\varepsilon)$. Then

By Fact 4.2 and (9) we have $\|\sum_{i\in G}\frac{w_{i}}{w_{g}}X_{i}-\mu\|_{2}\leq\delta_{2}$. Now by the triangle inequality we have

Using the fact that the variance is nonnegative we have

where in the last inequality we have used Fact 4.2 and (8). Hence altogether this implies that

As a corollary, we find that any set of weights in $\mathcal{C}_{\delta}$ immediately yields a good estimate for $\mu$:

Suppose that (8) and (9) hold. Let $w\in\mathcal{C}_{\delta}$ for $\delta=O(\varepsilon\log 1/\varepsilon)$. Then

Our main result in this section is an approximate separation oracle for $\mathcal{C}_{\delta}$ with $\delta=O(\varepsilon\log 1/\varepsilon)$.

Fix $\varepsilon>0$, and let $\delta=O(\varepsilon\log 1/\varepsilon).$ Suppose that (8) and (9) hold. Let $w^{\ast}$ denote the weights which are uniform on the uncorrupted points. Then there is a constant $c$ and an algorithm such that:

(Completeness) If $w=w^{\ast}$, then it outputs “YES”.

We remark that these two facts imply that for any $\tau>0$, the ellipsoid method with this separation oracle will output a $w^{\prime}$ such that $\|w-w^{\prime}\|_{\infty}<\varepsilon/(N\sqrt{d\log(N/\tau)})$, for some $w\in\mathcal{C}_{c\delta}$ in $\operatorname{poly}(d,1/\varepsilon,\log 1/\tau)$ steps.

The conditions that the separation oracle given here satisfy are slightly weaker than the traditional guarantees, given, for instance, in [GLS88]. However, the correctness of the ellipsoid algorithm with this separation oracle follows because outside $C_{c\delta}$, the separation oracle acts exactly as a separation oracle for $w^{*}$. Thus, as long as the algorithm continues to query points outside of $C_{c\delta}$, the action of the algorithm is equivalent to one with a separation oracle for $w^{*}$. Moreover, the behavior of the algorithm is such that it will never exclude $w^{*}$, even if queries are made within $C_{c\delta}$. From these two conditions, it is clear from the classical theory presented in [GLS88] that the ellipsoid method satisfies the guarantees given above.

The separation oracle is given in Algorithm 2. Next, we prove correctness for our approximate separation oracle:

Again, let $\Delta=\mu-\widehat{\mu}$, and let $M=\sum_{i=1}^{N}w_{i}Y_{i}Y_{i}^{T}-I$. By expanding out the formula for $M$, we get:

Suppose $w=w^{\ast}$. Then $\|M\|_{2}<\frac{c}{2}\delta$.

Recall that $w^{\ast}$ are the weights that are uniform on the uncorrupted points. Because $|E|\leq 2\varepsilon N$ we have that $w^{\ast}\in S_{N,\varepsilon}$. We can now use (8) to conclude that $w^{\ast}\in\mathcal{C}_{\delta_{1}}$. Now by Corollary 4.16 we have that $\|\Delta\|_{2}\leq O(\varepsilon\sqrt{\log 1/\varepsilon})$. Thus

Suppose that $w\not\in C_{c\delta}$. Then $|\lambda|>\frac{c}{2}\delta$.

Let us now split into two cases. If $\|\Delta\|_{2}\leq\sqrt{c\delta/10}$, then the first term above is at least $c\delta$ by definition and we can conclude that $|\lambda|>c\delta/2$. On the other hand, if $\|\Delta\|_{2}\geq\sqrt{c\delta/10}$, by Lemma 4.15, we have that

which for sufficiently small $\varepsilon$ also yields $|\lambda|>c\delta/2$. ∎

Suppose $\|\Delta\|_{2}\leq\sqrt{c\delta/10}$. By (11) we immediately have:

since $\lambda>c\delta/2$. On the other hand, if $\|\Delta\|_{2}\geq\sqrt{c\delta/10}$ then by (10) we have $\lambda=\Omega\left(\frac{\|\Delta\|_{2}^{2}}{\varepsilon}\right)$. Putting it all together we have:

where in the last line we used the fact that $\lambda>\Omega\left(\frac{\|\Delta\|_{2}^{2}}{\varepsilon}\right)$, and $\|\Delta\|_{2}^{2}\geq\Omega(\varepsilon^{2}\log 1/\varepsilon)$. This now completes the proof. ∎

3.3 The Full Algorithm

This separation oracle, along with the classical theory of convex optimization [GLS88], implies that we have shown the following:

Fix $\varepsilon,\tau>0$, and let $\delta=O(\varepsilon\sqrt{\log 1/\varepsilon})$. Let $X_{1},\ldots,X_{N}$ be an $\varepsilon$-corrupted set of points satisfying (8)-(9), for $\delta_{1}\leq\delta$ and $\delta_{2}\leq\delta\sqrt{\log 1/\varepsilon}$. Let $c$ be a sufficiently large constant. Then, there is an algorithm $\textsc{LearnApproxMean}(\varepsilon,\tau,X_{1},\ldots,X_{N})$ which runs in time $\operatorname{poly}(N,d,1/\varepsilon,\log 1/\tau)$, and outputs a set of weights $w^{\prime}\in S_{N,\varepsilon}$ such that there is a $w\in C_{c\delta}$ such that $\|w-w^{\prime}\|_{\infty}\leq\varepsilon/(N\sqrt{d\log(N/\tau)})$.

This algorithm, while an extremely powerful primitive, is technically not sufficient. However, given this, the full algorithm is not too difficult to state: simply run NaivePrune, then optimize over $C_{c\delta}$ using this separation oracle, and get some $w$ which is approximately in $C_{c\delta}$. Then, output $\sum_{i=1}^{N}w_{i}X_{i}$. For completeness, the pseudocode for the algorithm is given below. In the pseudocode, we assume that $\textsc{Ellipsoid}(\textsc{SeparationOracleUnknownMean},\varepsilon^{\prime})$ is a convex optimization routine, which given the SeparationOracleUnknownMean separation oracle and a target error $\varepsilon^{\prime}$, outputs a $w^{\prime}$ such that $\|w-w^{\prime}\|_{\infty}\leq\varepsilon^{\prime}$. From the classical theory of optimization, we know such a routine exists and runs in polynomial time.

Fix $\varepsilon,\tau>0$, and let $\delta=O(\varepsilon\sqrt{\log 1/\varepsilon})$. Let $X_{1},\ldots,X_{N}$ be an $\varepsilon$-corrupted set of samples, where

Let $\widehat{\mu}$ be the output of $\textsc{LearnMean}(\varepsilon,\tau,X_{1},\ldots,X_{N})$. Then with probability $1-\tau$, we have $\|\widehat{\mu}-\mu\|_{2}\leq\delta$.

By Fact 4.6, Lemma 4.3, and Lemma 4.4, we know that (7)-(9) hold with probability $1-\tau$, with $\delta_{1},\delta_{2}\leq\delta$. Condition on the event that this event holds. After NaivePrune, by Fact 4.14 we may assume that no uncorrupted points are removed, and all points satisfy $\|X_{i}-\mu\|_{2}\leq O(\sqrt{d\log(N/\tau)})$. Let $w^{\prime}$ be the output of the algorithm, and let $w\in C_{c\delta}$ be such that $\|w-w^{\prime}\|_{\infty}<\varepsilon/(N\sqrt{d\log(N/\tau)})$. By Corollary 4.16, we know that $\|\sum_{i=1}^{N}w_{i}X_{i}-\mu\|_{2}\leq O(\delta)$. Hence, we have

so the entire error is at most $O(\delta)$, as claimed. ∎

4 Finding the Covariance, Using a Separation Oracle

In this section, we consider the problem of approximating $\Sigma$ given $N$ samples from $\operatorname{\mathcal{N}}(0,\Sigma)$ in the full adversary model. Let $U_{i}=\Sigma^{-1/2}X_{i}$ such that if $X_{i}\sim\operatorname{\mathcal{N}}(0,\Sigma)$ then $U_{i}\sim\operatorname{\mathcal{N}}(0,I)$. Moreover let $Z_{i}=U_{i}^{\otimes 2}$. Our approach will parallel the one given earlier in Section 4.3. Again, we will work with a convex set

and our goal is to design an approximate separation oracle. Our results in this section will rely on the following deterministic conditions:

for all $w\in S_{N,\varepsilon}$, and all sets $T\subseteq G$ of size $|T|\leq 2\varepsilon N$. As before, by Fact 4.2, the renormalized weights over the uncorrupted points are in $S_{N,4\varepsilon}$. Hence, we can appeal to Fact 4.6, Corollary 4.8, Corollary 4.9, and Theorem 4.13 with $S_{N,4\varepsilon}$ instead of $S_{N,\varepsilon}$ to bound the probability that this event does not hold. Let $w^{\ast}$ be the set of weights which are uniform over the uncorrupted points; by (13) for $\delta\geq\Omega(\varepsilon\sqrt{\log 1/\varepsilon})$ we have that $w^{\ast}\in C_{\delta}$.

Let $\delta=O(\varepsilon\log 1/\varepsilon)$. Suppose that (13), (14), and 15 hold for $\delta_{1},\delta_{2}\leq O(\delta)$ and $\delta_{3}\leq O(\delta\log 1/\varepsilon)$. Then, there is a constant $c$ and an algorithm such that, given any input $w\in S_{N,\varepsilon}$ we have:

(Completeness) If $w=w^{\ast}$, the algorithm outputs “YES”.

As in the case of learning an unknown mean, by the classical theory of convex optimization this implies that we will find a point $w$ such that $\|w-w^{\prime}\|_{\infty}\leq\frac{\varepsilon}{\operatorname{poly}(N)}$ for some $w^{\prime}\in C_{c\delta}$, using polynomially many calls to this oracle. We make this more precise in the following subsubsection.

The pseudocode for the (approximate) separation oracle is given in Algorithm 4. Observe briefly that this algorithm does indeed run in polynomial time. Lines 2-7 require only taking top eigenvalues and eigenvectors, and so can be done in polynomial time. For any $\xi\in\{-1,+1\}$, line 8 can be run by sorting the samples by $w_{i}\left(\frac{\|Y_{i}\|^{2}}{\sqrt{d}}-\sqrt{d}\right)$ and seeing if there is a subset of the top $2\varepsilon N$ samples satisfying the desired condition, and line 9 can be executed similarly.

We now turn our attention to proving the correctness of this separation oracle. We require the following technical lemmata.

Fix $\delta<1$ and suppose that $M$ is symmetric. If $\|M-I\|_{F}\geq\delta$ then $\|M^{-1}-I\|_{F}\geq\frac{\delta}{2}$.

We will prove this lemma in the contrapositive, by showing that if $\|M^{-1}-I\|_{F}<\frac{\delta}{2}$ then $\|M-I\|_{F}<\delta$. Since the Frobenius norm is rotationally invariant, we may assume that $M^{-1}=\mbox{diag}(1+\nu_{1},\ldots,1+\nu_{d})$, where by assumption $\sum\nu_{i}^{2}<\delta^{2}/4$. By our assumption that $\delta<1$, we have $|\nu_{i}|\leq 1/2$ for all $i$. Thus

where we have used the inequality $|1-\frac{1}{1+x}|\leq|2x|$ which holds for all $|x|\leq 1/2$. This completes the proof. ∎

Let $N_{1},\ldots,N_{d}$ be the columns of $N$. Then

so the desired result follows by taking square roots of both sides. ∎

By the definition of $\|\cdot\|_{\operatorname{{\mathcal{S}}}}$, we have

By self duality of the Frobenius norm, we know that

since $I^{\flat}\in\operatorname{{\mathcal{S}}}^{\perp}$. The result now follows. ∎

Let us first prove completeness. Observe that by Theorem 4.12, we know that restricted to $\operatorname{{\mathcal{S}}}$, we have that $M_{4}=2I$. Therefore, by (15) we will not output a hyperplane in line 7. Moreover, by (14), we will not output a hyperplane in line 8. This proves completeness.

Thus it suffices to show soundness. Suppose that $w\not\in\mathcal{C}_{c\delta}$. We will make use of the following elementary fact:

Let $A=\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}$ and $B=\widehat{\Sigma}^{-1/2}\Sigma\widehat{\Sigma}^{-1/2}$. Then

In particular $A^{-1}=\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}$. Using this expression and the fact that all the matrices involved are symmetric, we can write

where in the third line we have used the fact that the trace of a product of matrices is preserved under cyclic shifts. ∎

Assume (13) holds with $\delta_{1}\leq O(\delta)$ and assume furthermore that $\|A-I\|_{F}\geq c\delta$. Then, if we let $\delta^{\prime}=\frac{(1-\varepsilon)c}{2}\delta=\Theta(\delta)$, we have

Let $A,B$ be as in Fact 4.28. Combining Lemma 4.25 and Fact 4.28 we have

We can rewrite (13) as the expression $\sum_{i\in G}w_{i}X_{i}X_{i}^{T}=w_{g}\Sigma^{1/2}(I+R)\Sigma^{1/2}$ where $R$ is symmetric and satisfies $\|R\|_{F}\leq\delta_{1}$. By the definition of $\widehat{\Sigma}$ we have that $\sum_{i=1}^{N}w_{i}Y_{i}Y_{i}^{T}=I$, and so

by applying Lemma 4.26. And putting it all together we have

It is easily verified that for $c>10$, we have that for all $\delta$, if $\|\widehat{\Sigma}^{-1/2}\Sigma\widehat{\Sigma}^{-1/2}-I\|_{F}\geq c\delta$, then

where $\delta^{\prime}=\frac{c(1-\varepsilon)}{2}\delta=\Theta(\delta)$. The desired result then follows from the Pythagorean theorem. ∎

Claim 4.29 tells us that if $w\not\in C_{c\delta},$ we know that one of the terms in (17) must be at least $\frac{1}{2}\delta^{\prime}$. We first show that if the first term is large, then the algorithm outputs a separating hyperplane:

Assume that (13)-(15) hold with $\delta_{1},\delta_{2}\leq O(\delta)$ and $\delta_{3}\leq O(\delta\log 1/\varepsilon)$. Moreover, suppose that

Then the algorithm outputs a hyperplane in line 7, and moreover, it is a separating hyperplane.

Let us first show that given these conditions, then the algorithm indeed outputs a hyperplane in line 7. Since $I^{\flat}\in S^{\perp}$, the first term is just equal to $\left\|\sum_{i\in E}w_{i}Z_{i}\right\|_{S}$. But this implies that there is some $M^{\flat}\in S$ such that $\|M^{\flat}\|_{2}=\|M\|_{F}=1$ and such that

The $w_{i}/w_{b}$ are a set of weights satisfying the conditions of Claim 4.24 and so this implies that

Let $\widetilde{\Sigma}=\widehat{\Sigma}^{-1}\Sigma$. By Theorem 4.12 and (15), we have that

since it is easily verified that $\delta\|\widetilde{\Sigma}\|^{2}\leq O(\|\widetilde{\Sigma}-I\|_{F})$ as long as $\|\widetilde{\Sigma}-I\|_{F}\geq\Omega(\delta)$, which it is by (17).

Equations 18 and 19 then together imply that

Now let us assume that the first term on the LHS is less than $\frac{1}{2}\delta^{\prime}$, such that the algorithm does not necessarily output a hyperplane in line 7. Thus, the second term on the LHS of Equation 16 is at least $\frac{1}{2}\delta^{\prime}$. We now show that this implies that this implies that the algorithm will output a separating hyperplane in line 9.

Assume that (13)-(15) hold. Moreover, suppose that

Then the algorithm outputs a hyperplane in line 9, and moreover, it is a separating hyperplane.

By the definition of $\operatorname{{\mathcal{S}}}^{\perp}$, the assumption implies that

which is equivalent to the condition that

for some $\xi\in\{-1,1\}$. In particular, the algorithm will output a hyperplane

in Step 9, where $S$ is some set of size at most $\varepsilon N$, and $\lambda=O(\delta^{\prime})$. Since it will not affect anything, for without loss of generality let us assume that $\xi=1$. The other case is symmetrical.

where $\|A\|_{F}=O\left(\delta\frac{N}{|T|}\right)$. Hence,

as long as $\delta^{\prime}\geq O(\delta)$. By self-duality of the Frobenius norm, using the test matrix $\frac{1}{\sqrt{d}}I$, this implies that

These two claims in conjunction directly imply the correctness of the theorem. ∎

As before, this separation oracle and the classical theory of convex optimization [GLS88] shows that we have demonstrated an algorithm FindApproxCovariance with the following properties:

Fix $\varepsilon,\tau>0$, and let $\delta=O(\varepsilon\log 1/\varepsilon)$. Let $c>0$ be a universal constant which is sufficiently large. Let $X_{1},\ldots,X_{N}$ be an $\varepsilon$-corrupted set of points satisfying (13-(15), for $\delta_{1},\delta_{2}\leq O(\delta)$ and $\delta_{3}\leq O(\delta\log 1/\varepsilon)$. Then $\textsc{FindApproxCovariance}(\varepsilon,\tau,X_{1},\ldots,X_{N})$ runs in time $\operatorname{poly}(N,d,1/\varepsilon,\log 1/\tau)$, and outputs a $u$ such that there is some $w\in C_{c\delta}$ such that $\|w-u\|_{\infty}\leq\varepsilon/(Nd\log(N/\tau))$.

As before, this is not quite sufficient to actually recover the covariance robustly. Naively, we would just like to output $\sum_{i=1}^{N}u_{i}X_{i}X_{i}^{T}$. However, this can run into issues if there are points $X_{i}$ such that $\|\Sigma^{-1/2}X_{i}\|_{2}$ is extremely large. We show here that we can postprocess the $u$ such that we can weed out these points. First, observe that we have the following lemma:

Assume $X_{1},\ldots,X_{N}$ satisfy (13). Let $w\in S_{N,\varepsilon}$. Then

This follows since by (13), we have that $\sum_{i\in G}w_{i}X_{i}X_{i}^{T}\succeq w_{g}(1-\delta_{1})\Sigma\succeq(1-O(\delta_{1}))\Sigma$. The lemma then follows since $\sum_{i\in E}w_{i}X_{i}X_{i}^{T}\succeq 0$ always. ∎

Let $u$ be such that there is $w\in C_{c\delta}$ such that $\|u-w\|_{\infty}\leq\varepsilon/(Nd\log(N/\tau))$. Then, $\sum_{i=1}^{N}u^{-}_{i}X_{i}X_{i}^{T}\preceq(1+O(\delta))\Sigma$.

By the definition of $C_{c\delta}$, we must have that $\sum_{i=1}^{N}w_{i}X_{i}X_{i}^{T}\preceq(1+c\delta)\Sigma$. Moreover, we must have $\widetilde{u}^{-}_{i}\leq w_{i}$ for every index $i\in[N]$. Thus we have that $\sum_{i=1}^{N}\widetilde{u}^{-}_{i}w_{i}X_{i}X_{i}^{T}\preceq(1+c\delta)\Sigma$, and hence $\sum_{i=1}^{N}u^{-}_{i}w_{i}X_{i}X_{i}^{T}\preceq(1+c\delta)\Sigma$, since $\sum_{i=1}^{N}u^{-}_{i}w_{i}X_{i}X_{i}^{T}\preceq(1+O(\varepsilon))\sum_{i=1}^{N}\widetilde{u}^{-}_{i}w_{i}X_{i}X_{i}^{T}$. ∎

We now give the full algorithm. The algorithm proceeds as follows: first run FindApproxCovariance to get some set of weights $u$ which is close to some element of $C_{c\delta}$. We then compute the empirical covariance $\Sigma_{1}=\sum_{i=1}^{N}u_{i}X_{i}X_{i}^{T}$ with the weights $u$, and remove any points which have $\|\Sigma_{1}^{-1/2}X_{i}\|_{2}^{2}$ which are too large. We shall show that this removes no good points, and removes all corrupted points which have $\|\Sigma^{-1/2}X_{i}\|_{2}^{2}$ which are absurdly large. We then rerun FindApproxCovariance with this pruned set of points, and output the empirical covariance with the output of this second run. Formally, we give the pseudocode for the algorithm in Algorithm 5.

We now show that this algorithm is correct.

Let $1/2\geq\varepsilon>0$, and $\tau>0$. Let $\delta=O(\varepsilon\log 1/\varepsilon)$. Let $X_{1},\ldots,X_{N}$ be a $\varepsilon$-corrupted set of samples from $\operatorname{\mathcal{N}}(0,\Sigma)$ where

Let $\widehat{\Sigma}$ be the output of $\textsc{LearnCovariance}(\varepsilon,\tau,X_{1},\ldots,X_{N})$. Then with probability $1-\tau$, $\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I\|_{F}\leq O(\delta)$.