Robustly Learning a Gaussian: Getting Optimal Error, Efficiently

TL;DR

This paper develops efficient algorithms for robustly estimating the parameters of high-dimensional Gaussian distributions with adversarial noise, achieving near-optimal error bounds with polynomial sample complexity.

Contribution

It introduces robust estimators that attain optimal or near-optimal error bounds in high-dimensional Gaussian parameter estimation under adversarial noise, with polynomial runtime.

Findings

Achieves estimation error $O(\varepsilon)$ in total variation distance.

Provides polynomial-time algorithms for mean-only and mean-and-covariance estimation.

Requires only polynomial samples, matching classical one-dimensional guarantees.

Abstract

We study the fundamental problem of learning the parameters of a high-dimensional Gaussian in the presence of noise -- where an $\varepsilon$-fraction of our samples were chosen by an adversary. We give robust estimators that achieve estimation error $O(\varepsilon)$ in the total variation distance, which is optimal up to a universal constant that is independent of the dimension. In the case where just the mean is unknown, our robustness guarantee is optimal up to a factor of $\sqrt{2}$ and the running time is polynomial in $d$ and $1/\epsilon$. When both the mean and covariance are unknown, the running time is polynomial in $d$ and quasipolynomial in $1/\varepsilon$. Moreover all of our algorithms require only a polynomial number of samples. Our work shows that the same sorts of error guarantees that were established over fifty years ago in the one-dimensional setting can also be achieved by efficient algorithms in high-dimensional settings.

Equations

$$\mathcal{D}=(1-\epsilon )\mathcal{N}(\mu ,{\sigma }^2)+\epsilon \mathcal{Z},$$

$$\hat{\mu }=\text{mbox}median(X_1,X_2,\dots ,X_n)\text{mbox}and\hat{\sigma }=\frac{\text{mbox}median(|X_i-\hat{\mu }|)}{{\Phi }^{-1}(3/4)},$$

$$d_{TV}(\mathcal{N}(\mu ,\Sigma ),\mathcal{N}(\hat{\mu },\hat{\Sigma }))\le C\epsilon$$

$$d_{TV}(\mathcal{N}(\mu ,I),\mathcal{N}(\hat{\mu },I))\le \left(\frac{1}{\sqrt{2}}+O\left(\frac{1}{\sqrt{\beta }}+{\epsilon }^2\right)\right)\epsilon .$$

$${\text{mbox}\mathbf{E}}_{X\sim \mathcal{N}(0,\Sigma )}[p(X)]=1\text{mbox}and{\text{mbox}\mathbf{E}}_{X\sim \mathcal{N}(0,\hat{\Sigma })}[p(X)]=1+C^{\prime }\epsilon .$$

$$d_{TV}(\mathcal{N}(0,\Sigma ),\mathcal{N}(0,\hat{\Sigma }))\le C\epsilon ,$$

$$\Sigma =\left[\begin{array}{cc}{\Sigma }_V& A^T\\ A& {\Sigma }_{V^{\perp }}\end{array}\right],$$

$$d_{TV}(\mathcal{N}(\mu ,\Sigma ),\mathcal{N}(\hat{\mu },\hat{\Sigma }))\le C\epsilon ,$$

$$\varphi (S^{\prime },G)=\frac{|G\setminus G^{\prime }|}{|S^{\prime }|},\psi (S^{\prime },G)=\frac{|E^{\prime }|}{|S^{\prime }|},\Delta (S^{\prime },G)=\psi (S^{\prime },G)+\varphi (S^{\prime },G)\log \frac{1}{\varphi (S^{\prime },G)}.$$

$$\left|{\text{mbox}\mathbf{E}}_G[p(X)]-{\text{mbox}\mathbf{E}}_{\mathcal{N}(\mu ,I)}[p(X)]\right|$$

$$\left|{\text{mbox}\mathbf{E}}_G[p^2(X)]-{\text{mbox}\mathbf{E}}_{\mathcal{N}(\mu ,I)}[p^2(X)]\right|$$

$$\underset{G}{\mathrm{Pr}}[p(X)\ge 0]$$

$$\left|{\text{mbox}\mathbf{E}}_G[p(X)]-{\text{mbox}\mathbf{E}}_{\mathcal{N}(0,\Sigma )}[p(X)]\right|$$

$$\left|{\text{mbox}\mathbf{E}}_G[p^2(X)]-{\text{mbox}\mathbf{E}}_{\mathcal{N}(0,\Sigma )}[p^2(X)]\right|$$

$$\underset{X\sim G}{\mathrm{Pr}}[p(X)\ge 0]$$

$$\left|{\text{mbox}\mathbf{E}}_G[p(X)]-{\text{mbox}\mathbf{E}}_{\mathcal{N}(0,\Sigma )}[p(X)]\right|$$

$$\underset{G}{\mathrm{Pr}}[p(X)\ge 0]$$

$$d_{\mathrm{TV}}\left(\mathcal{N}({\mu }_1,I),\mathcal{N}({\mu }_2,I)\right)=\left(\frac{1}{\sqrt{2\pi }}+o(1)\right)\epsilon .$$

$$\Vert {\Pi }_V(\mu -\tilde{\mu }){\Vert }_2=\frac{1+2\rho }{1-\rho }\left(\sqrt{\pi }+O\left(\frac{\gamma }{d}\right)\right)\epsilon .$$

$$\underset{X\sim S}{\mathrm{Pr}}[⟨v,X⟩>a]=\frac{|G|}{|S|}\underset{X\sim G}{\mathrm{Pr}}[⟨v,X⟩>a]+\frac{|E|}{|S|}\underset{X\sim E}{\mathrm{Pr}}[⟨v,X⟩>a].$$

$$\left|\underset{X\sim G}{\mathrm{Pr}}[⟨v,X⟩>a]-\underset{X\sim G_0}{\mathrm{Pr}}[⟨v,X⟩>a]\right|$$

$$\left|\underset{X\sim S}{\mathrm{Pr}}[⟨v,X⟩>a]-\frac{|G|}{|S|}\underset{X\sim G_0}{\mathrm{Pr}}[⟨v,X⟩>a]\right|\le \psi (S,G)+o(\epsilon ),$$

$$\left|\underset{X\sim S}{\mathrm{Pr}}[⟨v,X⟩<a]-\frac{|G|}{|S|}\underset{X\sim G_0}{\mathrm{Pr}}[⟨v,X⟩<a]\right|\le \psi (S,G)+o(\epsilon ).$$

$$\left|\underset{X\sim S}{\mathrm{Pr}}[⟨v,X⟩>{\Pi }_y\mu +a]-\frac{|G|}{|S|}\left(\frac{1}{2}-\frac{1}{\sqrt{2\pi }}a\right)\right|\le \psi (S,G)+O\left(\frac{\gamma \epsilon }{d}\right).$$

$$\underset{X\sim S}{\mathrm{Pr}}[⟨v,X⟩>{\Pi }_y\mu +\sqrt{\frac{\pi }{2}}a]$$

$$\le \frac{1}{2}-\frac{|G|}{|S|}\frac{1}{\sqrt{2\pi }}a+\frac{|E|}{2|S|}+O\left(\frac{\gamma \epsilon }{d}\right)+o(\epsilon ).$$

$$\mathcal{C}=\bigcap _{v\in \mathcal{F}}\{y\in {\mathbb{R}}^d:⟨v,y⟩\in \left[b_v-\beta ,b_v+\beta \right]\},$$

$$|⟨x-y,u⟩|$$

$$\le 2\beta +\rho \Vert x-y{\Vert }_2.$$

$$\min \Vert x-y{\Vert }_2$$

$$n=\Omega \left(\frac{(d\log (d/\gamma \epsilon \delta ){)}^6}{{\gamma }^2{\epsilon }^2}\right).$$

$$\mathrm{Pr}\left(\left|X^{T}AX-\text{mbox}\mathbf{E}[X^{T}AX]\right|>t\right)\le 2\exp \left(-c_0\cdot \min \left(\frac{t^2}{\Vert A{\Vert }_F^2},\frac{t}{\Vert A{\Vert }_2}\right)\right).$$

$$\Vert {\Sigma }^{1/2}M{\Sigma }^{1/2}-{\hat{\Sigma }}^{1/2}M{\hat{\Sigma }}^{1/2}{\Vert }_{S^1}$$

$$\Vert {\Sigma }^{1/2}M{\Sigma }^{1/2}-M{\Vert }_{S^1}$$

Full text

In the case where just the mean is unknown, our robustness guarantee is optimal up to a factor of $\sqrt{2}$ and the running time is polynomial in $d$ and $1/\varepsilon$. When both the mean and covariance are unknown, the running time is polynomial in $d$ and quasipolynomial in $1/\varepsilon$. Moreover all of our algorithms require only a polynomial number of samples. Our work shows that the same sorts of error guarantees that were established over fifty years ago in the one-dimensional setting can also be achieved by efficient algorithms in high-dimensional settings.

1 Introduction

1.1 Background

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:

Definition 1**.**

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

[TABLE]

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 challenging111None 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

[TABLE]

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 for222See 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

[TABLE]

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 problem333 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.

1.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:

Theorem 1**.**

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

[TABLE]

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 satisfies444More 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

[TABLE]

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:

Theorem 2**.**

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

[TABLE]

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

[TABLE]

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:

Corollary 1**.**

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

[TABLE]

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.

More Related Work

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.

1.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.

2 Preliminaries

In this section, we give various definitions and lemmata we will require throughout the paper. First, given a distribution $F$, we let $\mathop{\mbox{\bf E}}_{F}[f(X)]=\mathop{\mbox{\bf E}}_{X\sim F}[f(X)]$ denote the expectation of $f(X)$ under $F$. If $S$ is a finite set, we let $\mathop{\mbox{\bf E}}_{S}[f(X)]=\mathop{\mbox{\bf E}}_{X\sim\mathrm{unif}(S)}[f(X)]$ denote the expectation of $f(X)$ under the uniform distribution over points in $S$ (i.e., the empirical mean of $f$ under $S$). Given any subspace $V\subseteq\mathbb{R}^{d}$, we let $\Pi_{V}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ be the projection operator onto $V$. If $V=\mathrm{span}(v)$ is 1-dimensional, we will denote this projection as $\Pi_{v}$.

2.1 The Strong Contamination Model

Here we formally define the strong contamination model.

Definition 2**.**

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$. 2. 2.

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

[TABLE]

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.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.

2.2.1 Regularity Conditions for Unknown Mean

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

Definition 3**.**

Let $G$ be a multiset of points in $\mathbb{R}^{d}$ and $\eta,\delta>0$. We say that $G$ is $(\eta,\delta)$-good with respect to $\mathcal{N}(\mu,I)$ if the following hold:

• •

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

• •

For every affine function $L:\mathbb{R}^{d}\rightarrow\mathbb{R}$ we have $|\Pr_{G}(L(X)\geq 0)-\Pr_{\mathcal{N}(\mu,I)}(L(X)\geq 0)|\leq\eta/(d\log(d/\eta\delta))\;.$

• •

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

• •

We have that $\left\|\mathrm{Cov}_{G}[X]-I\right\|_{2}\leq\eta/d.$

• •

For any even degree-$2$* polynomial $p:\mathbb{R}^{d}\rightarrow\mathbb{R}$ we have that*

[TABLE]

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.2 Regularity Conditions for Unknown Covariance

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

Definition 4**.**

Let $G$ be a set of $n$ points of $\mathbb{R}^{d}$, and $\eta,\delta>0$. We say that $G$ is $(\eta,\delta)$-good with respect to $\mathcal{N}(0,\Sigma)$ if the following hold:

• •

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

• •

For any even degree-$2$* polynomial $p:\mathbb{R}^{d}\rightarrow\mathbb{R}$ we have*

[TABLE]

• •

For any even degree-$4$* polynomial $p:\mathbb{R}^{d}\rightarrow\mathbb{R}$ we have*

[TABLE]

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.

2.3 Bounds on the Total Variation Distance

We will require some simple bounds on the total variation distance between two Gaussians. These bounds are well-known. Roughly speaking, they say that the total variation distance between two Gaussians with identity covariance is governed by the $\ell_{2}$ norm between their means, and the total variation distance between two Gaussians with mean zero is governed by the Frobenius norm between their covariance matrices, provided that the matrices are close to the identity.

Lemma 1**.**

Let $\mu_{1},\mu_{2}\in\mathbb{R}^{d}$ be such that $\|\mu_{1}-\mu_{2}\|_{2}=\varepsilon$ for $\varepsilon<1$. Then

[TABLE]

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):

Lemma 2** (Cor. 2.14 in [DKK*+*16]).**

Let $\Sigma,\widehat{\Sigma}$ be such that $\|\Sigma-I\|_{F}\leq O(\varepsilon\log 1/\varepsilon)$, and $\|\Sigma-\widehat{\Sigma}\|_{F}\leq C\varepsilon$. Then $d_{\mathrm{TV}}(\mathcal{N}(0,\Sigma),\mathcal{N}(0,\widehat{\Sigma}))\leq O(\varepsilon)$.

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.

3 Robustly Learning the Mean in Low Dimensions

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

Theorem 3**.**

Fix $\mu\in\mathbb{R}^{d}$, and let $\varepsilon,\gamma,\delta>0$. Let $S_{0}=(G_{0},E_{0})$ be such that $G_{0}$ is a $(\gamma\varepsilon,\delta)$-good set with respect to $\mathcal{N}(\mu,I)$, and $|E_{0}|/|S_{0}|\leq\varepsilon$. Let $S=(G,E)$ be another set with $\Delta(S,S_{0})<\varepsilon$. Let $V\subseteq\mathbb{R}^{d}$ be a subspace. For all $0<\rho<1$, the algorithm $\textsc{LearnMeanLowD}(V,\gamma,\varepsilon,\delta,S,\rho)$ runs in time $\operatorname*{poly}(d,|S|,(1/\rho)^{O(\dim(V))},\log(\rho\varepsilon/(1-\rho)),\log(1/\rho))$ and returns a $\widetilde{\mu}$ so that

[TABLE]

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.

For simplicity, in the rest of this section, we will let $V=\mathbb{R}^{d}$, that is, we assume there is no projection. It should be clear that this can be done without loss of generality. Our algorithm proceeds as follows: First, we show that in one dimension, the median produces an estimate which is optimal, up to lower order terms, if the sample set is $(\gamma\varepsilon,\delta)$-good with respect to the underlying Gaussian. Then, we show that by using a net argument, we can produce a convex body in $\mathbb{R}^{d}$ with diameter at most $2R=2(\sqrt{\frac{\pi}{2}}+o(1))\varepsilon$ which must contain the true mean. Finally, we use an old result of Jung [Jun01] that such a set can be circumscribed by a ball of radius $\sqrt{2}R$ (see [BW41] for an English language version of the result). We use the center of the ball as our estimate $\widetilde{\mu}$.

3.1 Robustness of the Median

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:

Lemma 3**.**

Fix any $v\in\mathbb{R}^{d}$. Fix $\mu\in\mathbb{R}^{d}$, and let $\delta>0$. Let $S_{0}=(G_{0},E_{0})$ be so that $G_{0}$ be a $(\gamma\varepsilon,\delta)$-good set with respect to $\mathcal{N}(\mu,I)$, and $|E_{0}|/|S_{0}|\leq\varepsilon$. Let $S=(G,E)$ be another set with $\Delta(S,S_{0})<\varepsilon$. Let $b$ be the median of $S$ when projected onto $v$. Then, $|b-\Pi_{v}\mu|\leq\left(\sqrt{\frac{\pi}{2}}+O\left(\frac{\gamma}{d}\right)\right)\varepsilon$.

Proof.

For any $a\in\mathbb{R}$, we have

[TABLE]

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

[TABLE]

Hence, we have

[TABLE]

since by assumption $\Delta(S,G)\leq\varepsilon$. Similarly, we have that for all $a\in\mathbb{R}$,

[TABLE]

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

[TABLE]

Thus, for $|a|=O(\varepsilon)$ we have

[TABLE]

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)$. ∎

3.2 Finding a Minimum Radius Circumscribing Ball

For any $x\in\mathbb{R}^{d}$ and $r>0$, let $B(x,r)=\{y\in\mathbb{R}^{d}:\|x-y\|_{2}\leq r\}$ denote the closed ball of radius $r$ centered at $x$. The following classical result of Gale gives a bound on the radius of the circumscribing ball of any convex set in terms of its diameter:

Theorem 4** (see [Jun01, BW41]).**

Fix $R>0$. Let ${\cal C}\subseteq\mathbb{R}^{d}$ be a convex body so that for all $x,y\in{\cal C}$, we have $\|x-y\|_{2}\leq 2R$. Then ${\cal C}$ is contained within a ball of radius $R\sqrt{2}$.

The bound is asymptotically achieved for the standard simplex as we increase its dimension. The goal of this subsection is to show that the (approximately) minimum radius circumscribing ball can be found efficiently. We will assume we are given an approximate projection oracle for the convex body that given a point $y\in\mathbb{R}^{d}$, outputs a point which is almost the closest point in ${\cal C}$ to $x$:

Definition 5**.**

A $\rho$-projection oracle for a convex body ${\cal C}$ is a function $\mathcal{O}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, which, given a point $y\in\mathbb{R}^{d}$, outputs a point $x\in{\cal C}$ so that $\|x-y\|_{2}\leq\inf_{x^{\prime}\in{\cal C}}\|x^{\prime}-y\|_{2}+\rho$.

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.

Claim 1**.**

Fix $r>0$. Then, for any $\beta>0$, there is a $\beta$-net $\mathcal{F}$ for the sphere of radius $r$ around [math] in $\mathbb{R}^{d}$ of size $(r/\beta)^{O(d)}$. Moreover, this net can be constructed in time $\operatorname*{poly}(d,|\mathcal{F}|)$.

With this, we can show:

Lemma 4**.**

Fix $R,{\cal C}$ as in Theorem 4, and let $1>\rho>0$. Let $x\in{\cal C}$ be arbitrary. Let ${\cal O}$ be a $(\rho R/3)$-projection oracle for ${\cal C}$. Suppose a call to ${\cal O}$ runs in time $T$. Then, there is an algorithm $\textsc{CircumscribeNet}(R,\rho,{\cal O},x)$ which runs in time $\operatorname*{poly}((R/\rho)^{O(d)},T)$ and outputs a set $\mathcal{X}\subseteq\mathbb{R}^{d}$ so that $\mathcal{X}$ is a $(\rho R)$-net for ${\cal C}$, and moreover, $|\mathcal{X}|\leq(R/\rho)^{O(d)}$.

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.

Proof.

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. ∎

With this, we obtain:

Corollary 2**.**

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.

Proof.

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.3 The Full Low-Dimensional Algorithm

We now have all the tools to describe the full algorithm in low-dimensions. Let $S$ be our corrupted dataset as in Theorem 3. Fix $\rho>0$. We first produce a $\rho$-net $\mathcal{F}$ over the unit sphere in $\mathbb{R}^{d}$. By (a slight modification of) Claim 1, this net has size $(1/\rho)^{O(d)}$ and can be constructed in time $\operatorname*{poly}(d,|\mathcal{F}|)$. For each $v\in\mathcal{F}$, we project all points in $S$ onto $v$, and take the median of these points to produce $b_{v}$. We then construct the following set:

[TABLE]

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:

Claim 2**.**

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

Proof.

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

[TABLE]

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.

Claim 3**.**

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}))$.

Proof.

The projection problem may be stated as

[TABLE]

This is a convex minimization problem with linear constraints. By the classical theory of optimization [GLS88], finding a $\rho$-approximate $y$ can be done in $\operatorname*{poly}(d,\log(\mathrm{vol}({\cal C})/\rho^{\prime}))$ queries to a separation oracle for ${\cal C}$. Since the separation oracle must only consider the constraints induced by $\mathcal{F}$, this can be done in time $(1/\rho)^{O(d)}$. Since by Claim 2 we have $\mathrm{vol}({\cal C})\leq(2\beta/(1-\rho))^{O(d)}$, the desired runtime follows immediately. ∎

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.

Proof.

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. ∎

4 Robustly Learning the Mean in High Dimensions

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

Theorem 5**.**

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

[TABLE]

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:

Theorem 6**.**

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:

1. (i)

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

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.

4.1 Making Progress with Many Large Eigenvalues

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

Theorem 7**.**

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})$. 2. 2.

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.

4.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.

Lemma 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.

4.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,

Lemma 6**.**

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$.

Proof.

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

At this point, we conclude with the proof of Theorem 5. Within the subspace $V$, Lemma 3 guarantees that the mean is accurate up to $\ell_{2}$-error $\frac{\sqrt{\pi}+O(\gamma)}{1-\alpha}\varepsilon$. Within the subspace $V^{\perp}$, the contrapositive of the statement of Lemma 5 guarantees the mean is accurate up to $\ell_{2}$-error $\frac{\varepsilon}{\beta^{1/2}}$. The desired result follows from the Pythagorean theorem.

4.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:

Theorem 8**.**

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)$. 2. 2.

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.

5 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)]$.

5.1 Additional Preliminaries

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:

Lemma 7** (Hanson-Wright Inequality [LM00, Ver10]).**

Let $X\sim\mathcal{N}(0,I)\in\mathbb{R}^{d}$ and $A$ be a $d\times d$ matrix. Then for some absolute constant $c_{0}$, for every $t\geq 0$,

[TABLE]

We will also require the following lemmata:

Lemma 8** (Hölder’s inequality for Schatten norms).**

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}}$.

This implies the following corollary:

Corollary 3**.**

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$.

Proof.

We have

[TABLE]

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

[TABLE]

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

5.2 An Improvement Theorem

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

Theorem 9**.**

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,

[TABLE]

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.

Proof.

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

[TABLE]

Similarly, the samples have

[TABLE]

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

[TABLE]

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

Lemma 9**.**

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

[TABLE]

Proof.

Let $M$ be the symmetric matrix with $\|M\|_{F}=1$ such that $p(x)=x^{T}Mx$ for $x\in\mathbb{R}^{d}$. We can write $p(X)$ with $X\sim\mathcal{N}(0,I)$ as

[TABLE]

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:

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

To this end, consider how $\mathop{\mbox{\bf E}}[f(Z_{i,\lambda})]$ varies with $\lambda$. We can write $Z_{i,\lambda}=Z_{-i}+Z_{i}$, where $Z_{-i}=-\mathrm{tr}(M)+\sum_{j<i}b_{j}Y_{j}^{2}+\sum_{j>i}a_{j}Y_{j}^{2}$; and $Z_{i}=c_{i}Y_{i}^{2}$, where $c_{i}=((1-\lambda)a_{i}+\lambda b_{i})$. We assume for now that $c_{i}\neq 0$. Since only $Z_{i}$ depends on $\lambda$, we have

[TABLE]

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

Fact 1**.**

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:

[TABLE]

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

[TABLE]

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$,

[TABLE]

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

[TABLE]

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

[TABLE]

where $c=(b_{1},\dots,b_{i}-1,c_{i},a_{i+1},\dots,a_{d})$. Note that $\|c\|_{2}\leq 1+O(\|\Sigma-I\|_{F})\leq 2$. Also, $\mathop{\mbox{\bf E}}[Z_{i,\lambda}]+\mathrm{tr}(M)$ is the sum of coordinates of $c$, and so $|\mathop{\mbox{\bf E}}[Z_{i,\lambda}]|\leq\sum_{j=1}^{i}|b_{j}-a_{j}|\leq O(\|\Sigma-I\|_{F})\leq 1$. Putting this together, we obtain

[TABLE]

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

[TABLE]

Now we need to deal with the special case $c_{i}=0$. Note that since $f$ is Lipschitz, for any $\beta\in\mathbb{R}$, we have $|\mathop{\mbox{\bf E}}[f(Z_{-i}+\beta Y_{i}^{2})-f(Z_{-i})]|/|\beta|\leq|\mathop{\mbox{\bf E}}[Y_{i}^{2}]|\leq 1$. By considering the limit as $\beta$ tends to zero from above or below, we get that when $c_{i}=0$, the derivative still exists and $\left|\frac{d\mathop{\mbox{\bf E}}[f(Z_{i,\lambda})]}{d\lambda}\right|\leq 1$. Since the limit $2|(b_{i}-a_{i})|/C^{2}$ only does not apply at one point where $\mathop{\mbox{\bf E}}[f(Z_{i,\lambda})]$ is still continuous as a function of $\lambda$, this is not an issue. We still obtain that

[TABLE]

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

[TABLE]

as required. ∎

This completes the proof of the theorem. ∎

5.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:

Theorem 10**.**

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

[TABLE]

for all $p\in\mathrm{span}(W_{1}\cup W_{2})\cap\mathcal{P}_{2}$, with probability $1-\tau$.

In particular, this implies:

Corollary 4**.**

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)$. Let $S=\{X_{1},\ldots,X_{n}\}$ be so that $\Delta(S,G_{0})\leq\varepsilon$. Let $V$ be a subspace of $\mathbb{R}^{d}$. Then there is an algorithm LearnCovLowDim which given $S,\varepsilon,\xi,\tau,V$ runs in time $\operatorname*{poly}(|S|,2^{O(\dim(V)^{2})},\log 1/\tau)$ and returns a $\Sigma^{\prime}$ so that

[TABLE]

with probability $1-\tau$.

Proof.

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.

Proof of Theorem 10.

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}$,

[TABLE]

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. ∎

6 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:

Theorem 11**.**

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:

1. (i)

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

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:

Corollary 5**.**

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

[TABLE]

6.1 Technical Overview

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.

We overcome this hurdle in several steps. First, in Section 6.3, we show how to find a filter if there are many medium-sized eigenvalues of the empirical covariance matrix. This will proceed roughly in the same way that the filter for the unknown mean does. If no filter is created, then we know there are at most logarithmically large eigenvalues of the empirical covariance. In the subspace $V\subseteq\mathbb{R}^{d}$ spanned by their eigenvectors, we can then learn the covariance to high accuracy using our low-dimensional estimator.

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

[TABLE]

where here $\mathbb{R}^{d}$ is written as $V\oplus V^{\perp}$, this implies we have learned $\Sigma_{V}$ and $\Sigma_{V^{\perp}}$ to high accuracy. Thus, it remains to estimate the cross term $A$.

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.

6.2 Additional Preliminaries

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

6.2.1 The Agnostic Tournament

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]).

Theorem 12**.**

Fix $\varepsilon,\delta>0$. Let $D_{1},\ldots,D_{k},D$ be a set of distributions where $\min_{i}d_{\mathrm{TV}}(D_{i},D)=\gamma$. Set $n=\Omega\left(\frac{\log k+\log 1/\delta}{\varepsilon^{2}}\right)$. There is an algorithm Tournament which given oracles for evaluating the pdfs of $D_{1},\ldots,D_{k}$ along with $n$ independent samples $X_{1},\ldots X_{n}$ from $D$, outputs a $D_{i}$ so that $d_{\mathrm{TV}}(D_{i},D)\leq 3\gamma+\varepsilon$ with probability $1-\delta$. Moreover, the running time and number of oracle calls needed is at most $O(n^{2}/\varepsilon^{2})$.

Remark 1**.**

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)$.

6.2.2 The Fourth Moment Tensor of a Gaussian

As in [DKK*+*16], it will be crucial for us to understand the behavior of the fourth moment tensor of a Gaussian. Let $\otimes$ denote the Kronecker product on matrices. We will make crucial use of the following definition:

Definition 6**.**

For any matrix $M\in\mathbb{R}^{d\times d}$, let $M^{\flat}\in\mathbb{R}^{d^{2}}$ denote its canonical flattening into a vector in $\mathbb{R}^{d^{2}}$, and for any vector $v\in\mathbb{R}^{d^{2}}$, let $v^{\sharp}$ denote the unique matrix $M\in\mathbb{R}^{d\times d}$ so that $M^{\flat}=v$.

We will also require the following definition:

Definition 7**.**

Let \mathcal{S}_{\mathrm{sym}}=\{M^{\flat}\in\mathbb{R}^{d^{2}}:\mbox{M is symmetric}\}.

The following result was proven in [DKK*+*16]:

Theorem 13** (Theorem 4.15 in [DKK*+*16]).**

Let $X\sim\mathcal{N}(0,\Sigma)$. Let $M$ be the $d^{2}\times d^{2}$ matrix given by $M=\mathop{\mbox{\bf E}}[(X\otimes X)(X\otimes X)^{T}]$. Then, as an operator on $\mathcal{S}_{\mathrm{sym}}$, we have

[TABLE]

6.2.3 Polynomials in Gaussian Space

Here we review some basic facts about polynomials under Gaussian measure, which will be crucial for our algorithm for learning Gaussians with unknown covariance. We equip the set of polynomials over $\mathbb{R}^{d}$ with the Gaussian inner product, defined by $\langle f,g\rangle=\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,I)}[f(X)g(X)]$, and we let $\|f\|_{2}^{2}=\langle f,f\rangle$.

For any symmetric $M$ with $\|M\|_{F}=1$, define the degree-$2$ polynomial $p(x)=\frac{1}{\sqrt{2}}(x^{T}Mx-\mathrm{tr}(M))$. We call $p$ the polynomial associated to $M$. Observe that $p$ is even (i.e., has no degree-$1$ terms). We will use the following properties of such polynomials:

Lemma 10**.**

Let $M$ be symmetric, so that $\|M\|_{F}=1$. Let $p$ be its associated polynomial. Then, we have:

1. (i)

$\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,I)}[p(X)]=0$. 2. (ii)

More generally, for any positive definite matrix $\Sigma$, we have $\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[p(X)]=\langle M,\Sigma-I\rangle$. 3. (iii)

$\mathop{\mbox{\bf Var}}_{X\sim\mathcal{N}(0,I)}[p(X)]=\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,I)}[p^{2}(X)]=\langle p,p\rangle=1$. 4. (iv)

More generally, for any positive definite matrix $\Sigma$, we have

[TABLE]

Proof.

The first three properties are a straightforward calculation. We show the last one here. By definition, we have

[TABLE]

as claimed, where (a) follows from Theorem 13. ∎

Observe that Lemma 10(iv) implies that if we take the top eigenvector of the $d^{2}\times d^{2}$ matrix

[TABLE]

on the linear subspace

[TABLE]

then the associated polynomial maximizes $\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[p^{2}(X)]$, and so we can find these polynomials efficiently. More generally, if we take any linear subspace of degree two polynomials with associated matrix subspace $V^{\prime}$, so that $V^{\prime}\subseteq V$, then the top eigenvector of the same matrix restricted to $V^{\prime}$ allows us to find the polynomial in this subspace which maximizes $\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[p^{2}(X)]$ efficiently.

We have the following tail bound for degree-$2$ polynomials in Gaussian space: We will use $\Pi_{V}(x)$ and $\Pi_{V}(S)$ to denote projection to a subspace $V$, of a point $x$ and a set of points $S$, respectively. We will also need the following hypercontractivity theorem for low-degree polynomials in Gaussian space:

Theorem 14** (Hypercontractivity in Gaussian space, see e.g. [O’D14]).**

Let $p:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a degree $m$ polynomial, and let $q\geq 2$ be even. Then $\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,I)}[p(X)^{q}]^{1/q}\leq(\sqrt{q-1})^{m}\|p\|_{2}$.

We need the following definition:

Definition 8**.**

Let $\mathcal{P}_{k}$ denote the set of even degree-$k$ polynomials over $d$ variables satisfying $\mathop{\mbox{\bf Var}}_{X\sim\mathcal{N}(0,1)}[p(X)]=1$. Moreover, for any subspace $W\subseteq\mathbb{R}^{d}$, let $\mathcal{P}_{k}(W)$ denote the set of even polynomials over $d$ variables which only depend on the coordinates in $W$.

Then by the arguments above, we have that for any two matrices $\Sigma,\widehat{\Sigma}$,

[TABLE]

In particular, by Lemma 2, this implies that when $\|\Sigma-I\|_{2}$ is small, then learning a Gaussian with unknown covariance in total variation distance is equivalent to learning the expectation of every even degree-$2$ polynomial.

Theorem 14 implies the following concentration for degree-$4$ (more generally, low-degree) polynomials of Gaussians:

Corollary 6**.**

Let $p$ be a degree-$4$ polynomial. Then there is some $A,C\geq 0$ so that for all $t\geq C$, we have

[TABLE]

Proof.

Hypercontractivity in particular implies the following moment bound: for all $q\geq 2$ even, we have

[TABLE]

By a typical moment argument, and optimizing the choice of $q$, this gives the desired bound. ∎

Hermite polynomials

Hermite polynomials are what arise by Gram-Schmidt orthogonalization applied with respect to this inner product. For a vector of non-negative integers $a=(a_{1},\ldots,a_{d})$, we let $H_{a}(x):\mathbb{R}^{d}\rightarrow\mathbb{R}$ be the Hermite polynomial associated with multi-index $a$. It is well-known that the degree of $H_{a}$ is $|a|=\sum_{i=1}^{d}a_{i}$, and moreover, $\langle H_{a},H_{b}\rangle=\delta_{a,b}$. In particular, for any $r\geq 1$, the Hermite polynomials of degree at most $r$ form an orthonormal basis with respect to the Gaussian inner product for all polynomials with degree at most $r$.

Therefore, given any polynomial $p:\mathbb{R}^{d}\rightarrow\mathbb{R}$ with degree $r$, we may write it uniquely as

[TABLE]

We define the $k$th harmonic component of $p$ to be

[TABLE]

and we say $p$ is harmonic of degree $k$ if it equals its $k$th part.

6.3 Working with Many Large Eigenvalues of the Second and Fourth Moment

As in the unknown mean case, we will need a filter to detect if there are many directions of the empirical covariance which have too large an eigenvalue. Formally, we need:

Theorem 15**.**

Fix $\varepsilon,\delta>0$. Assume $\|\Sigma-I\|_{F}\leq\xi$, where $\xi=O(\varepsilon\log 1/\varepsilon)$. Suppose that $G_{0}$ is $(\varepsilon,\delta)$-good with respect to $\mathcal{N}(0,\Sigma)$. Let $S$ be a set so that $\Delta(S,G_{0})\leq\varepsilon$. Let $\widehat{\Sigma}=\mathop{\mbox{\bf E}}_{S}[XX^{T}]$. Then there is an algorithm FilterCovManyDeg2Eig and a universal constant $C$ such that the following guarantee holds:

If $\widehat{\Sigma}-I$ has more than $O(\log 1/\varepsilon)$ eigenvalues larger than $C\xi$, then the algorithm outputs a $S^{\prime}$ so that $\Delta(S^{\prime},G_{0})<\Delta(S,G_{0})$. 2. 2.

Otherwise, the algorithm outputs “OK”, and outputs an orthonormal basis $v_{1},\ldots,v_{k}$ for the subspace $V$ of vectors spanned by all eigenvectors of $\widehat{\Sigma}-I$ with eigenvalue larger than $C\xi$.

The filter developed here is almost identical to the one developed for unknown mean. Thus, for conciseness we describe and prove the theorem in Appendix D.1.

We will also need a subroutine to enforce the condition that not only does the fourth moment tensor have spectral norm which is at most $O(\varepsilon\log^{2}1/\varepsilon)$ (restricted to a certain subspace of polynomials), but there can only be at most $O(\operatorname*{poly}\log 1/\varepsilon)$ directions in which the eigenvalue is large. However, the techniques here are a bit more complicated, for a number of reasons. Intuitively, the main complication comes from the fact that we do not know what the fourth moment tensor looks like, whereas in the unknown mean case, we knew that the covariance was the identity by assumption. Our main result in this subsection is the following subroutine:

Theorem 16** (Filtering when there are many large eigenvalues).**

Fix $\varepsilon,\delta>0$. Assume $\|\Sigma-I\|_{F}\leq\xi$, where $\xi=O(\varepsilon\log 1/\varepsilon)$. Let $C$ be the universal constant in FilterCovManyDeg2Eig. Let $W\subseteq\mathbb{R}^{d}$ be a subspace, so that for all $v\in W$ with $\|v\|_{2}=1$, we have $v^{T}\mathop{\mbox{\bf E}}_{S}[XX^{T}]v\leq 1+C\xi$. Suppose that $G_{0}$ is $(\varepsilon,\delta)$-good with respect to $\mathcal{N}(0,\Sigma)$. Let $S$ be a set so that $\Delta(S,G_{0})\leq\varepsilon$. Let $k=O(\log^{4}1/\varepsilon)$. Then there is an algorithm FilterCovManyDeg4Eig and universal constants $C_{1},C_{2}$ such that the following guarantee holds:

If there exist $p_{1},\ldots,p_{k}\in\mathcal{P}_{2}(W)$ so that $\langle p_{j},p_{\ell}\rangle=\delta_{j\ell}$ for all $j,\ell$, and so that $\mathop{\mbox{\bf E}}_{S}[p_{j}^{2}(Y)]-1\geq C_{1}\varepsilon$ for all $j$, then the algorithm outputs an $S^{\prime}$ so that $\Delta(S^{\prime},G_{0})<\Delta(S,G_{0})$. 2. 2.

Otherwise, the algorithm outputs “OK”, and outputs an orthonormal basis $p_{1},\ldots,p_{k^{\prime}}$ for a subspace $V$ of degree-$2$* polynomials in $\mathcal{P}_{2}(W)$ with $k^{\prime}\leq k$ so that for all $p\in V^{\perp}\cap\mathcal{P}_{2}$, we have $\mathop{\mbox{\bf E}}_{S}[p^{2}(X)]-1\leq C_{2}\varepsilon$.*

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

Roughly, we will show that if there are many polynomials with large empirical variance, this implies that there is a degree-four polynomial whose value is much larger than it could be if $w$ were the set of uniform weights over the uncorrupted points. Moreover, we can explicitly construct this polynomial, and it has a certain low-rank structure which allows us to use the concentration bounds we have previously derived.

6.4 Stitching Together Two Subspaces

This section is dedicated to giving an algorithm which allows us to fully reconstruct the covariance matrix given that we know it up to small error on a low-dimensional subspace $V$ and on $W=V^{\perp}$.

Theorem 17**.**

Let $1>\xi>\eta>\varepsilon>0$, and let $\tau>0$. Let $\Sigma$ so that $\|\Sigma-I\|_{F}\leq\xi$. Suppose that $\mathbb{R}^{d}$ is written as $V\oplus W$ for orthogonal subspaces $V$ and $W$ with $\dim(V)=O(\log(1/\varepsilon))$. Suppose furthermore that

[TABLE]

with $\|\Sigma_{V}-I_{V}\|_{F},\|\Sigma_{W}-I_{W}\|_{F}=O(\eta)$. Let $S_{0}=(G_{0},E_{0})$ be an $\varepsilon$-corrupted set of samples from $\mathcal{N}(0,\Sigma)$, and let $S\subseteq S_{0}$ with $\Delta(S,G)\leq O(\varepsilon)$ of size $\operatorname*{poly}(d,1/\eta,\log 1/\delta)$.

Then, there exists a universal constant $C_{5}$ and an algorithm Stitching that given $V,W,\xi,\eta,\varepsilon,\tau$ and $S$ runs in polynomial time and with probability at least $1-\tau$ returns a matrix $\Sigma_{0}$ with $\|\Sigma_{0}-\Sigma\|_{F}=C_{5}\eta+O(\xi^{2})$.

In the latter, we will show the algorithm works when $\tau=2/3$. As usual the probability of success can be boosted by repeating it independently.555Observe the only randomness at this point is in the random choices made by the algorithm. Thus, one can just run this algorithm $O(\log 1/\delta)$ times to obtain $\Sigma_{0}^{(1)},\ldots,\Sigma_{0}^{(\ell)}$ and find any $\Sigma_{0}^{(j)}$ which is $O(\eta+\xi^{2})$ close to at least a $2/3+o(1)$ fraction of the other outputs. The basic idea of the proof is as follows. Since we already know good approximations to $\Sigma_{V}$ and $\Sigma_{W}$, it suffices to find an approximation to $A$. In order to do this, we note that if we take a sample $x$ from $G$ conditioned on its projection to $V$ being some vector $v$, we find that the distribution over $W$ is a Gaussian with mean approximately $Av$. Running our algorithm for approximating the mean of a noisy-Gaussian, we can then compute the mapping $v\rightarrow Av$, which will allow us to compute $A$.

There are three main technical obstacles to this approach. The first is that we cannot condition on $x_{V}$ taking a particular value, as we will likely see no samples from $X$ with exactly that projection. Instead, what we will do is given samples from $X$ we will reject them with probabilities depending on their projections to $V$ in such a way to approximate the conditioning we require. The second obstacle is that the errors in $X$ may well be concentrated around some particular projection to $V$. Therefore, some of these conditional distributions may have a much larger percentage of errors than $\varepsilon$. To circumvent this, we will show that by carefully choosing how we do our conditioning and by carefully picking the correct distribution over vectors $v$, that on average these errors are only $O(\varepsilon)$. Finally, we need to be able to reconstruct $A$ from a collection of noisy approximations to $Av$. We show that this can be done by computing these approximations at a suitably large random set of $v$’s, and finding the matrix $A$ that minimizes the average $\ell_{2}$ error between $Av$ and its approximation.

Our algorithm is given in Algorithm 9:

Before we prove Theorem 17, we will need the following definition.

Definition 9**.**

A function $f:\mathbb{R}^{d}\rightarrow R$ is a positive measure if $f\geq 0$ and $\int f\leq 1$. We write $|f|_{1}=\int f$, and we say $f_{1}\leq f_{2}$ if $f_{1}(x)\leq f_{2}(x)$ pointwise.

Proof of Theorem 17.

Throughout this proof, let $G=\mathcal{N}(0,\Sigma)$. It is clear that this algorithm has polynomial runtime and sample complexity. We have yet to show correctness. The first thing that we need to understand is the procedure of rejection sampling, where we reject $x$ except with probability $\exp(-\|x_{V}-v\|^{2}/2)$. Therefore, given a distribution $D$, we let the positive measure $D_{v}$ be what is obtained by sampling from $D$ and accepting a sample $x$ only with probability $\exp(-\|x_{V}-v\|^{2}/2)$. We need to understand the distributions $G_{v}$ and $X_{v}$.

Note that the pdf of $G$ is given by

[TABLE]

Therefore, the density of $G_{v}$ is

[TABLE]

Note that letting $\mu_{v}:=(\Sigma^{-1}+I_{V})^{-1}v=v/2+O(\delta\|v\|_{2})$ that this equals

[TABLE]

Note that this is a Gaussian with mean $\mu_{v}$ weighted by

[TABLE]

so long as $\|v\|_{2}\ll\delta^{-1/2}$. Therefore, if this condition holds, a random sample from $G$ is accepted by this procedure with probability $\Theta(2^{-\dim(V)}\exp(-\|v\|_{2}^{2}/4)$.

We also need to understand the fraction of samples from $X_{v}$ that are erroneous. In a slight abuse of notation, let $E$ also denote the distribution which is uniform over the points in $E$, and let $L$ be the distribution which is uniform over the points in $L$. Therefore, we define

[TABLE]

that is approximately the fraction of samples from $X_{v}$ that are errors (where subtractive errors are weighted more heavily).

Note that so long as $\|v\|_{2}<\sqrt{C\log(1/\varepsilon)}$ that so long as $|L_{v}|_{1}=o(1)$ that $T$ in Step 5 will have sufficiently many samples with high probability. Furthermore, if this is the case, we have (assuming our mean estimation algorithm succeeds) that $\|\tilde{\mu}-\pi_{W}(\mu_{v})\|_{2}\leq O(\varepsilon_{v}+\eta)$. Therefore, unless $\|\mu_{v}\|_{2}=\Omega(C\log(1/\varepsilon))$ (which can only happen in $\|v\|_{2}=\Omega(C\log(1/\varepsilon))$), we have with high probability that $\|a_{v}-\pi_{W}(\mu_{v})\|_{2}=O(\varepsilon_{v}+\eta)$. In fact, we assume that this holds for all $v\in S$ with $\|v\|_{2}\ll(C\log(1/\varepsilon))$.

In order to show that this is generally a good approximation, we need to know that $\varepsilon_{v}$ is not too large on average. In particular, we show:

Lemma 11**.**

If $v\sim N(0,2I_{V})$, then $\mathop{\mbox{\bf E}}_{v}[\varepsilon_{v}]=O(\varepsilon)$.

Proof.

Letting $F=E+\log(1/\varepsilon)L$, we have that

[TABLE]

∎

Note that $\pi_{W}(\mu_{v})=Mv$, where $M=\pi_{W}(\Sigma^{-1}+I_{V})^{-1}\pi_{V}^{T}$. Let $B$ be any $\dim(W)\times\dim(V)$ matrix. Note that

[TABLE]

Therefore, $\mathop{\mbox{\bf E}}_{v\sim N(0,2I_{V})}[\|Bv-Mv\|_{2}]=O(\|B-M\|_{F})$. Since $\|Bv-Mv\|_{2}^{2}$ is a degree-$2$ polynomial in $v$, by Corollary 4 in [Kan12],

[TABLE]

Combining these we find that

[TABLE]

Next we show that our choice of $S$ derandomizes this result.

Lemma 12**.**

With high probability over the choice of $S$, we have that

[TABLE]

for all $\dim(W)\times\dim(V)$ matrices $B$.

Proof.

This is equivalent to showing that for all matrices $U=B-M$ it holds

[TABLE]

By the standard scaling laws, it suffices to show this only for $U$ with $\|U\|_{F}=1$.

We also note that it suffices to show this only for $U$ in an $\varepsilon$-net for all such matrices. This is because if $\|U-U^{\prime}\|_{F}<\varepsilon$ and if $U^{\prime}$ satisfies the desired condition, then

[TABLE]

and with high probability $\mathop{\mbox{\bf E}}_{v\in_{u}S}[\|v\|_{2}]=O(\log(1/\varepsilon)).$

Note that such nets exist with size $\exp(\operatorname*{poly}(n/\varepsilon))$. Therefore, it suffices to show that this condition holds for each such $U$ with probability $\exp(-(n/\varepsilon)^{\Omega(C)})$.

As noted above, it suffices that

[TABLE]

and that

[TABLE]

The first follows because if $S=\{v_{1},\ldots,v_{m}\}$ then $\mathop{\mbox{\bf E}}_{v\in_{u}S}[\|Uv\|^{2}]$ is a degree-$2$ polynomial in the $v_{i}$ with mean $O(1)$ and variance $O(1/\sqrt{m})$, so by standard concentration results, is $O(1)$ with $1-\exp(-\Omega(|S|))$ probability. The latter follows from standard concentration bounds. This completes the proof. ∎

We also note that the random choice of $S$ has another nice property:

Lemma 13**.**

With $S$ and $a_{v}$ as above, with high probability we have that

[TABLE]

Proof.

Let $S=\{v_{1},\ldots,v_{m}\}$. Then $\|f(v_{i})-Mv_{i}\|_{2}$ are independent random variables with mean $\mathop{\mbox{\bf E}}[O(\varepsilon_{v}+\eta)]=O(\eta)$ and variance at most

[TABLE]

The result follows by Chernoff bounds. ∎

Now by Lemmas 12 and 13 with high probability over the choice of $S$ in Step 4, and the $a_{v}$ in Step 5 we have that for all $\dim(W)\times\dim(V)$-matrices $B$ that

[TABLE]

and have that

[TABLE]

Combining these statements, we have that for all $\dim(W)\times\dim(V)$-matrices $B$ that

[TABLE]

Note that by taking $B=M$, this quantity is $O(\eta)$. Therefore, the $B$ found in Step 15 satisfies

[TABLE]

and therefore, $\|B-M\|_{F}=O(\eta)$.

The rest of the proof is a simple computation of the matrices involved. In particular, recall that

[TABLE]

where the $O(\eta)$ denotes a matrix with Frobenius norm $O(\eta)$ and where $\|A\|_{F}=O(\delta)$. It is easy to see that

[TABLE]

Therefore,

[TABLE]

Hence,

[TABLE]

Therefore, $M=A/2+O(\eta+\delta^{2})$. Therefore, $A=2M+O(\eta+\delta^{2})=2B+O(\eta+\delta^{2})$. And finally, we conclude that

[TABLE]

∎

6.5 The Full High-Dimensional Algorithm

We now show how to prove Theorem 11, given the pieces we have. We first show that given enough samples from $\mathcal{N}(0,\Sigma)$, the empirical data set without corruptions satisfies the regularity conditions in Section 2.2.2 with high probability. For clarity of exposition, the proof of this lemma is deferred to Appendix D.3.

Lemma 14**.**

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

Finally, we require the following guarantee, which states that if there is a degree-$2$ polynomial whose expectation under $S$ and the truth differs by a lot (equivalently, if the empirical covariance differs from the true covariance in Frobenius norm substantially), then it must also have very large variance under $S$.

Lemma 15**.**

Fix $\varepsilon,\delta>0$. Assume $\|\Sigma-I\|_{F}\leq\xi$, where $\xi=O(\varepsilon\log 1/\varepsilon)$. Suppose that $G_{0}$ is $(\varepsilon,\delta)$-good with respect to $\mathcal{N}(0,\Sigma)$, and let $S\subseteq S_{0}$ be a set so that $\Delta(S,G_{0})\leq\varepsilon$. There is some absolute constant $C_{5}$ so that if $p\in\mathcal{P}_{2}$ is a polynomial so that $\left|\mathop{\mbox{\bf E}}_{S}[p(X)]-\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma)}[p(X)]\right|>C_{5}\sqrt{\xi\varepsilon}$, then $\mathop{\mbox{\bf E}}_{S}[p^{2}(X)]-1>C_{1}\xi$.

We defer the proof of this lemma to the Appendix.

We are now ready to present the full algorithm as Algorithm 10.

Proof of Theorem 11.

Condition on the events that neither LearnCovLowDim nor Stitching fail. This happens with probability at least $\operatorname*{poly}(\varepsilon,1/d,\delta)$. Observe that if we pass the “if” statement in Line 5, then by the guarantee of FilterCovManyDeg2Eig this is indeed an $S^{\prime}$ satisfying the desired properties. Otherwise, by the guarantees of FilterCovManyDeg2Eig, we have that $W$ satisfies the conditions needed by FilterCovManyDeg4Eig. Hence, if we pass the “if” statement in Line 11, then the guarantee of FilterCovManyDeg4Eig this is indeed a $S^{\prime}$ satisfying the desired properties. Otherwise, by Lemma 15, we know that for all polynomials $p\in\mathcal{P}_{2}$ over $W$ orthogonal to $U_{1}$, we have $|\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma)}-\mathop{\mbox{\bf E}}_{S}[XX^{T}]|\leq C_{5}\sqrt{\xi\varepsilon}$. Thus, $\Sigma_{W}$ satisfies the conditions needed by Stitching.

By Corollary 4, we know that $\Sigma_{V}$ satisfies the conditions for Stitching, and so the correctness of the algorithm follows from Theorem 17. ∎

6.6 The Barrier at Quasi-Polynomial

Here we explain why improving the running time from quasi-polynomial to polynomial in $1/\varepsilon$ will likely be rather difficult. Recall that our strategy is to project the problem onto lower dimensional subproblems and stitch together the answer. We need the dimension of the subspace to be large enough that we can find a polynomial $Q$ that is itself the sum of squares of $k$ orthogonal degree two polynomials $p_{i}$ so that the value of $Q$ on the corrupted points is considerably larger than the value on the uncorrupted points. More precisely, if we let $S=(G,E)$ denote our corrupted set of samples then we want $\mathop{\mbox{\bf E}}_{E}[Q(X)]$ to be larger than $Q(X)$ for all but a $\operatorname*{poly}(\varepsilon)$ fraction of $X\in G$. We then remove all points $X\in S$ with large $Q(X)$ and by the properties of $Q$ we are guaranteed that we throw out mostly corrupted points. It turns out that the most aggressive we could be is removing points where $Q(X)$ is more than $\sqrt{k}$ standard deviations away from its expectation under the true Gaussian. But since $Q$ is a degree-four polynomial and we want $Q(X)$ to be smaller than our cutoff for all but a $\operatorname*{poly}(\varepsilon)$ fraction of $X\in G$, we are forced to choose $\sqrt{k}=\Omega(\log 1/\varepsilon)$, which means that we need to reduce to $k=\Omega(\log^{2}1/\varepsilon)$ dimensional subproblems. Thus, if we solve low-dimensional subproblems in time exponential in the dimension, we naturally arrive at a quasi-polynomial running time. It seems that any approach for reducing the running time to polynomial would require fundamentally new ideas.

7 The General Algorithm

We now have all the tools to robustly learn the mean and covariance of an arbitrary high-dimensional Gaussian. We first show how to reduce the problem of robustly learning the covariance of $\mathcal{N}(\mu,\Sigma)$ to learning the covariance of $\mathcal{N}(0,\Sigma)$, by at most doubling error, a trick previously used in [DKK*+*16] and [LRV16]. Given an $\varepsilon$-corrupted set of samples $X_{1},\ldots,X_{2n}$ of size $2n$ from $\mathcal{N}(\mu,\Sigma)$, we may let $Y_{i}=(X_{i}-X_{n+i})/\sqrt{2}$. Then we see that if $X_{i}$ and $X_{n+i}$ are uncorrupted, then $Y_{i}\sim\mathcal{N}(0,\Sigma)$. Moreover, at most $2\varepsilon n$ of the $Y_{i}$ can be corrupted, since there are at most $2\varepsilon n$ corrupted $X_{i}$. Therefore, by doubling the error rate, we may assume that $\mu=0$. We may then apply the algorithm in Corollary 5 to obtain a $\widehat{\Sigma}$ so that with high probability, we have $\|\widehat{\Sigma}^{-1/2}\Sigma\widehat{\Sigma}^{-1/2}-I\|_{F}\leq O(\varepsilon)$ with polynomially many samples, and in $\operatorname*{poly}(d,(1/\varepsilon)^{O(\log^{4}1/\varepsilon)})$ time.

We may then take an additional set of $\varepsilon$-corrupted samples $\{X_{i}^{\prime},\ldots,X_{n}^{\prime}\}$, and let $Y_{i}^{\prime}=\widehat{\Sigma}^{-1/2}X_{i}^{\prime}$. Then, by our guarantee on $\widehat{\Sigma}$, we have that if $X_{i}^{\prime}$ is uncorrupted, then $Y_{i}^{\prime}\sim\mathcal{N}(0,\widetilde{\Sigma})$ where $\|\widetilde{\Sigma}-I\|_{F}\leq O(\varepsilon)$. We then run RecoverMeanNoisy with the $Y_{i}^{\prime}$ to obtain a $\widehat{\mu}$ so that $\|\widehat{\mu}-\widehat{\Sigma}^{-1/2}\mu\|_{2}\leq O(\varepsilon)$. This guarantees that $d_{\mathrm{TV}}(\mathcal{N}(\widehat{\mu},\widetilde{\Sigma}),\mathcal{N}(\widehat{\Sigma}^{-1/2}\mu,\widehat{\Sigma}))\leq O(\varepsilon)$, which in turn implies that $d_{\mathrm{TV}}(\mathcal{N}(\widehat{\mu},\widehat{\Sigma}),\mathcal{N}(\mu,\Sigma))\leq O(\varepsilon)$, as claimed.

Therefore, we have shown:

Theorem 18**.**

Fix $\varepsilon,\delta>0$. Given an $\varepsilon$-corrupted set of samples $S$ from $\mathcal{N}(\mu,\Sigma)$, where $n=\operatorname*{poly}(d,1/\varepsilon,\log 1/\delta)$, there is an algorithm RecoverGaussian which takes as input $S,\varepsilon,\delta$, and outputs a $\widehat{\mu},\widehat{\Sigma}$ so that

[TABLE]

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

Appendix A Lower Bounds on Agnostic Learning

In this section, we prove information theoretic lower bounds for robust distribution learning. In particular, we will prove that no algorithm (efficient or inefficient) can learn a general distribution from an $\varepsilon$-corrupted set of samples to total variation distance less than $\frac{\varepsilon}{1-\varepsilon}$. We note that our lower bound applies in more general settings than we consider in this paper in couple of ways:

Our construction is univariate; 2. 2.

Our construction works for any pair of distributions which are $\frac{\varepsilon}{1-\varepsilon}$-close, and not just for Gaussian distributions; and 3. 3.

Our construction holds for Huber’s $\varepsilon$-contamination model, which is weaker than the noise model studied in this paper.

In Huber’s $\varepsilon$-contamination model, data is drawn from a mixture distribution $(1-\varepsilon)P+\varepsilon Q$, where $P$ is some distribution that we wish to estimate and $Q$ is arbitrary (in particular, it might depend on $P$). We will show that any two distributions which are $\frac{\varepsilon}{1-\varepsilon}$-close can be made indistinguishable under this contamination model.

Lemma 16**.**

Let $p_{1}$ and $p_{2}$ be two distributions such that $d_{\mathrm{TV}}(p_{1},p_{2})=\frac{\varepsilon}{1-\varepsilon}$. Then there exist Huber $\varepsilon$-contaminations $p_{1^{\prime}}$ and $p_{2^{\prime}}$ of $p_{1}$ and $p_{2}$ respectively, such that $p_{1^{\prime}}=p_{2^{\prime}}$. Therefore, no algorithm can learn a distribution $p$ up to accuracy $<\frac{\varepsilon}{1-\varepsilon}$ in Huber’s $\varepsilon$-contamination model.

Proof.

Since $d_{\mathrm{TV}}(p_{1},p_{2})=\frac{\varepsilon}{1-\varepsilon}$, we can write

[TABLE]

where $p_{c},q_{1},q_{2}$ are distributions. Let $p_{1^{\prime}}$ be the Huber $\varepsilon$-contamination of $p_{1}$ with $q_{2}$, $p_{2^{\prime}}$ be the Huber $\varepsilon$-contamination of $p_{2}$ with $q_{1}$. In other words,

[TABLE]

and thus $p_{1^{\prime}}=p_{2^{\prime}}$ as desired. ∎

We note that a similar lower bound holds when $p_{1}$ and $p_{2}$ are required to be Gaussians, but weakened by a factor of $2$. This is due to the geometry of the space of Gaussian distributions, and we note that in the case where the variance is known, this is achieved by the median.

Lemma 17**.**

No algorithm can output an estimate for the mean of a unit-variance Gaussian at accuracy $<\left(\sqrt{\frac{\pi}{2}}-o(1)\right)\varepsilon$ with probability $>1/2$ in Huber’s $\varepsilon$-contamination model. Consequently, by Lemma 1, no algorithm can learn a Gaussian to total variation distance $<\left(\frac{1}{2}-o(1)\right)\varepsilon$ in Huber’s $\varepsilon$-contamination model.

Proof.

We will consider the distributions $p_{1}=\mathcal{N}(-\alpha,1)$ and $p_{2}=\mathcal{N}(\alpha,1)$, where $\alpha$ is to be specified later. We note that if $p_{1}$ and $p_{2}$ can be $\varepsilon$-corrupted into the same distribution, then the best estimate for the mean is to output [math] (by symmetry), and the result holds.

We will show that $p_{1}$ can be $\varepsilon$-corrupted into a distribution $f$, where $f(x)=\frac{\max\left\{p_{1}(x),p_{2}(x)\right\}}{\eta}$ and $\eta$ is a normalizing constant (the case of $p_{2}$ follows similarly). In other words, $f$ can be written as $(1-\varepsilon)p_{1}+\varepsilon q_{1}$ for some distribution $q_{1}$. Since $q_{1}$ is a distribution, it is non-negative, and thus we require that $(1-\varepsilon)p_{1}(x)\leq f(x)$ for all $x\in\mathbb{R}$. Note that $f(x)=\frac{\max\left\{p_{1}(x),p_{2}(x)\right\}}{\eta}\geq\frac{p_{1}(x)}{\eta}$.

We can compute $\eta$ as follows:

[TABLE]

Thus,

[TABLE]

Again, since we require that $(1-\varepsilon)p_{1}(x)\leq f(x)$, it suffices that $1-\varepsilon\leq 1-\sqrt{\frac{2}{\pi}}\alpha+O(\alpha^{2})$ and thus that $\sqrt{\frac{2}{\pi}}\alpha-O(\alpha^{2})\leq\varepsilon$. This corresponds to $\alpha\leq\left(\sqrt{\frac{\pi}{2}}-o(1)\right)\varepsilon$, and the proof is complete. ∎

Observe that this lower bound can be strengthened by a factor of two in the subtractive adversary model, where the adversary is able to both remove $\varepsilon N$ samples and add $\varepsilon N$ samples. In this infinite sample regime, this means that a distribution can be corrupted to any distribution which is $\varepsilon$-far in total variation distance. By noting that $\mathcal{N}(-\alpha,1)$ and $\mathcal{N}(\alpha,1)$ are both $\varepsilon$-far from $\mathcal{N}(0,1)$ for $\alpha=\left(\sqrt{2\pi}+o(1)\right)\varepsilon$, we obtain the following lemma. Note that this lower bound is also achieved by the median in this model.

Lemma 18**.**

No algorithm can output an estimate for the mean of a unit-variance Gaussian at accuracy $<\left(\sqrt{2\pi}+o(1)\right)\varepsilon$ with probability $>1/2$ in the subtractive adversary model. Consequently, by Lemma 1, no algorithm can learn a Gaussian to total variation distance $<\varepsilon$ in the subtractive adversary model.

Finally, we conclude by sketching a lower bound for mean estimation of sub-Gaussian distributions.

Lemma 19**.**

No algorithm can output an estimate for the mean of a sub-Gaussian distribution at accuracy $o(\varepsilon\log^{1/2}(1/\varepsilon))$ with probability $>1/2$ in the Huber’s $\varepsilon$-contamination model.

Proof.

We start with the distribution $q=\mathcal{N}(0,1)$. We construct $p_{1}$ by truncating the right tail of $q$ at the point $x_{r}=c_{1}\log^{1/2}(1/\varepsilon)$ for some constant $c_{1}$, and rescaling the rest of the distribution appropriately. Observe that, for an appropriate choice of $c_{1}$:

• •

$p_{1}$ is sub-Gaussian with a constant of $O\left(\frac{1}{1-\varepsilon}\right)$;

• •

The mean of $p_{1}$ is $-c_{2}\varepsilon\log^{1/2}(1/\varepsilon)$ for some constant $c_{2}$;

• •

$p_{1}$ can be corrupted in Huber’s $\varepsilon$-contamination model to be $q$.

We can similarly consider $p_{2}$, which is constructed by truncating the left tail of $q$ at the point $x_{l}=-c_{1}\log^{1/2}(1/\varepsilon)$. Since $p_{1}$ and $p_{2}$ are indistinguishable when they are both $\varepsilon$-corrupted to $q$, and the mean of all three distributions are separated by $\geq c_{2}\varepsilon\log^{1/2}(1/\varepsilon)$, the lemma follows. ∎

Appendix B Omitted Proofs from Section 2

B.1 Proof of Lemma 1

Proof.

Observe that by rotational and translational invariance, it suffices to consider the problem when $\mu_{1}=-\varepsilon e_{1}/2$ and $\mu_{2}=\varepsilon e_{1}/2$, where $e_{1}$ is the first standard basis vector. By the decomposability of TV distance, we have that the TV distance can in fact be written as a 1 dimensional integral:

[TABLE]

The value of the function $f(x)=e^{-(x-\varepsilon/2)^{2}/2}-e^{-(x+\varepsilon/2)^{2}/2}$ is negative when $x<0$ and positive when $x>0$, hence this integral becomes

[TABLE]

where $F(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-t^{2}/2}dt$ is the CDF of the standard normal Gaussian. By Taylor’s theorem, and since $F^{\prime\prime}(x)$ is bounded when $x\in[-1,1]$, we have

[TABLE]

which proves the claim. ∎

Appendix C Omitted Proofs from Section 4

C.1 Proof of Theorem 7

First, we require the following several claims about $p(x)$ under various distributions, including the true Gaussian and from choosing points uniformly at random from the sets of interest. Note that throughout this section we will use Lemma 3, which implies that $\|\Pi_{V^{\prime}}\mu-\tilde{\mu}\|_{2}<2\varepsilon$ when the parameters $\alpha,\gamma$ are chosen appropriately.

Claim 4**.**

$\mathop{\mbox{\bf E}}_{\mathcal{N}(\mu,I)}[p(X)]=\|\Pi_{V^{\prime}}\mu-\tilde{\mu}\|_{2}^{2}$.

Proof.

Letting $v_{1},\dots,v_{C_{1}\beta\log(1/\varepsilon)}$ be an orthonormal basis of $V^{\prime}$, we have

[TABLE]

∎

Claim 5**.**

For some absolute constant $c_{0}$ and all $t\geq\|\Pi_{V^{\prime}}\mu-\tilde{\mu}\|_{2}^{2}$,

[TABLE]

Proof.

This follows from Lemma 7, after re-centering the polynomial using Claim 4 and noting that the spectral norm and squared Frobenius norm of the corresponding $A$ matrix are at most $1$ and $C_{1}\beta\log(1/\varepsilon)$, respectively. ∎

Claim 6**.**

$\mathop{\mbox{\bf E}}_{S}[p(X)]\geq C_{1}\varepsilon\log(1/\varepsilon)$.

Proof.

Recall that $\hat{\mu}$ is the empirical mean of the point set.

[TABLE]

∎

Claim 7**.**

$\mathop{\mbox{\bf E}}_{G_{0}}[p(X)]\leq\|\mu-\tilde{\mu}\|_{2}^{2}+O(\gamma)\varepsilon$.

Proof.

[TABLE]

The inequalities follow from $(\gamma\varepsilon,\delta)$-goodness and Claim 4. ∎

Claim 8**.**

$\phi(S,G_{0})\mathop{\mbox{\bf E}}_{L}[p(X)]=(C_{1}\beta+O(\gamma)+o(1))\varepsilon$.

Proof.

[TABLE]

(4) and (5) follow from $G_{0}$ being $(\gamma\varepsilon,\delta)$-good, (6) is from Claim 5, and (7) is because

[TABLE]

∎

This gives:

Claim 9**.**

$\psi\mathop{\mbox{\bf E}}_{E}[p(X)]\geq C_{1}\varepsilon\log 1/\varepsilon-(C_{1}\beta+1+O(\gamma)+o(1))\varepsilon$.

Proof.

This immediately follows from Claims 6, 7, and 8. ∎

We now show that in the case that there are many large eigenvalues, there will be a $T$ satisfying the conditions of the filter.

Claim 10**.**

Suppose $\dim(V)\geq C_{1}\beta\log(1/\varepsilon)$. Then there is a $T$ satisfying the conditions in the algorithm.

Proof.

Suppose not. Then, we have

[TABLE]

By applying the contradiction assumption, we have

[TABLE]

and thus, we have

[TABLE]

which contradicts Claim 9. ∎

It now suffices to prove that if we construct a filter, then the invariant that $\Delta$ decreases is preserved. Formally, we show:

Claim 11**.**

Suppose $\dim(V)\geq C_{1}\beta\log(1/\varepsilon)$. Let $S^{\prime}$ be the set of points we return. Then $\Delta(S^{\prime},G)<\Delta(S,G)$.

Proof.

Let $T$ be the threshold we pick. If $T>C_{2}d\log(|S|/\delta)$ then the invariant is satisfied since we remove no good points, by $(\gamma\varepsilon,\delta)$-goodness. It suffices to show that in the other case, we remove $\log(1/\varepsilon)$ times many more bad points than good points. By definition we remove at least

[TABLE]

points. On the other hand, by $(\gamma\varepsilon,\delta)$-goodness, Claim 5, we know that we throw away at most

[TABLE]

By our choice of $T$, and since $\|\mu-\widetilde{\mu}\|_{2}\leq O(\delta)\ll 1$, the first term is upper bounded by

[TABLE]

and so

[TABLE]

Since we have maintained this invariant so far, in particular, we have thrown away more bad points than good points, and so $|G_{0}|\leq(1+\varepsilon)|S|$.

Therefore, we have $\log(1/\varepsilon)\cdot|G|\cdot\Pr_{G}[p(X)>T]\leq|S|\cdot\left(\exp\left(-\frac{c_{0}T}{2C_{3}}\right)+\frac{\gamma\varepsilon}{d\log|S|/\delta}\right)$, and hence, the invariant is satisfied. ∎

Claims 10 and 11 together imply the correctness of Theorem 7.

C.2 Proof of Lemma 5

We require the following basic fact of good sets (see, e.g., Fact 8.6 in [DKK*+*16]):

Fact 2**.**

Let $v\in\mathbb{R}^{d}$ be any unit vector. Let $G$ be $(\eta,\delta)$-good with respect to $\mathcal{N}(\mu,I)$. Then for any $T>0$,

[TABLE]

and

[TABLE]

Proof.

Let $\langle v,\mu-\widehat{\mu}\rangle=R>\frac{\varepsilon}{\beta^{1/2}}$. Observe that by direct calculation, we have $\mathop{\mbox{\bf E}}_{S}[\langle v,X-\mu\rangle^{2}]-\mathop{\mbox{\bf E}}_{S}[\langle v,X-\widehat{\mu}\rangle^{2}]\leq 3R^{2}$. Hence, it suffices to show that $\mathop{\mbox{\bf E}}_{S}[\langle v,X-\mu\rangle^{2}]\geq 1+\Omega(R^{2}/\varepsilon)-\left(\gamma+O\left(\frac{1}{\beta}\right)\right)\varepsilon.$

We write $S=(G_{0}\setminus L,E)$ where $|E|=\psi|S|$ and $|L|=\phi|S|$. First we consider the expectation of $v\cdot X$ over $X$ in $S$. This is

[TABLE]

where $\mu_{0}$ is the mean over $G_{0}$, $\mu_{L}$ the mean over $L$ and $\mu_{E}$ the mean over $E$. By $(\gamma\varepsilon,\delta)$-goodness, we have that $\left|\langle v,\mu_{0}-\mu\rangle\right|\leq\gamma\varepsilon$. We also have that

[TABLE]

Since we have

[TABLE]

by Fact 2 and the calculations done in the proof of Corollary 8.8 in [DKK*+*16], we have that $|\langle v,\mu_{L}-\mu\rangle|\leq O(\log|S|/|L|)\leq O(\log 1/\phi)$. Since by assumption we have $\langle v,\mu-\widehat{\mu}\rangle=R>\frac{\varepsilon}{\beta^{1/2}}$, this implies that $\langle v,\mu_{E}-\mu\rangle=\Omega(R/\varepsilon)$. In particular, this implies that

[TABLE]

Next we consider the expectation of $\langle v,X-\mu\rangle^{2}$. By $(\gamma\varepsilon,\delta)$-goodness, we have $\mathop{\mbox{\bf E}}_{G_{0}}[\langle v,X-\mu\rangle^{2}]=1+O(\varepsilon)$. We also have that

[TABLE]

Since

[TABLE]

this is at least $1+\Omega(R^{2}/\varepsilon)-\left(\gamma+O\left(\frac{1}{\beta}\right)\right)\varepsilon$.

∎

Appendix D Omitted Proofs from Section 6

D.1 Proof of Theorem 15

Our algorithm works as follows, just as for Algorithm 4. 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 11.

The proofs of the following claims are identical to the proofs of Claim 6-9, by applying the corresponding property of $(\gamma\varepsilon,\delta)$-goodness for this setting, and so we omit them.

Claim 12**.**

$\mathop{\mbox{\bf E}}_{S}[p(x)]\geq C\xi.$**

Claim 13**.**

$\mathop{\mbox{\bf E}}_{G_{0}}[p(X)]\leq\frac{\varepsilon}{(d\log 1/\varepsilon)^{2}}$.

Claim 14**.**

$\phi(S,G_{0})\mathop{\mbox{\bf E}}_{L}[p(X)]=(C+O(\gamma)+o(1))\varepsilon$.

Claim 15**.**

$\psi\mathop{\mbox{\bf E}}_{E}[p(X)]\geq C\xi-(C_{1}+1+O(\gamma)+o(1))\varepsilon$.

These claims imply just as before that there is a $T$ with the desired properties. Again, the proof is identical. Formally:

Claim 16**.**

Suppose $\dim(V)\geq C\log(1/\varepsilon)$. Then there is a $T$ satisfying the conditions in the algorithm.

Finally, we show the invariant that $\Delta$ decreases. This is almost identical to the proof of Claim 11, however, we need to slightly change our application of the Hanson-Wright inequality. Formally, we show:

Claim 17**.**

Suppose $\dim(V)\geq C_{1}\log 1/\varepsilon$. Let $S^{\prime}$ be the set of points we return. Then $\Delta(S^{\prime},G)<\Delta(S,G)$.

Proof.

Let $T$ be the threshold we pick. If $T>C_{2}d\log(|S|/\delta)$ then the invariant is satisfied since we remove no good points, by $(\gamma\varepsilon,\delta)$-goodness. It suffices to show that in the other case, we remove $\log 1/\varepsilon$ times many more bad points than good points. By definition we remove at least

[TABLE]

points. On the other hand, observe that if $v_{1},\ldots,v_{k}$ is an orthonormal basis for $V$, we have

[TABLE]

so that if $X\sim\mathcal{N}(0,\Sigma)$, we have $Y\sim\mathcal{N}(0,I)$. Let $M=\Sigma^{1/2}\left(\sum v_{i}v_{i}^{T}\right)\Sigma^{1/2}$. We have that $\|M\|_{F}\leq\sum_{i=1}^{n}\|\Sigma\|_{2}\leq(1+\xi)k$, and $\|M\|_{2}\leq\|\Sigma\|_{2}\leq(1+\xi)$. Since

[TABLE]

we have by Hanson-Wright that

[TABLE]

so by our choice of $T$, we have

[TABLE]

The remaining proof now proceeds identically to the proof of Claim 17. ∎

D.2 Proof of Theorem 16

Clearly, the only non-trivial condition to certify for Theorem 16 is that if we are in Case (1), the returned set satisfies the desired properties.

Our proof will roughly follow the same structure as the proof of Theorem 7. We will first show that the empirical average of the polynomial $Q$ can only be large because of the contribution of the points in $E$ (Claim 21). We will then show that this implies that there exists a threshold $T$ which the algorithm will find in this case (Claim 22). Finally, we will show that for any such $T$ we find, the returned set of points will indeed satisfy $\Delta(S^{\prime},G_{0})<\Delta(S,G_{0})$, which implies the correctness of the algorithm (Claim 23).

We first show the following claim:

Claim 18**.**

$\mathop{\mbox{\bf E}}_{S}[r_{i}]-\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma)}[r_{i}]\leq(4C+1)\xi$**

Proof.

Let us suppose that $p_{i}$ corresponds to the matrix $A_{i}$, given by $(A_{i})_{a,b}=\nabla_{a}\nabla_{b}p_{i}/\sqrt{2!}$, so that the $A_{i}$ are orthonormal with respect to the Frobenius norm. Then the constant harmonic part of $p_{i}^{2}$ corresponds to $\|A_{i}\|_{F}^{2}\leq 1$. The degree-$2$ harmonic part of $p_{i}^{2}$ corresponds to the matrix $2\sqrt{2}A_{i}^{2}$. This is because if we let $B_{i}$ be the matrix corresponding to $p_{i}^{2}$, we have

[TABLE]

where the second to last line follows since $\nabla_{a}\nabla_{b}p_{i}(X)$ is a constant, and the last line follows from explicit computation. In particular, this implies that the non-constant component of $r_{i}$ corresponds to matrix with trace norm at most $2\sqrt{2}\leq 4$. Therefore, $r_{i}$ can be written as $r_{i}(x)=\sum\alpha_{i}\langle v_{i},x\rangle^{2}+C_{0}$ for some constant $C_{0}$, where $\sum|\alpha_{i}|\leq 4$. Thus, by our assumption, we have $\mathop{\mbox{\bf E}}_{S}[r_{i}]\leq 1+C\xi$. The claim then follows since $\Sigma$ and $I$ are differ in Frobenius norm by at most $\xi$. ∎

Claim 19**.**

$\mathop{\mbox{\bf E}}_{S}[Q(X)]-\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[Q(X)]\geq(C_{1}-6)\xi k$**

Proof.

Observe that since $\|\Sigma-I\|_{F}\leq\xi$, in particular we have $\|\Sigma^{\otimes 2}-I^{\otimes 2}\|_{2}\leq\xi$, and hence by Lemma 10 we have that $\left|\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma)}[p^{2}(X)]-\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,I)}[p^{2}(X)]\right|\leq 2\xi$ for all $p\in\mathcal{P}_{2}$. Hence, in particular, we have

[TABLE]

By Claim 18, we have

[TABLE]

In particular, this implies that

[TABLE]

as claimed. ∎

We now show:

Claim 20**.**

There is some universal constant $B$ so that

[TABLE]

Before we prove this, we need the following lemma:

Lemma 20**.**

For any degree-$4$ polynomial $p$, and any $\Sigma$, if we let $F(y,p,\Sigma)$ denote the $y$th percentile of $p$ under $\Sigma$, then we have $F(1/4,p^{2},\Sigma),F(3/4,p^{2},\Sigma)=\Theta(\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[p^{2}(X)])$.

Proof.

Let $\mu^{\prime}=\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[p^{2}(X)]$. First, we note that $\Pr(p^{2}(X)>4\mu^{\prime})\leq 1/4,$ so $F(3/4,p^{2},\Sigma)\leq 4\mu^{\prime}$. On the other hand, by known anti-concentration bounds [CW01], we have that $\Pr(p^{2}(X)\leq\varepsilon\mu^{\prime})=\Pr(|p(X)|\leq\sqrt{\varepsilon\mu^{\prime}})=O(\varepsilon^{1/8})$. So, for $\varepsilon$ a sufficiently small constant, $\Pr(p^{2}(X)\leq\varepsilon\mu^{\prime})<1/4$, and therefore, $F(1/4,p^{2},\Sigma)\geq\varepsilon\mu^{\prime}$. Since $\varepsilon\mu^{\prime}\leq F(1/4,p^{2},\Sigma)\leq F(3/4,p^{2},\Sigma)\leq 4\mu^{\prime}$, this completes our proof. ∎

Proof of Claim 20.

Since $d_{\mathrm{TV}}(\mathcal{N}(0,\Sigma),\mathcal{N}(0,I))=O(\varepsilon\log 1/\varepsilon)=o(1)$, we have that $F(3/4,Q^{2},\Sigma)\leq F(1/4,Q^{2},\Sigma)$, so by the above lemma, we have that

[TABLE]

Hence, it suffices to bound $\|Q\|_{2}^{2}$.

Let us again suppose that $p_{i}$ corresponds to the matrix $A_{i}$, given by $(A_{i})_{a,b}=\nabla_{a}\nabla_{b}p_{i}/\sqrt{2!}$. Note that the $A_{i}$ are symmetric matrices that form an orthonormal set. We note that $q_{i}$ is the harmonic degree-$4$ polynomial corresponding to the rank-$4$ tensor

[TABLE]

By linearity, $Q$ corresponds to the rank-$4$ tensor $T=\sqrt{6}\sum_{i=1}^{k}\mathrm{Sym}(A_{i}\otimes A_{i})$. It thus suffices to show that $\|T\|_{2}=O(k)$, where here $\|\cdot\|_{2}$ denotes the square root of the sum of the squares of the entries of the tensor. In order to show this, we note that $\|T\|_{2}=\sup_{\|V\|_{2}=1}\langle V,T\rangle.$ Therefore, it suffices to show that for all $4$-tensors $V$ with $\|V\|_{2}\leq 1$, that $\langle V,T\rangle=O(\sqrt{k})$. We note that

[TABLE]

Note that $\mathrm{Sym}(V)$ is a symmetric $4$-tensor of $\ell_{2}$ norm at most $1$. Thinking of $\mathrm{Sym}(V)$ as a symmetric matrix over $2$-Tensors, we can write it as $\mathrm{Sym}(V)=\sum_{j=1}^{n^{2}}\lambda_{j}B_{j}\otimes B_{j}$, where the $B_{j}$’s are an orthonormal basis for the set of $2$-tensors and $\sum_{j=1}^{n^{2}}\lambda_{j}^{2}=\|\mathrm{Sym}(V)\|_{2}^{2}\leq 1$. Then, if $c_{i,j}:=\langle A_{i},B_{j}\rangle$, and $c_{j}=\sum_{i=1}^{k}c_{i,j}^{2}$, we have that

[TABLE]

Since the $A_{i}$ are orthonormal, $c_{j}\leq 1$ for all $j$. Furthermore,

[TABLE]

Therefore,

[TABLE]

∎

Let $\mu^{\prime}=\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[Q(X)]$. We now show that since $G_{0}$ is $\varepsilon$-good, then almost all of the difference in Claim 19 must be because of the points in $E$.

Claim 21**.**

$\mathop{\mbox{\bf E}}_{E}[Q(Y)]-\mu^{\prime}\geq\frac{C_{1}-6}{2}\cdot\frac{\xi k}{\phi}$**

Proof.

It suffices to show that

[TABLE]

that is, the good points do not contribute much to the difference. By $(\varepsilon,\delta)$-goodness of $G_{0}$ and Claim 20 we have $\mathop{\mbox{\bf E}}_{0}[Q(Y)]-\mu^{\prime}\leq O\left(\varepsilon\sqrt{k}\right)$. Moreover, we have

[TABLE]

where (a) follows from the boundedness condition of $(\varepsilon,\delta)$-goodness, (b) follows from the last condition of goodness, and (c) follows from hypercontractivity. This shows (9), which completes the proof. ∎

We now show that this implies that there must be a $T$ satisfying the conditions in Algorithm 8.

Claim 22**.**

If $\dim(V_{m})\geq k$, then Algorithm 8 returns a $T$ satisfying the conditions in the algorithm.

Proof.

Suppose not. By the assumption that $\|I-\Sigma\|_{F}\leq\xi$, we have $|\mu^{\prime}|=|\mathop{\mbox{\bf E}}_{X\sim\mathcal{N}(0,\Sigma)}[Q(X)]|\leq 2\xi\|Q\|_{2}\leq\sum_{i=1}^{k}\|p_{i}\|_{2}+\|r_{i}\|_{2}\leq 5\xi k$. Thus, we have

[TABLE]

But since we assume there is no $T$ satisfying the conditions in Algorithm 8, we have

[TABLE]

which is a contradiction, for our choice of $C_{1},C_{2},C_{3}$, and since we chose $k=O(\log^{4}1/\varepsilon)$. ∎

The final thing we must verify is that the number of good points we remove is much smaller than the number of bad points we remove. Formally, we show:

Claim 23**.**

If $\dim(V_{m})\geq k$, then Algorithm 8 returns a $S^{\prime}$ satisfying $\Delta(S^{\prime},G_{0})<\Delta(S,G_{0})$.

Proof.

Let $T$ be the threshold we pick. If $T>C_{3}d^{2}\sqrt{k}\log(|S|)$ then the invariant is satisfied since we remove no good points, by $(\varepsilon,\delta)$-goodness. It suffices to show that we remove $\log 1/\varepsilon$ times many more bad points than good points. Otherwise, by definition we remove a total of

[TABLE]

points. On the other hand, by $(\varepsilon,\delta)$-goodness, hypercontractivity, and Claim 20, we know that the total number of points we throw away is at most

[TABLE]

Since we have maintained this invariant so far, in particular, we have thrown away more bad points than good points, and so $|G_{0}|\leq(1+\varepsilon)|S|$. Moreover, since $T\geq 4B\sqrt{k}\log^{2}1/\varepsilon$, we have

[TABLE]

Therefore, we have $\log(1/\varepsilon)\cdot|G|\cdot\Pr_{G}[Q(X)>T]\leq|S|\cdot\left(\exp\left(-A\left(\frac{T}{4B\sqrt{k}}\right)^{1/2}\right)+\frac{\varepsilon^{2}}{d^{2}\log|G|/\delta}\right)$, and hence, the invariant is satisfied. ∎

Claims 22 and 23 together prove the theorem.

D.3 Proof of Lemma 6

By Lemma 8.16 in [DKK*+*16], the first two items hold together with probability $1-O(\delta)$ after taking $O(\operatorname*{poly}(d,1/\eta,\log 1/\delta))$ samples. Thus, it suffices to show that the last property holds with probability $1-O(\delta)$ given $O(\operatorname*{poly}(d,1/\eta,\log 1/\delta))$ samples. We may clearly WLOG take $\Sigma=I$. Moreover, these properties are clearly invariant under scaling, and hence, it suffices to prove them for degree-four polynomials with $\mathop{\mbox{\bf Var}}_{\mathcal{N}(0,I)}[p(X)]=1$. For any fixed degree-$4$ polynomial $p$, by hypercontractivity, since $P(X)=\frac{1}{n}\sum_{i=1}^{n}p(X_{i})$ has variance $\mathop{\mbox{\bf Var}}[P(X)]=\frac{1}{n}$, we have that

[TABLE]

In particular, if we take

[TABLE]

then

[TABLE]

Since there is a $1/3$-net over all degree-$4$ polynomials with unit variance of size $(1/3)^{O(d^{4})}$, by union bounding, we obtain that if we take $\Omega\left(\frac{d^{4}\log^{1/\delta}}{\eta^{2}}\right)$ samples, then the second to last property holds with probability $1-\delta$.

Finally, the same net technique used to prove Lemma 8.16, along with hypercontractivity, may be used to show the last property holds when given $\operatorname*{poly}(d,\log 1/\eta,\log 1/\delta)$ samples.

D.4 Proof of Lemma 15

Proof.

Let $\mu^{\prime}=\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma)}[p(X)]$, $s=\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma)}[p^{2}(X)]$, and let $\mathop{\mbox{\bf E}}_{S}[p(X)]-\mathop{\mbox{\bf E}}_{\mathcal{N}(0,\Sigma)}[p(X)]=K$. By $(\varepsilon,\delta)$-goodness and Lemma 10, we have

[TABLE]

Moreover, for some appropriate choice of $\beta_{1}$, we have

[TABLE]

where (a), (b), follows from the definition of goodness, and (c) follows from standard Guassian concentration bounds. In particular, (10) and (11) together imply that

[TABLE]

so this implies

[TABLE]

Expanding this yields that

[TABLE]

so since $|\mu^{\prime}|\leq\|\Sigma-I\|_{F}\leq\xi$, we have

[TABLE]

On the other hand, by Lemma 10, we have

[TABLE]

Also, by $(\varepsilon,\delta)$-goodness and Lemma 10, we have

[TABLE]

and, for $A$ as in Corollary 6, and $\beta$ appropriately chosen,

[TABLE]

where (a) follows from goodness, and (b) follows from Corollary 6 and since $\|\Sigma-I\|_{F}=o(1)$. Thus, (12), (13), (14), and (15) together imply that

[TABLE]

and so since $K\geq\Omega\left(\sqrt{\varepsilon\xi}\right)$, this gives the desired bound. ∎