Robustly Learning a Gaussian: Getting Optimal Error, Efficiently
Ilias Diakonikolas∗
CS, USC
diakonik@usc.edu

arXiv:1704.03866v2 [] 5 Nov 2017

Jerry Li §
EECS & CSAIL, MIT
jerryzli@mit.edu

Gautam Kamath†
EECS & CSAIL, MIT
g@csail.mit.edu
Ankur Moitra¶
Math & CSAIL, MIT
moitra@mit.edu

Daniel M. Kane‡
CSE & Math, UCSD
dakane@cs.ucsd.edu
Alistair Stewartk
CS, USC
alistais@usc.edu

November 7, 2017

Abstract

We study the fundamental problem of learning the parameters of a high-dimensional
Gaussian in the presence of noise — where an ε-fraction of our samples were chosen by
an adversary. We give robust estimators that achieve estimation error O(ε) 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 2 and the running time is polynomial in d and 1/ε. When both the mean
and covariance are unknown, the running time is polynomial in d and quasipolynomial
in 1/ε. 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.

∗

Supported by NSF CAREER Award CCF-1652862, a Sloan Research Fellowship, and a Google Faculty Research
Award.
†
Supported by NSF CCF-1551875, CCF-1617730, CCF-1650733, and ONR N00014-12-1-0999.
‡
Supported by NSF CAREER Award CCF-1553288 and a Sloan Research Fellowship.
§
Supported by NSF CAREER Award CCF-1453261, CCF-1565235, a Google Faculty Research Award, and an
NSF Graduate Research Fellowship.
¶
Supported by NSF CAREER Award CCF-1453261, CCF-1565235, a Packard Fellowship, a Sloan Research Fellowship, a grant from the MIT NEC Corporation, and a Google Faculty Research Award.
k
Supported by a USC startup grant.

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
D = (1 − ε)N (µ, σ 2 ) + εZ ,
where N (µ, σ 2 ) is a Gaussian of mean µ and variance σ 2 , and Z is an arbitrary distribution chosen
by an adversary.
Intuitively, among our samples, about a (1−ε) 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 challenging1 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
µ and σ 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 X1 , X2 , . . . , Xn , we can compute
µ
b = median(X1 , X2 , . . . , Xn ) and σ
b=

median(|Xi − µ
b|)
,
−1
Φ (3/4)

where Φ is the cumulative distribution of the standard Gaussian. (This scaling constant is needed
to ensure that σ
b is an unbiased estimator when there is no noise.) If n ≥ C logε21/δ , then with
probability at least 1 − δ we have that dT V (N (µ, σ 2 ), N (b
µ, σ
b2 )) ≤ Cε. In Huber’s contamination
model, this is the strongest type of error guarantee we could hope for2 and captures both the task
of learning the underlying parameters µ and σ 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,
1

None 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.
2
See Lemma 17.

1

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 µ) and scale (e.g., estimating σ 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 N (µ, Σ) 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 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
b ≤ Cε
dT V (N (µ, Σ), N (b
µ, Σ))

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 µ is unknown, our algorithm runs in time polynomial in
the dimension d and 1/ε. When both the mean and covariance are unknown, our algorithm runs
in time polynomial in d and quasi-polynomial in 1/ε. 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(ε log 1/ε) for this problem3 , again with respect to total variation distance. Concurrently, Lai, Rao and Vempala [LRV16] gave an algorithm which achieves estimation error roughly
O(ε1/2 log1/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 ε-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(ε log1/2 1/ε) 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.
3

We note that, as stated, the results in [DKK+ 16] give estimation error O(ε log3/2 1/ε). 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].

2

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 Σ = I is known and only the mean µ 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 dT V (N (µ, I), N (b
µ, I)) ≤
Cε. 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
1−ε
T
T b, and at least a 1−ε fraction of the samples satisfy
2 fraction of the samples satisfy u Xi ≥ u µ
2
T
T
u µ
b ≥ u Xi . However, if we have a candidate µ
b, 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 Θ(ε log1/2 1/ε). This is because along the
direction u, only when a point is Ω(log1/2 1/ε) 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 Θ(ε log1/2 1/ε). 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 = poly(d, 1/ε) samples from the strong contamination
model, where the underlying d-dimensional Gaussian is N (µ, I). Let ε ≤ ε0 , where ε0 is a positive
universal constant. For any β > 0, there is an algorithm to learn an estimate N (b
µ, I) that with
high probability satisfies



1
1
2
ε.
dT V (N (µ, I), N (b
µ, I)) ≤ √ + O √ + ε
β
2
Moreover, the algorithm runs in time poly(n, (1/ε)β ).
We prove an almost matching lower bound of 2ε + Ω(ε2 ) on the estimation error. Thus, our
√
robustness guarantee is optimal up to a factor√of 2, even among computationally inefficient robust
estimators. Interestingly, our extra factor of 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/ 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 ε/2 of the projection of the true mean µ. We can now confine µ to a slab of width ε, and by
taking the intersection of all such slabs we get a convex body that contains µ and has diameter of
3

√
at most ε. By covering
√ the body with a ball of radius ε/ 2, we are guaranteed that the center of
the ball is within ε/ 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(ε log1/2 1/ε)
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 µ = 0 is known and only the covariance Σ is
unknown. In the terminology of robust statistics, this is called robust estimation of scale. In this
b that satisfies4 kΣ − Σk
b F ≤ Cε. When Σ
b does not satisfy
case, we want to compute an estimate Σ
this condition, it can be shown (in Section 6.2.3) that there is a degree-two polynomial p(X), where
E

[p(X)] = 1 and

X∼N (0,Σ)

E

[p(X)] = 1 + C ′ ε .

b
X∼N (0,Σ)

It turns out that, even given the polynomial p(X), deciding whether or not the above conditions
b we can certainly compute E
approximately hold is challenging. Given p(X) and Σ,
b [p(X)].
X∼N (0,Σ)
But given only contaminated samples from N (0, Σ) and without knowing what Σ is, can we estimate
EX∼N (0,Σ) [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 ω(ε)), 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 = poly(d, 1/ε) samples from the strong contamination
model, where the underlying d-dimensional Gaussian is N (0, Σ). There is an algorithm to learn an
b that runs in time poly(n, (1/ε)O(log 4 1/ε) ) and with high probability satisfies
estimate N (0, Σ)
b ≤ Cε ,
dT V (N (0, Σ), N (0, Σ))

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

More precisely, to obtain O(ε) error guarantee with respect to the total variation distance, we need to robustly
approximate Σ within O(ε) 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 Σ is close to the identity matrix, in which case the
Frobenius error suffices.

4

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 Σ as


ΣV AT
,
Σ=
A ΣV ⊥

and have accurate estimates for ΣV and ΣV ⊥ , 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 µ and Σ 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 = poly(d, 1/ε) samples from the strong contamination
model, where the underlying d-dimensional Gaussian is N (µ, Σ). There is an algorithm to learn an
b that runs in time poly(n, (1/ε)O(log 4 1/ε) ) and with high probability satisfies
estimate N (b
µ, Σ)
b ≤ Cε ,
dT V (N (µ, Σ), N (b
µ, Σ))

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/ε, and the running time
depends polynomially on d and quasi-polynomially on 1/ε. 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/ε (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
5

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 EF [f (X)] = EX∼F [f (X)] denote the expectation of f (X) under F .
If S is a finite set, we let ES [f (X)] = EX∼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 ⊆ Rd , we let ΠV : Rd → Rd be the projection operator onto V . If V = span(v) is 1-dimensional,
we will denote this projection as Πv .

2.1

The Strong Contamination Model

Here we formally define the strong contamination model.
Definition 2. Fix ε > 0. We say a set of samples X1 , . . . , Xn was generated from the strong
contamination model on a distribution F , if it was generated via the following process:
1. We produce (1 − ε)n i.i.d. samples G from F .
2. An adversary is allowed to observe these samples and add εn points E arbitrarily.
We are then given the set of samples G ∪ E in random order. Also, we will say that the samples
X1 , . . . , Xn are ε-corrupted. Moreover given an ε-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 ′ ⊂ S, we will also write S ′ = (G′ , E ′ ), where G′ = S ′ ∩ G and E ′ = S ′ ∩ E denote
the set of uncorrupted points and corrupted points remaining in S ′ . L will denote G \ G′ , which is
the set of “lost” uncorrupted points.
Given a contaminated set S ′ = (G′ , E ′ ) and a set G so that G′ ⊆ G, define the following
quantities
φ(S ′ , G) =

|G \ G′ |
,
|S ′ |

ψ(S ′ , G) =

|E ′ |
,
|S ′ |

∆(S ′ , G) = ψ(S ′ , G) + φ(S ′ , G) log

1
.
φ(S ′ , G)

(1)

In particular, observe that if ∆(S ′ , G) < O(ε), then a simple calculation implies that φ(S ′ , G) ≤
O(ε/ log 1/ε). Equivalently, we have removed at most an O(ε/ log 1/ε) fraction of good points from
G. This is crucial, as if we throw out an ε-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 ε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 ε (for n sufficiently large) can be simulated by a
(1 + o(1))ε-corrupted set of samples, except with negligible failure probability. Hence, we lose
only an additive o(ε) 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.
6

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 Rd and η, δ > 0. We say that G is (η, δ)-good with
respect to N (µ, I) if the following hold:
p
• For all x ∈ G we have kx − µk2 ≤ O( d log(|G|/δ)).

• For every affine function L : Rd → R we have | PrG (L(X) ≥ 0) − PrN (µ,I) (L(X) ≥ 0)| ≤
η/(d log(d/ηδ)) .
• We have that k EG [X] − EN (µ,I) [X]k2 ≤ η.
• We have that kCovG [X] − Ik2 ≤ η/d.
• For any even degree-2 polynomial p : Rd → R we have that
E[p(X)] −
G

E[p2 (X)] −
G

E [p(X)] ≤ η

N (µ,I)

E [p2 (X)] ≤ η

N (µ,I)

E [p2 (X)]1/2 ,

N (µ,I)

E [p2 (X)], and

N (µ,I)

Pr[p(X) ≥ 0] ≤ Pr [p(X) ≥ 0] +
G

N (µ,I)

η
.
d log(|G|/δ)

It is easy to show (see Lemma 6) that given enough samples from N (µ, 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 Rd , and η, δ > 0. We say that G is (η, δ)-good with
respect to N (0, Σ) if the following hold:
• For all x ∈ G we have that xT Σ−1 x = O(d log(|G|/δ)).

7

• For any even degree-2 polynomial p : Rd → R we have
E[p(X)] −
G

E[p2 (X)] −
G

E [p(X)] ≤ η

N (0,Σ)

E [p2 (X)] ≤ η

N (0,Σ)

E [p2 (X)]1/2 ,

N (0,Σ)

E [p2 (X)] and

N (0,Σ)

Pr [p(X) ≥ 0] ≤ Pr [p(X) ≥ 0] +

X∼G

N (0,Σ)

η2
.
d log(|G|/δ)

• For any even degree-4 polynomial p : Rd → R we have
E[p(X)] −
G

E [p(X)] ≤ η Var [p(X)]1/2 ,

N (0,Σ)

N (0,Σ)

Pr[p(X) ≥ 0] ≤ Pr [p(X) ≥ 0] +
G

N (µ,I)

η2
.
2 log(1/ε)(d log(|G|/δ))2

As before, it is easy to show (see Lemma 14) that given enough samples from N (0, Σ), 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 ℓ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 µ1 , µ2 ∈ Rd be such that kµ1 − µ2 k2 = ε for ε < 1. Then


1
dTV (N (µ1 , I), N (µ2 , I)) = √ + o(1) ε .
2π
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):
b be such that kΣ−IkF ≤ O(ε log 1/ε), and kΣ− Σk
b F ≤
Lemma 2 (Cor. 2.14 in [DKK+ 16]). Let Σ, Σ
b ≤ O(ε).
Cε. Then dTV (N (0, Σ), N (0, Σ))

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:

8

Theorem 3. Fix µ ∈ Rd , and let ε, γ, δ > 0. Let S0 = (G0 , E0 ) be such that G0 is a (γε, δ)-good
set with respect to N (µ, I), and |E0 |/|S0 | ≤ ε. Let S = (G, E) be another set with ∆(S, S0 ) < ε.
Let V ⊆ Rd be a subspace. For all 0 < ρ < 1, the algorithm LearnMeanLowD(V, γ, ε, δ, S, ρ)
runs in time poly(d, |S|, (1/ρ)O(dim(V )) , log(ρε/(1 − ρ)), log(1/ρ)) and returns a µ
e so that
kΠV (µ − µ
e)k2 =

 γ 
1 + 2ρ √
ε.
π+O
1−ρ
d

√
In particular, as we let ρ, γ → 0, the √
parameter estimation error approaches πε (corresponding to
a total variation approximation of ε/ 2). In Lemma 17
p in the Appendix we show√that no algorithm
can achieve parameter estimation error better than π2 ε. Thus, we achieve a 2 approximation
to the optimal error.
For simplicity, in the rest of this section, we will let V = Rd , 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 (γε, δ)-good with respect to the underlying
d
Gaussian. Then, we show that
p π by using a net argument, we can produce a convex body in R with
diameter at most 2R = 2( 2 + o(1))ε 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 2R (see [BW41]
for an English language version of the result). We use the center of the ball as our estimate µ
e.

3.1

Robustness of the Median

First we show that if we projectponto one dimension, then the median of the corrupted data differs
from the true mean by at most π2 ε + o(ε). Our proof will rely only on the notion of a (γε, δ)-good
set with respect to N (µ, I) and thus it works even in the strong contamination model. Formally,
we show:
Lemma 3. Fix any v ∈ Rd . Fix µ ∈ Rd , and let δ > 0. Let S0 = (G0 , E0 ) be so that G0 be a
set 
with
(γε, δ)-good set with respect to N (µ, I), and |E0 |/|S0 | ≤ ε. Let S = (G, E) be another
pπ
γ
ε.
∆(S, S0 ) < ε. Let b be the median of S when projected onto v. Then, |b − Πv µ| ≤
2 +O d

Proof. For any a ∈ R, we have

Pr [hv, Xi > a] =

X∼S

Observe that we have

|E|
|S|

|G|
|E|
Pr [hv, Xi > a] +
Pr [hv, Xi > a] .
|S| X∼G
|S| X∼E

PrX∼E [hv, Xi > a] ≤ ψ(S, G). Moreover, by simple calculation we have
Pr [hv, Xi > a] − Pr [hv, Xi > a] ≤ 2φ(S, G) .

X∼G

X∼G0

Hence, we have
Pr [hv, Xi > a] −

X∼S

|G|
Pr [hv, Xi > a] ≤ ψ(S, G) + o(ε) ,
|S| X∼G0

since by assumption ∆(S, G) ≤ ε. Similarly, we have that for all a ∈ R,
Pr [hv, Xi < a] −

X∼S

|G|
Pr [hv, Xi < a] ≤ ψ(S, G) + o(ε) .
|S| X∼G0
9

For |a| = O(ε) we have that PrN (0,I) [X > a] = 12 − √12π a + O(ε3 ). Thus, by (γε, δ)-goodness of G0 ,
this implies that for |a| = O(ε), we have


 γε 
|G| 1
1
Pr [hv, Xi > Πy µ + a] −
.
(2)
− √ a ≤ ψ(S, G) + O
X∼S
|S| 2
d
2π
Thus, for |a| = O(ε) we have

r



 γε 
|G| 1
1
|E|
π
+ o(ε)
a] ≤
−√ a +
+O
Pr [hv, Xi > Πy µ +
X∼S
2
|S| 2
|S|
d
2π
 γε 
1 |G| 1
|E|
√ a+
≤ −
+ o(ε) .
+O
2
|S| 2π
2|S|
d


pπ 
p
In particular, we see that if a > π2 ε + O γε
d  + o(ε), then PrX∼S
p π hv, Xi > Πy µ + 2 ε < 1/2.
By symmetric logic, we also have that Pr
hv, Xi > Πy µ − 2 ε > 1/2. Thus, the median in
pX∼S
π
direction v differs from Πy µ by at most 2 ε + O γε
d + o(ε).

3.2

Finding a Minimum Radius Circumscribing Ball

For any x ∈ Rd and r > 0, let B(x, r) = {y ∈ Rd : kx − yk2 ≤ 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 C ⊆ Rd be a convex √
body so that for all x, y ∈ C,
we have kx − yk2 ≤ 2R. Then C is contained within a ball of radius R 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 ∈ Rd , outputs a point which is almost the closest point in C to x:
Definition 5. A ρ-projection oracle for a convex body C is a function O : Rd → Rd , which, given
a point y ∈ Rd , outputs a point x ∈ C so that kx − yk2 ≤ inf x′ ∈C kx′ − yk2 + ρ.
Our first step is to use such an oracle to construct a net for C. First, we need the following
well-known bound on the size of the net.
Claim 1. Fix r > 0. Then, for any β > 0, there is a β-net F for the sphere of radius r around 0
in Rd of size (r/β)O(d) . Moreover, this net can be constructed in time poly(d, |F|).
With this, we can show:
Lemma 4. Fix R, C as in Theorem 4, and let 1 > ρ > 0. Let x ∈ C be arbitrary. Let O be a
(ρR/3)-projection oracle for C. Suppose a call to O runs in time T . Then, there is an algorithm
CircumscribeNet(R, ρ, O, x) which runs in time poly((R/ρ)O(d) , T ) and outputs a set X ⊆ Rd
so that X is a (ρR)-net for C, and moreover, |X | ≤ (R/ρ)O(d) .
The algorithm is fairly straightforward. First, we observe that C is contained within B(x, 2R).
We then form a (ρ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 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.
10

Algorithm 1 Generating a net of C
1: function CircumscribeNet(R, ρ, O, x)
2:
Form an ρ/3-net F ′ of the sphere of radius 2 of size (1/ρ)O(d) as in Claim 1.
3:
Let F = R · F ′ + x.
4:
Let X ← ∅
5:
for each v ∈ F do
6:
Let uv ← O(v)
7:
if kv − uv k2 ≤ 2ρR/3 then
8:
Add uv to X
9:
end if
10:
end for
11:
return X
12: end function
Proof. The runtime bound follows from Claim 1. We now turn our attention to correctness. By
Claim 1, and rescaling and shifting, the set F is clearly a (ρR)/3-net for a ball B of radius 2R
containing C. We now claim that the set X is indeed a (ρR)/3-net for C. Fix y ∈ C. Since C ⊆ B,
this implies there is some v ∈ F so that ky − vk2 ≤ ρR/3. Thus, in Line 7, when processing v,
we must find some uv ∈ C so that kuv − vk2 ≤ 2ρR/3. The claim then follows from the triangle
inequality.
With this, we obtain:
Corollary 2. Fix R, C, ρ, O, x as in Lemma 4. Suppose a call to O runs in time T . Then, there is
O(d) , T ) and returns a point
an algorithm Circumscribe(R, ρ, O, x) which runs
√ in time poly((R/ρ)
yb so that C is contained within a ball of radius 2(1 + 2ρ)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 2(1 + ρ)R, and output any such point. The formal pseudocode for Circumscribe
is given in Algorithm 2.
Algorithm 2 Finding a circumscribing ball of small radius
1: function Circumscribe(R, ρ, O, x)
2:
Form an ρ/3-net F ′ of B(0, 2) of size (1/ρ)O(d) as in Claim 1.
3:
Let F = R · F ′ + x.
4:
Let X ← CircumscribeNet(R, ρ, O, x).
5:
for each v ∈ F do
√
6:
if for all u ∈ X , we have ku − vk2 ≤ 2(1 + ρ)R then
7:
return u
8:
end if
9:
end for
10: end function
Proof. The√runtime bound is immediate. By Theorem 4, there is some y ∈ B(x, √
2R) so that
′
C ⊆ B(y, R 2). Thus, by the triangle inequality, there is some y ∈ F so that C ⊆ B(y, 2(1+ρ)R).
Thus, the algorithm will output some point y ′′ ∈ F. By an additional
application of the triangle
√
′′
inequality, since X is a ρR-net for C, this implies that C ⊆ B(y , 2(1 + 2ρ)R), as claimed.
11

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 ρ > 0. We first produce a ρ-net F over the unit sphere in Rd .
By (a slight modification of) Claim 1, this net has size (1/ρ)O(d) and can be constructed in time
poly(d, |F|). For each v ∈ F, we project all points in S onto v, and take the median of these points
to produce bv . We then construct the following set:
\
C=
{y ∈ Rd : hv, yi ∈ [bv − β, bv + β]} ,
(3)
v∈F

pπ

γε 

where β = 2 ε + O d + o(ε) 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 C
has small diameter:
Claim 2. For all x, y ∈ C, we have kx − yk2 ≤ 2β/(1 − ρ).
Proof. Fix any x, y ∈ C. By definition of C, it follows that for all v ∈ F, we have |hx − y, vi| ≤ 2β.
For any u with kuk2 = 1, there is some v ∈ F with ku − vk2 ≤ ε, and so we have
|hx − y, ui| ≤ |hx − y, vi| + |hx − y, u − vi|
≤ 2β + ρkx − yk2 .

Taking the supremum over all unit vectors u and simplifying yields that kx − yk2 ≤ 2β/(1 − ρ), as
claimed.
The second property shows that we may find an α-projection oracle for C efficiently.
Claim 3. Fix ρ′ > 0. There is a ρ′ -projection oracle ProjOracle(y, ρ′ , C) for C which runs in
time poly((1/ρ)O(d) , log(γε/(1 − ρ)), log(1/ρ′ )).
Proof. The projection problem may be stated as
min kx − yk2 s.t.hv, yi ∈ [bv − β, bv + β], ∀v ∈ F .
This is a convex minimization problem with linear constraints. By the classical theory of optimization [GLS88], finding a ρ-approximate y can be done in poly(d, log(vol(C)/ρ′ )) queries to a
separation oracle for C. Since the separation oracle must only consider the constraints induced by
F, this can be done in time (1/ρ)O(d) . Since by Claim 2 we have vol(C) ≤ (2β/(1 − ρ))O(d) , the
desired runtime follows immediately.
We now finally describe LearnMeanLowD. Using convex optimization, we first find an arbitrary x ∈ C. By Lemma 3 we know µ ∈ C and so this step succeeds. After constructing 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.
C. By
√ 1+2ρBy Lemma 3, we know that µ√∈ 1+2ρ
Claim 2 and Corollary 2, the output y satisfies B(y, 2 1−ρ β). Thus, we have kµ−yk2 ≤ 2 1−ρ β,
as claimed.

12

Algorithm 3 Finding a circumscribing ball of small radius
1: function LearnMeanLowD(ε, δ, S, ρ)
2:
Form a ρ-net F of B(0, 1) of size (1/ρ)O(d) as in Claim 1.
3:
for each v ∈ F do
4:
Let bv be the median of S projected onto v.
5:
end for
6:
Form C as in Equation (3).
7:
Find an p
x ∈ C using convex
optimization.

+
o(ε)
8:
Let β = π2 ε + O γε
d
9:
Let R = β/(1 − ρ)
10:
Let O(·) = ProjOracle(·, (ρR)/3, C)
11:
return the output of Circumscribe(R, ρ, O, x)
12: end function

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 ε, γ, δ > 0, and let X1 , . . . , Xn be an ε-corrupted set of points from N (µ, I), where
kµk2 ≤ O(ε log 1/ε), and where


(d log(d/γεδ))6
.
n=Ω
γ 2 ε2
Then, for every α, β > 0, there is an algorithm RecoverMean(X1 , . . . , Xn , ε, δ, γ, α, β) which
runs in timepoly(d,
1/γ, 1/εβ , 1/α, log 1/δ) and outputs a µ
b so that with probability 1 − δ, we have
√
π+O(γ)
1
kb
µ − µk2 ≤
+ √β ε.
1−α
In particular, observe that Theorem 5, in conjunction with Lemma 1, gives us Theorem 1, if we set
γ = o(1). With this, we may state our primary algorithmic contribution:

Theorem 6. Fix ε, γ, α, δ, β > 0, and let S0 = (G0 , E0 ) be an ε-corrupted set of samples of size
n from N (µ, I), where kµk2 ≤ O(ε log 1/ε), and where n = poly(d, 1/(γε), log 1/δ). Suppose that
G0 is (γε, δ)-good with respect to N (µ, I). Let S ⊆ S0 be a set so that ∆(S, G0 ) ≤ ε. Then, there
exists an algorithm FilterMeanOpt that given S, ε, γ, α, β outputs one of two possible outcomes:
√

π+O(γ)
√1
+
ε.
(i) A µ
b, so that kb
µ − µk2 ≤
1−α
β
(ii) A set S ′ ⊂ S so that ∆(S ′ , G0 ) < ∆(S, G0 ).

Moreover, FilterMeanOpt runs in time poly(d, 1/γ, 1/εβ , 1/α, log 1/δ).
+
By
p first running the algorithm of [DKK 16] to obtain an estimate of the mean to error
O(ε log 1/ε), 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/ε) 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/ε) eigenvectors of the empirical covariance

13

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(log2 1/ε) 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 ∆ (which gives an additional logarithmic penalty to discarding good
points) and our guarantee that ∆ decreases, our filter can only afford to throw out an ε/ log(1/ε)
fraction of good points in total, since ∆ is initially ε. This is critical, as if we threw away an
ε-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(ε)-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 ε, γ, δ, α, β > 0, and let S0 = (G0 , E0 ) be an ε-corrupted set of samples of size
n from N (µ, I), where kµk2 ≤ O(ε log 1/ε), and where n = poly(d, 1/(γε), log 1/δ). Suppose that
b be the
G0 is (γε, δ)-good with respect to N (µ, I). Let S ⊆ S0 be a set so that ∆(S, G0 ) ≤ ε. Let Σ
sample covariance of S, let µ
b be the sample mean of S, and let V be the subspace of all eigenvectors
b − I with eigenvalue more than 1 ε. Then, there exists an algorithm FilterMeanManyEig
of Σ
β
that given S, ε, γ, δ, α, β outputs one of two possible outcomes:
1. If dim(V ) ≥ C1 β log(1/ε), then it outputs an S ′ so that ∆(S ′ , G0 ) < ∆(S, G0 ).

2. Otherwise, the algorithm outputs “OK”, and outputs an orthonormal basis for V .
b − I, and if there are too
Our algorithm works as follows: It finds all large eigenvalues of Σ
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
14

Algorithm 4 Filter if there are many large eigenvalues of the covariance
1: function FilterMeanManyEig(S, ε, γ, δ, α, β)
2:
Let C1 , C2 , C3 > 0 be sufficiently large constants.
b be the empirical mean and covariance of S, respectively.
3:
Let µ
b and Σ
b − I with eigenvalue more than
4:
Let V be the subspace of Rd spanned by eigenvectors of Σ
1
β ε.
5:
if dim(V ) ≥ C1 β log(1/ε) then
6:
Let V ′ be a subspace of V of dimension C1 β log(1/ε).
√
7:
Let µ̃ be an approximation to ΠV ′ (µ) with ℓ2 -error π+O(γ)
1−α ε, computed using
LearnMeanLowD(V, γ, ε, δ, S, γ).
8:
Let p(x) be the quadratic polynomial p(x) = kΠV ′ (x) − µ̃k22 − dim(V ′ ).
9:
Find a value T > 0 so that either:
(a) T > C2 d log(|S|/δ) and p(x) > T for at least one x ∈ S, or

(b) T > 2C3 log(1/ε)/c0 and PrS (p(x) > T ) > exp(−c0 T /(2C3 )) + γε/(d log(|S|/δ)).
10:
11:
12:
13:
14:

return S ′ = {x ∈ S : p(x) ≤ T }
else
return an orthonormal basis for V .
end if
end function

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.
Algorithm 5 Return a mean if there are few large eigenvalues of the covariance
1: function FilterMeanFewEig(S, ε, γ, δ, α, β, V )
√
π+O(γ)
2:
Let µ̃V be an approximation to ΠV (µ) with ℓ2 -error
1−α ε, computed using
LearnMeanLowD(V, γ, ε, δ, S, γ).
3:
Let µ̃V ⊥ be the empirical mean on V ⊥ , ΠV ⊥ µ̂.
4:
return µ̃V + µ̃V ⊥ .
5: end function
Lemma 5. Let µ, η, G0 , S be as in Theorem 7. Let µ
b be the sample mean of S, and let v be a unit
ε
vector. Suppose that hv, µ − µ
bi > β 1/2
. Then VarS [hv, Xi] > 1 + βε .
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 η, δ > 0. Let X1 , . . . , Xn be independent samples from N (µ, I), where n =
Ω((d log(d/ηδ))6 /η 2 ). Then, S = {X1 , . . . , Xn } is (η, δ)-good with respect to N (µ, I) with probability at least 1 − δ.
15

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
√
π+O(γ)
3 guarantees that the mean is accurate up to ℓ2 -error 1−α ε. Within the subspace V ⊥ , the
ε
contrapositive of the statement of Lemma 5 guarantees the mean is accurate up to ℓ2 -error β 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 χ, ε, δ > 0, and let X1 , . . . , Xn be an ε-corrupted set of points from N (µ, Σ),
where kΣ − Ik2 ≤ O(χ), kµk2 ≤ O(ε log 1/ε), and where n = poly(d, 1/χ, 1/ε, log 1/δ). For
any γ > 0, there is an algorithm RecoverMeanNoisy(X1 , . . . , Xn , ε, δ, γ, χ) which runs in time
poly(d, 1/χ, 1/ε, log 1/δ) and outputs a µ
b so that with probability 1 − δ, we have kb
µ − µk2 ≤
(C + γ)ε + O(χ).
This extension follows from two elementary observations:

1. 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(ε + α).
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(ε + α), we
would immediately get the desired bound.

5

Robustly Estimating the Mean of Degree Two Polynomials

In this section, we give robust estimates of E[p2 (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 E[p2 (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 ∼ N (0, I) ∈ Rd and A be a d × d
matrix. Then for some absolute constant c0 , for every t ≥ 0,


 2

t
t
T
T
,
.
Pr X AX − E[X AX] > t ≤ 2 exp −c0 · min
kAk2F kAk2
We will also require the following lemmata:

16

Lemma 8 (Hölder’s inequality for Schatten norms). Let A, B be matrices. Then, for all p, q so
that 1p + 1q = 1, we have kABkS 1 ≤ kAkS p kBkS q .
This implies the following corollary:
b M be so that kΣ − Σk
b F ≤ O(δ), and so that kM kF = 1. Then, we have
Corollary 3. Let Σ, Σ,
1/2
1/2
1/2
1/2
b
b
kΣ M Σ − Σ M Σ kS 1 ≤ 5δ.
Proof. We have

b 1/2 M Σ
b 1/2 kS 1 ≤ kΣ1/2 M Σ1/2 − M kS 1 + kM − Σ
b 1/2 M Σ
b 1/2 kS 1 .
kΣ1/2 M Σ1/2 − Σ

We will bound the first term on the RHS by 5δ/2; the second term is bounded symmetrically. We
have
kΣ1/2 M Σ1/2 − M kS 1 ≤ kΣ1/2 M Σ1/2 − M Σ1/2 kS 1 + kM Σ1/2 − M kS 1
kΣ1/2 − IkS 2 kM Σ1/2 kS 2 + kΣ1/2 − IkS 2 kM kS 2
≤ 5δ/2 ,

where the last line follows from Hö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 ε, δ, τ > 0. Let Σ be so that kΣ − IkF ≤ O(ε log 1/ε), and fix a p ∈ P2 , where P2
denotes the set of even degree-2 polynomials in d variables. Let G0 be an (ε, δ)-good set of samples
from N (0, Σ), and let S = {X1 , . . . , Xn } be so that ∆(S, G0 ) ≤ ε. Then, for any C > 0 there is an
algorithm LearnMeanChiSquared which, given p, X1 , . . . , Xn , and ε, outputs a µ
b so that with
probability 1 − τ over the randomness of the algorithm,
µ
b−

E

X∼N (0,Σ)

[p(X)] ≤ kΣ − IkF /C + O(log(C)ε) .

Moreover, the algorithm runs in time O(|S| + log(1/τ )/ε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 highdimensional 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 N (0, Σ), we use our corrupted data set to
estimate the mean of p(X) conditioned on the event that |p(X)| ≤ T . Then, to estimate the contribution of the mean from points X so that |p(X)| > T , we estimate this by EX∼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 ∼ 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
17

too much. Secondly, because we know that kΣ − IkF 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 ′ ) for X ′ ∼ N (0, I) explicitly and wish to use this to get
a better estimate for the mean of p(X) for X ∼ N (0, Σ) than might be given by the mean of the
ε-corrupted set of samples.
Algorithm 6 Approximating E[p(X)] for X ∼ N (0, Σ) with corrupted samples.
1: function LearnMeanChiSquared(X1 , . . . , Xn , p(x), ε, τ )
2:
Let T = O(log
 C).

x − T, for x ≥ T
3:
Let f (x) =
0,
for |x| ≤ T .


x + T, for x ≤ −T
P
4:
Compute α = ni=1 (p(Xi ) − f (p(Xi )))/n.
′ from X ′ ∼ N (0, I).
2
5:
Simulate m = O((ln
X1′ , . . . , Xm
Pmτ )/ε ) samples
′
6:
Return µ
b = α + i=1 f (p(Xi ))/n.
7: end function

It follows from (ε, δ)-goodness that | PrX ′ ∼N (0,Σ) [p(X ′ ) > t] − #{Xi : p(Xi ) > t}/n| ≤ 2ε for all
t. We need to express the expectation in terms that we can use this to bound. For Z = p(X ′ ), we
have that
Z T
Z T
Pr[Z < −t]dt .
Pr[Z > t]dt −
E[Z − f (Z)] = E[max{T, min{−T, Z}}] =
0

0

Similarly, the samples have
α=

n
X
i=1

(p(Xi ) − f (p(Xi ))/n =

Z

T
0

#{Xi : p(Xi ) > t}
dt −
n

Z

t
0

#{Xi : p(Xi ) < −t}
dt .
n

Thus, we have | E[Z − f (Z)] − α| ≤ 2T ε.
Since p ∈ P2 , we have E[p(X ′ )] = 1 for X ′ ∼ N (0, I). Thus, we have Var[f (p(X ′ ))] ≤
E[f (p(X ′ ))2 ] ≤ E[p(X ′ )2 ] = 1. It follows by
results that the empirical after
Pstandard concentration
′ )/n − E[p(X ′ )]| ≤ ε with probability 1 − τ .
taking m = O(ln(1 − τ )/ε2 ) samples has | m
f
(p(X
)
i
i=1
When this holds, we have
µ
b−

E

X∼N (0,Σ)

[p(X)] ≤ (2T + 1)ε +

E

[f (p(X))] −

X∼N (0,I)

E

[f (p(X))] .

X∼N (0,Σ)

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
E

[f (p(X))] −

X∼N (0,I)

E

[f (p(X ′ ))] ≤

X ′ ∼N (0,Σ)

18

kΣ − IkF
.
C

Proof. Let M be the symmetric matrix with kM kF = 1 such that p(x) = xT M x for x ∈ Rd . We
can write p(X) with X ∼ N (0, I) as
p(X) + tr(M ) = X T M X = X T O T DOX = Y T DY =

d
X

ai Yi2 ,

i=1

where O is orthogonal, D is diagonal, ai are the eigenvalues
of M , and Y ∼ N (0, I) hence Yi are
P 2
i.i.d. from N (0, 1). Since kM kF = 1, here we have
ai = 1. If instead we express p(X ′ ) in terms
of X ′ ∼ N (0, σ 2 ), we obtain:
′

p(X ) + tr(M ) = Y

′T

Σ

1/2

MΣ

1/2

′

Y =Y

′T

′T

′

′

′

T

′

O DOY =Y DY =

d
X

bi Yi2 ,

i=1

where O′ is orthogonal, D ′ is diagonal, bi are the eigenvalues of Σ1/2 M Σ1/2 , and Y ′ , Y ∼ N (0, I)
hence Yi are i.i.d. from N (0, 1). P
By Corollary 3, we have that i |bi − ai | ≤ (5/2)kΣ − IkF . Now consider the random variables
Zi,λ = −tr(M ) +
where

d
X
j=1

cj (Yi2 − 1) ,



 bj ,
cj = (1 − λ)ai + λbi


aj

for j < i
for j = i ,
for j > i

for 1 ≤ i ≤ d and 0 ≤ λ ≤ 1. Note that Zi,1 = Zi+1,0 , for 1 ≤ i ≤ d − 1. Note that, to prove the
lemma, it suffices to show that
X
2|bi − ai |/5C .
| E[f (Z1,0 )] − E[f (Zd,1 )]| ≤
i

To this end,
Pconsider how
P E[f (Zi,λ )] varies with λ. We can write Zi,λ = Z−i + Zi , where Z−i =
−tr(M ) + j<i bj Yj2 + j>i aj Yj2 ; and Zi = ci Yi2 , where ci = ((1 − λ)ai + λbi ). We assume for
now that ci 6= 0. Since only Zi depends on λ, we have
d E[f (Zi,λ )]
d
=
E[f (Z−i + Zi )]
dλ
dλ 

d
E [f (Z−i + Zi )]
= E
Z−i dλ Zi
 Z ∞

d
= E
f (Z−i + x)Pi (x)dx
Z−i dλ −∞

Z ∞
f (Z−i + x)dPi (x)/dλdx ,
= E
Z−i

−∞

where Pi (x) is the probability density function of the random variable Zi . Standard results about
the χ2 distribution give that:
Fact 1. Let Y1 , Y2 , Y3 ∼ N (0, 1). Then the probability density function of Y12 is

Y12

+

Y22

is

1 −x/2
,
2e

and of

Y12

+

Y22

+

Y32

is

√
√ x e−x/2 .
2π

19

√ 1 e−x/2 ,
2πx

of

This gives that Pi (x) =

√ 1
e−x/2ci .
2πxci

Now consider the derivative:

1
1
/dλ · e−x/2ci + √
· de−x/2ci /dλ
2πxci
2πxci
1
dci
x
1
dci −x/2ci
=− √
· 2 · e−x/2ci ·
·e
+√
·
3/2 dλ
dλ
2πxci 2ci
2 2πxci
√
(bi − ai ) x −x/2ci
bi − a i
e−x/2ci + √
e
=− √
3/2
5/2
2 2πxci
2πci

dPi (x)/dλ = d √

(3)

= ((bi − ai )/2ci )(Pi (x) − Pi (x)) ,
(3)

where Pi (x) is the distribution of Zi + Zi′ + Zi′′ , where Zi′ and Zi′′ are i.i.d. copies of Zi . We thus
have

Z ∞
d E[f (Zi,λ )]
f (Z−i + x)dPi (x)/dλdx
= E
Z−i
dλ
−∞
= ((bi − ai )/2ci ) E [E[f (Z−i + Zi )] − E[f (Z−i + Zi + Zi′ + Zi′′ )]]
Z−i

= ((bi − ai )/2ci ) E[f (Zi,λ ) − f (Zi,λ + Zi′ + Zi′′ )]

= ((bi − ai )/2ci ) E [f (Zi,λ ) − ′E ′′ [f (Zi,λ + Zi′ + Zi′′ )]] .
Zi,λ

Zi ,Zi

Since f has Lipschitz constant 1, |f (Zi,λ ) − EZi′ ,Zi′′ [f (Zi,λ + Zi′ + Zi′′ ]| ≤ E[Zi′ + Zi′′ ] = 2ci whatever
value Zi,λ takes.
Using the probability distribution of Zi′ + Zi′′ , in the case ci > 0, we have for all |z| ≤ T ,
Z ∞
′′
′
e−x/2ci /2ci · (z + x − T )dx
E ′′ [f (z + Zi + Zi )] =
′
Zi ,Zi

T −z
−(T −z)/2ci

=e

(z − T + 2ci + (T − z))

= 2ci e−(T −z)/2ci .

Since ci ≤ max ai , bi ≤ 1 + O(ε log(1/ε)) ≤ 2, for large enough T = O(log C), we have e−(T /2)/2ci ≤
1/C 2 T . We assume that this holds.
For −T ≤ z ≤ T /2, we have 2ci e−(T −z)/2ci ≤ 2ci /C 2 T . A similar argument when ci < 0 gives
that for −T /2 ≤ z ≤ T , we have 2|ci |e−(T −z)/2ci ≤ 2|ci |/C 2 T . Now we have enough to show that
d E[f (Zi,λ )]
≤ |(bi − ai )/2ci || E [f (Zi,λ ) − ′E ′′ [f (Zi,λ ) + Zi′ + Zi′′ )]|
Zi,λ
dλ
Zi ,Zi
≤ Pr [|Zi,λ − tr(M )| ≥ T /2]|(bi − ai )| + |(bi − ai )|/C 2 .
Zi,λ

By the Hanson-Wright inequality, we have that for any x,
√
Pr [|Zi,λ − E[Zi,λ ]| ≥ 2kck2 x + kck∞ x] ≤ 2e−x ,
Zi,λ

where c = (b1 , . . . , bi − 1, ci , ai+1 , . . . , ad ). Note that kck2 ≤ 1 + O(kΣ − IkF ) ≤ 2. Also, E[Zi,λ ] +
P
tr(M ) is the sum of coordinates of c, and so | E[Zi,λ ]| ≤ ij=1 |bj − aj | ≤ O(kΣ − IkF ) ≤ 1. Putting
this together, we obtain
Pr [|Zi,λ − tr(M )| ≥ T /2] ≤ 2 exp(−((T − 1)/8)2 ) ≤ 1/C 2 .

Zi,λ

20

Finally, we have that, assuming ci 6= 0,
d E[f (Zi,λ )]
≤ 2|(bi − ai )|/C 2 .
dλ
Now we need to deal with the special case ci = 0. Note that since f is Lipschitz, for any β ∈ R,
we have | E[f (Z−i + βYi2 ) − f (Z−i )]|/|β| ≤ | E[Yi2 ]| ≤ 1. By considering the limit as β tends to
d E[f (Zi,λ )]
≤ 1.
zero from above or below, we get that when ci = 0, the derivative still exists and
dλ

Since the limit 2|(bi − ai )|/C 2 only does not apply at one point where E[f (Zi,λ )] is still continuous
as a function of λ, this is not an issue. We still obtain that
| E[f (Zi,0 )] − E[f (Zi,1 )]| ≤ 2|(bi − ai )|/C 2 .

Recalling that Zi,1 = Zi+1,0 with Z1,0 = p(X ′ ) and Zd,1 = p(X) for X ′ ∼ N (0, 1) and X ∼ N (0, Σ),
we have by the Mean Value Theorem, that
X
| E[f (p(X ′ ))] − E[f (p(X)]| ≤
| E[f (Zi,0 )] − E[f (Zi,1 )]|
i

≤ 2|(bi − ai )|/C 2 ≤ O(kΣ − IkF /C 2 )

≤ kΣ − IkF /C ,
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 ε, τ > 0. Let Σ be so that kΣ−IkF ≤ O(ε log 1/ε). Let G0 be an (ε, δ)-good set of
samples from N (0, Σ), and let S = {X1 , . . . , Xn } be so that ∆(S, G0 ) ≤ ε. Let W1 be a subspace of
degree-2 polynomials, and let W2 be an orthogonal subspace of degree-2 polynomials, so that we have
b so that EX∼N (0,Σ) [p(X)] − E
aΣ
b [p(X)] ≤ ξ for all p ∈ W2 . Then there is an algorithm
X∼N (0,Σ)
b runs in time poly(d, |S|, 2O(dim(W1 )) , log 1/τ ),
LearnMeanPolyLowD which given ε, S, W1 , W2 , Σ
and returns a Σ′ so that
E

N (0,Σ)

[p(X)] −

E [p(X)] ≤ 4 (kΣ − IkF /C + O(log(C)ε) + ξ ,

N (0,Σ′ )

for all p ∈ span(W1 ∪ W2 ) ∩ P2 , with probability 1 − τ .
In particular, this implies:
Corollary 4. Fix ε, τ > 0. Let Σ be so that kΣ − IkF ≤ O(ε log 1/ε). Let G0 be an (ε, δ)good set of samples from N (0, Σ). Let S = {X1 , . . . , Xn } be so that ∆(S, G0 ) ≤ ε. Let V be a
subspace of Rd . Then there is an algorithm LearnCovLowDim which given S, ε, ξ, τ, V runs in
2
time poly(|S|, 2O(dim(V ) ) , log 1/τ ) and returns a Σ′ so that

kΠV Σ − Σ′ ΠV k2 ≤ 4 (kΣ − IkF /C + O(log(C)ε) ,
with probability 1 − τ .

21

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 W1 = W and W2 = ∅
(so that we may take ξ = 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 ∩ P2 , l earn
the mean over every polynomial in the net, and then find a Σ′ 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.
Algorithm 7 Filter if there are many large eigenvalues of the covariance
1: function LowDimCovLearning(S, ε, ξ, τ, W1 , W2 )
2:
Generate a 1/2-cover C for W1 ∩ P2 .
3:
Let τ ′ = 2−|C| τ
4:
for p ∈ C do
5:
Compute mp = LearnMeanChiSquared(S, p, ε, τ ′ ).
6:
Generate a linear constraint cp (Σ′ ): EN (0,Σ′ ) [p(X)] − mp ≤ kΣ − IkF /C + O(log(C))ε.
7:
end for
8:
Generate the convex constraint that EN (0,Σ′ ) [p(X)] − EN (0,Σ′ ) [p(X)] ≤ ξ for all p ∈ W2 .
9:
Using a convex program, return any matrix Σ′ which obeys cp (Σ′ ) for all p ∈ C.
10: end function
Observe that every constraint for each polynomial in W1 is indeed linear in Σ′ , by Lemma 10.
Moreover, the constraint for W2 has an explicit separation oracle, since it induces a norm, and
for any p ∈ W2 , we may explicitly compute EN (0,Σ′ ) [p(X)] − EN (0,Σ′ ) [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 ∈ C. By a union bound, this occurs with probability at least 1 − τ . Thus, in each p ∈ C, we
have that |mp − EX∼N (0,Σ) [p(X)]| ≤ β, where β = kΣ − IkF /C + O(log(C))ε. Let Σ′ be the matrix
we find. By the triangle inequality, we then have that for every p ∈ C, that | EN (0,Σ′ ) [p(X)] −
EN (0,Σ) [p(X)]| ≤ 2β. Hence, by the usual net arguments, we know that for every p ∈ V ∩ P2 ,
|

E

[p(X)] −

N (0,Σ′ )

E [p(X)]| ≤ 4β .

N (0,Σ)

Moreover, by triangle inequality, for every p ∈ W2 , we have EN (0,Σ′ ) [p(X)] − EN (0,Σ′ ) [p(X)] ≤ 2ξ.
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 highdimensions up to error O(ε). 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 ε, δ > 0, and let S0 = (G0 , E0 ) be an ε-corrupted set of samples of size n from
N (0, Σ), where kΣ − IkF ≤ ξ where ξ = O(ε log 1/ε), and where n = poly(d, 1/ε, log 1/δ). Suppose
22

that G0 is (ε, δ)-good with respect to N (0, Σ). Let S ⊆ S0 be a set so that ∆(S, G0 ) ≤ ε. Then, there
exists an algorithm ImproveCov that given S, ξ, ε, fails with probability at most poly(ε, 1/d, δ),
and otherwise outputs one of two possible outcomes:
b so that kΣ
b − ΣkF ≤ kΣ − IkF /2.
(i) A matrix Σ,

(ii) A set S ′ ⊂ S so that ∆(S ′ , G0 ) < ∆(S, G0 ).

Moreover, ImproveCov runs in time poly(d, (1/ε)O(log

4

1/ε) , log 1/δ).

By first applying the algorithm in [DKK+ 16] to produce an initial estimate for Σ, and then
iterating the above algorithm polynomially many times, this immediately yields:
Corollary 5. Fix ε, δ > 0, and let G0 be a set of i.i.d. samples from N (0, Σ), where n =
poly(d, 1/ε, log 1/δ). Let S be so that ∆(S, G0 ) ≤ ε. There is an universal constant C and an
b so that with probability 1 − δ, we have kΣ
b −1/2 ΣΣ
b −1/2 − IkF ≤ Cε. In
algorithm which outputs a Σ
particular, this implies that


b ≤ 2Cε .
dTV N (0, Σ), N (0, Σ)

6.1

Technical Overview

Our strategy for obtaining a high-dimensional estimate for the covariance based on solving lowdimensional subproblems will be substantially more challenging than it was for the unknown mean
case. The natural approach is to take the poly log(1/ε)-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(ε)-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 ⊆ Rd 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 polylogarithmic dimensional subspace V , and on V ⊥ . This does not in general imply that we have
learned the covariance in Frobenius norm. In block form, if we write


ΣV AT
,
Σ=
A ΣV ⊥
23

where here Rd is written as V ⊕ V ⊥ , this implies we have learned ΣV and ΣV ⊥ 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 ⊥ , 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/ε) 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 ε, δ > 0. Let D1 , . . . , Dk , D be a set of distributions where mini dTV (Di , D) = γ.
1/δ
. There is an algorithm Tournament which given oracles for evaluating
Set n = Ω log k+log
ε2
the pdfs of D1 , . . . , Dk along with n independent samples X1 , . . . Xn from D, outputs a Di so that
dTV (Di , D) ≤ 3γ + ε with probability 1 − δ. Moreover, the running time and number of oracle calls
needed is at most O(n2 /ε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 ε. Throughout the paper, we assume the
algorithm knows ε. However, if the algorithm is not given this information, and instead given an
η and asked to return something with error at most O(ε + η), we may simply grid over {η, (1 +
γ)η, (1 + γ)2 η, . . . , 1} (here γ is some arbitrary constant that governs a tradeoff between runtime
and accuracy), run our algorithm with ε 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(ε + (1 + γ)η).
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 ⊗ denote the Kronecker product on matrices. We will make crucial use of the
following definition:
2

Definition 6. For any matrix M ∈ Rd×d , let M ♭ ∈ Rd denote its canonical flattening into a vector
2
2
in Rd , and for any vector v ∈ Rd , let v ♯ denote the unique matrix M ∈ Rd×d so that M ♭ = v.
We will also require the following definition:
2

Definition 7. Let Ssym = {M ♭ ∈ Rd : M is symmetric}.
The following result was proven in [DKK+ 16]:

Theorem 13 (Theorem 4.15 in [DKK+ 16]). Let X ∼ N (0, Σ). Let M be the d2 × d2 matrix given
by M = E[(X ⊗ X)(X ⊗ X)T ]. Then, as an operator on Ssym , we have
   T
.
M = 2Σ⊗2 + Σ♭ Σ♭
24

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 Rd with the Gaussian inner product, defined by hf, gi = EX∼N (0,I) [f (X)g(X)], and we let
kf k22 = hf, f i.
For any symmetric M with kM kF = 1, define the degree-2 polynomial p(x) = √12 (xT M x −
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 kM kF = 1. Let p be its associated polynomial. Then,
we have:
(i) EX∼N (0,I) [p(X)] = 0.
(ii) More generally, for any positive definite matrix Σ, we have EX∼N (0,Σ) [p(X)] = hM, Σ − Ii.
(iii) VarX∼N (0,I) [p(X)] = EX∼N (0,I) [p2 (X)] = hp, pi = 1.
(iv) More generally, for any positive definite matrix Σ, we have
2
T
1
E
[p2 (X)] = M ♭ Σ⊗2 M ♭ +
hΣ − I, M ♭ i .
2
X∼N (0,Σ)
Proof. The first three properties are a straightforward calculation. We show the last one here. By
definition, we have
E

X∼N (0,Σ)

1
E
[(X T M X − tr(M ))2 ]
2 X∼N (0,Σ)
1
E
[(X T M X)2 − 2(X T M X)tr(M ) + tr(M )2 ]
=
2 X∼N (0,Σ)


1
T
♭
2
♭T
=
E
[M (X ⊗ X)(X ⊗ X) M ] − 2hM, Σitr(M ) + tr(M )
2 X∼N (0,Σ)


 T
T
(a) 1
M ♭ 2Σ⊗2 + Σ♭ Σ♭ M ♭ − 2hM, Σitr(M ) + tr(M )2
=
2

T
1 ♭
= M ♭ Σ⊗2 M ♭ +
hM , Σi2 − 2hM, Σitr(M ) + tr(M )2
2
2
1
♭
♭ T ⊗2
=M Σ M +
hΣ − I, M ♭ i ,
2

[p2 (X)] =

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

Observe that Lemma 10(iv) implies that if we take the top eigenvector of the d2 × d2 matrix
Σ⊗2 +
on the linear subspace

1  ♭   ♭ ⊤
M
M
2

V = {M ♭ : M is a symmetric d × d matrix} ,
then the associated polynomial maximizes EX∼N (0,Σ) [p2 (X)], and so we can find these polynomials
efficiently. More generally, if we take any linear subspace of degree two polynomials with associated
25

matrix subspace V ′ , so that V ′ ⊆ V , then the top eigenvector of the same matrix restricted to V ′
allows us to find the polynomial in this subspace which maximizes EX∼N (0,Σ) [p2 (X)] efficiently.
We have the following tail bound for degree-2 polynomials in Gaussian space: We will use ΠV (x)
and Π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 : Rd → R be a degree
√
m polynomial, and let q ≥ 2 be even. Then EX∼N (0,I) [p(X)q ]1/q ≤ ( q − 1)m kpk2 .
We need the following definition:
Definition 8. Let Pk denote the set of even degree-k polynomials over d variables satisfying
VarX∼N (0,1) [p(X)] = 1. Moreover, for any subspace W ⊆ Rd , let Pk (W ) denote the set of even
polynomials over d variables which only depend on the coordinates in W .
b
Then by the arguments above, we have that for any two matrices Σ, Σ,
!
b F = sup
kΣ − Σk
E
[p(X)] −
E
[p(X)] .
p∈P2

X∼N (0,Σ)

b
X∼N (0,Σ)

In particular, by Lemma 2, this implies that when kΣ − Ik2 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 ≥ 0 so that for all t ≥ C,
we have


Pr [|p(X) − E [p(X)]| ≥ tkpk2 ] ≤ exp −At1/2 .
N (0,I)

N (0,I)

Proof. Hypercontractivity in particular implies the following moment bound: for all q ≥ 2 even,
we have
E [(p(X) − E [p(X)])q ] ≤ (q − 1)qm/2 kp(X) − E [p(X)]kq2 .
N (0,I)

N (0,I)

N (0,I)

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 = (a1 , . . . , ad ),
we let Ha (x) : Rd → R be the Hermite
Pd polynomial associated with multi-index a. It is well-known
that the degree of Ha is |a| =
i=1 ai , and moreover, hHa , Hb i = δa,b . In particular, for any
r ≥ 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 : Rd → R with degree r, we may write it uniquely as
X
p(x) =
ca (p)Ha (x), where ca (p) = hp, Ha i .
|a|≤r

We define the kth harmonic component of p to be
X
p[k](x) =
ca (p)Ha (x) ,
|a|=k

and we say p is harmonic of degree k if it equals its kth part.
26

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 ε, δ > 0. Assume kΣ − IkF ≤ ξ, where ξ = O(ε log 1/ε). Suppose that G0
b = ES [XX T ].
is (ε, δ)-good with respect to N (0, Σ). Let S be a set so that ∆(S, G0 ) ≤ ε. Let Σ
Then there is an algorithm FilterCovManyDeg2Eig and a universal constant C such that the
following guarantee holds:
b − I has more than O(log 1/ε) eigenvalues larger than Cξ, then the algorithm outputs a
1. If Σ
′
S so that ∆(S ′ , G0 ) < ∆(S, G0 ).

2. Otherwise, the algorithm outputs “OK”, and outputs an orthonormal basis v1 , . . . , vk for the
b − I with eigenvalue larger than Cξ.
subspace V of vectors spanned by all eigenvectors of Σ

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(ε log2 1/ε) (restricted to a certain subspace of polynomials), but there can only be at most O(poly log 1/ε) 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 ε, δ > 0. Assume kΣ − IkF ≤
ξ, where ξ = O(ε log 1/ε). Let C be the universal constant in FilterCovManyDeg2Eig. Let
W ⊆ Rd be a subspace, so that for all v ∈ W with kvk2 = 1, we have v T ES [XX T ]v ≤ 1 + Cξ.
Suppose that G0 is (ε, δ)-good with respect to N (0, Σ). Let S be a set so that ∆(S, G0 ) ≤ ε. Let
k = O(log4 1/ε). Then there is an algorithm FilterCovManyDeg4Eig and universal constants
C1 , C2 such that the following guarantee holds:
1. If there exist p1 , . . . , pk ∈ P2 (W ) so that hpj , pℓ i = δjℓ for all j, ℓ, and so that ES [p2j (Y )] − 1 ≥
C1 ε for all j, then the algorithm outputs an S ′ so that ∆(S ′ , G0 ) < ∆(S, G0 ).
2. Otherwise, the algorithm outputs “OK”, and outputs an orthonormal basis p1 , . . . , pk′ for a
subspace V of degree-2 polynomials in P2 (W ) with k′ ≤ k so that for all p ∈ V ⊥ ∩ P2 , we
have ES [p2 (X)] − 1 ≤ C2 ε.
Moreover, FilterCovManyEig runs in time poly(d, 1/ε, log 1/δ).
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 ⊥ .
27

Algorithm 8 Filter if there are many large eigenvalues of the fourth moment tensor
1: function FilterCovManyEig(S, ε, ξ, δ, W )
b = ES [XX T ]
2:
Let Σ
3:
Let C1 , C2 , C3 be some universal constants sufficiently large
4:
Let A be the constant in Corollary 6
5:
Let B be the constant in Claim 20
6:
Let m = 0
7:
Let k = O(log4 1/ε)
8:
while there exists p ∈ P2 (W ) so that p ∈ V ⊥ and ES [p2 (X)] − 1 > C1 ξ do
9:
Let Vm+1 = span(Vm ∪ p)
10:
Let m ← m + 1
11:
end while
12:
Let p1 , . . . , pm be an orthonormal basis for Vm
13:
if m ≥ k then
[4]
14:
Let qi = (p2i ) be the 4th harmonic component of p2i
15:
Let ri = p2i − qi be the degree-2 component of p2i
P
16:
Let Q(x) = ki=1 qi
17:
Find a T so that either:
√
• T > C3 d2 k log(|S|) and p(X) > T for at least one x ∈ S ′ , OR
√
• T > 4A2 C2 B k log2 (1/ε) and
√
Pr [Q(X) > T ] > exp(−A(T /4B k)1/2 ) + ε2 /(d log(|S|/δ))2 .
X∈u S

18:
19:
20:
21:
22:

return the set S ′ = {X ∈ S : Q(X) ≤ T }
else
return “OK”, and output p1 , . . . , pm
end if
end function

Theorem 17. Let 1 > ξ > η > ε > 0, and let τ > 0. Let Σ so that kΣ − IkF ≤ ξ. Suppose that
Rd is written as V ⊕ W for orthogonal subspaces V and W with dim(V ) = O(log(1/ε)). Suppose
furthermore that


ΣV AT
,
Σ=
A ΣW
with kΣV − IV kF , kΣW − IW kF = O(η). Let S0 = (G0 , E0 ) be an ε-corrupted set of samples from
N (0, Σ), and let S ⊆ S0 with ∆(S, G) ≤ O(ε) of size poly(d, 1/η, log 1/δ).
Then, there exists a universal constant C5 and an algorithm Stitching that given V, W, ξ, η, ε, τ
and S runs in polynomial time and with probability at least 1 − τ returns a matrix Σ0 with kΣ0 −
ΣkF = C5 η + O(ξ 2 ).
In the latter, we will show the algorithm works when τ = 2/3. As usual the probability of
success can be boosted by repeating it independently.5 The basic idea of the proof is as follows.
Since we already know good approximations to ΣV and ΣW , it suffices to find an approximation to
5
Observe the only randomness at this point is in the random choices made by the algorithm. Thus, one can just
(j)
(ℓ)
(1)
run this algorithm O(log 1/δ) times to obtain Σ0 , . . . , Σ0 and find any Σ0 which is O(η + ξ 2 ) close to at least a
2/3 + o(1) fraction of the other outputs.

28

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 → 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 xV 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 ε. 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(ε). 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 ℓ2 error between Av and
its approximation.
Our algorithm is given in Algorithm 9:
Algorithm 9 Stitching the two subspaces together
1: function Stitching(V, W, δ, ε, τ, S)
2:
Given a vector x, let xV and xW be the projections onto V and W , respectively.
3:
Let C be a sufficiently large constant (where C may depend on the constants in the big-O
terms in the guarantee that dim(V ) = O(log(1/ε))).
4:
Generate a set V = {v1 , . . . , vm } of (n/ε)C independent random samples from N (0, 2IV ).
5:
for v ∈ V do
6:
For each sample x ∈ S, add xW to a new set T independently with probability
exp(−kxV − vk2 /2) .
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:

17:

Treat T as a collection of independent samples from a noisy Gaussian with covariance
matrix IW + O(η).
Set av equal to 0 if T did not contain enough samples for our algorithm or if kµ̃k2 >
C log(1/ε).
for ε ∈ {1, 1/2, 1/4, 1/8, . . . , η} do
Let µ
e be the output of RecoverMeanNoisy(T, ε, (ε/n)2C , o(1), O(η)).
end for
Run Tournament with the output hypotheses.
Set av = µ̃, where µ
e is the winning hypothesis.
end for
Use linear programming to find the dim(W ) × dim(V )-matrix B that minimizes the convex
function Ev∈u S [|av − Bv|].
return


IV 2B T
.
Σ0 =
2B IW
end function
Before we prove Theorem 17, we will need the following definition.

29

Definition
9. A function f : Rd → R is a positive measure if f ≥ 0 and
R
|f |1 = f , and we say f1 ≤ f2 if f1 (x) ≤ f2 (x) pointwise.

R

f ≤ 1. We write

Proof of Theorem 17. Throughout this proof, let G = N (0, Σ). 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(−kxV − vk2 /2). Therefore, given a distribution D, we let the positive measure
Dv be what is obtained by sampling from D and accepting a sample x only with probability
exp(−kxV − vk2 /2). We need to understand the distributions Gv and Xv .
Note that the pdf of G is given by
p
(2π)−n/2 det(Σ−1 ) exp(−xΣ−1 x/2)dx.

Therefore, the density of Gv is
p
(2π)−n/2 det(Σ−1 ) exp(−xΣ−1 x/2 − kxV − vk22 /2)dx
p
=(2π)−n/2 det(Σ−1 ) exp(−x(Σ−1 + IV )x/2 + x · v − kvk22 /2)dx.

Note that letting µv := (Σ−1 + IV )−1 v = v/2 + O(δkvk2 ) that this equals
p
=(2π)−n/2 det(Σ−1 ) exp(−(x − µv )(Σ−1 + IV )(x − µv )/2 − kvk22 /2 + µv · v/2)dx
p
=(2π)−n/2 det(Σ−1 ) exp(−(x − µv )(Σ−1 + IV )(x − µv )/2 − kvk22 (1 + O(δ))/4)dx.

Note that this is a Gaussian with mean µv weighted by
s
det(Σ−1 )
exp(−kvk22 (1 + O(δ))/4) = Θ(2− dim(V )/2 exp(−kvk22 /4)) ,
det(Σ−1 + IV )

so long as kvk2 ≪ δ−1/2 . Therefore, if this condition holds, a random sample from G is accepted
by this procedure with probability Θ(2− dim(V ) exp(−kvk22 /4).
We also need to understand the fraction of samples from Xv 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
εv := 2n/2 exp(kvk22 /4)(|Ev |1 + |Lv |1 log(1/ε)) ,
that is approximately the fraction of samples from Xv that are errors (where subtractive errors are
weighted more heavily).
p
Note that so long as kvk2 < C log(1/ε) that so long as |Lv |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 kµ̃ − πW (µv )k2 ≤ O(εv + η). Therefore,
unless kµv k2 = Ω(C log(1/ε)) (which can only happen in kvk2 = Ω(C log(1/ε))), we have with high
probability that kav − πW (µv )k2 = O(εv + η). In fact, we assume that this holds for all v ∈ S with
kvk2 ≪ (C log(1/ε)).
In order to show that this is generally a good approximation, we need to know that εv is not
too large on average. In particular, we show:
Lemma 11. If v ∼ N (0, 2IV ), then Ev [εv ] = O(ε).

30

Proof. Letting F = E + log(1/ε)L, we have that
Z
− dim(V )/2
E[εv ] = (2π)
exp(kvk22 /4) exp(−kvk22 /4)dFv dv
v
Z
Z
− dim(V )/2
2
= (2π)
exp(−kx − vk /2)dF (x)dv = dF = O(ε).
Note that πW (µv ) = M v, where M = πW (Σ−1 + IV )−1 πVT . Let B be any dim(W ) × dim(V )
matrix. Note that
E
[kBv − M vk22 ] = 2kB − M k2F .
v∼N (0,2IV )

Therefore, Ev∼N (0,2IV ) [kBv − M vk2 ] = O(kB − M kF ). Since kBv − M vk22 is a degree-2 polynomial
in v, by Corollary 4 in [Kan12],
#
"
r
E
[kBv − M vk42 ]/2 > 1/81.
Pr
[kBv − M vk2 > kB − M kF /2] >
Pr
kBv − M vk22 >
v∼N (0,2IV )

v∼N (0,2IV )

v∼N (0,2IV )

Combining these we find that
E

[kBv − M vk2 ] = Θ(kB − M kF ) .

v∼N (0,2IV )

Next we show that our choice of S derandomizes this result.
Lemma 12. With high probability over the choice of S, we have that
E [kBv − M vk2 ] = Θ(kB − M kF ).

v∈u S

for all dim(W ) × dim(V ) matrices B.
Proof. This is equivalent to showing that for all matrices U = B − M it holds
E [kU vk] = Θ(kU kF ).

v∈u S

By the standard scaling laws, it suffices to show this only for U with kU kF = 1.
We also note that it suffices to show this only for U in an ε-net for all such matrices. This is
because if kU − U ′ kF < ε and if U ′ satisfies the desired condition, then
E [kU vk2 ] = E [kU ′ vk2 ] + O( E [k(U − U ′ )vk2 ]) = Θ(1) + O(ε) E [kvk2 ] ,

v∈u S

v∈u S

v∈u S

v∈u S

and with high probability Ev∈u S [kvk2 ] = O(log(1/ε)).
Note that such nets exist with size exp(poly(n/ε)). Therefore, it suffices to show that this
condition holds for each such U with probability exp(−(n/ε)Ω(C) ).
As noted above, it suffices that
E [kU vk2 ] = O( E [kU vk2 ]) = O(1)

v∈u S

v∈u S

and that
Pr [kU vk2 > 1/4] ≥ Pr [kU vk2 ≥ 1/4] − 1/100.

v∈u S

v∈u S

The first follows because if S = {v1 , . . . , vm } then Ev∈u S [kU vk2 ] is a degree-2 polynomial in the
√
vi with mean O(1) and variance O(1/ m), so by standard concentration results, is O(1) with 1 −
exp(−Ω(|S|)) probability. The latter follows from standard concentration bounds. This completes
the proof.
31

We also note that the random choice of S has another nice property:
Lemma 13. With S and av as above, with high probability we have that
E [kav − M vk2 ] = O(η) .

v∈u S

Proof. Let S = {v1 , . . . , vm }. Then kf (vi ) − M vi k2 are independent random variables with mean
E[O(εv + η)] = O(η) and variance at most
E[kav k22 + kM vk22 ] ≤ O(log2 (1/ε)) + E[kvk22 ] = O(log2 (1/ε)).
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 av in
Step 5 we have that for all dim(W ) × dim(V )-matrices B that
E [kBv − M vk2 ] = Θ(kB − M kF ),

v∈u S

and have that
E [kav − M vk2 ] = O(η).

v∈u S

Combining these statements, we have that for all dim(W ) × dim(V )-matrices B that
E [kav − Bk2 ] = Θ(kB − M kF ) + O(η).

v∈u S

Note that by taking B = M , this quantity is O(η). Therefore, the B found in Step 15 satisfies
O(ε) = E [kav − Bk2 ] = Θ(kB − M kF ) + O(η) ,
v∈u S

and therefore, kB − M kF = O(η).
The rest of the proof is a simple computation of the matrices involved. In particular, recall that

 

IV AT
ΣV AT
+ O(η) ,
=
Σ=
A IW
A ΣW
where the O(η) denotes a matrix with Frobenius norm O(η) and where kAkF = O(δ). It is easy to
see that


I
−AT
+ O(η + δ2 ).
Σ−1 = V
−A IW

Therefore,

Σ

−1

Hence,
(Σ

−1


2IV −AT
+ O(η + δ2 ).
+ IV =
−A IW


−1

+ IV )




IV /2 AT /2
=
+ O(η + δ2 ).
A/2
IW

Therefore, M = A/2 + O(η + δ2 ). Therefore, A = 2M + O(η + δ2 ) = 2B + O(η + δ2 ). And finally,
we conclude that


IV 2B T
+ O(η + δ2 ) = Σ0 + O(η + δ2 ).
Σ=
2B IW

32

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 N (0, Σ), 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 η, δ > 0. Let X1 , . . . , Xn be independent samples from N (µ, I), where n =
poly(d, 1/η, log 1/δ). Then, S = {X1 , . . . , Xn } is (η, δ)-good with respect to N (µ, I) with probability
at least 1 − δ.
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 ε, δ > 0. Assume kΣ − IkF ≤ ξ, where ξ = O(ε log 1/ε). Suppose that G0 is
(ε, δ)-good with respect to N (0, Σ), and let S ⊆ S0 be a set so that ∆(S, G0 ) ≤ ε. There is some
√
absolute constant C5 so that if p ∈ P2 is a polynomial so that ES [p(X)] − EN (0,Σ) [p(X)] > C5 ξε,
then ES [p2 (X)] − 1 > C1 ξ.
We defer the proof of this lemma to the Appendix.
We are now ready to present the full algorithm as Algorithm 10.
Algorithm 10 Filter if there are many large eigenvalues of the covariance
1: function ImproveCov(S, ξ, ε, δ)
2:
Let C be the universal constant in FilterCovManyDeg2Eig
3:
Let τ = poly(ε, 1/d, δ).
4:
Run FilterCovManyDeg2Eig(S, ε, ξ)
5:
if FilterCovManyDeg2Eig outputs S ′ then
6:
return S ′
7:
else
8:
Let V be the subspace returned by FilterCovManyDeg2Eig
9:
Let W = V ⊥ .
10:
Run FilterCovManyDeg4Eig(S, ε, ξ, δ, W )
11:
if FilterCovManyDeg4Eig outputs S ′ then
12:
return S ′
13:
else
14:
Let U1 be the subspace of degree 2 polynomials over W it returns
15:
Let U2 be the perpendicular subspace of degree 2 polynomials over W
b = ES [XX T ]
16:
Let Σ
17:
Let ΣV = LearnCovLowDim(S, ε, ξ, τ, V )
b
18:
Let ΣW = LearnMeanPolyLowD(S, ε, ξ, τ, U1 , U2 , Σ)
′
19:
Take poly(n, 1/ε) fresh ε-corrupted samples S
20:
return Stitching(V, W, ΣV , ΣW , ξ, ε, S ′ ).
21:
end if
22:
end if
23: end function

33

Proof of Theorem 11. Condition on the events that neither LearnCovLowDim nor Stitching
fail. This happens with probability at least poly(ε, 1/d, δ). Observe that if we pass the “if”
statement in Line 5, then by the guarantee of FilterCovManyDeg2Eig this is indeed an S ′
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 ′ satisfying the desired properties. Otherwise, by Lemma 15, we know
√ that for all polynomials
p ∈ P2 over W orthogonal to U1 , we have | EN (0,Σ) − ES [XX T ]| ≤ C5 ξε. Thus, ΣW satisfies the
conditions needed by Stitching.
By Corollary 4, we know that Σ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/ε 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 pi 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 EE [Q(X)] to be larger than Q(X) for all but a poly(ε) fraction of X ∈ G.
We then remove all points X ∈ 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 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 poly(ε) fraction of X ∈ G, we are forced to choose k = Ω(log 1/ε),
which means that we need to reduce to k = Ω(log2 1/ε) 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 N (µ, Σ) to
learning the covariance of N (0, Σ), by at most doubling error, a trick previously used in [DKK+ 16]
and [LRV16]. Given an
√ ε-corrupted set of samples X1 , . . . , X2n of size 2n from N (µ, Σ), we may
let Yi = (Xi − Xn+i )/ 2. Then we see that if Xi and Xn+i are uncorrupted, then Yi ∼ N (0, Σ).
Moreover, at most 2εn of the Yi can be corrupted, since there are at most 2εn corrupted Xi .
Therefore, by doubling the error rate, we may assume that µ = 0. We may then apply the algorithm
b so that with high probability, we have kΣ
b −1/2 ΣΣ
b −1/2 − IkF ≤ O(ε)
in Corollary 5 to obtain a Σ
4
O(log
1/ε)
with polynomially many samples, and in poly(d, (1/ε)
) time.
b −1/2 X ′ .
We may then take an additional set of ε-corrupted samples {Xi′ , . . . , Xn′ }, and let Yi′ = Σ
i
′
′
b we have that if X is uncorrupted, then Y ∼ N (0, Σ)
e where
Then, by our guarantee on Σ,
i
i
e − IkF ≤ O(ε). We then run RecoverMeanNoisy with the Y ′ to obtain a µ
kΣ
b so that
i
−1/2
−1/2
b
e
b
b
kb
µ−Σ
µk2 ≤ O(ε). This guarantees that dTV (N (b
µ, Σ), N (Σ
µ, Σ)) ≤ O(ε), which in turn
b N (µ, Σ)) ≤ O(ε), as claimed.
implies that dTV (N (b
µ, Σ),
34

Therefore, we have shown:
Theorem 18. Fix ε, δ > 0. Given an ε-corrupted set of samples S from N (µ, Σ), where n =
poly(d, 1/ε, log 1/δ), there is an algorithm RecoverGaussian which takes as input S, ε, δ, and
b so that
outputs a µ
b, Σ
b ≤ O(ε) .
dTV (N (µ, Σ), N (b
µ, Σ))

Moreover, the algorithm runs in time poly(d, (1/ε)O(log

4

1/ε) , log 1/δ).

References
[BDLS17] S. Balakrishnan, S. S. Du, J. Li, and A. Singh. Computationally efficient robust sparse
estimation in high dimensions. In Proceedings of the 30th Annual Conference on Learning Theory, COLT ’17, 2017.
[Ber06]

T. Bernholt. Robust estimators are hard to compute. Technical report, University of
Dortmund, Germany, 2006.

[BNJ03]

D. M Blei, A. Y. Ng, and M. I. Jordan. Latent dirichlet allocation. Journal of machine
Learning research, 3(Jan):993–1022, 2003.

[BW41]

L. M. Blumenthal and G. E. Wahlin. On the spherical surface of smallest radius enclosing a bounded subset of n-dimensional euclidean space. Bull. Amer. Math. Soc.,
47:771–777, 1941.

[CSV17]

M. Charikar, J. Steinhardt, and G. Valiant. Learning from untrusted data. In Proc.
49th Annual ACM Symposium on Theory of Computing (STOC), pages 47–60. ACM
Press, 2017.

[CW01]

A. Carbery and J. Wright. Distributional and Lq norm inequalities for polynomials over
convex bodies in Rn . Mathematical Research Letters, 8(3):233–248, 2001.

[DBS17]

S. S. Du, S. Balakrishnan, and A. Singh. Computationally efficient robust estimation
of sparse functionals. CoRR, abs/1702.07709, 2017.

[DDS12]

C. Daskalakis, I. Diakonikolas, and R.A. Servedio. Learning Poisson Binomial Distributions. In Proceedings of the 44th Symposium on Theory of Computing, pages 709–728,
2012.

[DDS15]

A. De, I. Diakonikolas, and R. Servedio. Learning from satisfying assignments. In
Proceedings of the 26th Annual ACM-SIAM Symposium on Discrete Algorithms, SODA
2015, pages 478–497, 2015.

[DK14]

C. Daskalakis and G. Kamath. Faster and sample near-optimal algorithms for proper
learning mixtures of gaussians. In Proceedings of The 27th Conference on Learning
Theory, COLT 2014, pages 1183–1213, 2014.

[DKK+ 16] I. Diakonikolas, G. Kamath, D. M. Kane, J. Li, A. Moitra, and A. Stewart. Robust
estimators in high dimensions without the computational intractability. In Proceedings
of FOCS’16, 2016.

35

[DKK+ 17] I. Diakonikolas, G. Kamath, D. M. Kane, J. Li, A. Moitra, and A. Stewart. Being
robust (in high dimensions) can be practical. In Proceedings of the 34th International
Conference on Machine Learning, ICML ’17, pages 999–1008. JMLR, Inc., 2017.
[DKS16]

I. Diakonikolas, D. M. Kane, and A. Stewart. Robust learning of fixed-structure bayesian
networks. CoRR, abs/1606.07384, 2016.

[DKS17]

I. Diakonikolas, D. M. Kane, and A. Stewart. Statistical query lower bounds for robust
estimation of high-dimensional gaussians and gaussian mixtures. In Proceedings of
the 58th Annual IEEE Symposium on Foundations of Computer Science, FOCS ’17,
Washington, DC, USA, 2017. IEEE Computer Society.

[DL01]

L. Devroye and G. Lugosi. Combinatorial methods in density estimation. Springer Series
in Statistics, Springer, 2001.

[GCSR14] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin. Bayesian data analysis, volume 2.
Chapman & Hall/CRC Boca Raton, FL, USA, 2014.
[GLS88]

M. Grötschel, L. Lovász, and A. Schrijver. Geometric Algorithms and Combinatorial
Optimization, volume 2. Springer, 1988.

[Goo15]

R. Goodell. Personal communication, 2015.

[Ham01]

F. Hampel. Robust statistics: A brief introduction and overview. In First International
Symposium on Robust Statistics and Fuzzy Techniques in Geodesy and GIS. A. Carosio,
H. Kutterer (editors), Swiss Federal Institute of Technology Zurich (ETH), Institute of
Geodesy and Photogrammetry, IGP-Bericht, number 295, pages 13–17, 2001.

[HM13]

M. Hardt and A. Moitra. Algorithms and hardness for robust subspace recovery. In
COLT 2013, pages 354–375, 2013.

[HR09]

P.J. Huber and E. M. Ronchetti. Robust statistics. Wiley New York, 2009.

[HRRS86] F. R. Hampel, E. M. Ronchetti, P. J. Rousseeuw, and W. A. Stahel. Robust statistics.
The approach based on influence functions. Wiley New York, 1986.
[Hub64]

P. J. Huber. Robust estimation of a location parameter. The Annals of Mathematical
Statistics, 35(1):73–101, 1964.

[Jun01]

H. W. E. Jung. Über die kleinste kugel, die eine räumliche figur einschliesst. J. Reine
Angew. Math, 123:241–257, 1901.

[Kan12]

D. M. Kane. A structure theorem for poorly anticoncentrated gaussian chaoses and
applications to the study of polynomial threshold functions. In FOCS, pages 91–100,
2012.

[Kla07]

B. Klartag. A central limit theorem for convex sets.
168(1):91–131, 2007.

[Li17]

J. Li. Robust sparse estimation tasks in high dimensions. CoRR, abs/1702.05860, 2017.

[LM00]

B. Laurent and P. Massart. Adaptive estimation of a quadratic functional by model
selection. Ann. Statist., 28(5):1302–1338, 2000.
36

Inventiones mathematicae,

[LRV16]

K. A. Lai, A. B. Rao, and S. Vempala. Agnostic estimation of mean and covariance. In
Proceedings of FOCS’16, 2016.

[NJB+ 08] J. Novembre, T. Johnson, K. Bryc, Z. Kutalik, A. R. Boyko, A. Auton, A. Indap, K. S.
King, S. Bergmann, M. R. Nelson, et al. Genes mirror geography within europe. Nature,
456(7218):98–101, 2008.
[O’D14]

R. O’Donnell. Analysis of Boolean Functions. Cambridge University Press, 2014.

[OF96]

B. A. Olshausen and D. J. Field. Emergence of simple-cell receptive field properties by
learning a sparse code for natural images. Nature, 381(6583):607, 1996.

[RL05]

P. J. Rousseeuw and A. M. Leroy. Robust regression and outlier detection, volume 589.
John Wiley & Sons, 2005.

[Rou85]

P. Rousseeuw. Multivariate estimation with high breakdown point. Mathematical Statistics and Applications, pages 283–297, 1985.

[SCV18]

J. Steinhardt, M. Charikar, and G. Valiant. Resilience: A criterion for learning in
the presence of arbitrary outliers. In Proc. of the 9th Conference on Innovations in
Theoretical Computer Science, 2018. to appear.

[Tuk75]

J. W. Tukey. Mathematics and picturing of data. In Proceedings of ICM, volume 6,
pages 523–531, 1975.

[Ver10]

R. Vershynin. Introduction to the non-asymptotic analysis of random matrices, 2010.

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
ε
. We note that our lower
from an ε-corrupted set of samples to total variation distance less than 1−ε
bound applies in more general settings than we consider in this paper in couple of ways:
1. Our construction is univariate;
2. Our construction works for any pair of distributions which are
Gaussian distributions; and

ε
1−ε -close,

and not just for

3. Our construction holds for Huber’s ε-contamination model, which is weaker than the noise
model studied in this paper.
In Huber’s ε-contamination model, data is drawn from a mixture distribution (1 − ε)P + εQ,
where P is some distribution that we wish to estimate and Q is arbitrary (in particular, it might
ε
-close can be made indistindepend on P ). We will show that any two distributions which are 1−ε
guishable under this contamination model.
ε
Lemma 16. Let p1 and p2 be two distributions such that dTV (p1 , p2 ) = 1−ε
. Then there exist
Huber ε-contaminations p1′ and p2′ of p1 and p2 respectively, such that p1′ = p2′ . Therefore, no
ε
in Huber’s ε-contamination model.
algorithm can learn a distribution p up to accuracy < 1−ε

37

Proof. Since dTV (p1 , p2 ) =

ε
1−ε ,

we can write



ε
ε
p1 = 1 −
q1 ,
pc +
1−ε
1−ε


ε
ε
p2 = 1 −
q2 ,
pc +
1−ε
1−ε
where pc , q1 , q2 are distributions. Let p1′ be the Huber ε-contamination of p1 with q2 , p2′ be the
Huber ε-contamination of p2 with q1 . In other words,
p1′ = (1 − 2ε) pc + εq1 + εq2 ,

p2′ = (1 − 2ε) pc + εq2 + εq1 ,
and thus p1′ = p2′ as desired.

We note that a similar lower bound holds when p1 and p2 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.pNo algorithm
can output an estimate for the mean of a unit-variance Gaussian at

π
accuracy <
2 − o(1) ε with probability > 1/2 in Huber’s ε-contamination model. Consequently,

by Lemma 1, no algorithm can learn a Gaussian to total variation distance < 12 − o(1) ε in Huber’s
ε-contamination model.
Proof. We will consider the distributions p1 = N (−α, 1) and p2 = N (α, 1), where α is to be
specified later. We note that if p1 and p2 can be ε-corrupted into the same distribution, then the
best estimate for the mean is to output 0 (by symmetry), and the result holds.
2 (x)}
and
We will show that p1 can be ε-corrupted into a distribution f , where f (x) = max{p1 (x),p
η
η is a normalizing constant (the case of p2 follows similarly). In other words, f can be written as
(1 − ε)p1 + εq1 for some distribution q1 . Since q1 is a distribution, it is non-negative, and thus we
2 (x)}
require that (1 − ε)p1 (x) ≤ f (x) for all x ∈ R. Note that f (x) = max{p1 (x),p
≥ p1η(x) .
η
We can compute η as follows:
Z ∞
Z ∞
p2 (x)dx
max {p1 (x), p2 (x)} dx = 2
η=
0
−∞
r





α
α
1 1
2
+ erf √
α − O(α3 )
= 1 + erf √
=1+
=2
2 2
π
2
2
Thus,
p (x)
p1 (x)
q 1
=
f (x) ≥
≥
η
1 + π2 α − O(α3 )

1−

r

!
2
α + O(α2 ) p1 (x).
π

q
Again, since we require that (1 − ε)p1 (x) ≤ f (x), it suffices that 1 − ε ≤ 1 − π2 α + O(α2 ) and thus
q

pπ
that π2 α − O(α2 ) ≤ ε. This corresponds to α ≤
2 − o(1) ε, 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 εN samples and add εN samples. In this
38

infinite sample regime, this means that a distribution can be corrupted to any distribution which is
ε-far in total
√ variation
 distance. By noting that N (−α, 1) and N (α, 1) are both ε-far from N (0, 1)
for α =
2π + o(1) ε, 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
2π + o(1) ε with probability > 1/2 in the subtractive adversary model. Consequently,
accuracy <
by Lemma 1, no algorithm can learn a Gaussian to total variation distance < ε 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(ε log1/2 (1/ε)) with probability > 1/2 in the Huber’s ε-contamination model.
Proof. We start with the distribution q = N (0, 1). We construct p1 by truncating the right tail of
q at the point xr = c1 log1/2 (1/ε) for some constant c1 , and rescaling the rest of the distribution
appropriately. Observe that, for an appropriate choice of c1 :


1
;
• p1 is sub-Gaussian with a constant of O 1−ε
• The mean of p1 is −c2 ε log1/2 (1/ε) for some constant c2 ;

• p1 can be corrupted in Huber’s ε-contamination model to be q.
We can similarly consider p2 , which is constructed by truncating the left tail of q at the point
xl = −c1 log1/2 (1/ε). Since p1 and p2 are indistinguishable when they are both ε-corrupted to q,
and the mean of all three distributions are separated by ≥ c2 ε log1/2 (1/ε), the lemma follows.

B
B.1

Omitted Proofs from Section 2
Proof of Lemma 1

Proof. Observe that by rotational and translational invariance, it suffices to consider the problem
when µ1 = −εe1 /2 and µ2 = εe1 /2, where e1 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:
Z ∞
1
1
2
2
dTV (N (µ1 , I), N (µ2 , I)) = · √
e−(x−ε/2) /2 − e−(x+ε/2) /2 dx .
2
2π −∞
2 /2

The value of the function f (x) = e−(x−ε/2)
when x > 0, hence this integral becomes

2 /2

− e−(x+ε/2)

is negative when x < 0 and positive

Z ∞
1
2
2
e−(x−ε/2) /2 − e−(x+ε/2) /2 dx
dTV (N (µ1 , I), N (µ2 , I)) = √
2π 0
= F (ε/2) − F (−ε/2) ,

39

Rx
2
where F (x) = √12π −∞ e−t /2 dt is the CDF of the standard normal Gaussian. By Taylor’s theorem,
and since F ′′ (x) is bounded when x ∈ [−1, 1], we have
F (ε/2) − F (−ε/2) = F ′ (−ε/2)ε + O(ε3 )
1
2
= √ e−(ε/2) /2 ε + O(ε3 )
2π


1
= √ + o(1) ε ,
2π
which proves the claim.

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 kΠV ′ µ − µ̃k2 < 2ε when the
parameters α, γ are chosen appropriately.
Claim 4. EN (µ,I) [p(X)] = kΠV ′ µ − µ̃k22 .
Proof. Letting v1 , . . . , vC1 β log(1/ε) be an orthonormal basis of V ′ , we have
C1 β log(1/ε)

E [p(X)] =

N (µ,I)

X
i=1

C1 β log(1/ε)

E

N (viT µ,1)

2

[hvi , X − µ̃i − 1] =

X
i=1

(viT (µ − µ̃))2 = kΠV ′ µ − µ̃k22 .

Claim 5. For some absolute constant c0 and all t ≥ kΠV ′ µ − µ̃k22 ,



(t − kΠV ′ µ − µ̃k22 )2
2
Pr (p(X) > t) ≤ 2 exp −c0 · min
, t − kΠV ′ µ − µ̃k2
.
C1 β log(1/ε)
N (µ,I)
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
C1 β log(1/ε), respectively.
Claim 6. ES [p(X)] ≥ C1 ε log(1/ε).
Proof. Recall that µ̂ is the empirical mean of the point set.
E[p(X)] = E[kΠV ′ X − µ̃k22 ] − dim(V ′ )
S

S

≥ E[kΠV ′ X − ΠV ′ µ̂k22 ] − dim(V ′ )
S


1
≥ 1 + ε dim(V ′ ) − dim(V ′ ) = C1 ε log(1/ε) .
β
Claim 7. EG0 [p(X)] ≤ kµ − µ̃k22 + O(γ)ε.
40

Proof.
E [p(X)] ≤

G0

≤
≤

Z

0

Z

∞

Pr[p(X) ≥ t]dt
G0

O(d log(d/εδ))

Pr[p(X) ≥ t]dt
G0

0

Z

O(d log(d/εδ))

Pr [p(X) ≥ t]dt + O(γ)ε

N (µ,I)

0

≤ kµ − µ̃k22 + O(γ)ε

The inequalities follow from (γε, δ)-goodness and Claim 4.
Claim 8. φ(S, G0 ) EL [p(X)] = (C1 β + O(γ) + o(1))ε.
Proof.
φ(S, G0 ) E[p(X)] =
L

Z

∞

φ(S, G0 ) Pr[p(X) > t]dt
L

0

=

Z

O(d log(|G|/δ))

φ(S, G0 ) Pr[p(X) > t]dt

=

Z

(4)

L

0
C1 β log(1/ε)+4ε2

φ(S, G0 ) Pr[p(X) > t]dt

0

+

Z

L

O(d log(|G|/δ))
C1 β log(1/ε)+4ε2

φ(S, G0 ) Pr[p(X) > t]dt
L

2

≤ (C1 β log(1/ε) + 4ε )φ(S, G0 ) + 2

Z

O(d log(|G|/δ))

Pr[p(X) > t]dt

C1 β log(1/ε)+4ε2 G0

≤ (C1 β log(1/ε) + 4ε2 )φ(S, G0 )
Z O(d log(|G|/δ))
Pr [p(X) > t] +
+2
C1 β log(1/ε)+4ε2 N (µ,I)
2

γε
dt
d log(|G|/δ)

≤ (C1 β log(1/ε) + 4ε )φ(S, G0 ) + O(γε) + 8εc0 C1 β

≤ (C1 β + O(γ) + o(1))ε.

(5)
(6)
(7)

(4) and (5) follow from G0 being (γε, δ)-good, (6) is from Claim 5, and (7) is because
φ(S, G0 ) log(1/φ(S, G0 )) ≤ (1 + o(1))ε .
This gives:
Claim 9. ψ EE [p(X)] ≥ C1 ε log 1/ε − (C1 β + 1 + O(γ) + o(1))ε.
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 ) ≥ C1 β log(1/ε). Then there is a T satisfying the conditions in the
algorithm.
41

Proof. Suppose not. Then, we have
ψ E[p(X)] ≤
E

Z

2C3 log(1/ε)/c0

0

ψ Pr[p(X) ≥ t]dt +
E

≤ 2(C3 /c0 )ψ log 1/ε +
≤ 2(C3 /c0 )ψ log 1/ε +
≤ (2C3 /c0 )ε log 1/ε +

Z

Z

Z

C2 d log |S|/δ

Z

2(C3 /c0 ) log(1/ε)
C2 d log |S|/δ

∞
2C3 log(1/ε)/c0

ψ Pr[p(X) ≥ t]dt
E

ψ Pr[p(X) ≥ t]dt
E

Pr[p(X) ≥ t]dt

2(C3 /c0 ) log(1/ε) S
C2 d log |S|/δ

Pr[p(X) ≥ t]dt .

2(C3 /c0 ) log(1/ε) S

By applying the contradiction assumption, we have



c0 T
γε
Pr[p(X) ≥ t]dt ≤
exp −
dt
+
2C3
d log |S|/δ
2(C3 /c0 ) log(1/ε) S
2(C3 /c0 ) log(1/ε)


Z C2 d log |S|/δ
c0 T
exp −
≤
dt + C2 γε
2C3
2(C3 /c0 ) log(1/ε)


Z ∞
c0 T
dt + C2 γε
exp −
≤
2C3
2(C3 /c0 ) log(1/ε)
Z ∞
exp (−t) dt + C2 γε
≤ (C3 /c0 )

Z

Z

C2 d log |S|/δ

C2 d log |S|/δ

2 log(1/ε)
2

≤ (C3 /c0 ) · O(ε ) + C2 γε ,

and thus, we have
ψ E[p(X)] ≪ C1 ε log 1/ε − (C1 β + 1 + O(γ) + o(1))ε ,
E

which contradicts Claim 9.
It now suffices to prove that if we construct a filter, then the invariant that ∆ decreases is
preserved. Formally, we show:
Claim 11. Suppose dim(V ) ≥ C1 β log(1/ε). Let S ′ be the set of points we return. Then ∆(S ′ , G) <
∆(S, G).
Proof. Let T be the threshold we pick. If T > C2 d log(|S|/δ) then the invariant is satisfied since we
remove no good points, by (γε, δ)-goodness. It suffices to show that in the other case, we remove
log(1/ε) times many more bad points than good points. By definition we remove at least




γε
c0 T
+
.
|S| · exp −
2C3
d log |S|/δ
points. On the other hand, by (γε, δ)-goodness, Claim 5, we know that we throw away at most
|G| · Pr[p(X) > T ] ≤ |G0 | · Pr[p(X) > T ]
G
G0





γε
(T − kµ − µ̃k22 )2
2
, (T − kµ − µ̃k2 )
+
≤ |G0 | exp −c0 · min
.
C1 β log(1/ε)
2 log(1/ε)(d log(|S|/δ))
42

By our choice of T , and since kµ − µ
ek2 ≤ O(δ) ≪ 1, the first term is upper bounded by






(C3 − 1)2
(T − kµ − µ̃k22 )2
2
, (T − kµ − µ̃k2 )
≤ exp −c0 · min
T, (C3 − 1)T )
exp −c0 · min
C1 β log(1/ε)
C1 β
and so


exp −c0 · min






(C3 − 1)2
c0 T
T, (C3 − 1)T )
≤ exp −
C1 β
C3


c0 T 2
= exp −
2C3


c0 T
·ε
≤ exp −
2C3


1
c0 T
≪
.
exp −
2C3
log2 1/ε

Since we have maintained this invariant so far, in particular, we have thrown away more bad points
than good points, and so |G0 | ≤ (1 + ε)|S|.




γε
c0 T
Therefore, we have log(1/ε) · |G| · PrG [p(X) > T ] ≤ |S| · exp − 2C
+
d log |S|/δ , and hence,
3
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 ∈ Rd be any unit vector. Let G be (η, δ)-good with respect to N (µ, I). Then for any
T > 0,
Pr(|v · (X − µ)| > T ) ≤ 2 exp(−T 2 /2) + η/(d log(d/ηδ)) ,
G

and
Pr (|v · (X − µ)| > T ) ≤ 2 exp(−T 2 /2).

N (µ,I)

ε
Proof. Let hv, µ − µ
bi = R > β 1/2
. Observe that by direct calculation, we have ES [hv, X −
2
2
2
µi ] − ES [hv,
bi ] ≤ 3R . Hence, it suffices to show that ES [hv, X − µi2 ] ≥ 1 + Ω(R2 /ε) −
 X − µ

ε.
γ + O β1
We write S = (G0 \ L, E) where |E| = ψ|S| and |L| = φ|S|. First we consider the expectation
of v · X over X in S. This is
µ0 − φµL + ψµE
.
1−φ+ψ

where µ0 is the mean over G0 , µL the mean over L and µE the mean over E. By (γε, δ)-goodness,
we have that |hv, µ0 − µi| ≤ γε. We also have that
Z ∞
Pr [|hv, X − µi| ≥ t] dt
|hv, µL − µi| ≤
L
0
Z √
O(

d log(d/εδ))

Pr [|hv, X − µi| ≥ t] dt

=

L

0

43

Since we have



|G0 |
Pr [|hv, X − µi| > t] ,
Pr [|hv, (X − µ)i| > t] ≤ min 1,
L
|L| G0

(8)

by Fact 2 and the calculations done in the proof of Corollary 8.8 in [DKK+ 16], we have that
ε
,
|hv, µL − µi| ≤ O(log |S|/|L|) ≤ O(log 1/φ). Since by assumption we have hv, µ − µ
bi = R > β 1/2
this implies that hv, µE − µi = Ω(R/ε). In particular, this implies that


E hv, X − µi2 ≥ E [hv, X − µi]2
E
E
 2
R
≥Ω
.
ε2
Next we consider the expectation of hv, X −µi2 . By (γε, δ)-goodness, we have EG0 [hv, X −µi2 ] =
1 + O(ε). We also have that
Z ∞




2
φ Pr |hv, X − µi|2 ≥ t dt
φ E hv, X − µi ≤
L

0

L

5
≤ log(1/ε)φ +
β
≤

5
log(1/ε)φ +
β

Z

O(d log d/εδ)

(5/β) log(1/ε)

Z

O(d log d/εδ)



φ Pr |hv, X − µi|2 ≥ t dt
L



Pr |hv, X − µi|2 ≥ t dt

(5/β) log(1/ε) G0

Z

O(d log d/εδ)


5
Pr |hv, X − µi|2 ≥ t dt + γε
log(1/ε)φ +
β
(5/β) log(1/ε) N (µ,I)
  

1
= O
+ γ + o(1) ε.
β

=

Since
EG0 [hv, (X − µ)i2 ] − φ EL [hv, (X − µ)i2 ] + ψ EE [hv, (X − µ)i2 ]
,
S
1−φ+ψ

 
ε.
this is at least 1 + Ω(R2 /ε) − γ + O β1
E[hv, (X − µ)i2 ] =

D
D.1

Omitted Proofs from Section 6
Proof of Theorem 15

b − I, and
Our algorithm works as follows, just as for Algorithm 4. It finds all large eigenvalues of Σ
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 (γε, δ)-goodness for this setting, and so we omit them.
Claim 12. ES [p(x)] ≥ Cξ.
Claim 13. EG0 [p(X)] ≤

ε
.
(d log 1/ε)2

Claim 14. φ(S, G0 ) EL [p(X)] = (C + O(γ) + o(1))ε.
44

Algorithm 11 Filter if there are many large eigenvalues of the covariance
1: function FilterCovManyDeg2Eig(S, ε, ξ, δ)
b be the empirical second moment of S, respectively.
2:
Let Σ
b − I with eigenvalue more than
3:
Let V be the subspace of Rd spanned by eigenvectors of Σ
Cξ.
4:
if dim(V ) ≥ C1 log(1/ε) then
5:
Let V ′ be a subspace of V of dimension C1 log(1/ε).
6:
Let p(x) be the quadratic polynomial
p(x) = kΠV ′ (x)k22 − dim(V ′ ).
7:

Find a value T > 0 so that either:
• T > C2 d log(|S|/δ) and p(x) > T for at least one x ∈ S, or

• T > 2C3 log(1/ε)/c0 and PrS (p(x) > T ) > exp(−c0 T /(2C3 )) + ε/(d log(|S|/δ)).

8:

return

S ′ = {x ∈ S : p(x) ≤ T } .
else
return an orthonormal basis for V .
11:
end if
12: end function
9:
10:

Claim 15. ψ EE [p(X)] ≥ Cξ − (C1 + 1 + O(γ) + o(1))ε.
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 ) ≥ C log(1/ε). Then there is a T satisfying the conditions in the
algorithm.
Finally, we show the invariant that ∆ 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 ) ≥ C1 log 1/ε. Let S ′ be the set of points we return. Then ∆(S ′ , G) <
∆(S, G).
Proof. Let T be the threshold we pick. If T > C2 d log(|S|/δ) then the invariant is satisfied since we
remove no good points, by (γε, δ)-goodness. It suffices to show that in the other case, we remove
log 1/ε times many more bad points than good points. By definition we remove at least




c0 T
ε
|S| · exp −
+
2C3
d log |S|/δ
points. On the other hand, observe that if v1 , . . . , vk is an orthonormal basis for V , we have
X
kΠV xk22 =
hvi , Xi2

X
vi viT Σ1/2 X ,
= yΣ1/2
45

so that ifPX ∼ N (0, Σ), we have Y ∼ N (0, I). Let M = Σ1/2
kM kF ≤ ni=1 kΣk2 ≤ (1 + ξ)k, and kM k2 ≤ kΣk2 ≤ (1 + ξ). Since
E [kΠV xk22 ] − k =

N (0,Σ)

≤

k
X

i=1
k
X
i=1

P


vi viT Σ1/2 . We have that

viT (Σ − I)vi

kΣ − Ik2

≤ (1 + ξ)k ,
we have by Hanson-Wright that
Pr

N (0,Σ)







(T − (1 + ξ)k)2
, (T − kµ − (1 + ξ)kk22 )
,
kΠV xk22 − dim(V ) > T ≤ exp −c0 · min
C1 log(1/ε)

so by our choice of T , we have






(T − kµ − µ̃k22 )2
(C3 − 1)2
2
exp −c0 · min
, (T − kµ − µ̃k2 )
≤ exp −c0 · min
T, (C3 − 1)T )
.
C1 log(1/ε)
C1
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 ∆(S ′ , G0 ) < ∆(S, G0 ), which implies the correctness
of the algorithm (Claim 23).
We first show the following claim:
Claim 18. ES [ri ] − EN (0,Σ) [ri ] ≤ (4C + 1)ξ

√
Proof. Let us suppose that pi corresponds to the matrix Ai , given by (Ai )a,b = ∇a ∇b pi / 2!, so
that the Ai are orthonormal with respect to the Frobenius norm. Then the constant harmonic
2
2
2
part
√ of2 pi corresponds to kAi kF ≤ 1. The degree-2 harmonic part of p2i corresponds to the matrix
2 2Ai . This is because if we let Bi be the matrix corresponding to pi , we have


1
∇a ∇b p2i (X)
E
(Bi )a,b = √
2 X∼N (0,I)
2
=√
E
[∇b (pi (X)∇a pi (X))]
2 X∼N (0,I)
2
E
[∇a pi (X)∇b pi (X) + pi (X)∇a ∇b pi (X)]
=√
2 X∼N (0,I)
2
=√
E
[∇a pi (X)∇b pi (X)]
2 X∼N (0,I)
√
= 2 2h(Ai )a , (Ai )b i ,
46

where the second to last line follows since ∇a ∇b pi (X) is a constant, and the last line follows from
explicit computation. In particular,√this implies that the non-constant component of
Pri corresponds
to matrix with trace norm at most
2
2
≤
4.
Therefore,
r
can
be
written
as
r
(x)
=
αi hvi , xi2 +C0
i
i
P
for some constant C0 , where
|αi | ≤ 4. Thus, by our assumption, we have ES [ri ] ≤ 1 + Cξ. The
claim then follows since Σ and I are differ in Frobenius norm by at most ξ.
Claim 19. ES [Q(X)] − EX∼N (0,Σ) [Q(X)] ≥ (C1 − 6)ξk
Proof. Observe that since kΣ − IkF ≤ ξ, in particular we have kΣ⊗2 − I ⊗2 k2 ≤ ξ, and hence
by Lemma 10 we have that EN (0,Σ) [p2 (X)] − EX∼N (0,I) [p2 (X)] ≤ 2ξ for all p ∈ P2 . Hence, in
particular, we have
E[p2i (X)] − E [p2i (X)] ≥ (C1 − 2)ξ .
S

By Claim 18, we have
E
S

"

In particular, this implies that

N (0,Σ)

k
X
i=1

#

ri −

E [Q(X)] −
S

E

N (0,Σ)

"

k
X
i=1

#

ri ≤ (4C + 1)k .

E [Q(X)] ≥ (C1 − 4C − 5)ξk ,

N (0,I)

as claimed.
We now show:
Claim 20. There is some universal constant B so that
E

√
[Q2 (X)]1/2 ≤ B k .

X∼N (0,Σ)

Before we prove this, we need the following lemma:
Lemma 20. For any degree-4 polynomial p, and any Σ, if we let F (y, p, Σ) denote the yth percentile
of p under Σ, then we have F (1/4, p2 , Σ), F (3/4, p2 , Σ) = Θ(EX∼N (0,Σ) [p2 (X)]).
Proof. Let µ′ = EX∼N (0,Σ) [p2 (X)]. First, we note that Pr(p2 (X) > 4µ′ ) ≤ 1/4, so F (3/4, p2 , Σ) ≤
4µ′ . On the other hand,
anti-concentration bounds [CW01], we have that Pr(p2 (X) ≤
√ ′ by known
1/8
′
εµ ) = Pr(|p(X)| ≤ εµ ) = O(ε ). So, for ε a sufficiently small constant, Pr(p2 (X) ≤ εµ′ ) < 1/4,
and therefore, F (1/4, p2 , Σ) ≥ εµ′ . Since εµ′ ≤ F (1/4, p2 , Σ) ≤ F (3/4, p2 , Σ) ≤ 4µ′ , this completes
our proof.
Proof of Claim 20. Since dTV (N (0, Σ), N (0, I)) = O(ε log 1/ε) = o(1), we have that F (3/4, Q2 , Σ) ≤
F (1/4, Q2 , Σ), so by the above lemma, we have that
E

[Q2 (X)]1/2 ≤ O(

X∼N (0,Σ)

E

[Q2 (X)]1/2 ) .

X∼N (0,I)

Hence, it suffices to bound kQk22 .

47

√
Let us again suppose that pi corresponds to the matrix Ai , given by (Ai )a,b = ∇a ∇b pi / 2!.
Note that the Ai are symmetric matrices that form an orthonormal set. We note that qi is the
harmonic degree-4 polynomial corresponding to the rank-4 tensor
√
Ta,b,c,d = ∇a ∇b ∇c ∇d qi / 24
√
= ∇a ∇b ∇c ∇d (p2i )/ 24
√
= (2(∇a ∇b pi )(∇c ∇d pi ) + 2(∇a ∇c pi )(∇b ∇d pi ) + 2(∇a ∇d pi )(∇b ∇c pi ))/ 24
√
= (4(Ai )a,b (Ai )c,d + 4(Ai )a,c (Ai )b,d + 4(Ai )a,d (Ai )b,c )/ 24
√
= 6Sym(Ai ⊗ Ai ).

√ P
By linearity, Q corresponds to the rank-4 tensor T = 6 ki=1 Sym(Ai ⊗ Ai ). It thus suffices to
show that kT k2 = O(k), where here k · k2 denotes the square root of the sum of the squares of the
entries of the tensor. In order to show this, we note that kT k2 = supkV k2 =1 hV, T i. Therefore, it
√
suffices to show that for all 4-tensors V with kV k2 ≤ 1, that hV, T i = O( k). We note that
+
+
*
*
k
k
X
X
√
√
hV, T i = 6 V,
Sym(Ai ) ⊗ Ai = 6 Sym(V ),
Ai ⊗ Ai .
i=1

i=1

Note that Sym(V ) is a symmetric 4-tensor of ℓ2 norm at most 1. Thinking of Sym(V ) as a symmetric
Pn 2
matrix over 2-Tensors, we can write it as Sym(V ) =
j=1 λj Bj ⊗ Bj , where the Bj ’s are an
Pn 2 2
orthonormal basis for the set of 2-tensors and j=1 λj = kSym(V )k22 ≤ 1. Then, if ci,j := hAi , Bj i,
P
and cj = ki=1 c2i,j , we have that
2

2

k
n
n
X
√ X
√ X
c2i,j = 6
λj
λj cj .
hV, T i = 6
i=1

j=1

j=1

Since the Ai are orthonormal, cj ≤ 1 for all j. Furthermore,
2

2

n
X

cj =

c2i,j =

k
X
i=1

i=1 j=1

j=1

Therefore,

n
k X
X

kAi k2F = k.

v
v
u n2
u n2
n2
uX
√
√
√ X
√ uX
hV, T i = 6
λj cj ≤ 6t
λ2j t
c2j ≤ 6 · 1 · k = O( k) .
j=1

j=1

j=1

Let µ′ = EX∼N (0,Σ) [Q(X)]. We now show that since G0 is ε-good, then almost all of the
difference in Claim 19 must be because of the points in E.
Claim 21. EE [Q(Y )] − µ′ ≥

C1 −6
2

·

ξk
φ

Proof. It suffices to show that
E[Q(X)] − µ′ ≤
G

48

C1 − 6
ξk ,
2

(9)

that is, the good points do not contribute
 √ much to the difference. By (ε, δ)-goodness of G0 and
′
Claim 20 we have E0 [Q(Y )] − µ ≤ O ε k . Moreover, we have
φ(E[Q(X)] −
L

E

[Q(X)]) ≤

X∼N (0,Σ)

(a)

=

Z

∞

0

Z



φ Pr Q(X) − µ′ ≥ t dt

L
√ 2
O( kd log |S|/δ)



φ Pr Q(X) − µ′ ≥ t dt
L

0

Z
√
2
≤ φ · O( k log 1/ξ) +

√
≤ φ · O( k log2 1/ξ) +

√
≤ φ · O( k log2 1/ξ) +

(b)

√
O( kd2 log |S|/δ)

φ Pr
√
L
O( k log2 1/ξ)
Z O(√kd2 log |S|/δ)

Pr
√
G0
O( k log2 1/ξ)
Z O(√kd2 log |S|/δ)


Q(X) − µ′ ≥ t dt


Q(X) − µ′ ≥ t dt

Pr [Q(X) − µ′ ≥ t]dt + O(ξ)

√
O( k log2 1/ξ)

N (0,Σ)

√
≤ O(ξ k) + exp(−Ω(log 1/ξ))

(c)



≪ 5ξk ,

where (a) follows from the boundedness condition of (ε, δ)-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(Vm ) ≥ k, then Algorithm 8 returns a T satisfying the conditions in the algorithm.
Proof. Suppose not. By the assumption that kI − ΣkF ≤ ξ, we have |µ′ | = | EX∼N (0,Σ) [Q(X)]| ≤
P
2ξkQk2 ≤ ki=1 kpi k2 + kri k2 ≤ 5ξk. Thus, we have




C1 − 6
C1 − 6
− 5 ξk ≫
E[Q(X)] ≥
ξk .
E
2ψ
3ψ
But since we assume there is no T satisfying the conditions in Algorithm 8, we have
Z ∞
Pr[p(X) > t]dt
ψ E [p(X)] ≤ ψ
X∈E

0

E

√

2

≤ ψ · 4B k log (1/ξ) +
≤ 4Bψ ·

√

k log2 (1/ξ) +

≤ 4C2 Bψ ·

√

Z

√
O(d2 k log(|S|/δ))

ψ Pr[p(X)
√
E
4B k log2 1/ξ
Z O(d2 √k log(|S|/δ))
Pr[p(X)
√
S
4B k log2 1/ξ
Z O(d2 √k log(|S|/δ))

k log2 (1/ξ) +

√
4A2 C2 B k log2 1/ξ
√ Z ∞

√
≤ 4A2 C2 Bψ k log2 (1/ξ) + 4B k
√
e 2 ) ≪ 10ξk ,
≤ O( kξ) + O(ε

A2 log2 1/ξ

> t]dt

> t]dt


t
√
exp −A
4B k

1/2 !

ε2
+ 2
d log |G|/δ

e 2)
exp(−Ω(t1/2 ))dt + O(ε

which is a contradiction, for our choice of C1 , C2 , C3 , and since we chose k = O(log4 1/ε).
49

!

dt

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(Vm ) ≥ k, then Algorithm 8 returns a S ′ satisfying ∆(S ′ , G0 ) < ∆(S, G0 ).
√
Proof. Let T be the threshold we pick. If T > C3 d2 k log(|S|) then the invariant is satisfied since
we remove no good points, by (ε, δ)-goodness. It suffices to show that we remove log 1/ε times
many more bad points than good points. Otherwise, by definition we remove a total of
!
1/2 !

ε2
T
√
+ 2
|S| · exp −A
d log2 |S|/δ
2B k
points. On the other hand, by (ε, δ)-goodness, hypercontractivity, and Claim 20, we know that the
total number of points we throw away is at most
|G| · Pr[Q(X) > T ] ≤ |G0 | · Pr[Q(X) > T ]
G

G0

−1/2 !

!
ε2
+
≤ |G0 | exp −A
2 log(1/ε)(d log(|G|/δ))2
!

1/2 !
T
ε2
√
≤ |G0 | exp −A
.
+
2 log(1/ε)(d log(|G|/δ))2
B k


T
kQk2

Since we have maintained this invariant so far, in particular, we have√
thrown away more bad points
than good points, and so |G0 | ≤ (1 + ε)|S|. Moreover, since T ≥ 4B k log2 1/ε, we have
exp −A



T
√
B k

1/2 !

= exp −A
≤ exp −A
≤ exp −A
≤ exp −A









T
√
4B k
T
√
4B k
T
√
4B k
T
√
4B k

1/2 !2
1/2 !
1/2 !
1/2 !



· exp −A
· ε−1
·

!1/2 
√
2
4B k log 1/ε

√
4B k

1
.
log 1/ε




1/2 
Therefore, we have log(1/ε) · |G| · PrG [Q(X) > T ] ≤ |S| · exp −A 4BT√k
+

hence, the invariant is satisfied.

ε2
d2 log |G|/δ



, and

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(δ) after
taking O(poly(d, 1/η, log 1/δ)) samples. Thus, it suffices to show that the last property holds with
probability 1 − O(δ) given O(poly(d, 1/η, log 1/δ)) samples. We may clearly WLOG take Σ = I.
Moreover, these properties are clearly invariant under scaling, and hence, it suffices to prove them
50

for degree-four polynomials with VarN (0,I) [p(X)] = 1. For any fixed degree-4 polynomial p, by
P
hypercontractivity, since P (X) = n1 ni=1 p(Xi ) has variance Var[P (X)] = n1 , we have that
"
#
n


1X
Pr
p(Xi ) − E [p(X)] ≥ t ≤ exp −At1/2 n1/4 .
n
N (0,I)
i=1

In particular, if we take
n=Ω
then
Pr

"



log4 1/δ
η2



,

#
n
1X
p(Xi ) − E [p(X)] ≥ η ≤ O(δ) .
n
N (0,I)
i=1

4

Since there is a 1/3-net over all degree-4polynomials
with unit variance of size (1/3)O(d ) , by union

d4 log1/δ
samples, then the second to last property holds
bounding, we obtain that if we take Ω
η2
with probability 1 − δ.
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 poly(d, log 1/η, log 1/δ) samples.

D.4

Proof of Lemma 15

Proof. Let µ′ = EN (0,Σ) [p(X)], s = EN (0,Σ) [p2 (X)], and let ES [p(X)] − EN (0,Σ) [p(X)] = K. By
(ε, δ)-goodness and Lemma 10, we have
E [p(X)] − µ′ ≤ ε

E [p2 (X)]1/2 ≤ 2ε .

G0

N (0,Σ)

(10)

Moreover, for some appropriate choice of β1 , we have
Z ∞
′
φ Pr [|p(X) − µ| ≥ t] dt
φ E[p(X)] − φµ ≤
L

0

L

≤ β1 φ log 1/ε +
(a)

Z

≤ β1 φ log 1/ε +

(b)

≤ β1 φ log 1/ε +

(c)

O(d log |S|/δ)

β1 log 1/ε

|G0 |
|S|

|G0 |
|S|

Z

Z

φ Pr [|p(X) − µ| ≥ t] dt
L

O(d log |S|/δ)
β1 log 1/ε
O(d log |S|/δ)

β1 log 1/ε

Pr [|p(X) − µ| ≥ t] dt
G0

Pr [|p(X) − µ| ≥ t] dt + O(ε)

N (0,Σ)

≤ β1 φ log 1/ε + εO(β1 ) + O(ε)

= O(ε) ,

(11)

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
2K
K
≥ E[p(X) − µ′ ] ≥
,
E
ψ
2ψ

51

so this implies

K2
.
4ψ 2

E[(p(X) − µ′ )2 ] ≥
E

Expanding this yields that
E[(p(X) − µ′ )2 ] = E[p2 (X)] − 2µ′ E[p(X)] + (µ′ )2 ,
E

E

E

so since |µ′ | ≤ kΣ − IkF ≤ ξ, we have
E[p2 (X)] ≥
E

2K
K2
−
+ O(ξ 2 ) .
4ψ 2
ψ

(12)

On the other hand, by Lemma 10, we have
1
|s − 1| ≤ kΣ⊗2 − I ⊗2 k2 + kΣ − Ik2F
2
≤ ξ + O(ξ 2 ) .

(13)

Also, by (ε, δ)-goodness and Lemma 10, we have
E [p2 (X)] − s ≤ ε(1 + O(ξ)) ,

(14)

G0

and, for A as in Corollary 6, and β appropriately chosen,
Z ∞
2
φ(E[p (X)] − s) ≤
φ Pr[|p2 (X) − s| ≥ t]dt
L

0

L

2

≤ βAφ log 1/ε +

Z

O((d log |S|/δ)2 )

|G0 |
≤ βAφ log 1/ε +
|S|
2

(a)

≤ βAφ log2 1/ε +

(b)

≤ βAφ log2 1/ε +

φ Pr[|p2 (X) − s| ≥ t]dt
L

2A log2 1/ε

|G0 |
|S|
Z ∞

Z

O((d log |S|/δ)2 )

βA log 1/ε

Z

Pr[|p2 (X) − s| ≥ t]dt
G0

2

O((d log |S|/δ)2 )
2

βA log 1/ε
1/2 /2

e−At

2

Pr [|p2 (X) − s| ≥ t]dt + O(ε2 )

N (0,Σ)

dt + O(ε2 )

βA log 1/ε
2

≪ 2βAξ + O(ε) + O(ε ) ,

(15)

where (a) follows from goodness, and (b) follows from Corollary 6 and since kΣ − IkF = o(1). Thus,
(12), (13), (14), and (15) together imply that
E[p2 (X)] − 1 ≥ 1 − (ξ + O(ξ 2 )) − O(ε) − 2βAξO(ε) +
S

and so since K ≥ Ω

√ 
εξ , this gives the desired bound.

52

K2
− 2K + O(ξ 2 ) ,
4ψ

