Robustly Learning a Gaussian: Getting Optimal Error, Efficiently

Introduction

The most popular and widely used modeling assumption is that data is approximately Gaussian. This is a convenient simplification to make when modeling velocities of particles in an ideal gas [Goo15], measuring physical characteristics across a population (after controlling for gender), and even modeling fluctuations in a stock price on a logarithmic scale. However, real data is not actually Gaussian and is at best crudely approximated by a Gaussian (e.g., with heavier tails). What’s worse is that estimators designed under this assumption can perform poorly in practice and be heavily biased by just a few errant samples that do not fit the model.

For over fifty years, the field of robust statistics [HR09, HRRS86, RL05] has studied exactly this phenomenon — the sensitivity or insensitivity of estimators to small deviations in the model. Unsurprisingly, one of the central questions that shaped its development was the problem of learning the parameters of a one-dimensional Gaussian distribution when a small fraction of the samples are arbitrarily corrupted. More precisely, in 1964, Huber [Hub64] introduced the following model:

In Huber’s contamination model, we are given samples from a distribution

where $\mathcal{N}(\mu,\sigma^{2})$ is a Gaussian of mean $\mu$ and variance $\sigma^{2}$, and $\mathcal{Z}$ is an arbitrary distribution chosen by an adversary.

Intuitively, among our samples, about a $(1-\varepsilon)$ fraction will have been generated from a Gaussian and are called inliers, and the rest are called outliers or gross corruptions. We will work with an even more challengingNone of the results in our paper were previously known in Huber’s contamination model either. The reason we work with this stronger model is because we can — nothing in our analysis relies on the inliers and outliers being independent. model — called the strong contamination model (Definition 2) — where the adversary is allowed to look at the inliers and then decide on the outliers. The literature on robust statistics has given numerous explanations and empirical investigations [GCSR14, Ham01] into how such outliers might arise as the result of equipment failure, data being entered incorrectly, or even from a subpopulation that was not accounted for in a medical study. These types of errors are erratic and difficult to model, so instead our goal is to design a procedure that accurately estimates $\mu$ and $\sigma^{2}$ without making any assumptions about them.

In one dimension, the median and median absolute deviation are well-known robust estimators for the mean and variance respectively. In particular, given samples $X_{1},X_{2},\ldots,X_{n}$, we can compute

where $\Phi$ is the cumulative distribution of the standard Gaussian. (This scaling constant is needed to ensure that $\widehat{\sigma}$ is an unbiased estimator when there is no noise.) If $n\geq C\frac{\log 1/\delta}{\varepsilon^{2}}$, then with probability at least $1-\delta$ we have that $d_{TV}(\mathcal{N}(\mu,\sigma^{2}),\mathcal{N}(\widehat{\mu},\widehat{\sigma}^{2}))\leq C\varepsilon$. In Huber’s contamination model, this is the strongest type of error guarantee we could hope forSee Lemma 17. and captures both the task of learning the underlying parameters $\mu$ and $\sigma^{2}$, and finding the approximately best fit to the observed distribution within the family of one-dimensional Gaussians. In fact there are plentifully many other estimators — such as the trimmed mean, winsorized mean, Tukey’s biweight function, and the interquartile range — that achieve the same sorts of error guarantees, up to constant factors. The design of robust estimators for location (e.g., estimating $\mu$) and scale (e.g., estimating $\sigma^{2}$) is guided by certain overarching principles, such as the notion of the influence curve [HRRS86] or the notion of breakdown point [RL05]. In some cases, it is even possible to design robust estimators that are minimax optimal [Hub64].

These days, much of modern data analysis revolves around high-dimensional data — for example, when we model documents [BNJ03], images [OF96], and genomes [NJB+08] as vectors in a very high-dimensional space. The need for robust estimators is even more pressing in these applications, since it is infeasible to remove obvious outliers by inspection. However, adapting robust statistics to high-dimensional settings is fraught with challenges. The principles that guided the design of robust estimators in one dimension seem to inherently lead to high-dimensional estimators that are hard to compute [Ber06, HM13].

In this paper, we focus on the central problem of learning the parameters of a multivariate Gaussian $\mathcal{N}(\mu,\Sigma)$ in the strong contamination model. The textbook estimators for the mean and covariance – such as the Tukey median [Tuk75] and minimum volume enclosing ellipsoid [Rou85] – essentially search for directions where the projection of $\mathcal{D}$ is suitably non-Gaussian. However, trying to find a direction where the projection is non-Gaussian can be like looking for a needle in an exponentially-large haystack – these statistics are not efficiently computable, in general. Furthermore, a random projection will look Gaussian with high probability [Kla07].

In this paper, our main result is an efficiently computable estimator for a high-dimensional Gaussian that achieves error

in the strong contamination model, for a universal constant $C$ that is independent of the dimension. For a Gaussian distribution, we consider estimation in terms of total variation distance, which is equivalent to estimating the parameters under the natural measures. Our main idea is to use various regularity conditions satisfied by the inliers to make the problem of searching for non-Gaussian projections easier. When just the mean $\mu$ is unknown, our algorithm runs in time polynomial in the dimension $d$ and $1/\varepsilon$. When both the mean and covariance are unknown, our algorithm runs in time polynomial in $d$ and quasi-polynomial in $1/\varepsilon$. All of our algorithms achieve polynomial sample complexity.

Prior to our work, the best known algorithm of Diakonikolas et al. [DKK+16] achieved estimation error $O(\varepsilon\log 1/\varepsilon)$ for this problem We note that, as stated, the results in [DKK+16] give estimation error $O(\varepsilon\log^{3/2}1/\varepsilon)$. However, combining the techniques in [DKK+16] with the arguments in Section 7 of this paper gives the stated bound. This argument will be included in the full version of [DKK+16]., again with respect to total variation distance. Concurrently, Lai, Rao and Vempala [LRV16] gave an algorithm which achieves estimation error roughly $O(\varepsilon^{1/2}\log^{1/2}d)$. In fact, the algorithm of Diakonikolas et al. [DKK+16] works in a stronger model than what we consider here, where an adversary gets to look at the samples and then decides on an $\varepsilon$-fraction to move arbitrarily. Such errors are both additive and subtractive (because inliers are removed). Interestingly, Diakonikolas, Kane and Stewart [DKS17] proved that any Statistical Query learning algorithm that works in such an additive and subtractive model and achieves an error guarantee asymptotically better than $O(\varepsilon\log^{1/2}1/\varepsilon)$ must make a super-polynomial number of statistical queries. Our work shows a natural conclusion that in an additive only model it is possible to algorithmically achieve the same error guarantees as are possible in the one-dimensional case, up to a universal constant.

2 Our Results and Techniques

In what follows, we will explain both our work as well as prior work through the following lens:

At the core of any robust estimator is some procedure to certify that the estimates have not been moved too far away from the true parameters by a small number of corruptions.

First, we consider the subproblem where the covariance $\Sigma=I$ is known and only the mean $\mu$ is unknown. In the terminology of robust statistics, this is called robust estimation of location. If we could compute the Tukey median, we would have an estimate that satisfies $d_{TV}(\mathcal{N}(\mu,I),\mathcal{N}(\widehat{\mu},I))\leq C\varepsilon$. The way that the Tukey median guarantees that it is close to the true mean is that along every direction $u$ it is close to the median of the projection of the samples. More precisely, at least a $\frac{1-\varepsilon}{2}$ fraction of the samples satisfy $u^{T}X_{i}\geq u^{T}\widehat{\mu}$, and at least a $\frac{1-\varepsilon}{2}$ fraction of the samples satisfy $u^{T}\widehat{\mu}\geq u^{T}X_{i}$. However, if we have a candidate $\widehat{\mu}$, finding a direction $u$ that violates this condition is again like searching for a needle in an exponentially large haystack.

The approach of Diakonikolas et al. [DKK+16] was essentially a data-dependent way to search for appropriate directions $u$, by looking for directions where the empirical variance is larger than it should be (if there were no corruptions). However, because their approach considers only a single direction at a time, it naturally gets stuck at error $\Theta(\varepsilon\log^{1/2}1/\varepsilon)$. This is because along the direction $u$, only when a point is $\Omega(\log^{1/2}1/\varepsilon)$ away from most of the rest of the samples can we be relatively confident that it is an outlier. Thus, an adversary could safely place all the corruptions in the tails and move the mean by as much as $\Theta(\varepsilon\log^{1/2}1/\varepsilon)$. This would not affect the Tukey median by as much, but would affect an estimate based on the empirical mean (because the algorithm could find no other outliers to remove) by considerably more.

Our approach is to consider logarithmically many directions at once. Even though an inlier can be logarithmically many standard deviations away from the mean along a single direction $u$ with reasonable probability, it is unlikely to be that many standard deviations away simultaneously across many orthogonal directions. Essentially, this allows us to remove the influence of outliers on all but a logarithmic dimensional subspace. Combining this with an algorithm for robustly learning the mean in time exponential in the dimension (but polynomial in the number of samples), we obtain our first main result:

Suppose we are given a set of $n=\operatorname*{poly}(d,1/\varepsilon)$ samples from the strong contamination model, where the underlying $d$-dimensional Gaussian is $\mathcal{N}(\mu,I)$. Let $\varepsilon\leq\varepsilon_{0}$, where $\varepsilon_{0}$ is a positive universal constant. For any $\beta>0$, there is an algorithm to learn an estimate $\mathcal{N}(\widehat{\mu},I)$ that with high probability satisfies

Moreover, the algorithm runs in time $\operatorname*{poly}(n,(1/\varepsilon)^{\beta})$.

We prove an almost matching lower bound of $\frac{\varepsilon}{2}+\Omega(\varepsilon^{2})$ on the estimation error. Thus, our robustness guarantee is optimal up to a factor of $\sqrt{2}$, even among computationally inefficient robust estimators. Interestingly, our extra factor of $\sqrt{2}$ comes from the following geometric fact which we make crucial use of: Any convex body of diameter $D$ in any dimension can be covered by a ball of radius $D/\sqrt{2}$, and moreover such a ball can be (approximately) found in time exponential in the dimension. Suppose that along some direction $u$ we have an estimate $p$ that is guaranteed to be within $\varepsilon/2$ of the projection of the true mean $\mu$. We can now confine $\mu$ to a slab of width $\varepsilon$, and by taking the intersection of all such slabs we get a convex body that contains $\mu$ and has diameter of at most $\varepsilon$. By covering the body with a ball of radius $\varepsilon/\sqrt{2}$, we are guaranteed that the center of the ball is within $\varepsilon/\sqrt{2}$ of the true mean. This gives us a general way to combine one-dimensional robust estimates along a net of directions.

We note that, for general isotropic sub-Gaussian distributions, the bound of $O(\varepsilon\log^{1/2}1/\varepsilon)$ of [DKK+17] is optimal for robust mean estimation, even in one dimension. See Section A for a proof of this fact. However, our results can be seen to hold more generally than stated above – indeed, the same arguments work for a class of symmetric isotropic sub-Gaussian distributions which are sufficiently smooth near their mean. More precisely, we require that along any univariate projection, the mean is robustly estimated by the median.

We next consider the subproblem where the mean $\mu=0$ is known and only the covariance $\Sigma$ is unknown. In the terminology of robust statistics, this is called robust estimation of scale. In this case, we want to compute an estimate $\widehat{\Sigma}$ that satisfiesMore precisely, to obtain $O(\varepsilon)$ error guarantee with respect to the total variation distance, we need to robustly approximate $\Sigma$ within $O(\varepsilon)$ in Mahalanobis distance, which is a stronger metric than the Frobenius norm. As part of our approach, we are able to efficiently reduce to the case that $\Sigma$ is close to the identity matrix, in which case the Frobenius error suffices. $\|\Sigma-\widehat{\Sigma}\|_{F}\leq C\varepsilon$. When $\widehat{\Sigma}$ does not satisfy this condition, it can be shown (in Section 6.2.3) that there is a degree-two polynomial $p(X)$, where

It turns out that, even given the polynomial $p(X)$, deciding whether or not the above conditions approximately hold is challenging. Given $p(X)$ and $\widehat{\Sigma}$, we can certainly compute $\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\widehat{\Sigma})}[p(X)]$. But given only contaminated samples from $\mathcal{N}(0,\Sigma)$ and without knowing what $\Sigma$ is, can we estimate $\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[p(X)]$?

Often, univariate robust estimation problems are considered easy, with a simple recipe: Construct an unbiased estimator for the statistic for which each sample point has low influence. However, in our setting, it is highly non-trivial to construct such an estimator. The naive attempt in this case would be the median – this immediately fails since the distribution of $p(X)$ is asymmetric. Even if there were no noise, that would not necessarily be an unbiased estimator. So how can we dampen the influence of outliers, if there is no natural symmetry in the distribution? We construct a robust estimator crucially using the fact that $p(X)$ is the weighted sum of chi-squared random variables when there is no noise. The key structural fact we exploit is the following: Given two sums of chi-squared random variables, if the random variables are far in total variation distance, most of their difference must lie close to their means. We use this fact to show how, given a weak estimate of the mean (i.e., one which is only accurate up $\omega(\varepsilon)$), one can improve the estimate by a constant factor. Our result follows by an iterative application of this technique.

However, there is still a major complication in utilizing our low-dimensional estimator to obtain a high-dimensional estimator. In the unknown mean case, we knew the higher-order moments (since we assumed that the covariance is the identity). Here, we do not have control over the higher-order moments of the unknown Gaussian. Overcoming this difficulty requires several new techniques, which are quite complicated, and we defer the full details to Section 6. Our second main result is:

Suppose we are given a set of $n=\operatorname*{poly}(d,1/\varepsilon)$ samples from the strong contamination model, where the underlying $d$-dimensional Gaussian is $\mathcal{N}(0,\Sigma)$. There is an algorithm to learn an estimate $\mathcal{N}(0,\widehat{\Sigma})$ that runs in time $\operatorname*{poly}(n,(1/\varepsilon)^{O(\log^{4}1/\varepsilon)})$ and with high probability satisfies

for a universal constant $C$ that is independent of the dimension.

A key technical problem arises when we attempt to combine estimates for the covariance restricted to a subspace and its orthogonal complement. We refer to this as a stitching problem, where if we write $\Sigma$ as

and have accurate estimates for $\Sigma_{V}$ and $\Sigma_{V^{\perp}}$, we still need to accurately estimate $A$. Our algorithm utilizes an unexpected connection to the unknown mean case: We show that, under a carefully chosen projection scheme, we can simulate noisy samples from a Gaussian with identity covariance, where the mean of this distribution encodes the information needed to recover $A$. We defer the full details to Section 6.4.

It turns out that we can solve the general case when both $\mu$ and $\Sigma$ are unknown, by directly reducing to the previous subproblems, exactly as was done in [DKK+16] (with some caveats, addressed in Section 4.4). Since all of our error guarantees are optimal up to constant factors, there is only a constant factor loss in this reduction. Finally, we obtain the following corollary:

Suppose we are given a set of $n=\operatorname*{poly}(d,1/\varepsilon)$ samples from the strong contamination model, where the underlying $d$-dimensional Gaussian is $\mathcal{N}(\mu,\Sigma)$. There is an algorithm to learn an estimate $\mathcal{N}(\widehat{\mu},\widehat{\Sigma})$ that runs in time $\operatorname*{poly}(n,(1/\varepsilon)^{O(\log^{4}1/\varepsilon)})$ and with high probability satisfies

for a universal constant $C$ that is independent of the dimension.

This essentially settles the complexity of robustly learning a high-dimensional Gaussian. The sample complexity of our algorithm depends polynomially on $d$ and $1/\varepsilon$, and the running time depends polynomially on $d$ and quasi-polynomially on $1/\varepsilon$. Up to a constant factor, ours is the first high-dimensional algorithm that achieves the same error guarantees as in the one-dimensional case, where results were known for more than fifty years! It is an interesting open problem to reduce the running time to polynomial in $1/\varepsilon$ (while still being polynomial in $d$). As we explain in Section 6.6, this seems to require fundamentally new ideas.

In addition to the works mentioned above, there has been an exciting flurry of recent work on robust high-dimensional estimation. This includes studying graphical models in the presence of noise [DKS16], tolerating much more noise by allowing the algorithm to output a list of candidate hypotheses [CSV17], formulating general conditions under which robust estimation is possible [SCV18], developing robust algorithms under sparsity assumptions [Li17, DBS17, BDLS17] where the number of samples is sublinear in the dimension, and leveraging theoretical insights to give practical algorithms that can be applied to genomic data [DKK+17]. We note that, in comparison to all these other works, ours is the only to efficiently achieve the information-theoretically optimal error guarantee (up to constant factors). Despite all of this rapid progress, there are still many interesting theoretical and practical questions left to explore.

3 Organization

In Section 2, we go over preliminaries and notation that we will use throughout the paper. In Section 3, we describe an algorithm for robustly estimating the mean of a Gaussian in low-dimensional settings, and crucially apply it in the design of an algorithm for mean-estimation in high dimensions, described in Section 4. Similarly, in Section 5, we give an algorithm for robustly estimating the mean of degree-two polynomials in certain settings, which is applied in the context of our covariance-estimation algorithm in Section 6. Finally, we put these tools together and describe our general algorithm for robustly estimating a Gaussian in Section 7.

Preliminaries

Here we formally define the strong contamination model.

Fix $\varepsilon>0$. We say a set of samples $X_{1},\ldots,X_{n}$ was generated from the strong contamination model on a distribution $F$, if it was generated via the following process:

We produce $(1-\varepsilon)n$ i.i.d. samples $G$ from $F$.

An adversary is allowed to observe these samples and add $\varepsilon n$ points $E$ arbitrarily.

We are then given the set of samples $G\cup E$ in random order. Also, we will say that the samples $X_{1},\ldots,X_{n}$ are $\varepsilon$-corrupted. Moreover given an $\varepsilon$-corrupted set of samples $S$, we will write $S=(G,E)$ where $G$ is the set of uncorrupted points and $E$ is the set of corrupted points. Moreover, given a subset $S^{\prime}\subset S$, we will also write $S^{\prime}=(G^{\prime},E^{\prime})$, where $G^{\prime}=S^{\prime}\cap G$ and $E^{\prime}=S^{\prime}\cap E$ denote the set of uncorrupted points and corrupted points remaining in $S^{\prime}$. $L$ will denote $G\setminus G^{\prime}$, which is the set of “lost” uncorrupted points.

Given a contaminated set $S^{\prime}=(G^{\prime},E^{\prime})$ and a set $G$ so that $G^{\prime}\subseteq G$, define the following quantities

In particular, observe that if $\Delta(S^{\prime},G)<O(\varepsilon)$, then a simple calculation implies that $\phi(S^{\prime},G)\leq O(\varepsilon/\log 1/\varepsilon)$. Equivalently, we have removed at most an $O(\varepsilon/\log 1/\varepsilon)$ fraction of good points from $G$. This is crucial, as if we throw out an $\varepsilon$-fraction of good points then we essentially put ourselves in the subtractive model, and there our guarantees no longer hold.

There are two differences between the strong contamination model and Huber’s contamination model. First, the number of corrupted points is fixed to be $\varepsilon n$ instead of being a random variable. However, this difference is negligible. It follows from basic Chernoff bounds that $n$ samples from Huber’s contamination model with parameter $\varepsilon$ (for $n$ sufficiently large) can be simulated by a $(1+o(1))\varepsilon$-corrupted set of samples, except with negligible failure probability. Hence, we lose only an additive $o(\varepsilon)$ term when translating from Huber’s contamination model to the strong contamination model, which will not change any of the guarantees in our paper. The second difference is that the adversary is allowed to inspect the uncorrupted points before deciding on the corrupted points. This makes the model genuinely stronger since the samples we are given are no longer completely independent of each other.

2 Deterministic Regularity Conditions

In analyzing our algorithms, we only need certain deterministic regularity conditions to hold on the uncorrupted points. In this subsection, we formally state what these conditions are. It follows from known concentration bounds that these conditions all hold with high probability given a polynomial number of samples. Now with these regularity conditions defined once and for all, we will be able to streamline our proofs in the sense that each step in the analysis will only ever use one of these fixed set of conditions and will not use the randomness in the sampling procedure. We remark that some subroutines in our algorithm only need a subset of these conditions to hold, so we could improve the sample complexity by changing the regularity conditions we need at each step. However, since we will not be concerned with optimizing the sample complexity beyond showing that it is polynomial, we choose not to complicate our proofs in this manner.

In the unknown mean case, we will require the following condition:

For all $x\in G$ we have $\|x-\mu\|_{2}\leq O(\sqrt{d\log(|G|/\delta)})$.

We have that $\|\mathop{\mbox{\bf E}}_{G}[X]-\mathop{\mbox{\bf E}}_{\mathcal{N}(\mu,I)}[X]\|_{2}\leq\eta.$

It is easy to show (see Lemma 6) that given enough samples from $\mathcal{N}(\mu,I)$, the empirical data set will satisfy these conditions with high probability.

2.2 Regularity Conditions for Unknown Covariance

In the unknown covariance case, we will require the following condition:

For all $x\in G$ we have that $x^{T}\Sigma^{-1}x=O(d\log(|G|/\delta)).$

As before, it is easy to show (see Lemma 14) that given enough samples from $\mathcal{N}(0,\Sigma)$, the empirical data set will satisfy these conditions with high probability.

3 Bounds on the Total Variation Distance

For clarity of exposition we defer this calculation to the Appendix.

We also need to bound the total variation distance between two Gaussians with zero mean and different covariance matrices. The natural norm to use is the Mahalanobis distance. But in our setting, we will be able to use the more convenient Frobenius norm instead (because we effectively reduce to the case that the covariance matrices will be close to the identity):

These lemmata show that parameter estimation and approximation in total variation distance are essentially equivalent. Indeed, in this paper, we achieve both guarantees, but state our results in terms of total variation estimation.

Robustly Learning the Mean in Low Dimensions

This section is dedicated to the proof of the following theorem:

In particular, as we let $\rho,\gamma\rightarrow 0$, the parameter estimation error approaches $\sqrt{\pi}\varepsilon$ (corresponding to a total variation approximation of $\varepsilon/\sqrt{2}$). In Lemma 17 in the Appendix we show that no algorithm can achieve parameter estimation error better than $\sqrt{\frac{\pi}{2}}\varepsilon$. Thus, we achieve a $\sqrt{2}$ approximation to the optimal error.

First we show that if we project onto one dimension, then the median of the corrupted data differs from the true mean by at most $\sqrt{\frac{\pi}{2}}\varepsilon+o(\varepsilon)$. Our proof will rely only on the notion of a $(\gamma\varepsilon,\delta)$-good set with respect to $\mathcal{N}(\mu,I)$ and thus it works even in the strong contamination model. Formally, we show:

Observe that we have $\left|\frac{|E|}{|S|}\Pr_{X\sim E}[\langle v,X\rangle>a]\right|\leq\psi(S,G)$. Moreover, by simple calculation we have

For $|a|=O(\varepsilon)$ we have that $\Pr_{\mathcal{N}(0,I)}[X>a]=\frac{1}{2}-\frac{1}{\sqrt{2\pi}}a+O(\varepsilon^{3})$. Thus, by $(\gamma\varepsilon,\delta)$-goodness of $G_{0}$, this implies that for $|a|=O(\varepsilon)$, we have

In particular, we see that if $a>\sqrt{\frac{\pi}{2}}\varepsilon+O\left(\frac{\gamma\varepsilon}{d}\right)+o(\varepsilon)$, then $\Pr_{X\sim S}\left[\langle v,X\rangle>\Pi_{y}\mu+\sqrt{\frac{\pi}{2}}\varepsilon\right]<1/2$. By symmetric logic, we also have that $\Pr_{X\sim S}\left[\langle v,X\rangle>\Pi_{y}\mu-\sqrt{\frac{\pi}{2}}\varepsilon\right]>1/2$. Thus, the median in direction $v$ differs from $\Pi_{y}\mu$ by at most $\sqrt{\frac{\pi}{2}}\varepsilon+O\left(\frac{\gamma\varepsilon}{d}\right)+o(\varepsilon)$. ∎

2 Finding a Minimum Radius Circumscribing Ball

Our first step is to use such an oracle to construct a net for ${\cal C}$. First, we need the following well-known bound on the size of the net.

The algorithm is fairly straightforward. First, we observe that ${\cal C}$ is contained within $B(x,2R)$. We then form a $(\rho R)/3$-net of $B(x,2R)$ using Claim 1. We then iterate over every element $v$ of this net, and use our projection oracle to (approximately) find the closest point in ${\cal C}$ to $v$. If this point is too far away, we throw it out, otherwise, we add this projected point into the net. The formal pseudocode for CircumscribeNet is given in Algorithm 1.

The runtime bound follows from Claim 1. We now turn our attention to correctness. By Claim 1, and rescaling and shifting, the set $\mathcal{F}$ is clearly a $(\rho R)/3$-net for a ball $B$ of radius $2R$ containing ${\cal C}$. We now claim that the set $\mathcal{X}$ is indeed a $(\rho R)/3$-net for ${\cal C}$. Fix $y\in{\cal C}$. Since ${\cal C}\subseteq B$, this implies there is some $v\in\mathcal{F}$ so that $\|y-v\|_{2}\leq\rho R/3$. Thus, in Line 7, when processing $v$, we must find some $u_{v}\in{\cal C}$ so that $\|u_{v}-v\|_{2}\leq 2\rho R/3$. The claim then follows from the triangle inequality. ∎

Fix $R,{\cal C},\rho,{\cal O},x$ as in Lemma 4. Suppose a call to ${\cal O}$ runs in time $T$. Then, there is an algorithm $\textsc{Circumscribe}(R,\rho,{\cal O},x)$ which runs in time $\operatorname*{poly}((R/\rho)^{O(d)},T)$ and returns a point $\widehat{y}$ so that ${\cal C}$ is contained within a ball of radius $\sqrt{2}(1+2\rho)R$.

The algorithm at this point is very simple. Using the output of CircumscribeNet, we iterate over all points in a net over $B(x,2R)$, find an $x$ in this net so that the distance to all points in the net is at most $\sqrt{2}(1+\rho)R$, and output any such point. The formal pseudocode for Circumscribe is given in Algorithm 2.

The runtime bound is immediate. By Theorem 4, there is some $y\in B(x,2R)$ so that ${\cal C}\subseteq B(y,R\sqrt{2})$. Thus, by the triangle inequality, there is some $y^{\prime}\in\mathcal{F}$ so that ${\cal C}\subseteq B(y,\sqrt{2}(1+\rho)R)$. Thus, the algorithm will output some point $y^{\prime\prime}\in\mathcal{F}$. By an additional application of the triangle inequality, since $\mathcal{X}$ is a $\rho R$-net for ${\cal C}$, this implies that ${\cal C}\subseteq B(y^{\prime\prime},\sqrt{2}(1+2\rho)R)$, as claimed. ∎

3 The Full Low-Dimensional Algorithm

where $\beta=\sqrt{\frac{\pi}{2}}\varepsilon+O\left(\frac{\gamma\varepsilon}{d}\right)+o(\varepsilon)$ is as in Lemma 3. We now show two properties of this set, which in conjunction with the machinery above, allows us to prove Theorem 3. The first shows that ${\cal C}$ has small diameter:

For all $x,y\in{\cal C}$, we have $\|x-y\|_{2}\leq 2\beta/(1-\rho)$.

Fix any $x,y\in{\cal C}$. By definition of ${\cal C}$, it follows that for all $v\in\mathcal{F}$, we have $|\langle x-y,v\rangle|\leq 2\beta$. For any $u$ with $\|u\|_{2}=1$, there is some $v\in\mathcal{F}$ with $\|u-v\|_{2}\leq\varepsilon$, and so we have

Taking the supremum over all unit vectors $u$ and simplifying yields that $\|x-y\|_{2}\leq 2\beta/(1-\rho)$, as claimed. ∎

The second property shows that we may find an $\alpha$-projection oracle for ${\cal C}$ efficiently.

Fix $\rho^{\prime}>0$. There is a $\rho^{\prime}$-projection oracle $\textsc{ProjOracle}(y,\rho^{\prime},{\cal C})$ for ${\cal C}$ which runs in time $\operatorname*{poly}((1/\rho)^{O(d)},\log(\gamma\varepsilon/(1-\rho)),\log(1/\rho^{\prime}))$.

We now finally describe LearnMeanLowD. Using convex optimization, we first find an arbitrary $x\in{\cal C}$. By Lemma 3 we know $\mu\in{\cal C}$ and so this step succeeds. After constructing ${\cal C}$, we run Circumscribe with appropriate parameters, and return the outputted point. The formal pseudocode for LearnMeanLowD is given in Algorithm 3.

The runtime claim follows from the runtime claims for Circumscribe and ProjOracle. Thus, it suffices to prove correctness of this algorithm. By Lemma 3, we know that $\mu\in{\cal C}$. By Claim 2 and Corollary 2, the output $y$ satisfies $B(y,\sqrt{2}\frac{1+2\rho}{1-\rho}\beta)$. Thus, we have $\|\mu-y\|_{2}\leq\sqrt{2}\frac{1+2\rho}{1-\rho}\beta$, as claimed. ∎

Robustly Learning the Mean in High Dimensions

In this section, we prove the following theorem, which is our first main result:

Fix $\varepsilon,\gamma,\delta>0$, and let $X_{1},\ldots,X_{n}$ be an $\varepsilon$-corrupted set of points from $\mathcal{N}(\mu,I)$, where $\|\mu\|_{2}\leq O(\varepsilon\log 1/\varepsilon)$, and where

Then, for every $\alpha,\beta>0$, there is an algorithm $\textsc{RecoverMean}(X_{1},\ldots,X_{n},\varepsilon,\delta,\gamma,\alpha,\beta)$ which runs in time $\operatorname*{poly}(d,1/\gamma,1/\varepsilon^{\beta},1/\alpha,\log 1/\delta)$ and outputs a $\widehat{\mu}$ so that with probability $1-\delta$, we have $\|\widehat{\mu}-\mu\|_{2}\leq\left(\frac{\sqrt{\pi}+O(\gamma)}{1-\alpha}+\frac{1}{\sqrt{\beta}}\right)\varepsilon$.

In particular, observe that Theorem 5, in conjunction with Lemma 1, gives us Theorem 1, if we set $\gamma=o(1)$. With this, we may state our primary algorithmic contribution:

Fix $\varepsilon,\gamma,\alpha,\delta,\beta>0$, and let $S_{0}=(G_{0},E_{0})$ be an $\varepsilon$-corrupted set of samples of size $n$ from $\mathcal{N}(\mu,I)$, where $\|\mu\|_{2}\leq O(\varepsilon\log 1/\varepsilon)$, and where $n=\operatorname*{poly}(d,1/(\gamma\varepsilon),\log 1/\delta)$. Suppose that $G_{0}$ is $(\gamma\varepsilon,\delta)$-good with respect to $\mathcal{N}(\mu,I)$. Let $S\subseteq S_{0}$ be a set so that $\Delta(S,G_{0})\leq\varepsilon$. Then, there exists an algorithm FilterMeanOpt that given $S,\varepsilon,\gamma,\alpha,\beta$ outputs one of two possible outcomes:

A $\widehat{\mu}$, so that $\|\widehat{\mu}-\mu\|_{2}\leq\left(\frac{\sqrt{\pi}+O(\gamma)}{1-\alpha}+\frac{1}{\sqrt{\beta}}\right)\varepsilon$.

A set $S^{\prime}\subset S$ so that $\Delta(S^{\prime},G_{0})<\Delta(S,G_{0})$.

Moreover, FilterMeanOpt runs in time $\operatorname*{poly}(d,1/\gamma,1/\varepsilon^{\beta},1/\alpha,\log 1/\delta)$.

By first running the algorithm of [DKK+16] to obtain an estimate of the mean to error $O(\varepsilon\sqrt{\log 1/\varepsilon})$, then running FilterMeanOpt at most polynomially many times, we clearly recover the guarantee in Theorem 5. Thus, the rest of the section is dedicated to the proof of Theorem 6.

At a high level, the structure of the argument is as follows: We first show that if there is a subspace of eigenvectors of dimension at least $O(\log 1/\varepsilon)$ of the empirical covariance matrix with large associated eigenvalues, then we can produce a filter using a degree-2 polynomial (Section 4.1). Otherwise, we know that there are at most $O(\log 1/\varepsilon)$ eigenvectors of the empirical covariance with a large eigenvalue. We can learn the mean in this small dimensional subspace using our learning algorithm from the previous section, and then we can argue that the empirical mean on the remaining subspace is close to the true mean (Section 4.2).

This outline largely follows the structure of the filter arguments given in [DKK+16], however, the filtering algorithm we use here requires a couple of crucial new ideas. First, to produce the filter, instead of using a generic degree-2 polynomial over this subspace, we construct an explicit, structured, degree-2 polynomial which produces such a filter. Crucially, we can exploit the structure of this polynomial to obtain very tight tail bounds, e.g., via the Hanson-Wright inequality. This is critical to avoid a quasi-polynomial runtime. If instead we used arbitrary degree-$2$ polynomials in this subspace, it would need to be of dimension $O(\log^{2}1/\varepsilon)$ and the low-dimensional algorithm in the second step would take quasi-polynomial time.

Second, we must be careful to throw out far fewer good points than corrupted points. In particular, by our definition of $\Delta$ (which gives an additional logarithmic penalty to discarding good points) and our guarantee that $\Delta$ decreases, our filter can only afford to throw out an $\varepsilon/\log(1/\varepsilon)$ fraction of good points in total, since $\Delta$ is initially $\varepsilon$. This is critical, as if we threw away an $\varepsilon$-fraction of good points, then proving that the problem remains efficiently solvable becomes problematic. In particular, if these points were thrown away arbitrarily, then this becomes the full additive and subtractive model, for which a statistical query lower bound prevents us from getting an $O(\varepsilon)$-approximate answer in polynomial time [DKS17]. To avoid discarding too many good points, we exploit tight exponential tail bounds of Gaussians, and observe that by slightly increasing the threshold at which we filter away points, we decrease the fraction of good points thrown away dramatically.

We now give an algorithm for the case when there are many eigenvalues which are somewhat large. Formally, we show:

Fix $\varepsilon,\gamma,\delta,\alpha,\beta>0$, and let $S_{0}=(G_{0},E_{0})$ be an $\varepsilon$-corrupted set of samples of size $n$ from $\mathcal{N}(\mu,I)$, where $\|\mu\|_{2}\leq O(\varepsilon\log 1/\varepsilon)$, and where $n=\operatorname*{poly}(d,1/(\gamma\varepsilon),\log 1/\delta)$. Suppose that $G_{0}$ is $(\gamma\varepsilon,\delta)$-good with respect to $\mathcal{N}(\mu,I)$. Let $S\subseteq S_{0}$ be a set so that $\Delta(S,G_{0})\leq\varepsilon$. Let $\widehat{\Sigma}$ be the sample covariance of $S$, let $\widehat{\mu}$ be the sample mean of $S$, and let $V$ be the subspace of all eigenvectors of $\widehat{\Sigma}-I$ with eigenvalue more than $\frac{1}{\beta}\varepsilon$. Then, there exists an algorithm FilterMeanManyEig that given $S,\varepsilon,\gamma,\delta,\alpha,\beta$ outputs one of two possible outcomes:

If $\dim(V)\geq C_{1}\beta\log(1/\varepsilon)$, then it outputs an $S^{\prime}$ so that $\Delta(S^{\prime},G_{0})<\Delta(S,G_{0})$.

Otherwise, the algorithm outputs “OK”, and outputs an orthonormal basis for $V$.

Our algorithm works as follows: It finds all large eigenvalues of $\widehat{\Sigma}-I$, and if there are too many, produces an explicit degree-2 polynomial which, as we will argue, produces a valid filter. The formal pseudocode for our algorithm is in Algorithm 4.

For clarity of exposition, we defer the proof of Theorem 7 to Appendix C.

2 Returning an Estimate When There are Few Large Eigenvalues

At this point, we have run the filter of Algorithm 4 until there are few large eigenvalues. In the subspace with large eigenvalues, we again run the low dimensional algorithm to obtain an estimate for the mean in this subspace. Recall that Lemma 3 guarantees the accuracy of this estimator within this subspace. In the complement of this subspace, where the empirical covariance is very close to the identity, Lemma 5 (stated below) shows that the empirical mean is close to the true mean. This leads to a simple algorithm which outputs an estimate for the mean, described in Algorithm 5.

Let $\mu,\eta,G_{0},S$ be as in Theorem 7. Let $\widehat{\mu}$ be the sample mean of $S$, and let $v$ be a unit vector. Suppose that $\langle v,\mu-\widehat{\mu}\rangle>\frac{\varepsilon}{\beta^{1/2}}$. Then $\mathop{\mbox{\bf Var}}_{S}[\langle v,X\rangle]>1+\frac{\varepsilon}{\beta}$.

For clarity of exposition, we defer the proof of Lemma 5 to Appendix C.

3 The Full High-Dimensional Algorithm

We now have almost all the pieces needed to prove the full result. The last ingredient is the fact that, given enough samples, the good set condition is satisfied by the samples from the true distribution. Formally,

Fix $\eta,\delta>0$. Let $X_{1},\ldots,X_{n}$ be independent samples from $\mathcal{N}(\mu,I)$, where $n=\Omega((d\log(d/\eta\delta))^{6}/\eta^{2})$. Then, $S=\{X_{1},\ldots,X_{n}\}$ is $(\eta,\delta)$-good with respect to $\mathcal{N}(\mu,I)$ with probability at least $1-\delta$.

This follows from Lemmas 8.3 and 8.16 of [DKK+16]. ∎

4 An Extension, with Small Spectral Noise

For learning of arbitrary Gaussians, we will need a simple extension that allows us to learn the mean even in the presence of some spectral norm error in the covariance matrix. Since the algorithms and proofs are almost identical to the techniques above, we omit them for conciseness. Formally, we require:

Fix $\chi,\varepsilon,\delta>0$, and let $X_{1},\ldots,X_{n}$ be an $\varepsilon$-corrupted set of points from $\mathcal{N}(\mu,\Sigma)$, where $\|\Sigma-I\|_{2}\leq O(\chi)$, $\|\mu\|_{2}\leq O(\varepsilon\log 1/\varepsilon)$, and where $n=\operatorname*{poly}(d,1/\chi,1/\varepsilon,\log 1/\delta)$. For any $\gamma>0$, there is an algorithm $\textsc{RecoverMeanNoisy}(X_{1},\ldots,X_{n},\varepsilon,\delta,\gamma,\chi)$ which runs in time $\operatorname*{poly}(d,1/\chi,1/\varepsilon,\log 1/\delta)$ and outputs a $\widehat{\mu}$ so that with probability $1-\delta$, we have $\|\widehat{\mu}-\mu\|_{2}\leq(C+\gamma)\varepsilon+O(\chi)$.

This extension follows from two elementary observations:

For the learning in low dimensions, observe that the median is naturally robust to error in the covariance, and in general, by the same calculation we did, the error of the median becomes $O(\varepsilon+\alpha)$.

For the filter, observe that we only need concentration of squares of linear functions, and whatever error we have in this concentration goes directly into our error guarantee. Thus, by the same calculations that we had above, if we filtered for eigenvalues above $1+O(\varepsilon+\alpha)$, we would immediately get the desired bound.

Robustly Estimating the Mean of Degree Two Polynomials

In this section, we give robust estimates of $\mathop{\mbox{\bf E}}[p^{2}(X)]$ for degree-$2$ polynomials $p$ in subspaces of small dimension, which is an important prerequisite to learning the covariance in high-dimensions. A crucial ingredient in our algorithm is the following improvement theorem (stated and proved in the next section) which shows how to take any weak high-dimensional estimate for the covariance and use it to get an even better robust estimate for $\mathop{\mbox{\bf E}}[p^{2}(X)]$.

Here we give some additional preliminaries we require for the low-dimensional learning algorithm we present here. We will need the following well-known tail bound for degree-$2$ polynomials:

We will also require the following lemmata:

Let $A,B$ be matrices. Then, for all $p,q$ so that $\frac{1}{p}+\frac{1}{q}=1$, we have $\|AB\|_{S^{1}}\leq\|A\|_{S^{p}}\|B\|_{S^{q}}$.

Let $\Sigma,\widehat{\Sigma},M$ be so that $\|\Sigma-\widehat{\Sigma}\|_{F}\leq O(\delta)$, and so that $\|M\|_{F}=1$. Then, we have $\|\Sigma^{1/2}M\Sigma^{1/2}-\widehat{\Sigma}^{1/2}M\widehat{\Sigma}^{1/2}\|_{S^{1}}\leq 5\delta$.

We will bound the first term on the RHS by $5\delta/2$; the second term is bounded symmetrically. We have

where the last line follows from Hölder’s inequality for Schatten norms. ∎

2 An Improvement Theorem

Here we state and prove one of the main technical ingredients in our algorithm for robustly learning the covariance.

Fix $\varepsilon,\delta,\tau>0$. Let $\Sigma$ be so that $\|\Sigma-I\|_{F}\leq O(\varepsilon\log 1/\varepsilon)$, and fix a $p\in\mathcal{P}_{2}$, where $\mathcal{P}_{2}$ denotes the set of even degree-$2$ polynomials in $d$ variables. Let $G_{0}$ be an $(\varepsilon,\delta)$-good set of samples from $\mathcal{N}(0,\Sigma)$, and let $S=\{X_{1},\ldots,X_{n}\}$ be so that $\Delta(S,G_{0})\leq\varepsilon$. Then, for any $C>0$ there is an algorithm LearnMeanChiSquared which, given $p,X_{1},\ldots,X_{n}$, and $\varepsilon$, outputs a $\widehat{\mu}$ so that with probability $1-\tau$ over the randomness of the algorithm,

Moreover, the algorithm runs in time $O(|S|+\log(1/\tau)/\varepsilon^{2})$.

The way to think about how this result fits into the overall strategy is that robustly estimating the covariance is equivalent to robustly estimating the mean of every (normalized) degree-two polynomial $p$. The above theorem shows how a weak estimate in high-dimensions can be used to obtain stronger estimates in one dimension, which ultimately we will use to improve the high-dimensional estimate as well. The above theorem is the workhorse in our proof.

Our algorithm itself is simple, however, its correctness is quite non-trivial. We define some threshold $T$. Given our corrupted set of samples from $\mathcal{N}(0,\Sigma)$, we use our corrupted data set to estimate the mean of $p(X)$ conditioned on the event that $|p(X)|\leq T$. Then, to estimate the contribution of the mean from points $X$ so that $|p(X)|>T$, we estimate this by $\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,I)}[p(X)1_{|p(X)|>T}]$. In other words, we are replacing the contribution of the true tail by an estimate of the contribution of $p(X)$ when $X\sim\mathcal{N}(0,I)$ on this tail. The formal pseudocode is given in Algorithm 6.

Intuitively, this algorithm works because of two reasons. First, it is not hard to show that the influence of points $p(X)$ within the threshold $T$ on the estimator are bounded by at most $T$. Hence, the adversary cannot add corrupted points within this threshold and cause our estimator to deviate too much. Secondly, because we know that $\|\Sigma-I\|_{F}$ is small, by carefully utilizing smoothness properties of sums of chi-squared random variables, we are able to show that our estimate for the contribution of the tail is not too large. At a high level, this is because “most” of the distance between two chi-squared random variables must remain close to the means, so the difference in the tails is much smaller. Proving that this holds in a formal sense is the majority of the technical work of this section.

We know the distribution of $p(X^{\prime})$ for $X^{\prime}\sim N(0,I)$ explicitly and wish to use this to get a better estimate for the mean of $p(X)$ for $X\sim N(0,\Sigma)$ than might be given by the mean of the $\varepsilon$-corrupted set of samples.

It follows from $(\varepsilon,\delta)$-goodness that $|\Pr_{X^{\prime}\sim\mathcal{N}(0,\Sigma)}[p(X^{\prime})>t]-\#\{X_{i}:p(X_{i})>t\}/n|\leq 2\varepsilon$ for all $t$. We need to express the expectation in terms that we can use this to bound. For $Z=p(X^{\prime})$, we have that

Thus, we have $|\mathop{\mbox{\bf E}}[Z-f(Z)]-\alpha|\leq 2T\varepsilon$.

Since $p\in\mathcal{P}_{2}$, we have $\mathop{\mbox{\bf E}}[p(X^{\prime})]=1$ for $X^{\prime}\sim N(0,I)$. Thus, we have $\mathop{\mbox{\bf Var}}[f(p(X^{\prime}))]\leq\mathop{\mbox{\bf E}}[f(p(X^{\prime}))^{2}]\leq\mathop{\mbox{\bf E}}[p(X^{\prime})^{2}]=1$. It follows by standard concentration results that the empirical after taking $m=O(\ln(1-\tau)/\varepsilon^{2})$ samples has $|\sum_{i=1}^{m}f(p(X_{i})^{\prime})/n-\mathop{\mbox{\bf E}}[p(X^{\prime})]|\leq\varepsilon$ with probability $1-\tau$. When this holds, we have

To prove the correctness of the algorithm it remains to show that:

For any constant $C>0$, for $T=O(\log C)$, we obtain

where $O$ is orthogonal, $D$ is diagonal, $a_{i}$ are the eigenvalues of $M$, and $Y\sim\mathcal{N}(0,I)$ hence $Y_{i}$ are i.i.d. from $\mathcal{N}(0,1)$. Since $\|M\|_{F}=1$, here we have $\sum a_{i}^{2}=1$. If instead we express $p(X^{\prime})$ in terms of $X^{\prime}\sim N(0,\sigma^{2})$, we obtain:

where $O^{\prime}$ is orthogonal, $D^{\prime}$ is diagonal, $b_{i}$ are the eigenvalues of $\Sigma^{1/2}M\Sigma^{1/2}$, and $Y^{\prime},Y\sim\mathcal{N}(0,I)$ hence $Y_{i}$ are i.i.d. from $\mathcal{N}(0,1)$.

By Corollary 3, we have that $\sum_{i}|b_{i}-a_{i}|\leq(5/2)\|\Sigma-I\|_{F}$. Now consider the random variables

for $1\leq i\leq d$ and $0\leq\lambda\leq 1$. Note that $Z_{i,1}=Z_{i+1,0}$, for $1\leq i\leq d-1$. Note that, to prove the lemma, it suffices to show that

where $P_{i}(x)$ is the probability density function of the random variable $Z_{i}$. Standard results about the $\chi^{2}$ distribution give that:

Let $Y_{1},Y_{2},Y_{3}\sim\mathcal{N}(0,1)$. Then the probability density function of $Y_{1}^{2}$ is $\frac{1}{\sqrt{2\pi x}}e^{-x/2}$, of $Y_{1}^{2}+Y_{2}^{2}$ is $\frac{1}{2}e^{-x/2}$, and of $Y_{1}^{2}+Y_{2}^{2}+Y_{3}^{2}$ is $\frac{\sqrt{x}}{\sqrt{2\pi}}e^{-x/2}$.

This gives that $P_{i}(x)=\frac{1}{\sqrt{2\pi xc_{i}}}e^{-x/2c_{i}}.$ Now consider the derivative:

where $P^{(3)}_{i}(x)$ is the distribution of $Z_{i}+Z^{\prime}_{i}+Z^{\prime\prime}_{i}$, where $Z^{\prime}_{i}$ and $Z^{\prime\prime}_{i}$ are i.i.d. copies of $Z_{i}$. We thus have

Since $f$ has Lipschitz constant $1$, $|f(Z_{i,\lambda})-\mathop{\mbox{\bf E}}_{Z^{\prime}_{i},Z^{\prime\prime}_{i}}[f(Z_{i,\lambda}+Z^{\prime}_{i}+Z^{\prime\prime}_{i}]|\leq\mathop{\mbox{\bf E}}[Z^{\prime}_{i}+Z^{\prime\prime}_{i}]=2c_{i}$ whatever value $Z_{i,\lambda}$ takes.

Using the probability distribution of $Z^{\prime}_{i}+Z^{\prime\prime}_{i}$, in the case $c_{i}>0$, we have for all $|z|\leq T$,

Since $c_{i}\leq\max{a_{i},b_{i}}\leq 1+O(\varepsilon\log(1/\varepsilon))\leq 2$, for large enough $T=O(\log C)$, we have $e^{-(T/2)/2c_{i}}\leq 1/C^{2}T$. We assume that this holds.

For $-T\leq z\leq T/2$, we have $2c_{i}e^{-(T-z)/2c_{i}}\leq 2c_{i}/C^{2}T$. A similar argument when $c_{i}<0$ gives that for $-T/2\leq z\leq T$, we have $2|c_{i}|e^{-(T-z)/2c_{i}}\leq 2|c_{i}|/C^{2}T$. Now we have enough to show that

By the Hanson-Wright inequality, we have that for any $x$,

Finally, we have that, assuming $c_{i}\neq 0$,

Recalling that $Z_{i,1}=Z_{i+1,0}$ with $Z_{1,0}=p(X^{\prime})$ and $Z_{d,1}=p(X)$ for $X^{\prime}\sim N(0,1)$ and $X\sim N(0,\Sigma)$, we have by the Mean Value Theorem, that

This completes the proof of the theorem. ∎

3 Working in a Low-Dimensional Space of Degree-Two Polynomials

We now show that via similar techniques as before, we can patch our estimates together to find a matrix which agrees with the ground truth on all degree-two polynomials in a fixed subspace of low dimension. Formally, we show:

Fix $\varepsilon,\tau>0$. Let $\Sigma$ be so that $\|\Sigma-I\|_{F}\leq O(\varepsilon\log 1/\varepsilon)$. Let $G_{0}$ be an $(\varepsilon,\delta)$-good set of samples from $\mathcal{N}(0,\Sigma)$, and let $S=\{X_{1},\ldots,X_{n}\}$ be so that $\Delta(S,G_{0})\leq\varepsilon$. Let $W_{1}$ be a subspace of degree-$2$ polynomials, and let $W_{2}$ be an orthogonal subspace of degree-$2$ polynomials, so that we have a $\widehat{\Sigma}$ so that $\left|\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[p(X)]-\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\widehat{\Sigma})}[p(X)]\right|\leq\xi$ for all $p\in W_{2}$. Then there is an algorithm LearnMeanPolyLowD which given $\varepsilon,S,W_{1},W_{2},\widehat{\Sigma}$ runs in time $\operatorname*{poly}(d,|S|,2^{O(\dim(W_{1}))},\log 1/\tau)$, and returns a $\Sigma^{\prime}$ so that

Observe that the dimension of the space of degree-$2$ polynomials $W$ in $V$ is $O(\dim(V)^{2})$. Run the algorithm in Theorem 10 with the same parameters as before, with $W_{1}=W$ and $W_{2}=\emptyset$ (so that we may take $\xi=0$), and then the guarantee of that algorithm, along with Lemma 10, gives our desired guarantee. ∎

We now describe the algorithm for Theorem 10. Essentially, we do the same thing as we did for low-dimensional learning in the unknown mean case: we take a constant net over $V\cap\mathcal{P}_{2}$, l earn the mean over every polynomial in the net, and then find a $\Sigma^{\prime}$ which is close in each direction to the learned mean. Since we will not attempt to optimize the constant factor here, will will use a naive LP-based approach to find a point which is close to optimal. The formal pseudocode is given in Algorithm 7.

Observe that every constraint for each polynomial in $W_{1}$ is indeed linear in $\Sigma^{\prime}$, by Lemma 10. Moreover, the constraint for $W_{2}$ has an explicit separation oracle, since it induces a norm, and for any $p\in W_{2}$, we may explicitly compute $\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma^{\prime})}[p(X)]-\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma^{\prime})}[p(X)]$. Thus, we may use separating hyperplane techniques to solve this convex program in the claimed running time.

Let us condition on the event that LearnMeanChSquared succeeds for each $p\in\mathcal{C}$. By a union bound, this occurs with probability at least $1-\tau$. Thus, in each $p\in\mathcal{C}$, we have that $|m_{p}-\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[p(X)]|\leq\beta$, where $\beta=\|\Sigma-I\|_{F}/C+O(\log(C))\varepsilon$. Let $\Sigma^{\prime}$ be the matrix we find. By the triangle inequality, we then have that for every $p\in\mathcal{C}$, that $|\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma^{\prime})}[p(X)]-\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma)}[p(X)]|\leq 2\beta$. Hence, by the usual net arguments, we know that for every $p\in V\cap\mathcal{P}_{2}$,

Moreover, by triangle inequality, for every $p\in W_{2}$, we have $\left|\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma^{\prime})}[p(X)]-\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma^{\prime})}[p(X)]\right|\leq 2\xi$. The result then follows from the Pythagorean theorem. ∎

Robustly Learning the Covariance in High-Dimensions

In this section, we show how to robustly estimate the covariance of a mean-zero Gaussian in high-dimensions up to error $O(\varepsilon)$. We use our low-dimensional learning algorithm from the previous section as a crucial subroutine in what follows.

Our main algorithmic contribution is as follows:

Fix $\varepsilon,\delta>0$, and let $S_{0}=(G_{0},E_{0})$ be an $\varepsilon$-corrupted set of samples of size $n$ from $\mathcal{N}(0,\Sigma)$, where $\|\Sigma-I\|_{F}\leq\xi$ where $\xi=O(\varepsilon\log 1/\varepsilon)$, and where $n=\operatorname*{poly}(d,1/\varepsilon,\log 1/\delta)$. Suppose that $G_{0}$ is $(\varepsilon,\delta)$-good with respect to $\mathcal{N}(0,\Sigma)$. Let $S\subseteq S_{0}$ be a set so that $\Delta(S,G_{0})\leq\varepsilon$. Then, there exists an algorithm ImproveCov that given $S,\xi,\varepsilon$, fails with probability at most $\operatorname*{poly}(\varepsilon,1/d,\delta)$, and otherwise outputs one of two possible outcomes:

A matrix $\widehat{\Sigma}$, so that $\|\widehat{\Sigma}-\Sigma\|_{F}\leq\|\Sigma-I\|_{F}/2$.

A set $S^{\prime}\subset S$ so that $\Delta(S^{\prime},G_{0})<\Delta(S,G_{0})$.

Moreover, ImproveCov runs in time $\operatorname*{poly}(d,(1/\varepsilon)^{O(\log^{4}1/\varepsilon)},\log 1/\delta)$.

By first applying the algorithm in [DKK+16] to produce an initial estimate for $\Sigma$, and then iterating the above algorithm polynomially many times, this immediately yields:

Fix $\varepsilon,\delta>0$, and let $G_{0}$ be a set of $i.i.d.$ samples from $\mathcal{N}(0,\Sigma)$, where $n=\operatorname*{poly}(d,1/\varepsilon,\log 1/\delta)$. Let $S$ be so that $\Delta(S,G_{0})\leq\varepsilon$. There is an universal constant $C$ and an algorithm which outputs a $\widehat{\Sigma}$ so that with probability $1-\delta$, we have $\|\widehat{\Sigma}^{-1/2}\Sigma\widehat{\Sigma}^{-1/2}-I\|_{F}\leq C\varepsilon$. In particular, this implies that

Our strategy for obtaining a high-dimensional estimate for the covariance based on solving low-dimensional subproblems will be substantially more challenging than it was for the unknown mean case. The natural approach is to take the $\operatorname*{poly}\log(1/\varepsilon)$-dimensional subspace of degree-$2$ polynomials of largest empirical variance and construct a filter. However, this fails because, unlike in the mean case, we do not know the variance of these degree-$2$ polynomials to small error. For the unknown mean case, because we assumed that we knew the covariance was the identity (or spectrally close to the identity), this was not an issue. Now, the variance of our polynomials depends on the (unknown) covariance of the true Gaussian, which may be more than $O(\varepsilon)$-far from our current estimate. Indeed, it is not difficult to come up with counterexamples where there are many large eigenvalues of the empirical covariance matrix, but no filter can make progress.

Then, in Section 6.3, we show that if we restrict to the orthogonal subspace, i.e., the subspace where the empirical covariance matrix does not have large eigenvalues, we can indeed either produce a filter or improve our estimate of the covariance restricted to this subspace using our low-dimensional estimator. While the blueprint is similar to the filter for the unknown mean, the techniques are much more involved and subtle.

Supposing we have not yet created a filter, we have now estimated the covariance on a poly-logarithmic dimensional subspace $V$, and on $V^{\perp}$. This does not in general imply that we have learned the covariance in Frobenius norm. In block form, if we write

In Section 6.4, we show, given a polylogarithmically sized subspace $V$, and a good estimate of the covariance matrix on $V$ and $V^{\perp}$, how to fill in the entire covariance matrix. Roughly, we do this by randomly fixing directions in $V$, and performing rejection sampling based on the correlation in the direction in $V$, and showing that the problem reduces to one of robustly learning the mean of a Gaussian, which (conveniently) we have already solved. These steps together yield our overall algorithm ImproveCov. Finally, in Section 6.6 we explain why there is a natural barrier that makes reducing the running time from quasi-polynomial to polynomial (in $1/\varepsilon$) difficult.

2 Additional Preliminaries

Here we give some additional preliminaries we will require in this Section.

We also require the following classical result, which allows us to do agnostic hypothesis selection with corrupted samples (see e.g., [DL01, DDS12, DK14, DDS15]).

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$. However, if the algorithm is not given this information, and instead given an $\eta$ and asked to return something with error at most $O(\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(\varepsilon+(1+\gamma)\eta)$.