Being Robust (in High Dimensions) Can Be Practical∗

arXiv:1703.00893v4 [cs.LG] 13 Mar 2018

Ilias Diakonikolas†
CS, USC
diakonik@usc.edu
Jerry Li ¶
EECS & CSAIL, MIT
jerryzli@mit.edu

Gautam Kamath‡
EECS & CSAIL, MIT
g@csail.mit.edu
Ankur Moitrak
Math & CSAIL, MIT
moitra@mit.edu

Daniel M. Kane§
CSE & Math, UCSD
dakane@cs.ucsd.edu
Alistair Stewart∗∗
CS, USC
alistais@usc.edu

March 14, 2018

Abstract
Robust estimation is much more challenging in high dimensions than it is in one dimension: Most
techniques either lead to intractable optimization problems or estimators that can tolerate only a tiny
fraction of errors. Recent work in theoretical computer science has shown that, in appropriate distributional models, it is possible to robustly estimate the mean and covariance with polynomial time
algorithms that can tolerate a constant fraction of corruptions, independent of the dimension. However,
the sample and time complexity of these algorithms is prohibitively large for high-dimensional applications. In this work, we address both of these issues by establishing sample complexity bounds that are
optimal, up to logarithmic factors, as well as giving various refinements that allow the algorithms to
tolerate a much larger fraction of corruptions. Finally, we show on both synthetic and real data that our
algorithms have state-of-the-art performance and suddenly make high-dimensional robust estimation a
realistic possibility.

1

Introduction

Robust statistics was founded in the seminal works of [Tuk60] and [Hub64]. The overarching motto is
that any model (especially a parametric one) is only approximately valid, and that any estimator designed
for a particular distribution that is to be used in practice must also be stable in the presence of model
misspecification. The standard setup is to assume that the samples we are given come from a nice distribution,
but that an adversary has the power to arbitrarily corrupt a constant fraction of the observed data. After
several decades of work, the robust statistics community has discovered a myriad of estimators that are
provably robust. An important feature of this line of work is that it can tolerate a constant fraction of
corruptions independent of the dimension and that there are estimators for both the location (e.g., the
mean) and scale (e.g., the covariance). See [HR09] and [HRRS86] for further background.
It turns out that there are vast gaps in our understanding of robustness, when computational considerations are taken into account. In one dimension, robustness and computational efficiency are in perfect
∗A

version of this paper appeared in ICML 2017 [DKK+ 17].
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, a Google Faculty Research Award, and an NSF Fellowship.
k Supported by NSF CAREER Award CCF-1453261, a grant from the MIT NEC Corporation, and a Google Faculty Research
Award.
∗∗ Research supported by a USC startup grant.
Authors are in alphabetical order.
Code of our implementation is available at https://github.com/hoonose/robust-filter.
† Supported

1

harmony. The empirical mean and empirical variance are not robust, because a single corruption can arbitrarily bias these estimates, but alternatives such as the median and the interquartile range are straightforward
to compute and are provably robust.
But in high dimensions, there is a striking tension between robustness and computational efficiency.
Let us consider estimators for location. The Tukey median [Tuk60] is a natural generalization of the onedimensional median to high-dimensions. It is known that it behaves well (i.e., it needs few samples) when
estimating the mean for various symmetric distributions [DG92, CGR16]. However, it is hard to compute
in general [JP78, AK95] and the many heuristics for computing it degrade badly in the quality of their
approximation as the dimension scales [CEM+ 93, Cha04, MS10]. The same issues plague estimators for
scale. The minimum volume ellipsoid [Rou85] is a natural generalization of the one-dimensional interquartile
range and is provably robust in high-dimensions, but is also hard to compute. And once again, heuristics
for computing it [VAR09, RS98] work poorly in high dimensions.
The fact that robustness in high dimensions seems to come at such a steep price has long been a point of
consternation within robust statistics. In a 1997 retrospective on the development of robust statistics [Hub97],
Huber laments:
“It is one thing to design a theoretical algorithm whose purpose is to prove [large fractions of
corruptions can be tolerated] and quite another thing to design a practical version that can be
used not merely on small, but also on medium sized regression problems, with a 2000 by 50 matrix
or so. This last requirement would seem to exclude all of the recently proposed [techniques].”
The goal of this paper is to answer Huber’s call to action and design estimators for both the mean and
covariance that are highly practical, provably robust, and work in high-dimensions. Such estimators make
the promise of robust statistics – estimators that work in high-dimensions and guarantee that their output
has not been heavily biased by some small set of noisy samples – much closer to a reality.
First, we make some remarks to dispel some common misconceptions. There has been a considerable
amount of recent work on robust principal component analysis, much of it making use of semidefinite programming. Some of these works can tolerate a constant fraction of corruptions [CLMW11], however require
that the locations of the corruptions are evenly spread throughout the dataset so that no individual sample
is entirely corrupted. In contrast, the usual models in robust statistics are quite rigid in what they require
and they do this for good reason. A common scenario that is used to motivate robust statistical methods is
if two studies are mixed together, and one subpopulation does not fit the model. Then one wants estimators
that work without assuming anything at all about these outliers.
There have also been semidefinite programming methods proposed for robust principal component analysis with outliers [XCS10]. These methods assume that the uncorrupted matrix is rank r and that the
fraction of outliers is at most 1/r, which again degrades badly as the rank of the matrix increases. Moreover,
any method that uses semidefinite programming will have difficulty scaling to the sizes of the problems we
consider here. For sake of comparison – even with state-of-the-art interior point methods – it is not currently feasible to solve the types of semidefinite programs that have been proposed when the matrices have
dimension larger than a hundred.

1.1

Robustness in a Generative Model

Recent works in theoretical computer science have sought to circumvent the usual difficulties of designing
efficient and robust algorithms by instead working in a generative model. The starting point for our paper is
the work of [DKK+ 16] who gave an efficient algorithm for the problem of agnostically learning a Gaussian:
Given a polynomial number of samples from a high-dimensional Gaussian N (µ, Σ), where an
b that satisfy
adversary has arbitrarily corrupted an ε-fraction, find a set of parameters N 0 (b
µ, Σ)
0
∗
e
dT V (N , N ) ≤ O(ε)
.
Total variation distance is the natural metric to use to measure closeness of the parameters, since a
(1 − ε)-fraction of the observed samples came from a Gaussian. [DKK+ 16] gave an algorithm for the above
∗ We

use the notation Õ(·) to hide factors which are polylogarithmic in the argument – in particular, we note that this
bound does not depend on the dimension.

2

problem (note that the guarantees are dimension independent), whose running time and sample complexity
are polynomial in the dimension d and 1/ε. [LRV16] independently gave an algorithm for the unknown mean
√
e
log d), and in the unknown covariance case achieves guarantees in
case that achieves dT V (N , N 0 ) ≤ O(ε
a weaker metric that is not affine invariant. A crucial feature is that both algorithms work even when the
moments of the underlying distribution satisfy certain conditions, and thus are not necessarily brittle to the
modeling assumption that the inliers come from a Gaussian distribution.
A more conceptual way to view such work is as a proof-of-concept that the Tukey median and minimum
volume ellipsoid can be computed efficiently in a natural family of distributional models. This follows because
not only would these be good estimates for the mean and covariance in the above model, but in fact any
estimates that are good must also be close to them. Thus, these works fit into the emerging research direction
of circumventing worst-case lower bounds by going beyond worst-case analysis.
Since the dissemination of the aforementioned works [DKK+ 16, LRV16], there has been a flurry of research
activity on computationally efficient robust estimation in a variety of high-dimensional settings [DKS16,
DKS17, CSV17, DKK+ 17, Li17, DBS17, BDLS17, SCV18, DKK+ 18], including studying graphical distributional models [DKS16], understanding the computation-robustness tradeoff for statistical query algorithms [DKS17], tolerating much more noise by allowing the algorithm to output a list of candidate hypotheses [CSV17], and developing robust algorithms under sparsity assumptions [Li17, DBS17, BDLS17], where
the number of samples is sublinear in the dimension.

1.2

Our Results

Our goal in this work is to show that high-dimensional robust estimation can be highly practical. However,
there are two major obstacles to achieving this. First, the sample complexity and running time of the
algorithms in [DKK+ 16] is prohibitively large for high-dimensional applications. We just would not be able
to store as many samples as we would need, in order to compute accurate estimates, in high-dimensional
applications.
Our first main contribution is to show essentially tight bounds on the sample complexity of the filtering
based algorithm of [DKK+ 16]. Roughly speaking, we accomplish this with a new definition of the good set
which plugs into the existing analysis in a straightforward manner and shows that it is possible to estimate
2
e
e 2 /ε2 ) samples.
the mean with O(d/ε
) samples (when the covariance is known) and the covariance with O(d
Both of these bounds are information-theoretically optimal, up to logarithmic factors.
Our second main contribution is to vastly improve the fraction of adversarial corruptions that can be
tolerated in applications. The fraction of errors that the algorithms of [DKK+ 16] can tolerate is indeed a
constant that is independent of the dimension, but it is very small both in theory and in practice. This is
due to the fact that many of the steps in the algorithm are overly conservative. In fact, we found that a
naive implementation of the algorithm did not remove any outliers in many realistic scenarios. We combat
this by giving new ways to empirically tune the threshold for where to remove points from the sample set.
These optimizations dramatically improve the empirical performance.
Finally, we show that the same bounds on the error guarantee continue to work even when the underlying
distribution is sub-Gaussian. This theoretically confirms that the robustness guarantees of such algorithms
are in fact not overly brittle to the distributional assumptions. In fact, the filtering algorithm of [DKK+ 16]
is easily shown to be robust under much weaker distributional assumptions, while retaining near-optimal
sample and error guarantees. As an example, we show that it yields a near sample-optimal efficient estimator
for robustly estimating the mean of a distribution, under the assumption that its covariance is bounded.
Even in this regime, the filtering algorithm guarantees optimal error, up to a constant factor. Furthermore
we empirically corroborate this finding by showing that the algorithm works well on real world data, as we
describe below.
Now we come to the task of testing out our algorithms. To the best of our knowledge, there have been no
experimental evaluations of the performance of the myriad of approaches to robust estimation. It remains
mostly a mystery which ones perform well in high-dimensions, and which do not. To test out our algorithms,
we design a synthetic experiment where a (1 − ε)-fraction of the samples come from a Gaussian and the
rest are noise and sampled from another distribution (in many cases, Bernoulli). This gives us a baseline
to compare how well various algorithms recover µ and Σ, and how their performance degrades based on the
dimension. Our plots show a predictable and yet striking phenomenon: All earlier approaches have error

3

rates that scale polynomially with the dimension and ours is a constant that is almost indistinguishable from
the error that comes from sample noise alone. Moreover, our algorithms are able to scale to hundreds of
dimensions.
But are algorithms for agnostically learning a Gaussian unduly sensitive to the distributional assumptions
they make? We are able to give an intriguing visual demonstration of our techniques on real data. The
famous study of [NJB+ 08] showed that performing principal component analysis on a matrix of genetic
data recovers a map of Europe. More precisely, the top two singular vectors define a projection into the
plane and when the groups of individuals are color-coded with where they are from, we recover familiar
country boundaries that corresponds to the map of Europe. The conclusion from their study was that genes
mirror geography. Given that one of the most important applications of robust estimation ought to be in
exploratory data analysis, we ask: To what extent can we recover the map of Europe in the presence of noise?
We show that when a small number of corrupted samples are added to the dataset, the picture becomes
entirely distorted (and this continues to hold even for many other methods that have been proposed). In
contrast, when we run our algorithm, we are able to once again recover the map of Europe. Thus, even
when some fraction of the data has been corrupted (e.g., medical studies were pooled together even though
the subpopulations studied were different), it is still possible to perform principal component analysis and
recover qualitatively similar conclusions as if there were no noise at all!

2

Formal Framework

Notation. For a vector v, we will let kvk2 denote its Euclidean norm. If M is a matrix, we will let kM k2
denote its spectral norm and kM kF denote its Frobenius norm. We will write X ∈u S to denote that X is
drawn from the empirical distribution defined by S.
Robust Estimation. We consider the following powerful model of robust estimation that generalizes many
other existing models, including Huber’s contamination model:
Definition 2.1. Given ε > 0 and a distribution family D, the adversary operates as follows: The algorithm
specifies some number of samples m. The adversary generates m samples X1 , X2 , . . . , Xm from some (unknown) D ∈ D. It then draws m0 from an appropriate distribution. This distribution is allowed to depend
on X1 , X2 , . . . , Xm , but when marginalized over the m samples satisfies m0 ∼ Bin(ε, m). The adversary is
allowed to inspect the samples, removes m0 of them, and replaces them with arbitrary points. The set of m
points is then given to the algorithm.
In summary, the adversary is allowed to inspect the samples before corrupting them, both by adding
corrupted points and deleting uncorrupted points. In contrast, in Huber’s model the adversary is oblivious
to the samples and is only allowed to add corrupted points.
We remark that there are no computational restrictions on the adversary. The goal is to return the
b in D that are close to the true parameters in an appropriate metric. For the
parameters of a distribution D
case of the mean, our metric will be the Euclidean distance. For the covariance, we will use the Mahalanobis
b −1/2 − IkF . This is a strong affine invariant distance that implies corresponding
distance, i.e., kΣ−1/2 ΣΣ
bounds in total variation distance.
We will use the following terminology:
Definition 2.2. We say that a set of samples is ε-corrupted if it is generated by the process described in
Definition 2.1.

3

Nearly Sample-Optimal Efficient Robust Learning

In this section, we present near sample-optimal efficient robust estimators for the mean and the covariance
of high-dimensional distributions under various structural assumptions of varying strength. Our estimators
rely on the filtering technique introduced in [DKK+ 16].
We note that [DKK+ 16] gave two algorithmic techniques: the first one was a spectral technique to
iteratively remove outliers from the dataset (filtering), and the second one was a soft-outlier removal method
relying on convex programming. The filtering technique seemed amenable to practical implementation (as
4

it only uses simple eigenvalue computations), but the corresponding sample complexity bounds given in
[DKK+ 16] are polynomially worse than the information-theoretic minimum. On the other hand, the convex
programming technique of [DKK+ 16] achieved better sample complexity bounds (e.g., near sample-optimal
for robust mean estimation), but relied on the ellipsoid method, which seemed to preclude a practically
efficient implementation.
In this work, we achieve the best of both worlds: we provide a more careful analysis of the filter technique that yields sample-optimal bounds (up to logarithmic factors) for both the mean and the covariance.
Moreover, we show that the filtering technique easily extends to much weaker distributional assumptions
(e.g., under bounded second moments). Roughly speaking, the filtering technique follows a general iterative
recipe: (1) via spectral methods, find some univariate test which is violated by the corrupted points, (2)
find some concrete tail bound violated by the corrupted set of points, and (3) throw away all points which
violate this tail bound.
We start with sub-gaussian distributions. Recall that if P is sub-gaussian on Rd with mean vector µ and
parameter ν > 0, then for any unit vector v ∈ Rd we have that PrX∼P [|v · (X − µ)| ≥ t] ≤ exp(−t2 /2ν).
Theorem 3.1. Let G be a sub-gaussian distribution on Rd with parameter ν = Θ(1), mean µG , covariance
matrix I, and ε > 0. Let S be an ε-corrupted set of samples from G of size Ω((d/ε2 ) poly log(d/ε)). There
exists an efficient algorithm that, onp
input S and ε > 0, returns a mean vector µ
b so that with probability at
least 9/10 we have kb
µ − µG k2 = O(ε log(1/ε)).
[DKK+ 16] gave algorithms for robustly estimating the mean of a Gaussian distribution with known
covariance and for robustly estimating the mean of a binary product distribution. The main motivation for
considering these specific distribution families is that robustly estimating the mean within Euclidean distance
immediately implies total variation distance bounds for these families. The above theorem establishes that
these guarantees hold in a more general setting with near sample-optimal bounds. Under a bounded second
moment assumption, we show:
Theorem 3.2. Let P be a distribution on Rd with unknown mean vector µP and unknown covariance matrix
ΣP  σ 2 I. Let S be an ε-corrupted set of samples from P of size Θ((d/ε) log d). There exists
√ an efficient
algorithm that, on input S and ε > 0, with probability 9/10 outputs µ
b with kb
µ − µP k2 ≤ O( εσ).
A similar result on mean estimation under bounded second moments was concurrently shown in [SCV18].
The sample size above is optimal, up to a logarithmic factor, and the error guarantee is easily seen to be
the best possible up to a constant factor. The main difference between the filtering algorithm establishing
the above theorem and the filtering algorithm for the sub-gaussian case is how we choose the threshold for
the filter. Instead of looking for a violation of a concentration inequality, here we will choose a threshold at
random. In this case, randomly choosing a threshold weighted towards higher thresholds suffices to throw
out more corrupted samples than uncorrupted samples in expectation. Although it is possible to reject many
good samples this way, we show that the algorithm still only rejects a total of O(ε) samples with high
probability.
Finally, for robustly estimating the covariance of a Gaussian distribution, we have:
Theorem 3.3. Let G ∼ N (0, Σ) be a Gaussian in d dimensions, and let ε > 0. Let S be an ε-corrupted
set of samples from G of size Ω((d2 /ε2 ) poly log(d/ε)). There exists an efficient algorithm that, given S and
b so that with probability at least 9/10, it
ε, returns the parameters of a Gaussian distribution G0 ∼ N (0, Σ)
−1/2 b −1/2
holds kI − Σ
ΣΣ
kF = O(ε log(1/ε)).
We now provide a high-level description of the main ingredient which yields these improved sample
complexity bounds. The initial analysis of [DKK+ 16] established sample complexity bounds which were suboptimal by polynomial factors because it insisted that the set of good samples (i.e., before the corruption)
satisfied very tight tail bounds. To some degree such bounds are necessary, as when we perform our filtering
procedure, we need to ensure that not too many good samples are thrown away. However, the old analysis
required that fairly strong tail bounds hold uniformly. The idea for the improvement is as follows: If the
errors are sufficient to cause the variance of some polynomial p (linear in the unknown mean case or quadratic
in the unknown covariance case) to increase by more than ε, it must be the case that for some T , roughly
an ε/T 2 fraction of samples are error points with |p(x)| > T . As long as we can ensure that less than an
5

ε/T 2 fraction of our good sample points have |p(x)| > T , this will suffice for our filtering procedure to work.
For small values of T , these are much weaker tail bounds than were needed previously and can be achieved
with a smaller number of samples. For large values of T , these tail bounds are comparable to those used
in previous work [DKK+ 16] , but in such cases we can take advantage of the fact that |p(G)| > T only
with very small probability, again allowing us to reduce the sample complexity. The details are deferred to
Appendix A.

4

Filtering

We now describe the filtering technique more rigorously. We also describe some additional heuristics we
found useful in practice.

4.1

Robust Mean Estimation

We first consider mean estimation. The algorithms which achieve Theorems 3.1 and 3.2 both follow the
general recipe in Algorithm 1. We must specify three parameter functions:
• Thres(ε) is a threshold function—we terminate if the covariance has spectral norm bounded by Thres(ε).
• Tail(T, d, ε, δ, τ ) is an univariate tail bound, which would only be violated by a τ fraction of points if
they were uncorrupted, but is violated by many more of the current set of points.
• δ(ε, s) is a slack function, which we require for technical reasons.
Given these objects, our filter is fairly easy to state: first, we compute the empirical covariance. Then, we
check if the spectral norm of the empirical covariance exceeds Thres(ε). If it does not, we output the empirical
mean with the current set of data points. Otherwise, we project onto the top eigenvector of the empirical
covariance, and throw away all points which violate Tail(T, d, ε, δ, τ ), for some choice of slack function δ.
Algorithm 1 Filter-based algorithm template for robust mean estimation
1: Input: An ε-corrupted set of samples S, Thres(ε), Tail(T, d, ε, δ, τ ), δ(ε, s)
0
2: Compute the sample mean µS = EX∈u S 0 [X]
3: Compute the sample covariance matrix Σ
4: Compute approximations for the largest absolute eigenvalue of Σ, λ∗ := kΣk2 , and the associated unit
eigenvector v ∗ .
5: if kΣk2 ≤ Thres(ε) then
0
6:
return µS .
7: Let δ = δ(ε, kΣk2 ).
8: Find T > 0 such that
h
i
0
Pr 0 |v ∗ · (X − µS )| > T + δ > Tail(T, d, ε, δ, τ ).
X∈u S
0

9:

return {x ∈ S 0 : |v ∗ · (x − µS )| ≤ T + δ}.

Sub-gaussian case To
p concretely instantiate this algorithm for the subgaussian case, we take Thres(ε) =
O(ε log 1/ε), δ(ε, s) = 3 ε(s − 1), and
Tail(T, d, ε, δ, τ ) = 8 exp(−T 2 /2ν) + 8

ε
,
T 2 log(d log(d/ετ ))

where ν is the subgaussian parameter. See Section A.1 for details.
Second moment case To concretely instantiate this algorithm for the second moment case, we take
Thres(ε) = 9, δ = 0, and we take Tail to be a random rescaling of the largest deviation in the data set, in
the direction v ∗ . See Section A.2 for details.
6

4.2

Robust Covariance Estimation

Our algorithm for robust covariance follows the exact recipe outlined above, with one key difference—we
check for deviations in the empirical fourth moment tensor. Intuitively, just as in the robust mean setting,
we used degree-2 information to detect outliers for the mean (the degree-1 moment), here we use degree-4
information to detect outliers for the covariance (the degree-2 moment).
More concretely, this corresponds to finding a normalized degree-2 polynomial whose empirical variance
is too large. By then filtering along this polynomial, with an appropriate choice of Thres(ε), δ(ε, s), and Tail,
we achieve the desired bounds. See Section A.3 for the formal pseudocode and more details.

4.3

Better Univariate Tests

In the algorithms described above for robust mean estimation, after projecting onto one dimension, we center
the points at the empirical mean along this direction. This is theoretically sufficient, however, introduces
additional constant factors since the empirical mean along this direction may be corrupted. Instead, one can
use a robust estimate for the mean in one direction. Namely, it is well known that the median is a provably
robust estimator for the mean for symmetric distributions [HR09, HRRS86], and under certain models it is
in fact optimal in terms of its resilience to noise [DKW56, Mas90, Che98, DK14, DKK+ 17]. By centering
the points at the median instead of the mean, we are able to achieve better error in practice.

4.4

Adaptive Tail Bounding

In our empirical evaluation, we found that it was important to find an appropriate choice of Tail, to achieve
good error rates, especially for robust covariance estimation. Concretely, in this setting, our tail bound is
given by
Tail(T, d, ε, δ, τ ) = C1 exp(−C2 T ) + Tail2 (T, d, ε, δ, τ ) ,
for some function Tail2 , and constants C1 , C2 . We found that for reasonable settings, the term that dominated
was always the first term on the RHS, and that Tail2 is less significant. Thus, we focused on optimizing the
first term.
We found that depending on the setting, it was useful to change the constant C2 . In particular, in low
dimensions, we could be more stringent, and enforce a stronger tail bound (which corresponds to a higher
C2 ), but in higher dimensions, we must be more lax with the tail bound. To do this in a principled manner,
we introduced a heuristic we call adaptive tail bounding. Our goal is to find a choice of C2 which throws away
roughly an ε-fraction of points. The heuristic is fairly simple: we start with some initial guess for C2 . We
then run our filter with this C2 . If we throw away too many data points, we increase our C2 , and retry. If we
throw away too few, then we decrease our C2 and retry. Since increasing C2 strictly decreases the number
of points thrown away, and vice versa, we binary search over our choice of C2 until we reach something close
to our target accuracy. In our current implementation, we stop when the fraction of points we throw away
is between ε/2 and 3ε/2, or if we’ve binary searched for too long. We found that this heuristic drastically
improves our accuracy, and allows our algorithm to scale fairly smoothly from low to high dimension.

5

Experiments

We performed an empirical evaluation of the above algorithms on synthetic and real data sets with and
without synthetic noise. All experiments were done on a laptop computer with a 2.7 GHz Intel Core i5
CPU and 8 GB of RAM. The focus of this evaluation was on statistical accuracy, not time efficiency. In this
measure, our algorithm performs the best of all algorithms we tried. In all synthetic trials, our algorithm
consistently had the smallest error. In fact, in some of the synthetic benchmarks, our error was orders of
magnitude better than any other algorithms. In the semi-synthetic benchmark, our algorithm also (arguably)
performs the best, though there is no way to tell for sure, since there is no ground truth. We also note that
despite not optimizing our code for runtime, the runtime of our algorithm is always comparable, and in
many cases, better than the alternatives which provided comparable error. Code of our implementation is
available at https://github.com/hoonose/robust-filter.
7

1.5

excess `2 error

excess `2 error

0.15

1
0.5

0.1

0.05

0
100

200
300
dimension

100

400

Filtering
Sample mean w/ noise
RANSAC

200
300
dimension

400

LRVMean
Pruning
Geometric Median

Figure 1: Experiments with synthetic data for robust mean estimation: error is reported against dimension
(lower is better). The error is excess `2 error over the sample mean without noise (the benchmark). We
plot performance of our algorithm, LRVMean, empirical mean with noise, pruning, RANSAC, and geometric
median. On the left we report the errors achieved by all algorithms; however the latter four have much
larger error than our algorithm or LRVMean. On the right, we restrict our attention to only our algorithm
and LRVMean. Our algorithm has better error than all other algorithms.

100

excess Mahalanobis error

Skewed

100

excess Mahalanobis error

excess Mahalanobis error

excess Mahalanobis error

Isotropic
1.5
1
0.5
0
20

40
60
80
dimension

0.4

0.2

0
20

40
60
80
dimension

200

100

0
20

40
60
80
dimension

100

20

40
60
80
dimension

100

1

0.5

0

Filtering
Sample covariance w/ noise
RANSAC

LRVCov
Pruning

Figure 2: Experiments with synthetic data for robust covariance estimation: error is reported against dimension (lower is better). The error is excess Mahalanobis error over the sample covariance without noise (the
benchmark). We plot performance of our algorithm, LRVCov, empirical covariance with noise, pruning, and
RANSAC. We report two settings: one where the true covariance is isotropic (left column), and one where
the true covariance is very skewed (right column). In both, the latter three algorithms have substantially
larger error than ours or LRVCov. On the bottom, we restrict our attention to our algorithm and LRVCov.
The error achieved by LRVCov is quite good, but ours is better. In particular, our excess error is 4 orders of
magnitude smaller than LRVCov’s in high dimensions.

8

5.1

Synthetic Data

Experiments with synthetic data allow us to verify the error guarantees and the sample complexity rates
proven in Section 3 for unknown mean and unknown covariance. In both cases, the experiments validate the
accuracy and usefulness of our algorithm, almost exactly matching the best rate without noise.
Unknown mean The results of our synthetic mean experiment are shown in Figure 1. In the synthetic
mean experiment, we set ε = 0.1, and for dimension d = [100, 150, . . . , 400], we generate n = 10d
ε2 samples,
where a (1 − ε)-fraction come from N (µ, I), and an ε fraction come from a noise distribution. Our goal
is to produce an estimator which minimizes the `2 error the estimator has to the truth. As a baseline, we
compute the error that is achieved by only the uncorrupted sample points. This error will be used as the
gold standard for comparison, since in the presence of error, this is roughly the best one could do even if all
the noise points were identified exactly.†
On this data, we compared the performance of our Filter algorithm to that of (1) the empirical mean of
all the points, (2) a trivial pruning procedure, (3) the geometric median of the data, (4) a RANSAC-based
mean estimation algorithm, and (5) a recently proposed robust estimator for the mean due to [LRV16], which
we will call LRVMean. For (5), we use the implementation available in their Github.‡ In Figure 1, the x-axis
indicates the dimension of the experiment, and the y-axis measures the `2 error of our estimated mean minus
the `2 error of the empirical mean of the true samples from the Gaussian, i.e., the excess error induced over
the sampling error.
We tried various noise distributions, and found that the same qualitative pattern arose for all of them.
In the reported experiment, our noise distribution was a mixture of two binary product distributions, where
one had a couple of large coordinates (see Section B.1 for a detailed description). For all (nontrivial)
error distributions we tried, we observed that indeed the empirical mean, pruning, geometric median, and
RANSAC all have error which diverges as d grows, as the theory predicts. On the other hand, both our
algorithm and LRVMean have markedly smaller error as a function of dimension. Indeed, our algorithm’s
error is almost identical to that of the empirical mean of the uncorrupted sample points.
Unknown covariance The results of our synthetic covariance experiment are shown in Figure 2. Our
setup is similar to that for the synthetic mean. Since both our algorithm and LRVCov require access to
fourth moment objects, we ran into issues with limited memory on machines. Thus, we could not perform
experiments at as high a dimension as for the unknown mean setting, and we could not use as many
samples. We set ε = 0.05, and for dimension d = [10, 20, . . . , 100], we generate n = 0.5d
ε2 samples, where
a (1 − ε)-fraction come from N (0, Σ), and an ε fraction come from a noise distribution. We measure
distance in the natural affine invariant way, namely, the Mahalanobis distance induced by Σ to the identity:
b = kΣ−1/2 ΣΣ
b −1/2 − IkF . As explained above, this is the right affine-invariant metric for this problem.
err(Σ)
As before, we use the empirical error of only the uncorrupted data points as a benchmark.
On this corrupted data, we compared the performance of our Filter algorithm to that of (1) the empirical
covariance of all the points, (2) a trivial pruning procedure, (3) a RANSAC-based minimal volume ellipsoid
(MVE) algorithm, and (5) a recently proposed robust estimator for the covariance due to [LRV16], which
we will call LRVCov. For (5), we again obtained the implementation from their Github repository.
We tried various choices of Σ and noise distribution. Figure 2 shows two choices of Σ and noise. Again, the
x-axis indicates the dimension of the experiment and the y-axis indicates the estimator’s excess Mahalanobis
error over the sampling error. In the left figure, we set Σ = I, and our noise points are simply all located at
the all-zeros vector. In the right figure, we set Σ = I + 10e1 eT1 , where e1 is the first basis vector, and our
noise distribution is a somewhat more complicated distribution, which is similarly spiked, but in a different,
random, direction. We formally define this distribution in Section B.1. For all choices of Σ and noise we
tried, the qualitative behavior of our algorithm and LRVCov was unchanged. Namely, we seem to match the
empirical error without noise up to a very small slack, for all dimensions. On the other hand, the performance
of empirical mean, pruning, and RANSAC varies widely with the noise distribution. The performance of all
these algorithms degrades substantially with dimension, and their error gets worse as we increase the skew
of the underlying data. The performance of LRVCov is the most similar to ours, but again is worse by a large
† We

note that it is possible that an estimator may achieve slightly better error than this baseline.

‡ https://github.com/kal2000/AgnosticMean\AndCovarianceCode

9

constant factor. In particular, our excess risk was on the order of 10−4 for large d, for both experiments,
whereas the excess risk achieved by LRVCov was in all cases a constant between 0.1 and 2.
Discussion These experiments demonstrate that our statistical guarantees are in fact quite strong. In
particular, since our excess error is almost zero (and orders of magnitude smaller than other approaches),
this suggests that our sample complexity is indeed close to optimal, since we match the rate without noise,
and that the constants and logarithmic factors in the theoretical recovery guarantee are often small or
non-existent.
Original Data

Pruning Projection

0.2

-0.2
0.15

-0.1
0.1

0

0.05

0.1

0

-0.05

0.2

-0.1

0.3

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

-0.15

0.2

-0.2

The data projected onto the top two
directions of the original data set
without noise

-0.1

0

0.1

0.2

0.3

The data projected onto the top two directions
of the noisy data set after pruning

Filter Output

Filter Projection

-0.2

-0.2

-0.1

-0.1

0

0

0.1

0.1

0.2

0.2

0.3

0.3

The 0.1
filtered
set of
points
projected
onto-0.2the
0.15
0.05
0
-0.05
-0.1
-0.15
top two directions returned by the filter

The data
projected
onto
the -0.1
top two
0.1
0.05
0
-0.05
-0.15
directions returned by the filter

0.15

-0.2

Figure 3: Experiments with semi-synthetic data: given the real genetic data from [NJB+ 08], projected down
to 20-dimensions, and with added noise. The colors indicate the country of origin of the person, and match
the colors of the countries in the map of Europe in the center. Black points are added noise. The top left
plot is the original plot from [NJB+ 08]. We (mostly) recover Europe in the presence of noise whereas naive
methods do not.

5.2

Semi-synthetic Data

To demonstrate the efficacy of our method on real data, we revisit the famous study of [NJB+ 08]. In
this study, the authors investigated data collected as part of the Population Reference Sample (POPRES)
project. This dataset consists of the genotyping of thousands of individuals using the Affymetrix 500K single
nucleotide polymorphism (SNP) chip. The authors pruned the dataset to obtain the genetic data of over
1387 European individuals, annotated by their country of origin. Using principal components analysis, they
produce a two-dimensional summary of the genetic variation, which bears a striking resemblance to the map
of Europe.
10

Our experimental setup is as follows. While the original dataset is very high dimensional, we use a 20
dimensional version of the dataset as found in the authors’ GitHub§ . We first randomly rotate the data, as
then 20 dimensional data was diagonalized, and the high dimensional data does not follow such structure.
ε
fraction of points (so that they make up an ε-fraction of the final points).
We then add an additional 1−ε
These added points were discrete points, following a simple product distribution (see Section B.1 for full
details). We used a number of methods to obtain a covariance matrix for this dataset, and we projected the
data onto the top two singular vectors of this matrix. In Figure 3, we show the results when we compare
our techniques to pruning. In particular, our output was able to more or less reproduce the map of Europe,
whereas pruning fails to. In Section B.2, we also compare our result with a number of other techniques,
including those we tested against in the unknown covariance experiments, and other robust PCA techniques.
The only alternative algorithm which was able to produce meaningful output was LRVCov, which produced
output that was similar to ours, but which produced a map which was somewhat more skewed. We believe
that our algorithm produces the best picture.
In Figure 3, we also display the actual points which were output by our algorithm’s Filter. While it
manages to remove most of the noise points, it also seems to remove some of the true data points, particularly
those from Eastern Europe and Turkey. We attribute this to a lack of samples from these regions, and thus
one could consider them as outliers to a dataset consisting of Western European individuals. For instance,
Turkey had 4 data points, so it seems quite reasonable that any robust algorithm would naturally consider
these points outliers.
Discussion We view our experiments as a proof of concept demonstration that our techniques can be
useful in real world exploratory data analysis tasks, particularly those in high-dimensions. Our experiments
reveal that a minimal amount of noise can completely disrupt a data analyst’s ability to notice an interesting
phenomenon, thus limiting us to only very well-curated data sets. But with robust methods, this noise does
not interfere with scientific discovery, and we can still recover interesting patterns which otherwise would
have been obscured by noise.

Acknowledgments
We would like to thank Simon Du and Lili Su for helpful comments on a previous version of this work.

References
[AK95]

E. Amaldi and V. Kann. The complexity and approximability of finding maximum feasible
subsystems of linear relations. Theoretical Computer Science, 147:181–210, 1995.

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

[CEM+ 93] K. L. Clarkson, D. Eppstein, G. L. Miller, C. Sturtivant, and S.-H. Teng. Approximating
center points with iterated radon points. In Proceedings of the Ninth Annual Symposium on
Computational Geometry, SCG ’93, pages 91–98, New York, NY, USA, 1993. ACM.
[CGR16]

M. Chen, C. Gao, and Z. Ren. A general decision theory for huber’s ε-contamination model.
Electronic Journal of Statistics, 10(2):3752–3774, 2016.

[Cha04]

T. M. Chan. An optimal randomized algorithm for maximum tukey depth. In Proceedings of
the Fifteenth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 430–436,
2004.

[Che98]

Z. Chen. A note on bias robustness of the median. Statistics & probability letters, 38(4):363–368,
1998.

§ https://github.com/NovembreLab/Novembre_etal_2008_misc

11

[CLMW11] E. J. Candès, X. Li, Y. Ma, and J. Wright. Robust principal component analysis? J. ACM,
58(3):11, 2011.
[CSV17]

M. Charikar, J. Steinhardt, and G. Valiant. Learning from untrusted data. In Proceedings of
STOC’17, 2017.

[DBS17]

S. S. Du, S. Balakrishnan, and A. Singh. Computationally efficient robust estimation of sparse
functionals. In Proceedings of COLT’17, 2017.

[DG92]

D. L. Donoho and M. Gasko. Breakdown properties of location estimates based on halfspace
depth and projected outlyingness. Ann. Statist., 20(4):1803–1827, 12 1992.

[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.
Full version available at https://arxiv.org/pdf/1604.06443.pdf.
[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. Conference version available
at http://proceedings.mlr.press/v70/diakonikolas17a.html.
[DKK+ 18] I. Diakonikolas, G. Kamath, D. M. Kane, J. Li, A. Moitra, and A. Stewart. Robustly learning
a Gaussian: Getting optimal error, efficiently. In Proceedings of the 29th Annual ACM-SIAM
Symposium on Discrete Algorithms, SODA ’18, Philadelphia, PA, USA, 2018. SIAM.
[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.

[DKW56]

A. Dvoretzky, J. Kiefer, and J. Wolfowitz. Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator. Ann. Mathematical Statistics,
27(3):642–669, 1956.

[DL01]

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

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

[Hub97]

P. J. Huber. Robustness: Where are we now? Lecture Notes-Monograph Series, pages 487–498,
1997.

[JP78]

D. S. Johnson and F. P. Preparata. The densest hemisphere problem. Theoretical Computer
Science, 6:93–107, 1978.

[Li17]

J. Li. Robust sparse estimation tasks in high dimensions. In Proceedings of COLT’17, 2017.

12

[LRV16]

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

[Mas90]

P. Massart. The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality. Annals of Probability, 18(3):1269–1283, 1990.

[MS10]

G.L. Miller and D. Sheehy. Approximate centerpoints with proofs. Comput. Geom., 43(8):647–
654, 2010.

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

[Rou85]

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

[RS98]

P. J. Rousseeuw and A. Struyf. Computing location depth and regression depth in higher
dimensions. Statistics and Computing, 8(3):193–203, 1998.

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

[T+ 15]

J. A. Tropp et al. An introduction to matrix concentration inequalities. Foundations and Trends
in Machine Learning, 8(1-2):1–230, 2015.

[Tuk60]

J.W. Tukey. A survey of sampling from contaminated distributions. Contributions to probability
and statistics, 2:448–485, 1960.

[VAR09]

S. Van Aelst and P. Rousseeuw. Minimum volume ellipsoid. Wiley Interdisciplinary Reviews:
Computational Statistics, 1(1):71–82, 2009.

[Ver10]

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

[XCS10]

H. Xu, C. Caramanis, and S. Sanghavi. Robust pca via outlier pursuit. In Advances in Neural
Information Processing Systems, pages 2496–2504, 2010.

13

A

Omitted Details from Section 3

A.1

Robust Mean Estimation for Sub-Gaussian Distributions

In this section, we use our filter technique to give a near sample-optimal computationally efficient algorithm
to robustly estimate the mean of a sub-gaussian density with a known covariance matrix, thus proving
Theorem 3.1.
We emphasize that the algorithm and its analysis is essentially identical to the filtering algorithm given
in Section 8.1 of [DKK+ 16] for the case of a Gaussian N (µ, I). The only difference is a weaker definition
of the “good set of samples” (Definition A.4) and a simple concentration argument (Lemma A.5) showing
that a random set of uncorrupted samples of the appropriate size is good with high probability. Given these,
the analysis of this subsection follows straightforwardly from the analysis in Section 8.1 of [DKK+ 16] by
plugging in the modified parameters. For the sake of completeness, we provide the details below.
We start by formally defining sub-gaussian distributions:
Definition A.1. A distribution P on R with mean µ, is sub-gaussian with parameter ν > 0 if
EX∼P [exp(λ(X − µ))] ≤ exp(νλ2 /2)
for all λ ∈ R. A distribution P on Rd with mean vector µ is sub-gaussian with parameter ν > 0, if for all
unit vectors v, the one-dimensional random variable v · X, X ∼ P , is sub-gaussian with parameter ν.
We will use the following simple fact about the concentration of sub-gaussian random variables:
Fact A.2. If P is sub-gaussian on Rd with mean vector µ and parameter ν > 0, then for any unit vector
v ∈ Rd we have that PrX∼P [|v · (X − µ)| ≥ T ] ≤ exp(−T 2 /2ν).
The following theorem is a high probability version of Theorem 3.1:
Theorem A.3. Let G be a sub-gaussian distribution on Rd with parameter ν = Θ(1), mean µG , covariance
matrix I, and ε, τ > 0. Let S 0 be an ε-corrupted set of samples from G of size Ω((d/ε2 ) poly log(d/ετ )). There
exists an efficient algorithm that, on input
S 0 and ε > 0, returns a mean vector µ
b so that with probability at
p
G
least 1 − τ we have kb
µ − µ k2 = O(ε log(1/ε)).
P
1
Notation. We will denote µS = |S|
X∈S X and MS =
and modified sample covariance matrix of the set S.

1
|S|

P

X∈S (X

− µG )(X − µG )T for the sample mean

We start by defining our modified notion of good sample, i.e, a set of conditions on the uncorrupted set
of samples under which our algorithm will succeed.
Definition A.4. Let G be an identity covariance sub-gaussian in d dimensions with mean µG and covariance
matrix I and ε, τ > 0. We say that a multiset S of elements in Rd is (ε, τ )-good with respect to G if the
following conditions are satisfied:
p
(i) For all x ∈ S we have kx − µG k2 ≤ O( d log(|S|/τ )).
(ii) For every affine function L : Rd → R such that L(x) = v · (x − µG ) − T , kvk2 = 1, we have that
|PrX∈u S [L(X) ≥ 0] − PrX∼G [L(X) ≥ 0]| ≤ T 2 log dεlog( d ) .
(
ετ )
(iii) We have that kµS − µG k2 ≤ ε.
(iv) We have that kMS − Ik2 ≤ ε.
We show in the following subsection that a sufficiently large set of independent samples from G is (ε, τ )good (with respect to G) with high probability. Specifically, we prove:
Lemma A.5. Let G be sub-gaussian distribution with parameter ν = Θ(1) and with identity covariance,
and ε, τ > 0. If the multiset S is obtained by taking Ω((d/ε2 ) poly log(d/ετ )) independent samples from G,
it is (ε, τ )-good with respect to G with probability at least 1 − τ.
14

We require the following definition that quantifies the extent to which a multiset has been corrupted:
Definition A.6. Given finite multisets S and S 0 we let ∆(S, S 0 ) be the size of the symmetric difference of
S and S 0 divided by the cardinality of S.
The starting point of our algorithm will be a simple NaivePrune routine (Section 4.3.1 of [DKK+ 16])
that removes obvious outliers, i.e., points which are far from the mean. Then, we iterate the algorithm whose
performance guarantee is given by the following:
Proposition A.7. Let G be a sub-gaussian distribution on Rd with parameter ν = Θ(1), mean µG , covariance matrix I, ε > 0 be sufficiently small and τ > 0. Let S be an (ε, τ )-good
p set with respect to G. Let
S 0 be any multiset with ∆(S, S 0 ) ≤ 2ε and for any x, y ∈ S 0 , kx − yk2 ≤ O( d log(d/ετ )). There exists a
polynomial time algorithm Filter-Sub-Gaussian-Unknown-Mean that, given S 0 and ε > 0, returns one
of the following:
p
(i) A mean vector µ
b such that kb
µ − µG k2 = O(ε log(1/ε)).
def

d
(ii) A multiset S 00 ⊆ S 0 such that ∆(S, S 00 ) ≤ ∆(S, S 0 ) − ε/α, where α = d log(d/ετ ) log(d log( ετ
)).

We start by showing how Theorem A.3 follows easily from Proposition A.7.
Proof of Theorem A.3. By the definition of ∆(S, S 0 ), since S 0 has been obtained from S by corrupting an
ε-fraction of the points in S, we have that ∆(S, S 0 ) ≤ 2ε. By Lemma A.5, the set S of uncorrupted samples
is (ε, τ )-good with respect to G with probability at least 1p− τ. We henceforth condition on this event.
Since S is (ε, τ )-good, all x ∈ S have kx − µG k2 ≤ O( d log |S|/τ ). Thus, the NaivePrune procedure
does not remove from S 0 any member of S. Hence,
its output, S 00 , has ∆(S, S 00 ) ≤ ∆(S, S 0 ) and for any
p
00
x ∈ S , there p
is a y ∈ S with kxp
− yk2 ≤ O( d log |S|/τ ). By the triangle inequality, for any x, z ∈ S 00 ,
kx − zk2 ≤ O( d log |S|/τ ) = O( d log(d/ετ )).
Then, we iteratively apply the Filter-Sub-Gaussian-Unknown-Mean
procedure of Proposition A.7
p
until it terminates returning a mean vector µ with kb
µ − µG k2 = O(ε log(1/ε)). We claim that we need at
most O(α) iterations for this to happen. Indeed, the sequence of iterations results in a sequence of sets Si0 ,
so that ∆(S, Si0 ) ≤ ∆(S, S 0 ) − i · ε/α. Thus, if we do not output the empirical mean in the first 2α iterations,
in the next iteration there are no outliers left and the algorithm terminates outputting the sample mean of
the remaining set.
A.1.1

Algorithm Filter-Sub-Gaussian-Unknown-Mean: Proof of Proposition A.7

In this subsection, we describe the efficient algorithm establishing Proposition A.7 and prove its correctness.
0
Our algorithm calculates the empirical mean vector µS and empirical covariance matrix Σ. If the matrix Σ
0
has no large eigenvalues, it returns µS . Otherwise, it uses the eigenvector v ∗ corresponding to the maximum
0
magnitude eigenvalue of Σ and the mean vector µS to define a filter. Our efficient filtering procedure is
presented in detailed pseudocode below.
A.1.2

Proof of Correctness of Filter-Sub-Gaussian-Unknown-Mean

By definition, there exist disjoint multisets L, E, of points in Rd , where L ⊂ S, such that S 0 = (S \ L) ∪ E.
0
With this notation, we can write ∆(S, S 0 ) = |L|+|E|
|S| . Our assumption ∆(S, S ) ≤ 2ε is equivalent to |L|+|E| ≤
2ε · |S|, and the definition of S 0 directly implies that (1 − 2ε)|S| ≤ |S 0 | ≤ (1 + 2ε)|S|. Throughout the proof,
we assume that ε is a sufficiently small constant.
0
We define µG , µS , µS , µL , and µE to be the means of G, S, S 0 , L, and E, respectively.
Our analysis will make essential use of the following matrices:
• MS 0 denotes EX∈u S 0 [(X − µG )(X − µG )T ],
• MS denotes EX∈u S [(X − µG )(X − µG )T ],
• ML denotes EX∈u L [(X − µG )(X − µG )T ], and

15

Algorithm 2 Filter algorithm for a sub-gaussian with unknown mean and identity covariance
1: procedure Filter-Sub-Gaussian-Unknown-Mean(S 0 , ε, τ )
input: A multiset S 0 such that there exists an (ε, τ )-good S with ∆(S, S 0 ) ≤ 2ε
output: Multiset S 00 or mean vector µ
b satisfying Proposition A.7
0
2:
Compute the sample mean µS = EX∈u S 0 [X] and the sample covariance matrix Σ , i.e., Σ =
0
0
(Σi,j )1≤i,j≤d with Σi,j = EX∈u S 0 [(Xi − µSi )(Xj − µSj )].
3:
Compute approximations for the largest absolute eigenvalue of Σ − I, λ∗ := kΣ − Ik2 , and the
associated unit eigenvector v ∗ .
0
4:
if kΣ − Ik2 ≤ O(ε log(1/ε)), then return µS .
p
5:
Let δ := 3 εkΣ − Ik2 . Find T > 0 such that
h
i
0
ε
.
Pr 0 |v ∗ · (X − µS )| > T + δ > 8 exp(−T 2 /2ν) + 8 2
d
X∈u S
T log d log( ετ
)
0

return the multiset S 00 = {x ∈ S 0 : |v ∗ · (x − µS )| ≤ T + δ}.

6:

• ME denotes EX∈u E [(X − µG )(X − µG )T ].
Our analysis will hinge on proving the important claim that Σ − I is approximately (|E|/|S 0 |)ME . This
means two things for us. First, it means that if the positive errors align in some direction (causing ME to
have a large eigenvalue), there will be a large eigenvalue in Σ − I. Second, it says that any large eigenvalue
of Σ − I will correspond to an eigenvalue of ME , which will give an explicit direction in which many error
points are far from the empirical mean.
Useful Structural Lemmas. We begin by noting that we have concentration bounds on G and therefore,
on S due to its goodness.


d
G
2
Fact A.8. Let
 w ∈ R beGany unit
 vector, then2for any T > 0, εPrX∼G |w · (X − µ )| > T ≤ 2 exp(−T /2ν)
and PrX∈u S |w · (X − µ )| > T ≤ 2 exp(−T /2ν) + T 2 log d log( d ) .
(
ετ )
Proof. The first line is Fact A.2, and the former follows from it using the goodness of S.
By using the above fact, we obtain the following simple claim:
Claim A.9. Let w ∈ Rd be any unit vector, then for any T > 0, we have that:
0

0

Pr [|w · (X − µS )| > T + kµS − µG k2 ] ≤ 2 exp(−T 2 /2ν).

X∼G

and

0

0

Pr [|w · (X − µS )| > T + kµS − µG k2 ] ≤ 2 exp(−T 2 /2ν) +

X∈u S

0

T2

ε
.
d
log d log( ετ
)

0

Proof. This follows from Fact A.8 upon noting that |w · (X − µS )| > T + kµS − µG k2 only if |w · (X − µG )| >
T.
We can use the above facts to prove concentration bounds for L. In particular, we have the following
lemma:
Lemma A.10. We have that kML k2 = O (log(|S|/|L|) + ε|S|/|L|).
Proof. Since L ⊆ S, for any x ∈ Rd , we have that
|S| · Pr (X = x) ≥ |L| · Pr (X = x) .
X∈u S

X∈u L

(1)

Since ML is a symmetric matrix, we have kML k2 = maxkvk2 =1 |v T ML v|. So, to bound kML k2 it suffices to
bound |v T ML v| for unit vectors v. By definition of ML , for any v ∈ Rd we have that
|v T ML v| = EX∈u L [|v · (X − µG )|2 ].
16

For unit vectors v, the RHS is bounded from above as follows:
Z ∞




G 2
EX∈u L |v · (X − µ )| = 2
Pr |v · (X − µG )| > T T dT
X∈
L
u
0
Z √
O(

d log(d/ετ ))

Pr [|v · (X − µG )| > T ]T dT

=2
0

Z

√

O(

X∈u L
d log(d/ετ ))

≤2
0

Z 4√ν log(|S|/|L|)







|S|
G
min 1,
· Pr |v · (X − µ )| > T T dT
|L| X∈u S

T dT
0

√

O(

Z
+ (|S|/|L|)

√

4

d log(d/ετ ))



exp(−T 2 /2ν) +

ν log(|S|/|L|)


ε
 T dT
d
T 2 log d log( ετ )

 log(|S|/|L|) + ε · |S|/|L| ,
where the second line follows from the fact that kvk2 = 1, L ⊂ S, and S satisfies condition (i) of Definition A.4; the third line follows from (1); and the fourth line follows from Fact A.8.
As a corollary, we can relate the matrices MS 0 and ME , in spectral norm:
Corollary A.11. We have that MS 0 − I = (|E|/|S 0 |)ME + O(ε log(1/ε)), where the O(ε log(1/ε)) term
denotes a matrix of spectral norm O(ε log(1/ε)).
Proof. By definition, we have that |S 0 |MS 0 = |S|MS − |L|ML + |E|ME . Thus, we can write
MS 0 = (|S|/|S 0 |)MS − (|L|/|S 0 |)ML + (|E|/|S 0 |)ME
= I + O(ε) + O(ε log(1/ε)) + (|E|/|S 0 |)ME ,
where the second line uses the fact that 1 − 2ε ≤ |S|/|S 0 | ≤ 1 + 2ε, the goodness of S (condition (iv) in
Definition A.4), and Lemma A.10. Specifically, Lemma A.10 implies that (|L|/|S 0 |)kML k2 = O(ε log(1/ε)).
Therefore, we have that
MS 0 = I + (|E|/|S 0 |)ME + O(ε log(1/ε)) ,
as desired.
We now establish a similarly useful bound on the difference between the mean vectors:
p
p
0
0
E
G
Lemma A.12. We have that µS − µG = (|E|/|S
p |)(µ − µ ) + O(ε log(1/ε)), where the O(ε log(1/ε))
term denotes a vector with `2 -norm at most O(ε log(1/ε)).
Proof. By definition, we have that
0

|S 0 |(µS − µG ) = |S|(µS − µG ) − |L|(µL − µG ) + |E|(µE − µG ).
Since S is a good set, by condition (iii) of Definition A.4, we have kµS −µG k2 = O(ε). Since 1−2ε ≤ |S|/|S 0 | ≤
1 + 2ε, it follows that (|S|/|S 0 |)kµS − µG k2 = O(ε). Using the valid inequality
kML k2 ≥ kµL − µG k22 and

p
p
L
G
Lemma A.10, we obtain that kµ − µ k2 ≤ O
log(|S|/|L|) + ε|S|/|L| . Therefore,


p
p
p
(|L|/|S 0 |)kµL − µG k2 ≤ O (|L|/|S|) log(|S|/|L|) + ε|L|/|S| = O(ε log(1/ε)) .
In summary,

p
0
µS − µG = (|E|/|S 0 |)(µE − µG ) + O(ε log(1/ε)) ,

as desired. This completes the proof of the lemma.

17

By combining the above, we can conclude that Σ−I is approximately proportional to ME . More formally,
we obtain the following corollary:
Corollary A.13. We have Σ − I = (|E|/|S 0 |)ME + O(ε log(1/ε)) + O(|E|/|S 0 |)2 kME k2 , where the additive
terms denote matrices of appropriately bounded spectral norm.
0

0

Proof. By definition, we can write Σ − I = MS 0 − I − (µS − µG )(µS − µG )T . Using Corollary A.11 and
Lemma A.12, we obtain:
Σ − I = (|E|/|S 0 |)ME + O(ε log(1/ε)) + O((|E|/|S 0 |)2 kµE − µG k22 ) + O(ε2 log(1/ε))
= (|E|/|S 0 |)ME + O(ε log(1/ε)) + O(|E|/|S 0 |)2 kME k2 ,
where the second line follows from the valid inequality kME k2 ≥ kµE − µG k22 . This completes the proof.
0

Case of Small Spectral Norm. We are now ready to analyze the case that the mean vector µS is returned
def
by the algorithm in Step 4. In this case, we have that λ∗ = kΣ − Ik2 = O(ε log(1/ε)). Hence, Corollary
A.13 yields that
(|E|/|S 0 |)kME k2 ≤ λ∗ + O(ε log(1/ε)) + O(|E|/|S 0 |)2 kME k2 ,
which in turns implies that
(|E|/|S 0 |)kME k2 = O(ε log(1/ε)) .
On the other hand, since kME k2 ≥ kµE − µG k22 , Lemma A.12 gives that
p
p
p
0
kµS − µG k2 ≤ (|E|/|S 0 |) kME k2 + O(ε log(1/ε)) = O(ε log(1/ε)).
This proves part (i) of Proposition A.7.
Case of Large Spectral Norm. We next show the correctness of the algorithm when it returns a filter in
Step 5.
def
We start by proving that if λ∗ = kΣ − Ik2 > Cε log(1/ε), for a sufficiently large universal constant C,
then a value T satisfying the condition in Step 5 exists. We first note that that kME k2 is appropriately
large. Indeed, by Corollary A.13 and the assumption that λ∗ > Cε log(1/ε) we deduce that
(|E|/|S 0 |)kME k2 = Ω(λ∗ ) .

(2)

Moreover, using the inequality kME k2 ≥ kµE − µG k22 and Lemma A.12 as above, we get that
p
p
0
kµS − µG k2 ≤ (|E|/|S 0 |) kME k2 + O(ε log(1/ε)) ≤ δ/2 ,

(3)

p
def √
where we used the fact that δ = ελ∗ > C 0 ε log(1/ε).
Suppose for the sake of contradiction that for all T > 0 we have that
h
i
0
ε
.
Pr 0 |v ∗ · (X − µS )| > T + δ ≤ 8 exp(−T 2 /2ν) + 8 2
d
X∈u S
T log d log( ετ
)
Using (3), we obtain that for all T > 0 we have that
Pr

X∈u

S0

 ∗

|v · (X − µG )| > T + δ/2 ≤ 8 exp(−T 2 /2ν) + 8

ε
.
d
T 2 log d log( ετ
)

(4)

Since E ⊆ S 0 , for all x ∈ Rd we have that |S 0 | PrX∈u S 0 [X = x] ≥ |E| PrY ∈u E [Y = x]. This fact combined
with (4) implies that for all T > 0
!
 ∗

ε
 .
Pr |v · (X − µG )| > T + δ/2  (|S 0 |/|E|) exp(−T 2 /2ν) + 2
(5)
d
X∈u E
T log d log( ετ
)

18

We now have the following sequence of inequalities:
Z ∞




kME k2 = EX∈u E |v ∗ · (X − µG )|2 = 2
Pr |v ∗ · (X − µG )| > T T dT
X∈u E
0
Z O(√d log(d/ετ ))


Pr |v ∗ · (X − µG )| > T T dT
=2
X∈u E
0


Z O(√d log(d/ετ ))


|S 0 |
min 1,
· Pr 0 |v ∗ · (X − µG )| > T T dT
≤2
|E| X∈u S
0
Z O(√d log(d/ετ )) 
Z 4√ν log(|S 0 |/|E|)+δ
0
T dT + (|S |/|E|) √
exp(−T 2 /2ν) +

0
0

4

0

ν log(|S |/|E|)+δ


ε

T dT
d
)
T 2 log d log( ετ

0

2

 log(|S |/|E|) + δ + O(1) + ε · |S |/|E|
 log(|S 0 |/|E|) + ελ∗ + ε · |S 0 |/|E| .
Rearranging the above, we get that
(|E|/|S 0 |)kME k2  (|E|/|S 0 |) log(|S 0 |/|E|) + (|E|/|S 0 |)ελ∗ + ε = O(ε log(1/ε) + ε2 λ∗ ).
Combined with (2), we obtain λ∗ = O(ε log(1/ε)), which is a contradiction if C is sufficiently large. Therefore,
it must be the case that for some value of T the condition in Step 5 is satisfied.
The following claim completes the proof:
def

d
Claim A.14. Fix α = d log(d/ετ ) log(d log( ετ
)). We have that ∆(S, S 00 ) ≤ ∆(S, S 0 ) − 2ε/α .

Proof. Recall that S 0 = (S \ L) ∪ E, with E and L disjoint multisets such that L ⊂ S. We can similarly write
S 00 = (S \ L0 ) ∪ E 0 , with L0 ⊇ L and E 0 ⊂ E. Since
∆(S, S 0 ) − ∆(S, S 00 ) =

|E \ E 0 | − |L0 \ L|
,
|S|

it suffices to show that |E \ E 0 | ≥ |L0 \ L| + ε|S|/α. Note that |L0 \ L| is the number of points rejected by the
filter that lie in S ∩ S 0 . Note that the fraction of elements of S that are removed to produce S 00 (i.e., satisfy
0
|v ∗ · (x p
− µS )| > T + δ) is at most 2 exp(−T 2 /2ν) + ε/α. This follows from Claim A.9 and the fact that
T = O( d log(d/ετ )).
Hence, it holds that |L0 \ L| ≤ (2 exp(−T 2 /2ν) + ε/α)|S|. On the other hand, Step 5 of the algorithm
ensures that the fraction of elements of S 0 that are rejected by the filter is at least 8 exp(−T 2 /2ν) + 8ε/α).
Note that |E \ E 0 | is the number of points rejected by the filter that lie in S 0 \ S. Therefore, we can write:
|E \ E 0 | ≥ (8 exp(−T 2 /2ν) + 8ε/α)|S 0 | − (2 exp(−T 2 /2ν) + ε/α)|S|
≥ (8 exp(−T 2 /2ν) + 8ε/α)|S|/2 − (2 exp(−T 2 /2ν) + ε/α)|S|
≥ (2 exp(−T 2 /2ν) + 3ε/α)|S|
≥ |L0 \ L| + 2ε|S|/α ,
where the second line uses the fact that |S 0 | ≥ |S|/2 and the last line uses the fact that |L0 \ L|/|S| ≤
2 exp(−T 2 /2ν) + ε/α. Noting that log(d/ετ ) ≥ 1, this completes the proof of the claim.
A.1.3

Proof of Lemma A.5

Proof. Let N = Ω((d/ε2 ) poly log(d/ετ ))pbe the number of samples drawn from G. For (i), the probability
that a coordinate of a sample is at least 2ν log(N d/3τ ) is at most p
τ /3dN by Fact A.2. By a union bound,
the probability that
all
coordinates
of
all
samples
are
smaller
than
2ν log(N d/3τ ) is at least 1 − τ /3. In
p
p
this case, kxk2 ≤ 2νd log(N d/3τ ) = O( dν log(N ν/τ )).
After translating by µG , we note that (iii) follows immediately from Lemmas 4.3 of [DKK+ 16] and (iv)
follows from Theorem 5.50 of [Ver10], as long as N = Ω(ν 4 d log(1/τ )/ε2 ), with probability at least 1 − τ /3.
It remains to show that, conditioned on (i), (ii) holds with probability at least 1 − τ /3.
19

To simplify some expressions, let δ := ε/(log(d log d/ετ )) and R = C
that for all unit vectors v and all 0 ≤ T ≤ R that

p

d log(|S|/τ ). We need to show

Pr [|v · (X − µG )| > T ] − Pr [|v · (X − µG ) > T ≥ 0] ≤

X∈u S

X∼G

δ
.
T2

(6)

Firstly, we show that for all unit vectors v and T > 0
Pr [|v · (X − µG )| > T ] − Pr [|v · (X − µG )| > T ≥ 0] ≤

X∈u S

X∼G

δ
10ν ln(1/δ)

with probability at least 1 − τ /6. Since the VC-dimension of the set of all halfspaces is d + 1, this follows
2
from the VC inequality [DL01], since
pwe have more than Ω(d/(δ/(10ν log(1/δ)) ) samples. We thus only
need to consider the case when T ≥ 10ν ln(1/δ).
p
Lemma A.15. For any fixed unit vector v and T > 10ν ln(1/δ), except with probability exp(−N δ/(6Cν)),
we have that
δ
Pr [|v · (X − µG )| > T ] ≤
,
X∈u S
CT 2
where C = 8.
Proof. Let E be the event that |v · (X − µG )| > T . Since G is sub-gaussian, Fact A.2 yields that PrG [E] =
PrY ∼G [|v · (X − µG )| > T ] ≤ exp(−T 2 /(2ν)). Note that, thanks to our assumption on T , we have that
T ≤ exp(T 2 /(4ν))/2C, and therefore T 2 PrG [E] ≤ exp(−T 2 /(4ν))/2C ≤ δ/2C.
Consider ES [exp(t2 /(3ν) · N PrS [E])]. Each individual sample Xi for 1 ≤ i ≤ N , is an independent copy
of Y ∼ G, and hence:
"
#

 2

 2 X
n
T
T
ES exp
· N Pr[E]
= ES exp
·
1Xi ∈E )
S
3ν
3ν
i=1
"
#
 2 X
N
n
Y
T
=
EXi exp
·
1Xi ∈E )
3ν
i=1
i=1
N

 2
T
Pr[G] + 1
= exp
3ν G



N
(a)
T2
≤ exp
+1
6ν
(b)

≤ (1 + δ 5/3 )N

(c)

≤ exp(N δ 5/3 ) ,

where (a) follows from sub-gaussianity, (b) follows from our choice of T , and (c) comes from the fact that
1 + x ≤ ex for all x.
Hence, by Markov’s inequality, we have




δ
δN
5/3
Pr Pr[E] ≥
≤ exp N δ
−
S
CT 2
3C
= exp(N δ(δ 2/3 − 1/(3C))) .
Thus, if δ is a sufficiently small constant and C is sufficiently large, this yields the desired bound.
Now let C be a 1/2-cover in Euclidean distance
for the set of unit vectors of size 2O(d) . By a union bound,
p
0
0
for all v ∈ C and T a power of 2 between 4ν ln(1/δ) and R, we have that
Pr [|v 0 · (X − µG )| > T 0 ] ≤

X∈u S

20

δ
8T 2

except with probability
2O(d) log(R) exp(−N δ/6Cν) = exp (O(d) + log log R − N δ/6Cν) ≤ τ /6 .
p
However, for any unit vector v and 4ν ln(1/δ) ≤ T ≤ R, there is a v 0 ∈ C and such a T 0 such that for all
x ∈ Rd , we have |v · (X − µG )| ≥ |v 0 · (X − µG )|/2, and so |v 0 · (X − µG )| > 2T 0 implies |v 0 · (X − µG )| > T.
Then, by a union bound, (6) holds simultaneously for all unit vectors v and all 0 ≤ T ≤ R, with
probability a least 1 − τ /3. This completes the proof.

A.2

Robust Mean Estimation Under Second Moment Assumptions

In this section, we use our filtering technique to give a near sample-optimal computationally efficient algorithm to robustly estimate the mean of a density with a second moment assumption. We show:
Theorem A.16. Let P be a distribution on Rd with unknown mean vector µP and unknown covariance
matrix ΣP  I. Let S be an ε-corrupted set of samples from P of size Θ((d/ε)
√ log d). Then there exists an
algorithm that given S, with probability 2/3, outputs µ
b with kb
µ − µP k2 ≤ O( ε) in time poly(d/ε).
Note that Theorem 3.2 follows straightforwardly from the above (divide every sample by σ, run the
algorithm of Theorem A.16, and multiply its output by σ).
As usual in our filtering framework, the algorithm will iteratively look at the top eigenvalue and eigenvector of the sample covariance matrix and return the sample mean if this eigenvalue is small (Algorithm
3). The main difference between this and the filter algorithm for the sub-gaussian case is how we choose
the threshold for the filter. Instead of looking for a violation of a concentration inequality, here we will
choose a threshold at random (with a bias towards higher thresholds). The reason is that, in this setting,
the variance in the direction we look for a filter in needs to be a constant multiple larger – instead of the
typical Ω̃(ε) relative for the sub-gaussian case. Therefore, randomly choosing a threshold weighted towards
higher thresholds suffices to throw out more corrupted samples than uncorrupted samples in expectation.
Although it is possible to reject many good samples this way, the algorithm still only rejects a total of O(ε)
samples with high probability.
We would like our good set of samples to have mean close to that of P and bounded variance in all
directions. This motivates the following definition:
P
Definition A.17. We call a set S ε-good for a distribution
√ P withS mean µ and covariance ΣP  I if the
S
S
S
P
mean µ and covariance Σ of S satisfy kµ − µ k2 ≤ ε and kΣ k2 ≤ 2.

However, since we have no assumptions about higher moments, it may be be possible for outliers to affect
our sample covariance too much. Fortunately, such outliers have small probability and do not contribute too
much to the mean, so we will later reclassify them as errors.
Lemma A.18. Let S be N = Θ((d/ε) log d) samples drawn from P . Then, with probability at least 9/10, a
random X ∈u S satisfies
√
(i) kES [X] − µP k2 ≤ ε/3,
h
p i
(ii) PrS kX − µP k2 ≥ 80 d/ε ≤ ε/160,
(iii)
(iv)

h
ES (X − µP ) · 1kX−µP k

√

2 ≤80

i
d/ε

h
ES (X − µP )(X − µP )T · 1kX−µP k

√

≤

ε/3, and

2

√

2 ≤80

i
d/ε

≤ 3/2.
2

Proof. For (i), note that
ES [kE[X] − µP k22 ] =

X

2
ES [(E[X]i − µP
i ) ] ≤ d/N ≤ ε/360 ,

i

and so by Markov’s inequality, with probability at least 39/40, we have kE[X] − µP k22 ≤ ε/9.
21

For (ii), similarly to (i), note that
E[kY − µP k22 ] =

X 

2
E (Yi − µP
≤d,
i )
i

p
for Y ∼ P . By Markov’s inequality, Pr[kY − µP k2 ≥ 80 d/ε] ≤ ε/160 with probability at least 39/40.
For (iii), let ν = EX∼P [X · 1kX−µP k ≤80√d/ε ] be the true mean of the distribution when we condition on
2
p
the event that kX − µP k2 ≤ 80 d/ε. By the same argument as (i), we know that
h
i
√
EX∈u S X · 1kX−µP k ≤80√d/ε − ν ≤ ε/9 ,
2

2

with probability at least 39/40. Thus it suffices to show that ν − µP · 1kX−µP k

√

2 ≤80

so, it suffices to show that for all unit vectors v ∈ Rd , we have
D
E
√
v, ν − µP · 1kX−µP k ≤80√d/ε < ε/10 .

d/ε 2

≤

√

ε/10. To do

2

Observe that for any such v, we have
D
E
i
h
v, µP · 1kX−µP k ≤80√d/ε − ν = EX∼P v, X − µP · 1kX−µP k ≤80√d/ε
2
2
r
(a)
p
≤ EX∼P [hv, X − µP i2 ] Pr [kX − µP k2 ≥ 80 d/ε]
X∼P
r
h
p i
(b)
= v T ΣP v · Pr kX − µP k2 ≥ 80 d/ε
X∼P

(c)

≤

√

ε/10 ,

where (a) follows from Cauchy-Schwarz, and (b) follows from the definition of the covariance, and (c) follows
from the assumption that ΣP  I and from Markov’s inequality.
For (iv), we require the following Matrix Chernoff bound:
Lemma A.19 (Part of Theorem 5.1.1 of [T+ 15]). Consider
a sequence of d×d positive semi-definite random
P
matrices Xk with kXk k2 ≤ L for all k. Let µmax = k k E[Xk ]k2 . Then, for θ > 0,
#
"
X
Xk
≤ (eθ − 1)µmax /θ + L log(d)/θ ,
E
k

2

and for any δ > 0,
"
Pr

#
X

max

≥ (1 + δ)µ

Xk

k

max

≤ d(eδ /(1 + δ)1+δ )µ

/L

.

2

We apply this lemma with Xk = (xk − µP )(xk − µP )T 1kx

√

P
k −µ k2 ≤80

d/ε

for {x1 , . . . , xN } = S. Note that

kXk k2 ≤ (80)2 d/ε = L and that µmax ≤ N kΣP k2 ≤ N .
Suppose that µmax ≤ N/80. Then, taking θ = 1, we have
E[

X

] ≤ (e − 1)N/80 + O(d log(d)/ε) .

Xk

k

2

P
By Markov’s inequality, except with probability 39/40, we have k k Xk k2 ≤ N + O(d log(d)/ε) ≤ 3N/2,
for N a sufficiently high multiple of d log(d)/ε.
Suppose that µmax ≥ N/80, then we take δ = 1/2 and obtain
#
"
X
max
Pr
Xk ≥ 3µ
2 ≤ d(e3/2 /(5/2)3/2 )N ε/20d .
k

2

22

P
For N a sufficiently high multiple of d log(d)/ε, wePget that Pr[k k Xk k2 ≥ 3µmax /2] ≤ 1/40. Since
µmax ≤ N , we have with probability at least 39/40, k k Xk k2 ≤ 3N/2.
P
Noting that k k Xk k2 /N = kE[1kX−µP k ≤80√d/ε (X − µP )(X − µP )T ]k2 , we obtain (iv). By a union
2

bound, (i)-(iv) all hold simultaneously with probability at least 9/10.
Now we can get a 2ε-corrupted good set from an ε-corrupted set of samples satisfying Lemma A.18, by
reclassifying outliers as errors:
Lemma A.20. Let S = R ∪ E \ L, where R is a set of N = Θ(d log d/ε) samples drawn from P and E and
L are disjoint sets with |E|, |L| ≤ ε. Then, with probability 9/10, we can also write S = G ∪ E 0 \ L0 , where
G ⊆ R is ε-good, L0 ⊆ L and E 0 ⊆ E 0 has |E 0 | ≤ 2ε|S|.
p
Proof. Let G = {x ∈ R : kxk2 ≤ 80 d/ε}. Condition on the event that R satisfies Lemma A.18. By Lemma
A.18, this occurs with probability at least 9/10.
Since R satisfies (ii) of Lemma A.18, |G| − |R| ≤ ε|R|/160 ≤ ε|S|. Thus, E 0 = E ∪ (R \ G) has
|E 0 | ≤ 3ε/2. Note that (iv) of Lemma A.18 for R in terms of G is exactly |G|kΣG k2 /|R| ≤ 3/2, and so
kΣG k2 ≤ 3|R|/(2|G|) ≤ 2.
√
It remains to check that kµG − µP k2 ≤ ε. We have
i
h
|G| · µG − |G| · µP 2 = |R| · EX∼u R (X − µP ) · 1kX−µP k ≤80√d/ε
2
2
√
≤ |R| · ε/3 ,
where the last line follows from (iii) of Lemma A.18. Since we argued above that |R|/|G| ≥ 2/3, dividing
this expression by |G| yields the desired claim.

Algorithm 3 Filter under second moment assumptions
1: function FilterUnder2ndMoment(S)
2:
Compute µS , ΣS , the mean and covariance matrix of S.
3:
Find the eigenvector v ∗ with highest eigenvalue λ∗ of ΣS .
4:
if λ∗ ≤ 9 then
5:
return µS
6:
else
7:
Draw Z from the distribution on [0, 1] with probability density function 2x.
8:
Let T = Z max{|v ∗ · x − µS | : x ∈ S}.
9:
Return the set S 0 = {x ∈ S : |v ∗ · (X − µS )| < T }.
An iteration of FilterUnder2ndMoment may throw out more samples from G than corrupted samples.
However, in expectation, we throw out many more corrupted samples than from the good set:
Proposition A.21. If we run FilterUnder2ndMoment on a set S = G ∪ E \ L for some ε-good√ set G
and disjoint E, L with |E| ≤ 2ε|S|, |L| ≤ 9ε|S|, then either it returns µS with kµS − µP k2 ≤ O( ε), or
else it returns a set S 0 ⊂ S with S 0 = G ∪ E 0 \ L0 for disjoint E 0 and L0 . In the latter case we have
EZ [|E 0 | + 2|L0 |] ≤ |E| + 2|L|.
For D ∈ {G, E, L, S}, let µD be the mean of D and MD be the matrix EX∈u D [(X − µS )(X − µS )T ].
p
Lemma A.22. If G is an ε-good set with x ≤ 40 d/ε for x ∈ S ∪ G, then kMG k2 ≤ 2kµG − µS k22 + 2 .
Proof. For any unit vector v, we have
v T MG v = EX∈u G [(v · (X − µS ))2 ]
= EX∈u G [(v · (X − µG ) + v · (µP − µG ))2 ]
= v T ΣG v + (v · (µG − µS ))2
≤ 2 + 2kµG − µS k22 .

23

Lemma A.23. We have that |L|kML k2 ≤ 2|G|(1 + kµG − µS k22 ) .
Proof. Since L ⊆ G, for any unit vector v, we have
|L|v T ML v = |L|EX∈u L [(v · (X − µS ))2 ]
≤ |G|EX∈u G [(v · (X − µS ))2 ]
≤ 2|G|(1 + kµG − µS k22 ) .

Lemma A.24. kµG − µS k2 ≤

p
√
2εkMS k2 + 12 ε.

Proof. We have that |E|ME ≤ |S|MS + |L|ML and so
|E|kME k2 ≤ |S|kMS k2 + 2|G|(1 + kµG − µS k22 ) .
By Cauchy Schwarz, we have that kME k2 ≥ kµE − µS k22 , and so
q
p
|E|kµE − µS k2 ≤ |S|kMS k2 + 2|G|(1 + kµG − µS k22 ) .
By Cauchy-Schwarz and Lemma A.23, we have that
q
p
p
|L|kµL − µS k2 ≤ |L|kML k2 ≤ 2|G|(1 + kµG − µS k22 ) .
Since |S|µS = |G|µG + |E|µE − |L|µL and |S| = |G| + |E| − |L|, we get
|G|(µG − µS ) = |E|(µE − µS ) − |L|(µE − µS ) .
Substituting into this, we obtain
q
q
|G|kµG − µS k2 ≤ |E||S|kMS k2 + 2|E||G|(1 + kµG − µS k22 ) + 2|L||G|(1 + kµG − µS k22 ) .
√

√

√
x + y, we have
p
p
p
|G|kµG − µS k2 ≤ |E||S|kMS k2 + ( 2|E||G| + 2|L||G|)(1 + kµG − µS k2 ) .

Since for x, y > 0,

x+y ≤

Since ||G| − |S|| ≤ ε|S| and |E| ≤ 2ε|S|, |L| ≤ 9ε|S|, we have
p
√
kµG − µS k2 ≤ 2εkMS k2 + (6 ε)(1 + kµG − µS k2 ) .
√
Moving the kµG − µS k2 terms to the LHS, using 6 ε ≤ 1/2, gives
p
√
kµG − µS k2 ≤ 2εkMS k2 + 12 ε .

Since λ∗ = kMS k2 , the correctness if we return the empirical mean is immediate.
√
Corollary A.25. If λ∗ ≤ 9, we have that kµG − µS k2 = O( ε).
From now on, we assume λ∗ > 9. In this case we have kµG − µS k22 ≤ O(ελ∗ ). Using Lemma A.22, we
have
kMG k2 ≤ 2 + O(ελ∗ ) ≤ 2 + λ∗ /5
for sufficiently small ε. Thus, we have that
v ∗T MS v ∗ ≥ 4v ∗T MG v ∗ .

(7)

Now we can show that in expectation, we throw out many more corrupted points from E than from G\L:
24

Lemma A.26. Let S 0 = G ∪ E 0 \ L0 for disjoint E 0 , L0 be the set of samples returned by the iteration. Then
we have EZ [|E 0 | + 2|L0 |] ≤ |E| + 2|L|.
Proof. Let a = maxx∈S |v ∗ · x − µS |. Firstly, we look at the expected number of samples we reject:


EZ [|S 0 |] − |S| = EZ |S| Pr [|X − µS | ≥ aZ]
X∈u S

1

Z
= |S|

Pr

X∈u S
Z0 a

 ∗

|v · (X − µS )| ≥ ax 2xdx



Pr |v ∗ · (X − µS )| ≥ T (2T /a)dT
0 X∈u S


= |S|EX∈u S (v ∗ · (X − µS ))2 /a

= |S|

= (|S|/a) · v ∗T MS v ∗ .
Next, we look at the expected number of false positive samples we reject, i.e., those in L0 \ L.




0
S
EZ [|L |] − |L| = EZ (|G| − |L|) Pr
|X − µ | ≥ T
X∈u G\L


∗
S
≤ EZ |G| Pr [|v · (X − µ )| ≥ aZ]
X∈u G

Z
= |G|
= |G|

1

Pr [|v ∗ · (X − µS )| ≥ ax]2x dx

X∈u G
Z0 a

Pr [|v ∗ · (X − µS )| ≥ T ](2T /a) dT

X∈u G
Z0 ∞

Pr [|v ∗ · (X − µS )| ≥ T ](2T /a) dT


= |G|EX∈u G (v ∗ · (X − µS ))2 /a
≤ |G|

0

X∈u G

= (|G|/a) · v ∗T MG v ∗ .
Using (7), we have |S|v ∗T MS v ∗ ≥ 4|G|v ∗T MG v ∗ and so EZ [S 0 ] − S ≥ 3(EZ [L0 ] − L). Now consider that
|S 0 | = |G| + |E 0 | − |L0 | = |S| − |E| + |E 0 | + |L| − |L0 |, and thus |S 0 | − |S| = |E| − |E 0 | + |L0 | − |L|. This yields
that |E| − EZ [|E 0 |] ≥ 2(EZ [L0 ] − L), which can be rearranged to EZ [|E 0 | + 2|L0 |] ≤ |E| + 2|L|.
∗
S
P
Proof
√ of Proposition A.21. If λ ≤ 9, then we return the mean in Step 5, and by Corollary A.25, kµ −µ k2 ≤
O( ε).
If λ∗ > 9, then we return S 0 . Since at least one element of S has |v ∗ · X| = maxx∈S |v ∗ · X|, whatever
value of Z is drawn, we still remove at least one element, and so have S 0 ⊂ S. By Lemma A.26, we have
EZ [|E 0 | + 2|L0 |] ≤ |E| + 2|L|.

Proof of Theorem A.16. Our input is a set S of N = Θ((d/ε) log d) ε-corrupted samples so that with probability 9/10, S is a 2ε-corrupted set of ε-good samples for P by Lemmas A.18 and A.20. We have a set
S = G ∪ E 0 \ L, where G0 is an ε-good set, |E| ≤ 2ε, and |L| ≤ ε. Then, we iteratively apply FilterUnder2ndMoment until it outputs an approximation to the mean. Since each iteration removes a sample,
this must happen within N iterations. The algorithm takes at most poly(N, d) = poly(d, 1/ε) time.
As long as we can
√ show that the conditions of Proposition A.21 hold in each iteration, it ensures that
kµS − µP k2 ≤ O( ε). However, the condition that |L| ≤ 9ε|S| need not hold in general. Although in
expectation we reject many more samples in E than G, it is possible that we are unlucky and reject many
samples in G, which could make L large in the next iteration. Thus, we need a bound on the probability
that we ever have |L| > 9ε.
We analyze the following procedure: We iteratively run FilterUnder2ndMoment starting with a set
Si ∪ Ei \ Li of samples with S0 = S and producing a set Si+1 = G ∪ Ei+1 \ Li+1 . We stop if we output
an approximation to the mean or if |Li+1 | ≥ 13ε|S|. Since we do now always satisfy the conditions of
Proposition A.21, this gives that EZ [|Ei+1 | + |Li+1 |] = |Ei | + 2|Li |. This expectation is conditioned on the
25

state of the algorithm after previous iterations, which is determined by Si . Thus, if we consider the random
variables Xi = |Ei | + 2|Li |, then we have E[Xi+1 |Si ] ≤ Xi , i.e., the sequence Xi is a sub-martingale with
respect to Xi . Using the convention that Si+1 = Si , if we stop in less than i iterations, and recalling that we
always stop in N iterations, the algorithm fails if and only if |LN | > 9ε|S|. By a simple induction or standard
results on sub-martingales, we have E[XN ] ≤ X0 . Now X0 = |E0 | + 2|L0 | ≤ 3ε|S|. Thus, E[XN ] ≤ 3ε|S|. By
Markov’s inequality, except with probability 1/6, we have XN ≤ 18ε|S|. In this case, |LN | ≤ XN /2 ≤ 9ε|S|.
Therefore, the probability that we ever have |Li | > 9ε is at most 1/6.
By a union bound, the probability that the uncorrupted samples satisfy Lemma A.18 and Proposition
A.21 applies to every iteration is at least 9/10−1/6
≥ 2/3. Thus, with at least 2/3 probability, the algorithm
√
outputs a vector µ
b with kb
µ − µP k2 ≤ O( ε).

A.3

Robust Covariance Estimation

In this subsection, we give a near sample-optimal efficient robust estimator for the covariance of a zero-mean
Gaussian density, thus proving Theorem 3.3. Our algorithm is essentially identical to the filtering algorithm
given in Section 8.2 of [DKK+ 16]. As in Section A.1 the only difference is a weaker definition of the “good
set of samples” (Definition A.27) and a concentration argument (Lemma A.28) showing that a random set
of uncorrupted samples of the appropriate size is good with high probability. Given these, the analysis of
this subsection follows straightforwardly from the analysis in Section 8.2 of [DKK+ 16] by plugging in the
modified parameters.
The algorithm Filter-Gaussian-Unknown-Covariance to robustly estimate the covariance of a mean
0 Gaussian in [DKK+ 16] is as follows:
Algorithm 4 Filter algorithm for a Gaussian with unknown covariance matrix.
1: procedure Filter-Gaussian-Unknown-Covariance(S 0 , ε, τ )
input: A multiset S 0 such that there exists an (ε, τ )-good set S with ∆(S, S 0 ) ≤ 2ε
output: Either a set S 00 with ∆(S, S 00 ) < ∆(S, S 0 ) or the parameters of a Gaussian G0 with dT V (G, G0 ) =
O(ε log(1/ε)).
Let C > 0 be a sufficiently large universal constant.
2:
Let Σ0 be the matrix EX∈u S 0 [XX T ] and let G0 be the mean 0 Gaussian with covariance matrix Σ0 .
3:
if there is any x ∈ S 0 so that xT (Σ0 )−1 x ≥ Cd log(|S 0 |/τ ) then
4:
return S 00 = S 0 − {x : xT (Σ0 )−1 x ≥ Cd log(|S 0 |/τ )}.
5:
Compute an approximate eigendecomposition of Σ0 and use it to compute Σ0−1/2
6:
Let x(1) , . . . , x(|S 0 |) be the elements of S 0 .
⊗2
7:
For i = 1, . . . , |S 0 |, let y(i) = Σ0−1/2 x(i) and z(i) = y(i)
.
0
P
|S |
[
[T
0
T
8:
Let TS 0 = −I I + (1/|S |) i=1 z(i) z(i) .
9:
Approximate the top eigenvalue λ∗ and corresponding unit eigenvector v ∗ of TS 0 ..
10:
Let p∗ (x) = √12 ((Σ0−1/2 x)T v ∗] (Σ0−1/2 x) − tr(v ∗] ))
11:
12:
13:
14:

if λ∗ ≤ (1 + Cε log2 (1/ε))QG0 (p∗ ) then
return G0
Let µ be the median value of p∗ (X) over X ∈ S 0 .
Find a T ≥ C 0 so that
Pr (|p∗ (X) − µ| ≥ T + 4/3) ≥ Tail(T, d, ε, τ )

X∈u S 0

15:

return S 00 = {X ∈ S 0 : |p∗ (X) − µ| < T }.

In [DKK+ 16], we take Tail(T, d, ε, τ ) = 12 exp(−T ) + 3ε/(d log(N/τ ))2 , where N = Θ((d log(d/ετ ))6 /ε2 )
is the number of samples we took there.
To get a near sample-optimal algorithms, we will need a weaker definition of a good set. To use this, we
will need to weaken the tail bound in the algorithm to Tail(T, d, ε, τ ) = ε/(T 2 log2 (T )), when T ≥ 10 log(1/ε).
For T ≤ 10 log(1/ε), we take Tail(T, d, ε, τ ) = 1 so that we always choose T ≥ 10 log(1/ε). It is easy to show
26

that the integrals of this tail bound used in the proofs of Lemma 8.19 and Claim 8.22 of [DKK+ 16] have
similar bounds. Thus, our analysis here will sketch that these tail bounds hold for a set of Ω(d2 log5 (d/ετ )/ε2 )
samples from the Guassian.
Firstly, we state the new, weaker, definition of a good set:
Definition A.27. Let G be a Gaussian in Rd with mean 0 and covariance Σ. Let ε > 0 be sufficiently small.
We say that a multiset S of points in Rd is ε-good with respect to G if the following hold:
√
1. For all x ∈ S, xT Σ−1 x < d + O( d log(d/ε)).
2. We have that kΣ−1/2 Cov(S)Σ−1/2 − IkF = O(ε).
3. For all even degree-2 polynomials p, we have that Var(p(S)) = Var(p(G))(1 + O(ε)).
4. For p an even degree-2 polynomial with E[p(G)] = 0 and Var(p(G)) = 1, and for any T > 10 log(1/ε)
we have that
Pr (|p(x)| > T ) ≤ ε/(T 2 log2 (T )).
x∈u S

It is easy to see that the algorithm and analysis of [DKK+ 16] can be pushed through using the above
weaker definition. That is, if S is a good set, then G can be recovered to Õ(ε) error from an ε-corrupted
version of S. Our main task will be to show that random sets of the appropriate size are good with high
probability.
Proposition A.28. Let N be a sufficiently large constant multiple of d2 log5 (d/ε)/ε2 . Then a set S of N
independent samples from G is ε-good with respect to G with high probability.
Proof. First, note that it suffices to prove this when G = N (0, I).
Condition 1 follows by standard concentration bounds on kxk22 .
Condition 2 follows by estimating the entry-wise error between Cov(S) and I.
Condition 3 is slightly more involved. Let {pi } be an orthonormal basis for the set of even, degree-2,
mean-0 polynomials with respect to G. Define the matrix Mi,j = Ex∈u S [pi (x)pj (x)] − δi,j . This condition is
equivalent to kM k2 = O(ε). Thus, it suffices to show that for every v with kvk2 = 1 that v T M v = O(ε). It
O(d2 )
actually
. For each v, let
P suffices to consider a cover of such v’s. Note that this cover will be of size 2
pv = i vi pi . We need to show that Var(pv (S)) = 1 + O(ε). We can show this happens with probability
2
1 − 2−Ω(d ) , and thus it holds for all v in our cover by a union bound.
Condition 4 is substantially the most difficult of these conditions to prove. Naively, we would want to
find a cover of all possible p and all possible T , and bound the probability that the desired condition fails.
Unfortunately, the best a priori bound on Pr(|p(G)| > T ) are on the order of exp(−T ). As our cover would
2
need to be of size 2d or so, to make this work with T = d, we would require on the order of d3 samples in
order to make this argument work.
However, we will note that this argument is sufficient to cover the case of T < 10 log(1/ε) log2 (d/ε).
Fortunately, most such polynomials p satisfy much better tail bounds. Note that any even, mean zero
polynomial p can be written in the form p(x) = xT Ax−tr(A) for some matrix A. We call A the associated matrix to p. We note by the Hanson-Wright inequality that Pr(|p(G)| > T ) = exp(−Ω(min((T /kAkF )2 , T /kAk2 ))).
Therefore, the tail bounds above are only as bad as described when A has a single large eigenvalue. To take
advantage of this, we will need to break p into parts based on the size of its eigenvalues. We begin with a
definition:
Definition A.29. Let Pk be the set of even, mean-0, degree-2 polynomials, so that the associated matrix A
satisfies:
1. rank(A) ≤ k
√
2. kAk2 ≤ 1/ k.
√
√
Note that for p ∈ Pk that |p(x)| ≤ |x|2 / k + k.
Importantly, any polynomial can be written in terms of these sets.

27

Lemma A.30. Let p be an even, degree-2 polynomial with E[p(G)] = 0, Var(p(G)) = 1. Then if t = blog2 (d)c,
it is possible to write p = 2(p1 + p2 + . . . + p2t + pd ) where pk ∈ Pk .
Proof. Let A be the associated matrix to p. Note that kAkF = Var p = 1. Let Ak be the matrix corresponding
to the top k eigenvalues of A. We now let p1 be the polynomial associated to A1 /2, p2 be associated to
(A2 − A1 )/2, p4 be associated to (A4 − A2 )/2, and so on. It is clear that p = 2(p1 + p2 + . . . + p2t + pd ).
It is also clear that the matrix
associated to pk has rank at most k. If the matrix associated to pk had an
√
k,
it
would
need to be the case that the k/2nd largest eigenvalue of A had size at
eigenvalue
more
than
1/
√
least 2/ k. This is impossible since the sum of the squares of the eigenvalues of A is at most 1.
This completes our proof.
We will also need covers of each of these sets Pk .
Lemma A.31. For each k, there exists a set Ck ⊂ Pk so that
1. For each p ∈ Pk there exists a q ∈ Ck so that kp(G) − q(G)k2 ≤ (ε/d)2 .
2. |Ck | = 2O(dk log(d/ε)) .
√
Pk
Proof. We note that any such p is associated to a matrix A of the form A = i=1 λi vi viT , for λi ∈ [0, 1/ k]
Pk
and vi orthonormal. It suffices to let q correspond to the matrix A0 = i=1 µi wi wiT for with |λi −µi | < (ε/d)3
and |vi − wi | < (ε/d)3 for all i. It is easy to let µi and wi range over covers of the interval and the sphere
with appropriate errors. This gives a set of possible q’s of size 2O(dk log(d/ε)) as desired. Unfortunately, some
of these q will not be in Pk as they will have eigenvalues that are too large. However, this is easily fixed by
replacing each such q by the closest element of Pk . This completes our proof.
We next will show that these covers are sufficient to express any polynomial.
Lemma A.32. Let p be an even degree-2 polynomial with E[p(G)] = 0 and Var(p(G)) = 1. It is possible to
write p as a sum of O(log(d)) elements of some Ck plus another polynomial of L2 norm at most ε/d.
Proof. Combining the above two lemmas we have that any such p can be written as
t

p = (q1 + p1 ) + (q2 + p2 ) + . . . (q2t + p2t ) + (qd + pd ) = q1 + q2 + . . . + q 2 + q d + p0 ,
where qk above is in Ck and kpk (G)k2 < (ε/d)2 . Thus, p0 = p1 + p2 + . . . + p2t + pd has kp0 (G)k2 ≤ (ε/d).
This completes the proof.
p
The key observation now is that if |p(x)| ≥ T for kxk2 ≤ d/ε, then writing p = q1 + q2 + q4 + . . . + qd + p0
as above, it must be the case that |qk (x)| > (T −1)/(2 log(d)) for some k. Therefore, to prove our main result,
it suffices to show that, with high probability over the choice of S, for any T ≥ 10 log(1/ε) log2 (d/ε) and
any q ∈ Ck for some k, that Prx∈u S (|q(x)| > T /(2 log(d))) < ε/(2T 2 log2 (T ) log(d)). Equivalently, it suffices
2
to show that for T ≥ 10 log(1/ε) log(d/ε) it holds Prx∈u S (|q(x)| > T /(2 log(d))) < ε/(2T 2 log2 (T ) logp
(d)).
Note that this holds automatically for T > (d/ε), as p(x) cannot possibly be that large for kxk2 ≤ d/ε.
Furthermore, note that losing a constant factor in the probability, it suffices to show this only for T a power
of 2.
√
Therefore, it suffices to show for every k ≤ d, every q ∈ Ck and every d/ kε  T  log(1/ε) log(d/ε)
that with probability at least 1 − 2−Ω(dk log(d/ε)) over the choice of S we have that Prx∈u S (|q(x)| > T ) 
ε/(T 2 log4 (d/ε)). However, by the Hanson-Wright inequality, we have that
√
Pr(|q(G)| > T ) = exp(−Ω(min(T 2 , T k))) < (ε/(T 2 log4 (d/ε)))2 .
Therefore, by Chernoff bounds, the probability that more than a ε/(T 2 log4 (d/ε))-fraction of the elements
of S satisfy this property is at most
√
√
exp(−Ω(min(T 2 , T k))|S|ε/(T 2 log4 (d/ε))) = exp(−Ω(|S|ε/(log4 (d/ε)) min(1, k/T )))
≤ exp(−Ω(|S|ε2 /(log4 (d/ε))k/d))
≤ exp(−Ω(dk log(d/ε))) ,
28

as desired.
This completes our proof.

29

B
B.1

Omitted Details from Section 5
Full description of the distributions for experiments

Here we formally describe the distributions we used in our experiments. In all settings, our goal was to find
noise distributions so that noise points were not “obvious” outliers, in the sense that there is no obvious
pointwise pruning process which could throw away the noise points, which still gave the algorithms we tested
the most difficulty. We again remark that while other algorithms had varying performances depending on
the noise distribution, it seemed that the performance of ours was more or less unaffected by it.
Distribution for the synthetic mean experiment Our uncorrupted points were generated by N (µ, I),
where µ is the all-ones vector. Our noise distribution is given as
N=

1
1
Π1 + Π2 ,
2
2

where Π1 is the product distribution over the hypercube where every coordinate is 0 or 1 with probability
1/2, and Π2 is a product distribution where the first coordinate is ether 0 or 12 with equal probability, the
second coordinate is −2 or 0 with equal probability, and all remaining coordinates are zero.
Distribution for the synthetic covariance experiment For the isotropic synthetic covariance experiment, our uncorrupted points were generated by N (0, I), and the noise points were all zeros. For the skewed
synthetic covariance experiment, our uncorrupted points were generated by N (0, I + 100e1 eT1 ), where e1 is
the first unit vector, and our noise points were generated as follows: we took a fixed random rotation of
points of the form Yi ∼ Π, where Π is a product distribution whose first d/2 coordinates are each uniformly
selected from {−0.5, 0, 0.5}, and whose next d/2 − 1 coordinates are each 0.8 × Ai , where for each coordinate
i, Ai is an independent random integer between −2 and 2, and whose last coordinate is a uniformly random
integer between [−100, 100].
Setup for the semi-synthetic geographic experiment We took the 20 dimensional data from [NJB+ 08],
which was diagonalized, and randomly rotated it. This was to simulate the higher dimensional case, since
the singular vectors that [NJB+ 08] obtained did not seem to be sparse or analytically sparse. Our noise was
distributed as Π, where Π is a product distribution whose first d/2 coordinates are each uniformly random
integers between 0 and 2 and whose last d/2 coordinates are each uniformly randomly either 2 or 3, all scaled
by a factor of 1/24.

B.2

Comparison with other robust PCA methods on semi-synthetic data

In addition to comparing our results with simple pruning techniques, as we did in Figure 3 in the main
text, we also compared our algorithm with implementations of other robust PCA techniques from the literature with accessible implementations. In particular, we compared our technique with RANSAC-based
techniques, LRVCov, two SDPs ([CLMW11, XCS10]) for variants of robust PCA, and an algorithm proposed
by [CLMW11] to speed up their SDP based on alternating descent. For the SDPs, since black box methods were too slow to run on the full data set (as [CLMW11] mentions, black-box solvers for the SDPs are
impractical above perhaps 100 data points), we subsample the data, and run the SDP on the subsampled
data. For each of these methods, we ran the algorithm on the true data points plus noise, where the noise
was generated as described above. We then take the estimate of the covariance it outputs, and project the
data points onto the top two singular values of this matrix, and plot the results in Figure 4.
Similar results occurred for most noise patterns we tried. We found that only our algorithm and LRVCov
were able to reasonably reconstruct Europe, in the presence of this noise. It is hard to judge qualitatively
which of the two maps generated is preferable, but it seems that ours stretches the picture somewhat less
than LRVCov.

30

Original Data

Filter Projection

-0.2

-0.2

-0.1

-0.1

0

0

0.1

0.1

0.2

0.2

0.3

0.3

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

0.15

0.1

0.05

RANSAC Projection

0.2

0

-0.05

-0.1

-0.15

-0.2

0.1

0.15

0.2

0.1

0.15

0.2

LRV Projection
-0.2

0.15
-0.1
0.1

0.05

0

0

0.1

-0.05
0.2
-0.1
0.3
-0.15
-0.2

-0.1

0

0.1

0.2

0.3

-0.15

-0.1

-0.05

0.2

0

0.05

CLMW SDP Projection

CLMW ADMM Projection
0.3

0.15

0.2
0.1

0.05

0.1

0

0

-0.05

-0.1
-0.1

-0.2
-0.15
-0.2

-0.1

0

0.1

0.2

-0.15

0.3

-0.1

-0.05

0

0.05

XCS Projection
0.3

0.2

0.1

0

-0.1

-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

Figure 4: Comparison with other robust methods on the Europe semi-synthetic data. From left to right,
top to bottom: the original projection without noise, what our algorithm recovers, RANSAC, LRVCov, the
ADMM method proposed by [CLMW11], the SDP proposed by [XCS10] with subsampling, and the SDP
proposed by [CLMW11] with subsampling.
31

