arXiv:1711.11581v2 [] 23 Dec 2017

Outlier-robust moment-estimation
via sum-of-squares
Pravesh K. Kothari∗

David Steurer†

December 27, 2017

Abstract
We develop eﬃcient algorithms for estimating low-degree moments of unknown distributions in the presence of adversarial outliers. The guarantees of our algorithms improve in many cases signiﬁcantly over the best previous ones, obtained in recent works
of [DKK+ 16, LRV16, CSV17]. We also show that the guarantees of our algorithms match
information-theoretic lower-bounds for the class of distributions we consider. These improved
guarantees allow us to give improved algorithms for independent component analysis and
learning mixtures of Gaussians in the presence of outliers.
Our algorithms are based on a standard sum-of-squares relaxation of the following
conceptually-simple optimization problem: Among all distributions whose moments are
bounded in the same way as for the unknown distribution, ﬁnd the one that is closest in
statistical distance to the empirical distribution of the adversarially-corrupted sample.

∗ Princeton

University and Institute for Advanced Study
Zürich. Much of this work was done while at Cornell University and the Institute for Advanced Study.
Supported by a Microsoft Research Fellowship, a Alfred P. Sloan Fellowship, NSF awards, and the Simons Collaboration
for Algorithms and Geometry.
† ETH

Contents
1 Introduction
1.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 Robust identifiability of low-degree moments

1
5
10

3 Preliminaries
13
3.1 Pseudo-distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Sum-of-squares proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4 Moment estimation algorithm
16
4.1 Putting things together . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5 Certifiably subgaussian distributions

22

6 Applications
25
6.1 Outlier-Robust Independent Component Analysis . . . . . . . . . . . . . . . . . . . . 25
6.2 Outlier-Robust Mixtures of Gaussians . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
7 Information Theoretic Lower Bounds

31

References

32

A Sum-of-squares toolkit

36

1 Introduction
We consider the problem of outlier-robust parameter estimation: We are given independent draws
x1 , . . . , x n from an unknown distribution D over d and the goal is to estimate parameters of the
distribution, e.g., its mean or its covariance matrix. Furthermore, we require that the estimator is
outlier-robust: even if an ε-fraction of the draws are corrupted by an adversary, the estimation error
should be small.
These kind of estimators have been studied extensively in statistics (under the term robust
statistics) [Tuk75, MMY06, Hub11, HRRS11]. However, many robust estimators coming out of this
research eﬀort are computationally eﬃcient only for low-dimensional distributions (because the
running time is say exponential in the dimension d) [Ber06].
A recent line of research developed the ﬁrst robust estimators for basic parameter estimation
problems (e.g., estimating the mean and covariance matrix of a Gaussian distribution) that are
computationally eﬃcient for the high-dimensional case, i.e., the running time is only polynomial
in the dimension d [DKK+ 16, LRV16, CJN17, CSV17].
Our work continues this line of research. We design eﬃcient algorithms to estimate lowdegree moments of distributions. Our estimators succeed under signiﬁcantly weaker assumptions
about the unknown distribution D, even for the most basic tasks of estimating the mean and
covariance matrix of D. For example, in order to estimate the mean of D (in an appropriate
norm) our algorithms do not need to assume that D is Gaussian or has a covariance matrix with
small spectral norm in contrast to assumptions of previous works. Similary, our algorithms for
estimating covariance matrices work, unlike previous algorithms, for non-Gaussian distributions
and distributions that are not (aﬃne transformations of) product distributions.
Besides these qualitative diﬀerences, our algorithms also oﬀer quantitative improvements. In
particular, for the class of distributions we consider, the guarantees of our algorithms—concretely,
the asymptotic behavior of the estimation error as a function of the fraction ε of corruptions—match,
for the ﬁrst time in this generality, information-theoretic lower bounds.
Outlier-robust method of moments Our techniques for robust estimation of mean vectors and
covariance matrices extend in a natural way to higher-order moment tensors. This fact allows
us to turn many non-outlier-robust algorithms in a black-box way into outlier-robust algorithms.
The reason is that for many parameter estimation problems the best known algorithms in terms
of provable guarantees are based on the method of moments, which means that they don’t require
direct access to a sample from the distribution but instead only to its low-degree moments (e.g.
[Pea94, KMV10, MV10]). (Often, a key ingredient of these algorithms in the high-dimensional
setting is tensor decomposition [MR05, AFH+12, AGH+ 14, HK13b, BCMV14, BKS15, GM15, MSS16a].)
If there were no outliers, we could run these kinds of algorithms on the empirical moments of the
observed sample. However, in the presence of outliers, this approach fails dramatically because
even a single outlier can have a huge eﬀect on the empirical moments. Instead, we apply method-ofmoment-based algorithms on the output of our outlier-robust moment estimators. Following this
strategy, we obtain new outlier-robust algorithms for independent component analysis (even if the
underlying unknown linear transformation is ill-conditioned) and mixtures of spherical Gaussians
1

(even if the number of components of the mixture is large and the means are not separated).
Estimation algorithms from identifiability proofs Our algorithms and their analysis follow a
recent paradigm for computationally-eﬃcient provable parameter estimation that has been developed in the context of the sum-of-squares method. We say that a parameter estimation problem
satisﬁes identiﬁability if it is information-theoretically possible to recover the desired parameter
from the available data (disregarding computational eﬃciency). The key idea of this paradigm
is that a proof of identiﬁability can be turned into an eﬃcient estimation algorithm if the proof
is captured by a low-complexity proof system like sum-of-squares. Many estimation algorithms
based on convex relaxations, in particular sum-of-squares relaxations, can be viewed as following
this paradigm (e.g., compressed sensing and matrix completion [CRTV05, CR09, Gro11, Rec11]).
Moreover, this paradigm has been used, with a varying degree of explicitness, in order to design
a number of recent algorithms based on sum-of-squares for unsupervised learning, inverse, and
estimation problems like overcomplete tensor decomposition, sparse dictionary learning, tensor
completion, tensor principal component analysis [HSS15, BKS15, BM16, PS17, MSS16a].
In many of these settings, including ours, the proof of identiﬁability takes a particular form:
given a (corrupted) sample X ⊆ d from a distribution D and a candidate parameter θ̂ that is close
to true parameter θ of D, we want to be able to eﬃciently certify that the candidate parameter θ̂
is indeed close to θ. (This notion of certiﬁcate is similar to the one for NP.) Following the above
paradigm, if the certiﬁcation in this identiﬁability proof can be accomplished by a low-degree
sum-of-squares proof, we can derive an eﬃcient estimation algorithm (that computes an estimate
θ̂ just given the sample X).
Next we describe the form of our identiﬁability proofs for outlier-robust estimation.
Identifiability in the presence of outliers Suppose we are given an ε-corrupted sample X 
{x1 , . . . , x n } ⊆ d of an unknown distribution D and our goal is to estimate the mean µ ∈ d of D.
To hope for the mean to be identiﬁable, we need to make some assumptions about D. Otherwise,
D could for example have probability mass 1 − ε on 0 and probability mass ε on µ/ε. Then, an
adversary can erase all information about the mean µ in an ε-corrupted sample from D (because
only an ε fraction of the draws from D carry information about µ). Therefore, we will assume that
D belongs to some class of distributions C (known to the algorithm). Furthermore, we will assume
that the class C is deﬁned by conditions on low-degree moments so that if we take a large enough
(non-corrupted) sample X 0 from D, the uniform distribution over X 0 also belongs to the class C
with high probability. (We describe the conditions that deﬁne C in a paragraph below.)
With this setup in place, we can describe our robust identiﬁability proof: it consists of a
(multi-)set of vectors X ′  {x1′ , . . . , x ′n } ⊆ d that satisﬁes two conditions:
1. x ′i  x i for all but an ε fraction of the indices i ∈ [n],
2. the uniform distribution over X ′ is in C.
Note that given X and X ′, we can eﬃciently check the above conditions, assuming that the conditions on the low-degree moments that deﬁne C are eﬃciently checkable (which they will be).
2

Also note that the above notion of proof is complete in the following sense: if X is indeed an
ε-corruption of a typical1 sample X 0 from a distribution D ∈ C, then there exists a set X ′ that
satisﬁes the above conditions, namely the uncorrupted sample X 0 .
In Section 2, we show that the above notion of proof is also sound: if X is indeed an ε-corruption
of a (typical) sample X 0 from a distribution D ∈ C and X ′ satisﬁes the above conditions, then the
Í
empirical mean µ′ : n1 ni1 x ′i of the uniform distribution over X ′ is close to µ. We can rephrase
this soundness as the following concise mathematical statement, which we prove in Section 2: if
D and D ′ are two distributions in C that have small statistical distance, then their means are close
to each other (and their higher-order moments are close to each other as well).
Furthermore, the above soundness is captured by a low-degree sum-of-squares proof (using
for example the sum-of-squares version of Hölder’s inequality); this fact is the basis of our eﬃcient
algorithm for outlier-robust mean estimation (see Section 4).
Outlier-robustness and (certifiable) subgaussianity To motivate the aforementioned conditions
we impose on the distributions, consider the following scenario: we are given an ε-corrupted
sample of a distribution D over  and our goal is to estimate its variance σ2 up to a constant factors.
To hope for the variance to be identiﬁable, we need to rule out for example that D outputs 0
with probability 1 − ε and a Gaussian N(0, σ2 /ε 2 ) with probability ε. Otherwise, an adversary can
remove all information about σ by changing an ε fraction of the draws in a sample
p from D. If we look
k
1/k
at the low-degree moments of this distribution D, we see that ( D(x) x ) ≈ k/ε 1−2/k ( D(x) x 2 )1/2 .
So for large enough k (i.e., k ≈ log(1/ε)), the ratio between the L k norm of D and the L 2 norm of
√
D exceeds that of a Gaussian by a factor of roughly 1/ ε. In order to rule out this example, we
impose the following condition on the low-degree moments of D, and we show that this condition
is enough to estimate the variance of a distribution D over  up to constant factors given an
ε-corrupted sample,



D(x)

(x − µ D )

k

 1/k

6

√

Ck ·



2

D(x)

(x − µ D )

 1/2

for C > 0 and even k ∈  with Ck · ε 1−2/k ≪ 1 . (1.1)

Here, µ D is the mean of the distribution D.
In the high-dimensional setting, a natural idea is to impose the condition Eq. (1.1) for every
direction u ∈ d ,
 1/2


 1/k √
2
k
.
(1.2)
hx − µ D , ui
hx − µ D , ui
6 Ck ·
D(x)

D(x)

In Section 2, we show that this condition is indeed enough to ensure identiﬁability and that it is
information-theoretically possible to estimate the covariance matrix ΣD of D up to constant factors
(in the Löwner order sense) assuming again that Ck · ε 1−2/k ≪ 1. Unfortunately, condition Eq. (1.2)
is unlikely to be enough to guarantee an eﬃcient estimation algorithm. The reason is that Eq. (1.2)
might hold for the low-degree moments of some distribution D but every proof of this fact requires
1The sample X 0 of D should be typical in the sense that the empirical low-degree moments of the sample X 0 are
close to the (population) low-degree moments of the distribution D. If the sample is large enough (polynomial in the
dimension), this condition is satisﬁed with high probability.

3

exponential size. (This phenomenon is related to the fact that ﬁnding a vector u ∈ d that violates
Eq. (1.2) is an NP-hard problem in general.)
To overcome this source of intractability, we require that inequality Eq. (1.2) is not only true but
also has a low-degree sum-of-squares proof. With this additional condition, we give a polynomialtime algorithm to estimate the covariance matrix ΣD of D up to constant factors (in the Löwner
order sense) assuming again that Ck · ε 1−2/k ≪ 1.
Definition 1.1 (Certiﬁable subgaussianity of low-degree moments). A distribution2 D over d
with mean µ is (k, ℓ)-certiﬁably subgaussian with parameter C > 0 if there for every positive integer
k ′ 6 k/2, there exists a degree-ℓ sum-of-squares proof3 of the degree-2k ′ polynomial inequality
over the unit sphere,

∀u ∈ d−1 .

′

D(x)



hx − µ, ui 2k 6 C · k ′

D(x)

hx − µ, ui 2

 k′

(1.3)

In Section 5, we show that a wide range of distributions are certiﬁably subgaussian and that
many operations on distributions preserve this property. In particular, (aﬃne transformations
of) products of scalar-valued subgaussian distributions and mixtures thereof satisfy this property. In a subsequent work, Kothari and Steinhardt [KS17] show that certiﬁable subgaussianity
holds for distributions that satisfy a Poincaré inequality (which includes all strongly log-concave
distributions).
Since any valid polynomial inequality over the sphere has a sum-of-squares proof if we allow
the degree to be large enough and the inequality has positive slack, we say that D is (k, ∞)-certiﬁably
subgaussian with parameter C > 0 if the inequalities Eq. (1.3) hold for all k ′ 6 k/2. In most cases
we consider, a sum-of-squares proof of degree ℓ  k is enough. In this case, we just say that D
is k-certiﬁably subgaussian. Observe that k-certiﬁable subgaussianity only restricts moments up
to degree k of the underlying distribution. In this sense, it is a much weaker assumption than the
usual notion of subgaussianity, which imposes Gaussian-like upper bounds on all moments.
Sum-of-squares and quantifier alternation Taking together the above discussion of robust identiﬁability proofs and certiﬁable subgaussianity, our approach to estimate the low-degree moments of a certiﬁable subgaussian distribution is the following: given an ε-corrupted sample
X  {x1 , . . . , x n } ⊆ d from D, we want to ﬁnd a (multi-)set of vectors X ′  {x1′ , . . . , x ′n } ⊆ d
such that the uniform distributions over X and X ′ are ε-close in statistical distance and the uniform
distribution over X ′ is certiﬁably subgaussian.
It is straightforward to formulate these conditions as a system E of polynomial equations
over the reals such that x1′ , . . . , x ′n are (some of the) variables (see Section 4). Our outlier-robust
2We emphasize that our notion of certiﬁable subgaussianity is a property only of the low-degree moments of a
distribution and, in this way, signiﬁcantly less restrictive then the usual notion of subgaussianity (which restricts all
moments of a distribution).
3 In the special case of proving that p(u) > 0 holds for every vector u ∈  d−1 , a degree-ℓ sum-of-squares proof
consists of a polynomial q(u) with deg q 6 ℓ − 2 and polynomials r1 (u), . . . , r t (u) with deg r τ 6 ℓ/2 such that p 
q · (kuk 2 − 1) + r12 + · · · + r t2 . See Section 3 for a more general deﬁnition.

4

moment estimation proceeds by solving a standard sum-of-squares relaxation of this system E.
The solution to this relaxation can be viewed as a pseudo-distribution D ′ over vectors x1′ , . . . , x ′n that
satisﬁes the system of equations E. (This pseudo-distribution behaves in many ways like a classical
probability distribution over vectors x ′1, . . . , x ′n that satisfy the equations E. See Section 3 for a
deﬁnition.) Our moment estimation algorithm simply outputs the expected empirical moments
˜ D′ (x′ ,...,x′ ) 1 Ín (x ′ )⊗r with respect to the pseudo-distribution D ′.
i1 i
n n
1
We remark that previous work on computationally-eﬃciently outlier-robust estimation also
used convex optimization techniques albeit in diﬀerent ways. For example, Diakonikolas et al.
solve an implicitly-deﬁned convex optimization problem using a customized separation oracle
[DKK+ 16]. (Their optimization problem is implicit in the sense that it is deﬁned in terms of the
uncorrupted sample which we do not observe.)
An unusual feature of the aforementioned system of equations E is that it also includes variables
for the sum-of-squares proof of the inequality Eq. (1.2) because we want to restrict the search to
those sets X ′ such that the uniform distribution X ′ is certiﬁably subgaussian. It is interesting to
note that in this way we can use sum-of-squares as an approach to solve ∃ ∀-problems as opposed
to just the usual ∃-problems. (The current problem is an ∃ ∀-problem in the sense that we want to
ﬁnd X ′ such that for all vectors u the inequality Eq. (1.2) holds for the uniform distribution over
X ′.)
We remark that the idea of using sum-of-squares to solve problems with quantiﬁer alternation
also plays a role in control theory (where the goal is ﬁnd a dynamical system together with an
associated Lyapunov functions, which can be viewed as sum-of-squares proof of the fact that the
dynamical system behaves nicely in an appropriate sense). However, to the best of our knowledge,
this work is the ﬁrst that uses this idea for the design of computationally-eﬃcient algorithms with
provable guarantees. We remark that in a concurrent and independent work, Hopkins and Li use
similar ideas to learn mixtures of well-separated spherical Gaussians [HL17]. In a subsequent
paper, Kothari and Steinhardt use those ideas for clustering [KS17].

1.1 Results
Without any assumptions about the underlying distribution, the best known eﬃcient algorithms
for robust mean estimation incur an estimation error that depends on the spectral norm of the
√
covariance matrix Σ of the underlying distribution and is proportional to ε (where ε > 0 is the
fraction of outliers) [DKK+ 17a, SCV17]. Concretely, given an ε-corrupted sample of suﬃciently
larger polynomial size from a distribution D, they compute an estimate µ̂ for the mean µ of D
√
such that with high probability k µ̂ − µk 6 O( ε) · kΣk 1/2 . Furthermore, this bound is optimal
for general distributions in the sense that up to constant factors no better bound is possible
information-theoretically in terms of ε and the spectral norm of Σ.
In the following theorem, we show that better bounds for the mean estimation error are possible
for large classes of distributions. Concretely, we assume (certiﬁable) bounds on higher-order
moments of the distribution (degree 4 and higher). These higher-order moment assumptions
√
allow us to improve the estimation error as a function of ε (instead of a ε bound as for the
unconditional mean estimation before we obtain an ε 1−1/k if we assume a bound on the degree-k
5

moments). Furthermore, we also obtain multiplicative approximations for the covariance matrix
(in the Löwner order sense) regardless of the spectral norm of the covariance. (Note that our notion
of certiﬁable subgaussianity does not restrict the covariance matrix in any way.)
Theorem 1.2 (Robust mean and covariance estimation under certiﬁable subgaussianity). For every
C > 0 and even k ∈ , there exists a polynomial-time algorithm that given a (corrupted) sample S ⊆ d
outputs a mean-estimate µ̂ ∈ d and a covariance-estimate Σ̂ ∈ d×d with the following guarantee: there
exists n 0 6 (C + d)O(k) such that if S is an ε-corrupted sample with size |S| > n 0 of a k-certiﬁably
C-subgaussian distribution D over d with mean µ ∈ d and covariance Σ ∈ d×d , then with high
probability
kµ − µ̂k 6 O(Ck)1/2 · ε 1−1/k · kΣk 1/2
kΣ−1/2 (µ − µ̂)k 6 O(Ck)1/2 · ε 1−1/k

(1 − δ)Σ  Σ̂  (1 + δ)Σ for δ 6 O(Ck) · ε

(1.4)
(1.5)
1−2/k

.

(1.6)

For the last two bounds, we assume in addition Ck · ε 1−2/k 6 Ω(1).4
Note that the second guarantee for the mean estimation error µ− µ̂ is stronger because kµ− µ̂k 6
kΣk 1/2 · kΣ−1/2 (µ − µ̂)k. We remark that kΣ−1/2 (µ − µ̂)k is the Mahalanobis distance between µ̂ and
D.
In general, our results provide information theoretically tight way to utilize higher moment
information and give optimal dependence of the error on the fraction of outliers in the input sample
for every (k, ℓ)-certiﬁably O(1)-subgaussian distribution. This, in particular, improves on
pthe mean
1/2
1/2
estimation algorithm of [LRV16] by improving on their error bound of O(ε )kΣk
log (d) to
O(ε 3/4 kΣk 1/2 under their assumption of bounded 4th moments (whenever “certiﬁed” by SoS).
Frobenius vs spectral norm for covariance-matrix estimation Previous work for robust covariance estimation [LRV16, DKK+ 17b] work with Frobenius norms for measuring the estimation error
Σ − Σ̂ and obtain in this way bounds that can be stronger than ours. However, it turns out that
assuming only k-certiﬁable subgaussianity makes it information-theoretically impossible to obtain
dimension-free bounds in Frobenius norm and that we have to work with spectral norms instead.
In this sense, the assumptions we make about distributions are substantially weaker compared to
previous works.
Concretely, [LRV16] show a bound of Σ − Σ̂ F 6 Õ(ε 1/2 )kΣk assuming a 4-th moment bound
and that the distribution is an aﬃne transformation of a product distribution. [DKK+ 17b] show a
bound Σ − Σ̂ F 6 O(ε)kΣk assuming the distribution is Gaussian.5
However, even for simple k-certiﬁably subgaussian distributions with parameter√C 6 O(1), the
information-theoretically optimal error in terms of Frobenius norm is Σ − Σ̂ F 6 O( d·ε 1−2/k )·kΣk.
4This notation means that we require Ck · ε1−2/k 6 c0 for some absolute constant c0 > 0 (that could in principle be
extracted from the proof).
5 In the Gaussian case, [DKK+ 17b] establish the stronger bound Σ−1/2 (Σ − Σ̂)Σ−1/2 F 6 O(ε). This norm can be
viewed as the Frobenius norm in a transformed space. Our bounds are for the spectral norm in the same transformed
space.

6

Concretely, consider the mixture of N(0, I) and N(0, ε −2/k I) with weights 1 − ε and ε, respectively.
Then, it’s easy to conﬁrm that both the mixture and the standard gaussian N(0, I) are k-certiﬁably
2-subgaussian (using Lemma 5.3 and Lemma 5.7) and at most ε far in total variation distance. Thus,
given only k-certiﬁable 2-subgaussianity, it is information theoretically impossible to decide which
of the above two distributions generated a given ε-corrupted sample.√Finally, the diﬀerence of their
covariance equals ε(ε −2/k − 1)I which has Frobenius norm Ω(ε 1−2/k ) d for any ε < 1. For this case,
our algorithms guarantee the signiﬁcantly stronger6 spectral-norm bound kΣ− Σ̂k 6 O(ε 1−2/k )· kΣk.
Multiplicative vs. additive estimation error Another beneﬁt of our covariance estimation algorithm is that provide a multiplicative approximation guarantee, that is, the quadratic form of the
estimated covariance at any vector u is within (1 ± δ) of the quadratic form of the true covariance.
This strong guarantee comes in handy, for example, in whitening or computing an isotropic transformation of the data—a widely used primitive in algorithm design. Indeed, this ability to use the
estimated covariance to whiten the data is crucial in our outlier-robust algorithm for independent
component analysis (see Section 6.1). The Frobenius norm error guarantees, in general, do not
imply good multiplicative approximations and thus cannot be used for this application.
We note that in the special case of gaussian distributions, the results of Diakonikolas et al.
[DKK+ 16] allow recovering the mean and covariance in ﬁxed polynomial time with better dependence of the error on the fraction of outliers that grows as Õ(ε)kΣk 1/2 . Our results when
applied for the special case of gaussian mean and covariance estimation will require a running
√
time of d O( log (1/ε)) to achieve a similar error guarantee. Their algorithm for gaussian covariance
estimation also provides multiplicative error guarantees.
Robust Estimation of Higher Moments Our techniques for robustly estimating mean and covariance of a distribution extend in a direct way to robustly estimating higher-order moment tensors.
In order to measure the estimation error, we use a variant of the injective tensor norm (with respect
to a transformed space where the distribution is isotropic), which generalizes the norm we use for
the estimation error of the covariance matrix. This error bound allows us to estimate for every
direction u ∈ d , the low-degree moments of the distribution in direction u with small error
compared to the second moment in direction u.
The approaches in previous works fact inherent obstacles in generalizing to the problem of
estimating the higher moments with multiplicative (i.e. in every direction u) error guarantees.
This type of error is in fact crucial in applications for learning latent variable models such as
mixtures of Gaussians and independent component analysis.
In fact, our guarantees are in some technical way stronger, which is crucial for our applications
of higher-order moment estimates. Unlike spectral norms, injective norms are NP-hard to compute
(even approximately, under standard complexity assumptions). For this reason, it is not clear how to
make use of an injective-norm guarantee when processing moment-estimates further. Fortunately,
6 Similarly to [DKK+ 17b], we also work with norms in a transformed space (in which the distribution is isotropic)
and obtain the stronger bound kΣ−1/2 (Σ − Σ̂)Σ−1/2 k 6 O(ε1−2/k ).

7

it turns out that our algorithm not only guarantees an injective-norm bound for the error but
also a good certiﬁcate for this bound, in form of a low-degree sum-of-squares proof. It turns out
that this kind of certiﬁcate is precisely what we need for our applications—in particular, recent
tensor decomposition algorithms based on sum-of-squares [MSS16a] can tolerate errors with small
injective norm if that is certiﬁed by a low-degree sum-of-squares proof.
Theorem 1.3 (Robust Higher Moment Estimation). For every C > 0 and even k ∈ , there exists
a polynomial-time algorithm that given a (corrupted) sample S ⊆ d outputs a moment-estimates M̂2 ∈
k
2
d , . . . , M̂ k ∈ d with the following guarantee: there exists n 0 6 (C + d)O(k) such that if S is an
ε-corrupted sample with size |S| > n 0 of a k-certiﬁably C-subgaussian distribution D over d with
k
2
moment-tensors M2 ∈ d , . . . , M k ∈ d such that Ck · ε 1−2/k 6 Ω(1), then with high probability for
every r 6 k/2,

∀u ∈ d . hM r − M̂ r , u ⊗r i 2 6 δ r · hM2 , u ⊗2i r

r

for δ r 6 O (Ck)r/2 · ε 1− k

(1.7)

Furthermore, there exist degree-k sum-of-squares proofs of the above polynomial inequalities in u.
Information-theoretic optimality We show in Section 7 that the error guarantees in our robust
moment-estimation algorithms are tight in their dependence on both k and ε. For example, we
show that there are√two k-certiﬁably O(1)-subgaussian distributions with statistical distance ε but
means that are Ω( kε 1−1/k ) apart. A similar statement holds for higher-order moments. The
distributions are just mixtures of two one-dimensional Gaussians.
Application: independent component analysis As an immediate application of our robust
moment estimation algorithm, we get an algorithm for Outlier Robust Independent Component
Analysis. Independent component analysis (also known as blind source separation) is a fundamental problem in signal processing, machine learning and theoretical computer science with
applications to diverse areas including neuroscience. Lathauwer et. al. [DLCC07], following
up on a long line of work gave algorithms for ICA based on 4th order tensor decomposition. A
noise-tolerant version of this algorithm was developed in [MSS16a]. There is also a line of work in
theoretical computer science on designing eﬃcient algorithms for ICA [GVX13, VX15].
In the ICA problem, we are given a non-singular7 mixing matrix A ∈ d×d with condition
number κ and a product distribution on d . The observations come from the model {Ax} that is,
the observed samples are linear transformations (using the mixing matrix) of independent draws
of the product random variable x. The goal is to recover columns of A up to small relative error
in Euclidean norm (up to signs and permutations) from samples. It turns out that information
theoretic recovery of A is possible whenever at most one source is non-gaussian. A widely used
convention in this regard is the 4th moment assumption: for each i, [x 4i ] , 3 [x 2i ]. It turns out
that we can assume [x 2i ]  1 without loss of generality so this condition reduces to asserting
[x 4i ] , 3.
7We could also consider the case that A is rectangular and its columns are linearly independent. Essentially the same
algorithm and analysis would go through in this case. We focus on the quadratic case for notational simplicity.

8

Outlier robust version of ICA was considered as an application of the outlier-robust mean and
covariance estimation problems in [LRV16]. Theypsketched an algorithm with the guarantee that
the relative error in the columns of A is at most ε log (d) poly(κ).
1
. Here,
In particular, this guarantee is meaningful only if the fraction of outliers ε ≪ √
log (d) poly(κ)

we improve upon their result by giving an outlier-robust ICA algorithm that recovers columns of
A up to an error that is independent of both dimension d and the condition number κ of the mixing
matrix A.
Our algorithm directly follows by applying 4th order tensor decomposition. However, a crucial
step in the algorithm involves “whitening” the 4th moments by using an estimate of the covariance
matrix. Here, the multiplicative guarantees obtained in estimating the covariance matrix are crucial
- estimates with respect to the Frobenius norm error do not give such whitening transformation
in general. This whitening step essentially allows us to pretend that the mixing matrix A is
well-conditioned leading to no dependence on the condition number in the error.
Theorem 1.4 (Robust independent component analysis). For every C > 1 and even k ∈ , there
exists a polynomial-time algorithm that given a (corrupted) sample S ⊆ d outputs component estimates
â1 , . . . , â d ∈ d with the following guarantees: Suppose A ∈ d×d is a non-singular matrix with condition
number κ and columns a1 , . . . , a d ∈ d . Suppose x is a centered random vector with d independent
√
coordinates such that every coordinate i ∈ [d] satisﬁes [x2i ]  1, [x4i ] − 3  γ , 0, and [xki ]1/k 6 Ck.
Then, if S is an ε-corrupted sample of size |S| > n 0 from the distribution {Ax}, where n 0 6 (C + κ + d)O(k) ,
the component estimates satisfy with high probability
max min hA−1 â i , A−1 a π(i) i 2 > 1 − δ

for δ < (1 +

π∈S d i∈[d]

1
|γ| )

· O(C 2 k 2 ) · ε 1−4/k .

(1.8)

The quantity hA−1 â i , A−1 a π(i) i is closely related to the Mahalanobis distance between â i and
a π(i) with respect to the distribution {Ax}
Application: learning mixtures of Gaussians As yet another immediate application of our
robust moment estimation algorithm, we get an outlier-robust algorithm for learning mixtures
of spherical Gaussians. Our algorithm works under the assumption that the means are linearly
independent (and that the size of the sample grows with their condition number). In return, our
algorithm does not require the means of the Gaussians to be well-separated. Our algorithm can
be viewed as an outlier-robust version of tensor-decomposition based algorithms for mixtures of
Gaussians [HK13a, BCMV14].
Theorem 1.5 (Robust estimation of mixtures of spherical Gaussians (See Theorem 6.9 for detailed
error bounds)). Let D be mixtures of N(µ i , I) for i 6 q with uniform8 mixture weights. Assume that µ i s
Í
is Ω(1).
are linearly independent and, further, assume that κ, the smallest non-zero eigenvalue of q1 i µ i µ⊤
i
Given an ε-corrupted sample of size n > n 0  Ω((d log (d))k/2 /ε 2 ), for every k > 4, there’s a
poly(n)d O(k) time algorithm that recovers µ̂1 , µ̂2 , . . . µ̂ q so that there’s a permutation π : [q] → [q]
8While our algorithm generalizes naturally to arbitrary mixture weights, we restrict to this situation for simplicity

9

satisfying
max k(
i

1Õ
−1/2
µ i µ⊤
(µ̂ i − µ π(i) )k 6 O(qk)ε 1/3−1/k .
i )
q
i

[DKK+ 16]

gave an outlier-robust algorithm that learns mixtures of q
Diakonikolas et. al.
√
gaussians with error ≈ q ε in each of the recovered means. Their algorithm is polynomial in the
dimension but has an exponential dependence on number of components q in the running time.
Under the additional assumption that the means are linearly independent, our algorithm (say for
k  8) recovers similar error guarantees as theirs but runs in time polynomial in both q and d. The
key diﬀerence is the power of our algorithm to recover a multiplicative approximation to the 4th
moment tensor which allows us to apply blackbox tensor decomposition based methods and run
in ﬁxed polynomial time [HK13b].

2 Robust identifiability of low-degree moments
In this section, we show that low-degree moments of certiﬁably subgaussian distributions are
identiﬁable in the presence of outliers. As we will formalize in Section 4, these proofs of identiﬁability are captured by the sum-of-squares proof system at low-degree. This fact is the basis of our
polynomial-time algorithms for robust moment estimation (also Section 4).
We recall the basic strategy for the identiﬁability proofs as described in the introduction: Let
D0 be a certiﬁably subgaussian distribution whose moments we aim to estimate. If we take a large
enough sample X ⊆ d from D0 , then with high probability, the empirical distribution D (uniform
over the sample) has approximately the same low-degree moments as the (population) D0 . In
this case, the low-degree moments of D are also certiﬁably subgaussian9 (for essentially the same
parameters C as D0 ). Our input consists of an ε-corruption Y ⊆ d of the sample X. In order to
show that the low-degree moments of D0 (or equivalently D) are identiﬁable from Y, we analyze
the following (a-priori ineﬃcient) estimation procedure:
1. ﬁnd a certiﬁably subgaussian distribution D ′ that has statistical distance at most ε from the
uniform distribution over Y,
2. output the low-degree moments of D ′.
We can typically ﬁnd a distribution D ′ as above, because with high-probability D satisﬁes
the conditions. What remains to prove is that no matter what distribution D ′ satisfying those
conditions we choose, our moment estimates have low error. To this end, we show several concrete
and quantiative instances of the following general mathematical statement:
If two distributions D and D ′, whose low-degree moments are (certiﬁably) subgaussian, have
small statistical distance, then their low-degree moments are close.
Our ﬁrst bound of this kind, controls the distance between the means of D and D ′ in terms of
their covariances Σ and Σ′. (Later bounds will also relate Σ and Σ′ .)
9 Roughly speaking, certiﬁable subgaussianity is inherited by the empirical distribution from the population distribution because it is a property of the low-degree moments of distributions. We prove this claim in Section 5.

10

Lemma 2.1 (Robust identiﬁability of mean). Let k ∈  be even. Let D, D ′ be two distributions over d
with means µ, µ′ ∈ d and covariances Σ, Σ′ ∈ d×d . Suppose kD − D ′ k TV 6 ε < 0.9 and that D and D ′
are (k, ℓ)-certiﬁably subgaussian with parameter C > 0. Then,

∀u ∈ d . |hu, µ − µ′i| 6 δ · hu, (Σ + Σ′)ui 1/2

where δ  O

√



Ck · ε 1−1/k .

(2.1)

In particular, kµ − µ′ k 6 δkΣ + Σ′ k 1/2 and k(Σ + Σ′ )−1/2 (µ − µ′)k 6 δ. (Here, d is endowed with the
standard Euclidean norm and d×d with the corresponding operator norm.)
Proof. Consider a coupling between the distributions D(x) and D ′(x ′ ) such that {x  x ′ } > 1 − ε.
By Hölder’s inequality, we have for every vector u ∈ d ,



|hu, µ − µ′i|  | hu, x − x ′i 1x ,x′ | 6

At the same time,



′ k

hu, x − x i

 1/k



6





′

′

hu, x − x ′i k

′

hu, x − µ − (x − µ ) + (µ − µ )i
hu, x − µi k

 1/k

+



k

 1/k

hu, x ′ − µ′i k

 1/k 

 1/k

k/(k−1)

 1−1/k

1x ,x′
|¨¨¨¨¨¨¨¨¨¨¨¨¨¨{z¨¨¨¨¨¨¨¨¨¨¨¨¨¨}

(2.2)

6 ε 1−1/k

(2.3)
+ |hu, µ − µ′i|

(2.4)

(using triangle inequality for L k -norm)

√
 1/2
 1/2 
6 Ck ·
+ |hu, µ − µ′i|
hu, x − µi 2
+
hu, x ′ − µ′i 2

(using certiﬁable subgaussianity)
√
6 O( Ck ) · hu, (Σ + Σ′)ui 1/2 + |hu, µ − µ′i| .

(2.5)

(2.6)

Combining the two bounds yields


 √
|hu, µ − µ′i| 6 ε 1−1/k · O( Ck ) · hu, (Σ + Σ′)ui 1/2 + |hu, µ − µ′i| .

Since we
that ε < 0.9, we can conclude that |hu, µ − µ′i| 6 δ · hu, (Σ + Σ′ )ui 1/2 for
 assume
√ 

δ 6 O ε 1−1/k · Ck as desired.

Next we show that if two certiﬁably subgaussian distributions have small statistical distance,
then their (raw) second moments are close in the Löwner order sense.
Lemma 2.2 (Robust identiﬁability of second moment). Let D, D ′ be two distributions over d with
second-moment matrices M, M ′ ∈ d×d . Suppose kD − D ′ k TV 6 ε < 0.9 and that D and D ′ are
(k, ℓ)-certiﬁably subgaussian with parameters C > 0. Then,



− δ · (M + M ′)  M − M ′  δ · (M + M ′) for δ 6 O Ck · ε 1−2/k





In particular, 1 − O(δ) M ′  M  1 + O(δ) M ′ if δ < 0.9.
11



(2.7)

Proof. Consider a coupling between the distributions D(x) and D ′(x ′ ) such that {x  x ′ } > 1 − ε.
Then, for every vector u ∈ d ,
|hu, (M − M ′)ui|  | hu, xi 2 − hu, x ′i 2 |

6 ε 1−2/k ·
6ε

1−2/k

·





hu, xi 2 − hu, x ′i 2
hu, xi

k

6 ε 1−2/k · O(Ck) ·

 2/k

+



 k/2  2/k
′ k

hu, x i

hu, xi 2 +

(using Hölder inequality)

 2/k 

hu, x ′i 2

 O(Ck · ε 1−2/k ) · hu, (M + M ′)ui ,



(by L k/2 triangle inequ.)
(by subgaussianity)



which implies the desired bound −δ · (M + M ′)  M − M ′  δ · (M + M ′) for δ 6 O Ck · ε 1−1/k . For
the third inequality, we use Lemma 5.2 (moment bounds for shifts of subgaussian distributions).
For δ 6 1, we can now rearrange this bound to get the following relationship between M and
′
M,
′
′
1+δ
1−δ
(2.8)
1+δ · M  M  1−δ · M ,





which implies the desired bound 1 − O(δ) M ′  M  1 + O(δ) M ′ for δ < 0.9.



Next, we combine the previous bounds for the ﬁrst and second moments in order to establish
the identiﬁability of the covariance matrix.
Corollary 2.3 (Covariance identiﬁability). Under the same conditions as the previous Lemma 2.2, the
covariance matrices Σ and Σ′ of the distributions D and D ′ satisfy,
− δ · (Σ + Σ′)  Σ − Σ′  δ · (Σ + Σ′)

for δ 6 O(Ck · ε 1−2/k ) .

(2.9)

As before, in the case that δ 6 0.9, the above bounds imply a strong multiplicative approximation of the covariance so that (1 − δ′) · Σ′  Σ  (1 + δ′)Σ′ for some δ′ 6 O(δ).
Proof of Corollary 2.3. Let Dµ , Dµ′ be the distributions D, D ′ shifted by µ so that Dµ has expectation
0 and Dµ′ has expectation µ′ − µ. Note that Dµ − Dµ′

TV

 kD − D ′ k TV 6 ε. By Lemma 2.2

(second-moment identiﬁability applied to Dµ , Dµ′ ), every vector u ∈ d satisﬁes,

|hu, (Σ − Σ′ − (µ − µ′)(µ − µ′)T )ui| 6 δ 1 · hu, (Σ + Σ′ + (µ − µ′)(µ − µ′)T )ui .

(2.10)

Here, δ1 6 O(Ck · ε 1−1/k ). By Lemma 2.1 (mean identiﬁability), hu, (µ − µ′)(µ − µ′)T ui 6 δ 2 · hu, (Σ +
Σ′)ui, where δ 2 6 O(Ck · ε 1−1/k ). Combining the inequalities gives the desired bound.

Taking together Corollary 2.3 and Lemma 2.1, we get a mean estimation error that depends
only on the covariance of D (as opposed to the covariances of both D and D ′).
Corollary 2.4 (Stronger mean identiﬁability). Under the same conditions as the previous Lemma 2.2
and the additional condition Ck · ε 1−2k 6 Ω(1), the means µ and µ′ of the distributions D and D ′ satisfy,

√
(2.11)
∀u ∈ d . |hu, µ − µ′i| 6 δ · hu, Σui 1/2 where δ  O Ck · ε 1−1/k ,
12

where Σ ∈ d×d is the covariance of D. In particular, kµ − µ′ k 6 δkΣk 1/2 and kΣ−1/2 (µ − µ′)k 6 δ. (Here,
d is endowed with the standard Euclidean norm and d×d with the corresponding operator norm.)
Next we bound the distances of higher order moments for certiﬁably subgaussian distributions
that are close in statistical distance. Compared to ﬁrst and second order moments, distances for
higher order moments are less standard. Our choice for the distance is informed by the applications
of independent component analysis and learning mixtures of Gaussians. We view the order-2r
moments of D and D ′ as two polynomials p and p ′ of degree 2r and we bound their diﬀerence
p(u) − p ′(u) relative to the variance of the distribution in direction u.
Lemma 2.5 (Robust identiﬁability of higher moments). Let r ∈  and let D, D ′ be two distributions
2
r
over d with second moments M2 , M2′ ∈ d and r-th moments M r , M r′ ∈ d . Suppose kD − D ′ k TV 6
ε < 0.9 and that D and D ′ are (k, ℓ)-certiﬁably subgaussian with parameters C > 0 and 2r 6 k. Then, for
every u ∈ d ,
r

|hM r − M r′ , u ⊗r i| 6 δ · hM2 + M2′ , u ⊗2 i r/2

for δ 6 O(Ck)r/2 · ε 1− k

(2.12)

Proof. Consider a coupling between the distributions D(x) and D ′(x ′ ) such that {x  x ′ } > 1 − ε.
Then, for every vector u ∈ d ,
|hu, (M r − M r′ )ui|  | hu, xi r − hu, x ′i r |

6ε

1−r/k

·

6 ε 1−r/k ·





r

′ r k/r

(hu, xi − hu, x i )
hu, xi k

6 ε 1−r/k · O(C r k r ) ·

 r/k



+



 r/k

hu, x ′i k

hu, xi 2

 r/2

+

(using Hölder’s inequality)

 r/k 

(by L k/r triangle inequ.)

hu, x ′i 2

r

 r/2 

(by subgaussianity)

which implies the desired bound for δ 6 δ 6 O(Ck)r/2 · ε 1− k . For the third inequality, we rely on
Lemma 5.2 (moment bounds for shifts of subgaussian distributions).


3 Preliminaries
In this section, we deﬁne pseudo-distributions and sum-of-squares proofs. See the lecture notes
[BS16] for more details and the appendix in [MSS16b] for proofs of the propositions appearing
here.
Let x  (x1 , x2 , . . . , x n ) be a tuple of n indeterminates and let [x] be the set of polynomials
with real coeﬃcients and indeterminates x1 , . . . , x n . We say that a polynomial p ∈ [x] is a
sum-of-squares (sos) if there are polynomials q 1 , . . . , q r such that p  q 12 + · · · + q 2r .

3.1 Pseudo-distributions
Pseudo-distributions are generalizations of probability distributions. We can represent a discrete
(i.e., ﬁnitely supported) probability distribution over n by its probability mass function D : n →
13

 such that D > 0 and

Í

D(x)  1. Similarly, we can describe a pseudo-distribution by
its mass function. Here, we relax the constraint D > 0 and only require that D passes certain
low-degree non-negativity tests.
Concretely, a level-ℓ pseudo-distribution is a ﬁnitely-supported function D : n →  such that
Í
Í
2
x D(x)  1 and
x D(x) f (x) > 0 for every polynomial f of degree at most ℓ/2. (Here, the
summations are over the support of D.) A straightforward polynomial-interpolation argument
shows that every level-∞-pseudo distribution satisﬁes D > 0 and is thus an actual probability
distribution. We deﬁne the pseudo-expectation of a function f on d with respect to a pseudodistribution D, denoted ˜ D(x) f (x), as
x∈supp(D)

˜ D(x) f (x) 

Õ
x

D(x) f (x) .

(3.1)

The degree-ℓ moment tensor of a pseudo-distribution D is the tensor D(x) (1, x1 , x2 , . . . , x n )⊗ℓ .
In particular, the moment tensor has an entry corresponding to the pseudo-expectation of all
monomials of degree at most ℓ in x. The set of all degree-ℓ moment tensors of probability
distribution is a convex set. Similarly, the set of all degree-ℓ moment tensors of degree d pseudodistributions is also convex. Key to the algorithmic utility of pseudo-distributions is the fact that
while there can be no eﬃcient separation oracle for the convex set of all degree-ℓ moment tensors
of an actual probability distribution, there’s a separation oracle running in time n O(ℓ) for the convex
set of the degree-ℓ moment tensors of all level-ℓ pseudodistributions.
Fact 3.1 ([Sho87, Par00, Nes00, Las01]). For any n, ℓ ∈ , the following set has a n O(ℓ) -time weak
separation oracle (in the sense of [GLS81]):



˜ D(x) (1, x1 , x2 , . . . , x n )⊗d | degree-d pseudo-distribution D over n

(3.2)

.

This fact, together with the equivalence of weak separation and optimization [GLS81] allows
us to eﬃciently optimize over pseudo-distributions (approximately)—this algorithm is referred to
as the sum-of-squares algorithm.
The level-ℓ sum-of-squares algorithm optimizes over the space of all level-ℓ pseudo-distributions
that satisfy a given set of polynomial constraints—we formally deﬁne this next.
Definition 3.2 (Constrained pseudo-distributions). Let D be a level-ℓ pseudo-distribution over n .
Let A  { f1 > 0, f2 > 0, . . . , f m > 0} be a system of m polynomial inequality constraints. We say
that D satisfies the system of constraints A at degree r, denoted D r A, if for every S ⊆ [m] and every
Í
sum-of-squares polynomial h with deg h + i∈S max{deg f i , r},
˜Dh·

We write D

Ö

fi > 0 .

i∈S

A (without specifying the degree) if D

0

A holds. Furthermore, we say that D
−n ℓ

Î

r

A

· kh k · i∈S k f i k,
holds approximately if the above inequalities are satisﬁed up to an error of 2
where k·k denotes the Euclidean norm10 of the coﬃcients of a polynomial in the monomial basis.
ℓ

10The choice of norm is not important here because the factor 2−n swamps the eﬀects of choosing another norm.

14

We remark that if D is an actual (discrete) probability distribution, then we have D A if and
only if D is supported on solutions to the constraints A.
We say that a system A of polynomial constraints is explicitly bounded if it contains a constraint
of the form {kx k 2 6 M}. The following fact is a consequence of Fact 3.1 and [GLS81],
Fact 3.3 (Eﬃcient Optimization over Pseudo-distributions). There exists an (n + m)O(ℓ) -time algorithm
that, given any explicitly bounded and satisfiable system11 A of m polynomial constraints in n variables,
outputs a level-ℓ pseudo-distribution that satisfies A approximately.

3.2 Sum-of-squares proofs
Let f1 , f2 , . . . , f r and 1 be multivariate polynomials in x. A sum-of-squares proof that the constraints
{ f1 > 0, . . . , f m > 0} imply the constraint {1 > 0} consists of polynomials (p S )S⊆[m] such that

Õ

1

S⊆[m]

p S · Πi∈S f i .

(3.3)

We say that this proof has degree ℓ if for every set S ⊆ [m], the polynomial p S Πi∈S f i has degree at
most ℓ. If there is a degree ℓ SoS proof that { f i > 0 | i 6 r} implies {1 > 0}, we write:
{ f i > 0 | i 6 r}

ℓ

{1 > 0} .

(3.4)

Sum-of-squares proofs satisfy the following inference rules. For all polynomials f , 1 : n → 
and for all functions F : n → m , G : n → k , H : p → n such that each of the coordinates of
the outputs are polynomials of the inputs, we have:
{ f > 0, 1 > 0} A
,
A ℓ { f + 1 > 0}

A
A

ℓ

ℓ

B, B

A

ℓ·ℓ′

ℓ′

ℓ

{ f > 0}, A
A

ℓ+ℓ′

ℓ′

{1 > 0}

{ f · 1 > 0}

C

(transitivity)

C

{F > 0}
{F(H) > 0}

ℓ

(addition and multiplication)

{G > 0}

ℓ·deg(H)

{G(H) > 0}

(substitution)

.

Low-degree sum-of-squares proofs are sound and complete if we take low-level pseudodistributions as models.
Concretely, sum-of-squares proofs allow us to deduce properties of pseudo-distributions that
satisfy some constraints.
Fact 3.4 (Soundness). If D
proof A

r′

B, then D

r·r ′+r ′

r

B.

A for a level-ℓ pseudo-distribution D and there exists a sum-of-squares

11Here, we assume that the bitcomplexity of the constraints in A is (n + m)O(1) .

15

If the pseudo-distribution D satisﬁes A only approximately, soundness continues to hold if
we require an upper bound on the bit-complexity of the sum-of-squares A r ′ B (number of bits
required to write down the proof).
In our applications, the bit complexity of all sum of squares proofs will be n O(ℓ) (assuming that
all numbers in the input have bit complexity n O(1) ). This bound suﬃces in order to argue about
pseudo-distributions that satisfy polynomial constraints approximately.
The following fact shows that every property of low-level pseudo-distributions can be derived
by low-degree sum-of-squares proofs.
Fact 3.5 (Completeness). Suppose d > r ′ > r and A is a collection of polynomial constraints with degree
Í
at most r, and A ⊢ { ni1 x 2i 6 B} for some finite B.

Let {1 > 0} be a polynomial constraint. If every degree-d pseudo-distribution that satisfies D
satisfies D r ′ {1 > 0}, then for every ε > 0, there is a sum-of-squares proof A d {1 > −ε}.

r

A also

4 Moment estimation algorithm
Our robust moment estimation algorithm is a direct translation of the robust identiﬁability results
from Section Section 2.
Let Y  { y1 , . . . , yn } ⊆ d be an ε-corrupted sample of a distribution D over d . We consider
the following system A : AY,ε of quadratic equations in scalar-valued variables w 1 , . . . , w n and
vector-valued variables x1′ , . . . , x ′n ,

AY,ε :







Ín

i1

w i  (1 − ε) · n 



w 2i  w i



 ∀i ∈ [n]. w i · (yi − x ′)  0

i
∀i ∈ [n].








(4.1)

The following fact shows that the solutions to AY,ε correspond to all ε-corruptions of the
(already corrupted) sample Y. In particular, if Y  { y1 , . . . , yn } is an ε-corruption of a sample
X  {x1 , . . . , x n } of D, then one solution to A consists of the uncorrupted sample X so that
x1′  x1 , . . . , x ′n  x n .
Fact 4.1. An assignment to the variables x1′ , . . . , x ′n can be extended to a solution for AY,ε if and only if
x ′i  yi holds for all but at most an ε-fraction of the indices i ∈ [n].
Note that in order to extend such an assignment, we can choose w i  1 for (1 − ε)n indices
i ∈ [n] such that x i  yi and w i  0 for the remaining indices.
Furthermore, we consider the following system of polynomial equations B : BC,k in vectorvalued variables x1′ , . . . , x ′n , vector-valued variables p1 , . . . , p k/2 , and matrix-valued variables
Q 1 , . . . , Q k/2 . The idea is that p1 , . . . , p k/2 and Q 1 , . . . , Q k/2 are degree-ℓ sum-of-squares proofs of
the inequalities Eq. (1.3) (C-subgaussianity of moments up to degree k) for the uniform distribution

16

over x1′ , . . . , x ′n .

BC,k,ℓ :




 ∀k ′ ∈ 6 k/2 .






1
n

n
Õ
i1

hx ′i ,

ui

2k ′



′

 Ck ·

1
n

n
Õ
i1

hx ′i ,

ui

2

 k′




+ hp k′ , (1, u)⊗ℓ i · (1 − ku k 2) 


− kQ k′ · (1, u)⊗ℓ/2 k 2






(4.2)

We emphasize that in the above system of equations the vector-valued variable u  (u 1 , . . . , u d ) is
not part of the solution (i.e., we are not seeking an assignment for it). Instead, each equation in BC,k,ℓ
is about two polynomials with indeterminates u 1 , . . . , u d , whose coeﬃcients are polynomials in the
variables x1′ , . . . , x ′n , p1 , . . . , p k/2 , and Q 1 , . . . , Q k/2 . In particular, we could eliminate the variable
u from the equations in BC,k,ℓ by explicitly writing down the expressions for the coeﬃcients of the
polynomials in u. We refrain from doing so because these expressions would be rather complicated
and would obfuscate the meaning of the equations.
The following fact shows that the solutions to the system BC,k correspond to certiﬁably subgaussian distributions.
Fact 4.2. An assignment to the vector-valued variables x1′ , . . . , x ′n can be extended to a solution for BC,k,ℓ
if and only if the uniform distribution over x1′ , . . . , x ′n is (k, ℓ)-certifiably subgaussian with parameter C.
The following algorithm is the key ingredient of our eﬃcient outlier-robust estimators.
Algorithm 4.3 (Algorithm for moment estimation via sum-of-squares).
Given: ε-corrupted sample Y  { y1 , . . . , yn } ⊆ d of a (k, ℓ)-certiﬁably 0.9C-subgaussian
distribution D0 over d
Estimate: moments M1 , . . . , M k with M r 

D0 (x)

x ⊗r

Operation:
1. ﬁnd a level-ℓ pseudo-distribution D ′ that satisﬁes AY,ε and BC,k,ℓ
 Í

2. output moment-estimates M̂1 , . . . , M̂ k with M̂ r  ˜ D′(x′ ) 1 n (x ′ )⊗r
n

i1

i

The following lemma is the sum-of-squares analog of Lemma 2.1.

Lemma 4.4. Let D be the uniform distribution over X ⊆ d with mean µ ∈ d and covariance Σ ∈ d×d .
Suppose D is (k, ℓ)-certifiably subgaussian with parameter C > 0 and that Ck · ε 1−2/k 6 Ω(1).12 Let Y
Í
be an ε-corruption of X with ε < 0.01. Let µ′ , Σ′ be the polynomials µ′(X ′)  n1 ni1 x ′i and Σ′(X ′) 
Í
n
1
′
′ T
′
d
′
′ T
′
i1 (x i )(x i ) − (µ )(µ ) in the vector-valued variables x 1 , . . . , x n . Then, for every vector u ∈  ,
n
where δ  O

√

AY,ε ∪ BC,k,ℓ



X′
ℓ



hu, µ − µ′(X ′)i k 6 δ k hu, (Σ + Σ′(X ′))ui k/2 ,

(4.3)

Ck · ε 1−1/k . Furthermore, the bit complexity of the sum-of-squares proof is polynomial.

12 This notation means that we require Ckε1−2/k 6 γ for some absolute constant γ > 0, e.g., γ  10−10 .

17

Í

Proof. Choose r1 , . . . , rn ∈ {0, 1} such that i ri  (1 − ε)n and ri (x i − yi )  0 for all i ∈ [n]. (Such
a choice exists because Y is an ε-corruption of X.)
Letting z i  ri · w i , we claim that
w
AY,ε 2

(

n
Õ
i1

z i > (1 − 2ε)n

)

.

Í

wi

Indeed, since {w 2i  w i } 2 {(1 − z i ) 6 (1 − w i ) + (1 − ri )}, we have A Y,ε 2 { ni1 (1 − z i ) 6 2εn}
By the sum-of-squares version of Hölder’s inequality, we have for every vector u ∈ d , using
that A Y,ε k {(1 − z i )k−1  (1 − z i )}.

n

w,x′
ℓ

AY,ε ∪ BC,k,ℓ

hu, µ − µ′i k 

6

 Õ
1
n

 Õ
i

1
n

 Õ
i

6

1
n

i



For brevity, we will use
the fact that

a,b,c
k

AY,ε ∪ BC,k,ℓ

i∈[n]

c)k

(a + b +

w,x′
ℓ

n



i∈[n]



hu, x i −

x ′i i k

k

 Õ


1
n

i

(2ε)

(4.4)

(1 − z i )k

k−1

o

 k−1

.

to denote the average over indices i ∈ [n]. At the same time, using

6 3k (a k + b k + c k ) (similar to Lemma A.2),


hu, x i − µ − (x ′i − µ′) + (µ − µ′)i k

 Õ
1
n

i

k

hu, x i − x ′i i k

hu, x i − x ′i i k

i∈[n]

6 3k

hu, x i − x ′i i(1 − z i )

6 3 Ck



hu, x i − µi k + 3k

i∈[n]

hu, x i − µi


2 k/2

 Õ
1
n

i





hu, x ′i − µ′i k + 3k hu, µ − µ′i k

k

+ 3 Ck

i∈[n]

hu, x ′i − µ′i 2

(using certiﬁable subgaussianity)

6 O(Ck)

k/2

′

· hu, (Σ + Σ )ui

Combining the two bounds yields

k/2

k

′ k

+ 3 hu, µ − µ i .

o

 k/2



+ 3k hu, µ − µ′i k

(4.5)



hu, µ − µ′i k 6 (2ε)k−1 · O(Ck)k/2 · hu, (Σ + Σ′ )ui k/2 + 3k hu, µ − µ′i k .
Since we assume that ε < 0.01, we can conclude that
AY,ε ∪ BC,k,ℓ

w,x′
ℓ


√ 
for δ 1 6 O ε 1−1/k · Ck as desired.

hu, µ − µ′i k 6 δ 1k · hu, (Σ + Σ′)ui k/2 ,

(4.6)


18

The following lemma is the sum-of-squares analog of Lemma 2.2.
Lemma 4.5 (Robust identiﬁability of second moment). Let D be the uniform distribution over X ⊆ d
with second moment M ∈ d×d . Suppose D is (k, ℓ)-certifiably subgaussian with parameter C > 0 with
k divisible by 4 and that Ck · ε 1−2/k 6 Ω(1). Let Y be an ε-corruption of X. Let M ′ be the quadratic
Í
polynomial n1 i x ′i x ′iT in vector-valued variables x1′ , x2′ , . . . , x ′n . Then, for every u ∈ d ,
A Y,ε ∪ BC,k,ℓ

X′
ℓ





hu, (M − M ′)ui k/2 6 δ k/2 hu, M + M ′)ui k/2

Furthermore, if δ 6 0.0001,
AY,ε ∪ BC,k,ℓ

X′
ℓ

n

1 − δ′ 6

1
hu,Mui hu,

M ′ui 6 1 + δ′

for δ 6 O Ck · ε 1−2/k

o



(4.7)

for δ′  100δ .

Furthermore, the bit complexity of the sum-of-squares proof is polynomial.
Proof. By an argument analogous to the one in the proof of Lemma 4.4, we can obtain:
AY,ε ∪ BC,k,ℓ



w,x′
ℓ



hu, (M − M ′)ui k/2 6 δ k/2 · hu, (M + M ′))ui k/2 ,

(4.8)

for δ 6 O ε 1−2/k · Ck . Lemma A.4 together with (4.8) implies for some δ′ 6 100δ,
A Y,ε ∪ BC,k,ℓ

X′
ℓ

n

1 − δ′ 6

1
hu,Mui hu,

o

M ′(X ′)ui 6 1 + δ′ .

The following lemma is the sum-of-squares analog of Corollary 2.3.



Corollary 4.6 (Covariance estimation). Let D be the uniform distribution over X ⊆ d with covariance
Σ ∈ d×d . Suppose D is (k, ℓ)-certifiably subgaussian with parameter C > 0 with k divisible by 4 and
Í
that Ck · ε 1−2/k 6 Ω(1). Let Y be an ε-corruption of X. Let Σ′  n1 ni1 (x ′i − µ′)(x ′i − µ′)T be the formal
Í
covariance of the uniform distribution over x1′ , . . . , x ′n , where µ′  n1 i x ′i .
Then, for every u ∈ d ,
AY,ε ∪ BC,k,ℓ

X′
ℓ

1−δ 6

′
′
1
hu,Σui hu, Σ (X )ui

6 1 + δ for some δ 6 O(Ck) · ε 1−2/k .

(4.9)

Furthermore, the bit complexity of the sum-of-squares proof is polynomial.
Proof. Let Dµ be the distributions D shifted by µ so that Dµ has expectation 0. Then, Dµ is (k, ℓ)w,X ′ 
w i (z i − (x i − µ))
2
′ ,u
Í
Í
X
′
′
using Lemma 5.2, for every k ′ 6 k/2, BC,k,ℓ ℓ n1 i hx ′i , ui 2k 6 (2Ck ′)k /2 ( n1 i hx ′i ,
Thus, applying Lemma 4.5, we obtain for every vector u ∈ d ,

certiﬁably 2C-subgaussian. Further, for z i  x ′i − µ i AY,ε

AY,ε ∪ BC,k,ℓ

ui 2)k/2 .

w,X ′
ℓ

hu, (Σ − Σ′ − (µ − µ′)(µ − µ′)T )ui k/2 6



 0 . Further,

(4.10)
k/2
δ2

· hu, (Σ + Σ′ + (µ − µ′)(µ − µ′)T )ui k/2



6 (2δ2 )k/2 hu, (Σ + Σ′)ui k/2 + hu, (µ − µ′)(µ − µ′)T ui k/2 .
19

(4.11)
(4.12)

Here, δ 2 6 O(Ck · ε 1−2/k ). Further,
X ′ ,w
k



hu, (Σ − Σ′)ui k/2 6 2k/2 hu, (Σ − Σ′ − (µ − µ′)(µ − µ′)T )ui k/2 + 2k/2 hu, (µ − µ′)(µ − µ′)T ui k/2 .

By Lemma 2.1 (mean estimation), hu, (µ − µ′)(µ − µ′)T ui k/2 6 δ 21 · hu, (Σ + Σ′ )ui k/2 , where
δ1 6 O(Ck · ε 1−1/k ). Combining the above inequalities, we obtain:
w,X ′
ℓ

AY,ε ∪ BC,k,ℓ

hu, (Σ − Σ′)ui k/2 6 O(δ 2 )k/2 hu, (Σ + Σ′)ui k/2 .

(4.13)

Further,
w,X ′
k

hu, (Σ + Σ′ )ui k/2  (hu, 2Σui − hu, (Σ − Σ′ )ui)k/2

6 2k/2 hu, 2Σui + 2k/2 hu, (Σ − Σ′)ui k/2 .

Using this in conjunction with (4.13), we obtain:
AY,ε ∪ BC,k,ℓ

w,X ′
ℓ

Now, using Lemma A.3 with f 

hu, (Σ − Σ′)ui k/2 6 O(δ 2 )k/2 hu, Σui k/2.

hu,(Σ−Σ′ )ui
O(δ2 ) k/2 hu,Σui

completes the proof.



The following is the analog of the Corollary 2.4.
Corollary 4.7 (Stronger Mean Estimation). Let D be the uniform distribution over X ⊆ d with mean,
covariance µ ∈ d and Σ ∈ d×d . Suppose D is (k, ℓ)-certifiably subgaussian with parameter C > 0 with
k divisible by 4 and that Ck · ε 1−2/k 6 Ω(1). Let Y be an ε-corruption of X. Let µ′(X ′) be the linear
Í
polynomial n1 i x ′i in the vector-valued variables x1′ , x2′ , . . . , x ′n .

∀u ∈ d .AY,ε ∪ BC,k,ℓ

w,X ′
ℓ

n

± hu, µ − µ′i 6 δ · hu, Σui 1/2

w,X ′

where δ  O

√

Ck · ε 1−1/k

o

. (4.14)

In particular, AY,ε ∪ BC,k,ℓ ℓ kµ − µ′ k 6 δkΣk 1/2 and kΣ−1/2 (µ − µ′)k 6 δ. Furthermore, the bit
complexity of the sum-of-squares proof is polynomial.
√
Proof. Combining Lemma 4.4 and Corollary 4.6 we have for δ 1  O( Ckε 1−1/k ) and δ 2 
O(Ck)ε 1−2/k ,
w,X ′
ℓ
′ k

AY,ε ∪ BC,k,ℓ

hu, µ − µ′(X )i 6 δ 1k hu, (Σ + Σ′)ui k/2
2

6 (2δ1 )



hu, Σui

k/2

6 (6δ1 )2 hu, Σui k/2.

′

+ hu, Σ ui

k/2



Observe that hu, Σui k/2 is a constant, i.e., does not depend on any variable in the SoS proof.
hu,µ−µ′ (X ′ )i

Using Lemma A.3 with f  (6δ )2 hu,Σui1/2 now completes the proof.
1

Finally, the following lemma is the sum-of-squares analog of Lemma 2.5. The proof is similar
to the ones presented before, so we omit it here.
20

Lemma 4.8. Let D be the uniform distribution over X ⊆ d with second moments M2 ∈ d×d . Suppose D
is (k, ℓ)-certifiably subgaussian with parameter C > 0 for even k divisible by r. Let Y be ε-corruption sample
Í
be a r-order tensor with entries that are r-degree polynomial in vector-valued
of X. Let M r′  n1 ni1 x ′⊗r
i
′
′
variables x1 , . . . , x n . Then, the following inequality has a sum of squares proof of degree ℓ in variables of
AY,ε ∪ BC,k,ℓ and u ∈ d ,
u,w,X ′
ℓ







−δhu, M2 + M2′ ui k/2 6 hM r − M r′ , u ⊗r i k/r 6 δhu, M2 + M2′ ui k/2 ,
(4.15)
where δ  O (C r k r ) · ε 1−r/k . Consequently,
AY,ε ∪ BC,k,ℓ



AY,ε ∪ BC,k,ℓ ∪ ku k22 6 1

u,w,X ′
ℓ



−δhu, M2 ui r/2 6 hM r − M r′ , u ⊗r i 6 δhu, M2 ui r/2 .

(4.16)

Furthermore, the bit complexity of the sum-of-squares proof is polynomial.
Remark 4.9. We point out a subtle technical diﬀerence with respect to the two results above. In
Lemma 4.4 and Corollary 4.6, we do not force u to be a variable in the sum of squares proof system.
On the other hand, in Lemma 4.8, the resulting polynomial inequality has a sum of squares proof
in u (in addition to the other variables). This is important because, the inequalities are degree 2 (in
u) for Corollary 4.7 and Corollary 4.6 while higher degree in general. While unconstrained (in u)
inequalities of degree 6 2 always have a SoS proof, this is not true for higher degrees. The lemma
above relies on Proposition A.5 to obtain this stronger conclusion.

4.1 Putting things together
We can now complete the proof of Theorem 1.2.
Proof of Theorem 1.2. We just put together the implications of the above results to analyze Algorithm
4.3 here. Let S′ be the true sample used to produce the ε-corrupted sample S. Then, US′ - the
uniform distribution on S′ is (k, ℓ)-certiﬁably 2C-subgaussian by an application of Lemma 5.5.
Let ζ be the pseudo-distribution of degree at least ℓ over variables w, y satisfying AY,ε ∪ BC,k,ℓ
Í
Í
output by Algorithm 4.3. Deﬁne µ′  n1 ni1 yi and Σ′  n1 ni1 (yi − µ̂)(yi − µ̂)⊤ be the polynomials
deﬁned in Corollary 4.7 and Corollary 4.6 respectively.
We round the pseudo-distribution simply by setting µ̂  ˜ ζ µ′.
From Corollary 4.7, we have that:



AY,ε ∪ BC,k,ℓ ⊢ℓ hu, µ − µ′i 2 6 δ 2 hu, Σui ,
Thus, for every u ∈ d , ˜ ζ hu, µ − µ′i 2 6 δ 2 hu, Σui
By Cauchy-Schwarz inequality for pseudo-distributions, we have:

2

6 hu, µ − ˜ ζ µ′i 2 6 δ2 hu, Σui.
2
Using that ˜ µ′  µ̂ thus yields hu, µ − µ̂i 6 δ 2 hu, Σui.
˜ ζ hu, µ − µ′i

21

(4.17)

Using u → Σ−1/2 u thus immediately yields:



−1/2

hu, Σ

(µ − µ̂)i

2

6 δ2 ku k 2.

Taking square roots thus implies hu, Σ−1/2 (µ − µ̂)i 6 δku k or equivalently, kΣ−1/2 (µ − µ̂)k 6 δ.
On the other hand, taking square roots in (4.17) gives: |hu, µ − µ̂i|  kµ − µ̂k 6 δkΣk 1/2 as
required.
Next, we analyze the covariance estimation procedure from Algorithm 4.3. As before, our
rounding is simple: Σ̂  ˜ ζ Σ′ .
Applying Corollary 4.6 and using that degree ℓ pseudo-expectations must respect degree ℓ
SoS proofs (Fact 3.4), we get that for every u ∈ d , whenever δ < 1, (1 − δ)Σ  Σ̂  (1 + δ)Σ as
required.

The proof of Theorem 1.3 is entirely analogous. That the inequality in the conclusion has a
sum-of-squares proof in u follows from Proposition A.5.

5 Certifiably subgaussian distributions
In this section, we give suﬃcient conditions for distributions to be certifably subgaussian. As
concrete consequences, we derive that (aﬃne transformations of) products of scalar-valued subgaussian distributions are certiﬁably subgaussian and that certain mixture models are certiﬁably
subgaussian.
For the beneﬁt of the reader, we restate here the deﬁnition of certiﬁable subgaussianity.
Definition (Certiﬁable subgaussianity). A distribution D over d with mean µ is (k, ℓ)-certifiably
subgaussian with parameter C > 0 if there for every positive integer k ′ 6 k/2, there exists a degree-ℓ
sum-of-squares proof of the polynomial inequality

∀u ∈ 

d−1

.

D(x)

hx − µ, ui

2k ′



6 C·k

′

D(x)

hx − µ, ui

2

 k′

(5.1)

We ﬁrst note some basic invariance properties of certiﬁable subgaussianity.
Lemma 5.1 (Invariance Under Linear Transformations). Suppose x ∈ d is a random variable with a
(2k, ℓ)-certifiably C-subgaussian distribution. Then, for every matrix A ∈ d×d , the random variable Ax
is (2k, ℓ)-certifiably C-subgaussian.
Proof. Fix any v ∈ d . Invoking certiﬁable subgaussianity of x, we get that there’s a degree ℓ-SoS
proof of the inequality:
′

[hAx, vi 2k ] 

′

′

[hx, A⊤ vi 2k ] 6 (Ck ′ [hx, A⊤ vi 2 ])k  (Ck ′ [hAx, vi 2])k

′



22

Lemma 5.2 (Shifts of certiﬁable subgaussian distributions). Let k ∈  be even and let D be a mean
zero, (k, ℓ)-certifiably C-subgaussian distribution over d . Then, D satisfies for every number k ′ 6 k/2
and vector s ∈ d ,
′
′
(5.2)
hx + s, ui 2)k
∀u ∈ d−1 .
hx + s, ui 2k 6 (Ck ′
D(x)

D(x)

Furthermore, the above inequality Eq. (5.2) has a degree-ℓ sum-of-squares proof.
Next, we construct examples of various classes of distributions that are certiﬁably subgaussian.
We start by noting that the standard gaussian distribution is k-certiﬁably O(1)-subgaussian for
every k.
Lemma 5.3 (Certiﬁable Subgaussianity of Standard Gaussian). Let x be the standard gaussian random
variable with mean 0 and covariance Id on d . Then, x is (2k, 2k)certifiably 1-subgaussian for every k.

Proof. By rotational symmetry of the standard gaussian and the fact that [x 2k
]  Γk , [hx, ui 2k ] 
i
Γk [hx, ui 2]k for some Γk < (2k)k .

Next, we prove a similar fact about uniform mixtures of arbitrary mean gaussians.
Lemma 5.4. Let µ1 , µ2 , . . . µ q be arbitrary. Let D be the mixture of N(µ i , I) with mixture weights
each component. Then, D is (k, k)-certifiably C-subgaussian with C  O(q).

1
q

for

Proof. We have using Lemma 5.2 and certiﬁable C-subgaussianity of the gaussian distribution for
C  O(1),
hx, ui 2k 

1Õ
hx, ui 2k
q
x∼N(µ i ,I)


k
1Õ
2
hx, ui
6 (2Ck)
q
x∼N(µ i ,I)
i

k

i

1
1Õ
 (2Ck) (hu, µ i i 2 + 1)k 6 q k−1 (Ck)k
(hu, µ i i 2 + 1)
q
q
k

hx, ui 2 

On the other hand,
subgaussian.

1
q

Í

i

i (hu,

!k

.

µ i i 2 + 1). Thus, we obtain that D is (k, k)-certiﬁably (2Cq)



Lemma 5.5 (Sampling Preserves Certiﬁable Subgaussianity). Let D be a (k, ℓ)-certifiably subgaussian
distribution on d . Let S be an i.i.d. sample from D of size n > n 0  Ω((d log (d/δ))k/2 /ε 2 ). Then, with
probability at least 1 − δ over the draw of S according to D, the uniform distribution on S, US , is (k, ℓ)
certifiably (C + ε)-subgaussian.
We will use the following matrix concentration inequality in the proof.
Fact 5.6 (Matrix Rosenthal Inequality). Fix p > 1.5. For a finite sequence Pk for k > 1 of independent,
2p
random psd matrices that satisfy kPk k2p < ∞,
k

Õ
k

2p
Pk k2p

! 1/2p

© Õ
6 ­k
« k

1/2

p

Pk k2p + 4p − 2
23

Õ
k

2p
kPk k2p

! 1/4p

2

ª
® .
¬

Proof. Without loss of generality, assume that M2  D xx ⊤ is full-rank. If not, we just work in the
range-space of D xx ⊤ instead, as every sample x would live in this subspace.
−1/2
Now, consider the aﬃne transformation x → M2 x which makes second moment is now
I. By Lemma 5.1, it is enough to show certiﬁable subgaussianity of the distribution after this
transformation. Thus, without loss of generality, we assume that the D is isotropic, that is, M2  I.
Now, by sum-of-squares version of the Cauchy-Schwarz inequality
u
2k ′

±(Ck ′)

k ′ /2

(

k′

US

hx, ui 2) −

US

hx, ui

*

2k ′

 ± (Ck ′)

k ′ /2

′ k ′ /2

′

where k(Ck ′)k /2

US

x ⊗2

 ⊗k′

6 k(Ck )
−

US

′





x ⊗2

US

US

x

⊗2

 ⊗k′

 ⊗k′

−

−

US

x

⊗2k ′

,u

′

US

⊗2k ′

+

x ⊗2k k∞ ku k22k
′

′

(5.3)

′

x ⊗2k k∞ is the spectral norm of the canonical d k × d k ﬂattening



′

⊗k ′

′

− US x ⊗2k . Let T(D) be the corresponding
of the 2k ′ -order tensor T(US )  (Ck ′)k /2 US x ⊗2
tensor for D.
By standard matrix concentration, whenever n 0 ≫ (d log d/ε 2 ), US hx, ui 2 > 1− ε. By standard
argument combining Fact 5.6 with Matrix-Chebyshev inequality, we can obtain that with probability at least 1 − δ over the draw of the sample S of size n > n 0  Ω((d log (d/δ))k/2 /ε 2 ), for any
k ′ 6 k/2,
kT(D) − T(US )k∞ 6 ε.
(5.4)
Thus, along with (5.3), we obtain that:
u
2k ′

′ k ′ /2

(Ck )

2 k′

(

US

hx, ui ) −

Rearranging thus yields that:
u
2k ′



′

US



hx, ui

(Ck ′)k /2 + ε) (

US

2k ′

6

′
εku k22k

US

hx, ui

2

 k′

.

′

′

hx, ui 2)k −

< 2ε



US

hx, ui 2k > 0

completing the proof.


Lemma 5.7 (Certiﬁable Subgaussianity of Mixtures). Let D1 and D2 be two (k, ℓ)-certifiably Csubgaussian distributions with covariances Σ1 , Σ2 respectively. Suppose further that ∆Σ1  Σ2  Σ1 .
Let D be a mixture of D1 and D2 with mixture weights 1 − λ and λ respectively such that λ < ∆−k/2 .
Then, D is (k, ℓ)-certifiably 2C-subgaussian.
Proof. Assume that D1 and D2 are mean zero. We have for any u,
⊢2 D1 hx, ui 2 6 D hx, ui 2 6 (1 − λ + λ∆) D1 hx, ui 2.
2k ′ 6 (1 − λ) Ck
On the other hand, for any k ′ 6 k/2, ⊢2k′
D hx, ui
λ

D2

Ck ′

D2 hx,

 k′
ui 2

6 (1 − λ +

′
λ∆k )

Ck ′

Thus, D is (k, ℓ)-certiﬁably 2C-subgaussian.

D1 hx,

24

 k′
ui 2 .

D1 hx,

ui 2

 k′

+



The proof of this lemma follows immediately from the following proposition.
Proposition 5.8 (Invariance Under Additive Shifts). Let k ∈  be even and let X be a mean zero,
scalar random variable satisfying ( X k )1/k 6 (C X 2 )1/2 . Then, every s ∈ R satisfies ( (X + s)k )1/k 6
2(C (X + s)2 )1/2 .


k 1/k


2 1/2

Proof. By the triangle inequality for L k norms, we have

(X + s)k

 1/k

6

Xk

 1/k

+ |s |. Therefore,

6 C X
+ |s |. On the other hand (X + s)2  X 2 + s 2 using the fact that
(X + s)
X  0. Thus, 2( (X + s)2 )1/2 > ( X 2 )1/2 + |s |. Combining the last two inequalties yields the
lemma.

Lemma 5.9 (Certiﬁable Subgaussianity of Product Distributions). Let D  D1 ⊗ D2 ⊗ · · · ⊗ Dd be
a product distribution on d such that each Di has ean zero, variance 1 and certifiably subgaussian with
constant C i . Then, D is (2k, 2k)-certifiably subgaussian with constant C  max i C i .

Proof. Let 1 6 C  maxi C i . Observe that in that case, Di x 2 j 6 C 2 j · N(0,1) x 2 j for every j 6 2k.
Let D ′ be the distribution obtained by the following process: 1) sample x ∼ D and 2) output x
and −x with uniform probability. By linearity, D hx, ui 2k  D ′ hx, ui 2k .
Call a multi-set S of [d] even if every individual element in S appears even number of times.
Let S be an even multi-set of size |S|  2k. Then, D ′ x S  D x S 6 C |S| · N(0,I) x S using the
sub-gaussianity of each x i .
Now, let u ∈ d . Then, D hx, ui 2  D ′ hx, ui 2  ku k22. We have:
D

hx, ui 2k 
′

Õ

S: even multi-set

[x S ] · u S .

Since S is even, we can upper bound the term corresponding to S above by C |S| N(0,I) x S u S . Thus,
Í
2k 6 C |S|
S S
2k
2k
D ′ hx, ui
N(0,I) x u . Thus, the RHS is at most σ · x∼N(0,I) hx, ui .
S: even multi-set
By rotational symmetry, x∼N(0,I) hx, ui 2k  Γk · ( x∼N(0,I) hx, ui 2)k  Γk · ( D ′ hx, ui 2)k , where
Γk 6 (2k)k .
Thus, D ′ hx, ui 2k 6 (2k)k C 2k ( D ′ hx, ui 2)k . And thus, D hx, ui 2k 6 (2k)k C 2k ( D hx, ui 2)k .
Showing that D is (2k, 2k) certiﬁably C-subgaussian.


6 Applications
6.1 Outlier-Robust Independent Component Analysis
In this section, we apply our robust moment estimation algorithm to give an eﬃcient algorithm for
outlier-robust independent component analysis. We recall the deﬁnition of subgaussian random
variables ﬁrst.
Definition 6.1 (Subgaussian Distributions). A distribution D on  is said to be K-subgaussian if
√
for every p > 1, x∼D [x p ]1/p 6 K( z∼N(0,1) [z p ])1/p 6 K · p. A distribution D on n is said to be
K-subgaussian if for every unit vector v ∈ n , hx, vi is K-subgaussian random variable on .
25

Next, we describe the ICA problem:
Definition 6.2 (Fully Determined Independent Component Analysis). Let A be an unknown mixing matrix in d×d of full rank. Let x be an unknown product distribution with C-subgaussian
marginals satisfying [x 2i ]  1 and [x 4i ] − 1 , 0 for every i. Given samples from y  {Ax} for i.i.d.
samples x, estimate a matrix Â such that there’s a permutation π and some signing σ ∈ {±1} d on
columns of Â satisfying kA i − σi · Â π(i) k2 6 δkA i k for the smallest possible δ.
Remark 6.3.
1. The assumption that [x 2i ]  1 is without the loss of generality. Since we only
see samples from {Ax} one can’t uniquely recover the scalings on the column A i and [x 2i ]
simultaneously for any i.
2. The goal is to recover a matrix A whose columns are close to that of A up to permutation and
signs. This is again, necessary, since the order of the columns cannot be uniquely recovered
from samples of {Ax}.
3. Information theoretically, columns of A can be recovered (up to permutation and scaling) if
at most one x i is has gaussian distribution.
Our main result is a outlier-robust algorithm for ICA. Our algorithm is based on 4th order
tensor decomposition. We will use noise-tolerant generalization of the FOOBI algorithm due to
Cardoso [DLCC07] presented in [MSS16a].
Theorem 6.4 (Outlier-Robust ICA). For every C > 1 and even k ∈ , there exists a (dn)O(k) that
given a (corrupted) sample S ⊆ d of size n outputs component estimates â1 , . . . , â d ∈ d with the
following guarantees: Suppose A ∈ d×d is a non-singular matrix with condition number κ and columns
a1 , . . . , a d ∈ d . Suppose x is a centered random vector with d independent
coordinates such that every
√
2
4
k 1/k
coordinate i ∈ [d] satisfies [xi ]  1, [xi ]−3  γ , 0, and [xi ] 6 Ck. Then, if S is an ε-corrupted
sample of size |S| > n 0 from the distribution {Ax}, where n 0 6 (C + κ + d)O(k) , the component estimates
satisfy with high probability
max min hA−1 â i , A−1 a π(i) i 2 > 1 − δ
π∈S d i∈[d]

for δ < (1 +

1
|γ| )

· O(C 2 k 2 ) · ε 1−4/k .

(6.1)

Our algorithm is simple. It ﬁrst estimates the 2nd and 4th moments of the observed distribution
from an ε-corrupted sample using Algorithm 4.3 and then applies a blackbox tensor decomposition
algorithm to a “whitened” version of 4th moments.
Tensor Decomposition We state the tensor decomposition algorithm (implicit in [MSS16a], analog of Theorem 1.1 for 4th order tensor decomposition) we use before proceeding.
It is convenient to deﬁne the following relaxation of the injective tensor norm in describing the
guarantees of this algorithm.
Definition 6.5 (Sum of Squares Norms). Given a tensor T ∈ (d )⊗2t for any 2t, the degree 2t sum
of squares norm of T, kT ksos2t is deﬁned by
kT ksos2t

 2t



sup
µ(u): deg 2t pseudo-dist over sphere

26



˜ µ(u) hT, u ⊗2t i .

 ⊗4

Fact 6.6 ([MSS16a]). There exists a polynomial time algorithm that given a symmetric 4-tensor T ∈ d ,
outputs a set of vectors {c 1′ , c 2′ , . . . , c ′n } ⊆ d , such that for every orthonormal set {c 1 , c 2 , . . . , c n } ∈ d ,
there’s a permutation π : [n] → [n] satisfying

Õ
T
−
c i⊗4 ksos4 .
kT ksos4
n

max kc i − c ′π(i) k 2 6 O(1) · k
i

i1

2 j−1

Without loss of generality, we can assume that the distribution x satisﬁes [x i ]  0 for every
j > 0. This can be accomplished by negating each sample independently with probability 1/2. The
fraction of corruptions in the modiﬁed sample is still ε.
Algorithm 6.7 (Algorithm for Outlier-Robust ICA via Sum of Squares).
Given: ε-corrupted sample from {Ax}: Y  { y1 , . . . , yn } ⊆ d .
Estimate: Mixing Matrix A ∈ d×d .
Operation:
1. Apply Algorithm 4.3 to estimate the covariance Σ̂ and 4th moments M̂4 of sample
Y.
2. Deﬁne the “whitened” estimated 4th moment tensor M̂4′ so that for every u ∈ d ,
h M̂4′ , u ⊗4 i  h M̂4 , (Σ−1/2 u)⊗4 i.

Let σ1 , σ2 , . . . , σd be the columns of Σ̂−1/2 . Then, the above is equivalent to setting:
Í
M̂4′  j1 , j2 ,..., j4 σ j1 ⊗ σ j2 ⊗ σ j3 ⊗ σ j4 · (M̂4 )(j1 , j2 , j3 , j4 ).

3. Let T̂  M̂4′ − 3I ⊗ I.

4. Apply Tensor Decomposition from Theorem 6.6 to T̂ to recover components
â1 , â2 , . . . , â d .
5. Output Σ̂1/2 Â where Â  (â1 , . . . , â d ).
We now analyze the algorithm. The following lemma shows that the components of the tensor
T that we decompose are close to the columns of A.
Lemma 6.8. Let T  M̂4′ − 3I ⊗ I for the whitenened 4th moments M̂4′ . Let a1 , a2 , . . . , a d be the columns
of A and let b i  (AA⊤ )−1/2 a i be the whitened versions. Then,
k T̂ − γ

Õ
i

1/4

b i⊗4 ksos4 6 δ 4 (2 + 3δ 2 )1/4 .

Proof. First, let’s understand why the tensor decomposition should produce the columns of A if
all the moment estimates were exactly correct. We will then account for the estimation errors.
For the distribution {Ax}, we can compute M2  Σ  [(Ax)(Ax)⊤ ]  AA⊤ using [xx ⊤ ]  I.
Í
Í
Thus, for b i  Σ−1/2 a i , i b i b ⊤
 i Σ−1/2 a i a ⊤
Σ−1/2  I.
i
i
27

Similarly, we have
hM4 , u ⊗4 i 

[hAx, ui 4] 

Õ
i

[x 4i ]ha i , ui 4 + 3

Õ

Õ

4

 (γ − 3)ha i , ui + 3(

i, j

i

[x 2i x 2j ]ha i , ui 2 ha j , ui 2

ha i , ui 2)2 .

The whitened and shifted 4th moment T then, by deﬁnition, satisﬁes
hT, u ⊗4i  hM4′ − 3I ⊗ I, u ⊗4 i  hM4 , (Σ−1/2 u)⊗4 i − 3ku k24

Õ

 γhb i , ui 4 + 3(

i

Í

hb i , ui 2)2 − 3ku k24  γhb i , ui 4 .

Í

As a result, ⊢u,4 hT − γ i b i⊗4 , u ⊗4 i  0 and thus, kT − γ i b i⊗4 ksos4  0.
Now, let Σ̂, M̂4 , M̂4′ , T̂, b̂ i stand for their estimated (via Algorithm 4.3) counterparts. To complete
the analysis, we now understand the error kT − T̂ ksos4 . The proof is then complete by using:
k T̂ − γ

Õ
i

b i⊗4 ksos4 6 k T̂ − T ksos4 + kT − γ

Õ
i

b i⊗4 ksos4  k T̂ − T ksos4 .

First, we check that our algorithm (Algorithm 4.3) can indeed be applied to this setting. Since
x is a product distribution with ith marginal C i -subgaussian. By Lemma 5.9, x is (k, k)-certiﬁably
C-subgaussian. By Lemma 5.1, {Ax} is (k, k)-certiﬁably subgaussian. 13 Finally, applying Lemma
5.5 shows that with probability at least 1 − 1/n 2 , the sample the facts above all hold for sample Y
of size m > Ω̃(d k/2 /ε 2 ).
We will use Σ, M4 , M4′ , T, b i  Σ−1/2 a i for the covariance, 4th moment, whitenened 4th moment,
shifted, whitenened 4th moment and orthogonalized columns of A respectively.
By Corollary 4.6, for every u, hu, Σ̂ui  (1 ± δ 2 )hu, Σui for some δ2 < O(Ck)ε 1−2/k . Thus:
(1 − δ 2 )I 
By Theorem 1.3, for δ 4  O(C 2 k 2 )ε 1−4/k ,
u,k
4

Thus,
u,k
4



Õ
i



b̂ i b̂ ⊤
i  (1 + δ 2 )I.

(6.2)



± hM4 − M̂4 , u ⊗4 i 6 δ 4 · hΣ, u ⊗2 i 2 + hΣ̂, u ⊗2 i 2 .

−1/2

± hM4 − M̂4 , (M̂2







−1/2

u)⊗2 i 2 + h M̂2 , (M̂2

6 δ4 (2 + 3δ2 )ku k24

Or,
u,k
4

−1/2

u)⊗4 i 6 δ4 · hM2 , (M̂2

(6.3)

u)⊗2 i 2



± hM4′ − M̂4′ , u ⊗4 i 6 δ 4 (2 + 3δ 2 )ku k24

13Negating samples with probability 1/2 preserves 2nd and 4th moments and doesn’t aﬀect certiﬁable subgaussianity.

28

Since M4′ − M̂4′  T − T̂, this gives:
u,k
4



± hT − T ′ , u ⊗4 i 6 δ 4 (2 + 3δ 2 )ku k24

(6.4)



The proof of Theorem 1.4 then follows easily.
Proof of Theorem 1.4. We are now ready to apply Fact 6.6 to obtain that decomposing the tensor
T recovers components c 1 , c 2 , . . . , c d such that there’s a permutation π : [d] → [d] such that
1
kc i − b π(i) k22 6 O(1) · |γ|
· δ 4 (2 + 3δ 2 ).
Applying the reverse whitening transform â i  Σ̂1/2 c i then recovers the set of vectors

a1 , a2 , . . . , a d satisfying the requirements of Theorem 1.4.

6.2 Outlier-Robust Mixtures of Gaussians
In this section, we describe our outlier-robust algorithm for learning mixtures of spherical gaussians. As before, this algorithm can be seen as a direct analog of the non-robust algorithm that
learns mixtures of spherical gaussians with linearly independent means by decomposing the 3rd
order moment tensor of the observed sample.
Theorem 6.9 (Outlier-Robust Mixtures of Gaussians). Let D be mixtures of N(µ i , I) for i 6 q with
uniform14 mixture weights. Assume that µ i s are linearly independent and, further, let κ be the smallest
Í
non-zero eigenvalue of 1q i µ i µ⊤
.
i
Given an ε-corrupted sample of size n > n 0  Ω((κd log (d))k/2 ), for every k > 4, there’s a poly(n)d O(k)
time algorithm that recovers µ̂1 , µ̂2 , . . . µ̂ q so that there’s a permutation π : [q] → [q] satisfying
max k(
i

p
p
√
1Õ
µ i µ⊤
)−1/2 (µ̂ i − µ π(i) )k 6 O( qk)ε 1−1/k / κ + O(qk)ε 3/2−2/k /κ + O( qk)ε 1/3−1/k .
i
q
i

We will rely on 3rd order tensor decomposition algorithm (Theorem 1.1) from [MSS16b] here.

 ⊗3

Fact 6.10 ([MSS16a]). There exists a polynomial time algorithm that given a symmetric 3-tensor T ∈ d ,
outputs a set of vectors {c 1′ , c 2′ , . . . , c ′n } ⊆ d , such that for every orthonormal set {c 1 , c 2 , . . . , c n } ∈ d ,
there’s a permutation π : [n] → [n] satisfying

Õ
T
−
c i⊗3 ksos4 .
6 O(1) · k
kT ksos4
n

max kc i −
i

c ′π(i) k 2

i1

14While our algorithm generalizes naturally to arbitrary mixture weights, we restrict to this situation for simplicity

29

Algorithm 6.11 (Algorithm for Outlier-Robust Mixture of Spherical Gaussians via Sum of
Squares).
Given: ε-corrupted sample X from
Estimate: µ1 , µ2 , . . . , µ q .

Í

1
i q N(µ i , I)

for µ1 , µ2 , . . . , µ q .

Operation:
1. Apply Algorithm 4.3 to estimate the ﬁrst three moments, M̂1 , M̂2 , and M̂3 of sample
X0.
2. Deﬁne the “whitened” version of the estimated 3rd moment tensor T̂ so that for
every u ∈ d , hT̂, u ⊗3 i  h M̂3 − 3M̂1 ⊗ I, ((M̂2 − I)−1/2 u)⊗3 i.

Concretely, set M̂1′  (M̂2 −I)−1/2 M̂1 . Let σ1 , σ2 , . . . , σd be the columns of (M̂2 −I)−1/2 .
Í
′
Set T̂  j1 , j2 , j3 σ j1 ⊗ σ j2 ⊗ σ j3 · (M̂3 − 3M̂1 ⊗ I)(j1 , j2 , j3 ).

3. Apply Tensor Decomposition from Theorem 6.6 to T̂ to recover components
ν̂1 , ν̂2 , . . . , ν̂q .
4. Output µ̂ i  (M̂2 − I)1/2 ν̂i for i 6 q as the estimated means.
The following is the main technical lemma required in the proof of 6.9.
Lemma 6.12. Let Y be an ε-corrupted sample of size n > n 0  Ω((κd log d)k /ε 2 ) from D0 

Í

1
i q N(µ i , I).

Let M̂ i be the estimated ith moment of D0 given by Algorithm 4.3 when run on Y. Let M̂ i′ be the whitening
of the estimated ith moment of D0 as in Algorithm 6.11.
Then, with high probability over the draw of the ε-corrupted sample Y, for T̂  M̂3′ − 3M̂1′ ⊗ I and
νi  (M2 − I)−1/2 µ i for i 6 q and δ i  O(qk)i/2 ε 1−i/k for i  1, 2, 3,
k T̂ −

√
1 Õ ⊗3
1/3
1/2
νi ksos4 6 δ1 / κ + δ1 δ 2 /κ + 2δ 3 .
q
i

Proof. As before, we will ﬁrst analyze the tensor T assuming all estimates of the moments are
exactly correct and then account for the estimation errors. We will write M1 , M2 , M3 for the true
ﬁrst three moments of the input mixture of gaussian and Σ  M2 − I. We will let M̂1 , M̂2 , M̂3 , Σ̂
be the corresponding quantities estimated by our algorithm. We will write M1′ , M3′ for the true
whitened moments and M̂1′ , M̂3′ , the estimated counterparts.
Í
By direct computtion, hM3 , u ⊗3 i  1q i (hµ i , ui 3 + 3hµ i , ui). Thus, hM3 − 3M1 ⊗ I, u ⊗3 i 
1
q

Í

i hµ i ,

ui 3.

Now, by deﬁnition, hM3′ − 3M1′ ⊗ I, u ⊗3i  hM3 − 3M1 ⊗ I, (Σ−1/2 u)⊗3 i.
Í
As a result, for T  M3′ − 3M1′ ⊗ I, we have: kT − 1q i (Σ−1/2 µ i )⊗3 ksos4  0.

We now account for the estimation errors and bound kT − T̂ ksos4 . First, we use Lemma 5.4 and
Lemma 5.5 to conclude that with high probability, the uniform distribution on the input sample is
k-certiﬁably O(q)-subgaussian.
30

Now (via triangle inequality),
kT − T̂ ksos4 6 kM3′ − M̂3′ ksos4 + 3kM1′ ⊗ I − M̂1′ ⊗ I ksos4 .

(6.5)

|hM1′ − M̂1′ , ui| 6 |hΣ−1/2 (M1 − M̂1 ), ui| + |h(Σ−1/2 − Σ̂−1/2 )M̂1 , ui|

(6.6)

Now,

p

From Theorem 1.2, for δ 1  O( qk)ε 1−1/k , δ 2  O(qk)ε 1−2/k , hM1 − M̂1 , Σ−1/2 ui 6
δ1 hu, M2Σ−1 ui 1/2 .
Í
Í
Using that M2  I + 1q i µ i µ⊤
, M2 Σ−1  ( 1q i µ i µ⊤
)−1 + I. Thus, the ﬁrst term in (6.6) can be
i
i
√
upper bounded by δ 1 ku k2 + δ 1 hu, Σ−1 , ui 1/2 6 δ 1 (1 + 1/ κ)ku k2 .
From Theorem 1.2, for any u, hu, (M2 − M̂2 )ui 6 δ 2 hu, M2 ui. Writing M2  Σ + I for Σ  1q µ i µ⊤
,
i

we obtain that hu, Σ̂ui  (1 ± O(δ2 /κ))hu, Σui.
1/2 √
Thus, h(Σ−1/2 − Σ̂−1/2 )M̂1 , ui 6 O(δ2 / κ) · hΣ−1/2 M̂1 , ui. Writing hΣ−1/2 M̂1 , ui  hΣ−1/2 (M̂1 −
Í
M1 )ui + hΣ−1/2 M1 ui and using that 1q i hΣ−1/2 µ i , ui 6 √1q ku k2, we ﬁnally obtain anupper bound


√
√
√
1/2 √
1/2
on (6.6) of: δ1 (1 + 1/ κ) + O(δ2 / κ)δ 1 1 + 1/ κ + √1q ku k2 6 O(δ 2 δ1 /κ + δ1 / κ)ku k2 .
For the ﬁrst term of (6.5), we apply Theorem 1.3 to get that for δ 3  O(qk)3/2 ε 1−3/k ,
u
6
u

′

h M̂3 − M3 , u ⊗3 i 2 6 δ 23 hM2 , u ⊗2 i 3 6 2δ3 hΣ, u ⊗2i 3 ku k23 .

As a result, 6 h M̂3 − M3′ , u ⊗3 i 2 6 2δ23 ku k23 giving us an upper bound on the ﬁrst term in the RHS
of (6.5).
This completes the proof.


We can now complete the proof of Theorem 6.9.
Proof of Theorem 6.9. Applying Fact 6.10 to T̂ as in Lemma 6.12 yields that the resulting components
√
1/3
1/2
satisfy min π:[q]→[q] max i kνi − ν̂π(i) k2 6 O(1)δ 1 / κ+ δ1 δ 2 /κ+2δ 3 . The theorem now follows. 

7 Information Theoretic Lower Bounds
In this section, we give simple examples that show that our robust moment estimation results
achieve sharp information theoretic limits for (certiﬁably) subgaussian distributions. In both
results below, for any choice of k, we will show that there are two distributions that 1) both satisfy
(k, k)-certiﬁable C-subgaussianity and 2) are at most ε apart in total variation distance but have lowdegree moments that are exactly as far apart as the estimation error in our algorithm. In fact, these
distributions are really simple and in fact just one dimensional and thus the certiﬁability follows
from generic resuls about sum of squares representations for univariate non-negative polynomials.
We begin by showing a lower bound on the error incurred in robust estimation of mean.
Lemma 7.1. There exist two distributions D1 and D2 on  such that:
1. D1 , D2 are both 2-subgaussian in all moments up to k.
31

2. dTV (D1 , D2 ) 6 ε.
3. k

D1

x−

D2

√

xk >

kε 1−1/k .

Proof. Let D1 be the standard gaussian distribution N(0, 1). Let D2 be the distribution
that with
√
−1/k
. Observe
probability 1 − ε outputs a sample from N(0, 1) and with probability ε, outputs kε
that by construction, dTV (D1 , D2 ) 6 ε.√
The mean of D2 is easily seen to be kε 1−1/k . The variance of D2 can be computed to be between
1 and (1 + kε 1−2/k ).
It is standard to check that D1 is C-subgaussian in all moments for C  1 < 2. For D2 , we
√ k
observe that D2 (x − D2 x)2k 6 D2 x 2k 6 (1− ε)k k/2 + ε k ε −1 < 2k k/2 6 2k k/2 ( D2 (x − D2 x)2 )k/2 .
Thus, D2 is 2-subgaussian.
This completes the proof.


Next, we show that our robust estimation error for the covariance and higher moments are
also tight up to constant factors. As before, we will construct two zero-mean, 2-subgaussian
distributions on  that are at most ε apart in total variation but the variances diﬀer by at least
k/2ε 1−2/k . The construction is very similar to the one above for mean estimation.
Lemma 7.2. There exist two mean zero distributions D1 and D2 on  such that:
1. D1 , D2 are both 2-subgaussian in all moments up to k.
2. dTV (D1 , D2 ) 6 ε.
3. For any ε < 1/2, |

D1

4. For any ε < 2−k/2r , |

x2 −
D1

D2

x 2r −

x 2 | > 2k ε 1−2/k .
D2

x 2r | >

kr
2

ε 1−2r/k .

Proof. Let D1 be the standard gaussian distribution N(0, 1). Let D2 be the distribution that with
probability
1 − ε outputs a sample from N(0, 1) and with probability ε, outputs a uniform element
√
−1/k
. Observe that by construction, dTV (D1 , D2 ) 6 ε.
of ± kε
The mean of D1 , D2 is easily seen to be 0. The same computation as in the proof of Lemma 7.1
establishes that D1 , D2 are both 2-subgaussian.
The variance of D1 is 1 and for D2 can be computed to be at least (1 − ε)(1 + kε 1−2/k − kε 2−2/k ) >
1 + 2k ε 1−2/k .
The 2rth moment of D1 is ∼ k r and for D2 can be computed to as: k r (1 − ε + ε 1−2r/k ) >
r
k (1 + 0.5ε 1−2r/k ) for ε < 2−k/2r .
This completes the proof.



32

References
[AFH+ 12]

Anima Anandkumar, Dean P. Foster, Daniel J. Hsu, Sham Kakade, and Yi-Kai Liu, A
spectral algorithm for latent dirichlet allocation, NIPS, 2012, pp. 926–934. 1

[AGH+ 14] Animashree Anandkumar, Rong Ge, Daniel J. Hsu, Sham M. Kakade, and Matus
Telgarsky, Tensor decompositions for learning latent variable models, Journal of Machine
Learning Research 15 (2014), no. 1, 2773–2832. 1
[BCMV14] Aditya Bhaskara, Moses Charikar, Ankur Moitra, and Aravindan Vĳayaraghavan,
Smoothed analysis of tensor decompositions, STOC, ACM, 2014, pp. 594–603. 1, 9
[Ber06]

Thorsten Bernholt, Robust estimators are hard to compute, Tech. report, Technical Report/Universität Dortmund, SFB 475 Komplexitätsreduktion in Multivariaten Datenstrukturen, 2006. 1

[BKS15]

Boaz Barak, Jonathan A. Kelner, and David Steurer, Dictionary learning and tensor
decomposition via the sum-of-squares method, STOC, ACM, 2015, pp. 143–151. 1, 2, 36

[BM16]

Boaz Barak and Ankur Moitra, Noisy tensor completion via the sum-of-squares hierarchy,
COLT, JMLR Workshop and Conference Proceedings, vol. 49, JMLR.org, 2016, pp. 417–
445. 2

[BS16]

Boaz Barak and David Steurer, Proofs, beliefs, and algorithms through the lens of sum-ofsquares, 2016, Lecture notes in preparation, available on http://sumofsquares.org.
13

[CJN17]

Yeshwanth Cherapanamjeri, Prateek Jain, and Praneeth Netrapalli, Thresholding based
outlier robust PCA, COLT, Proceedings of Machine Learning Research, vol. 65, PMLR,
2017, pp. 593–628. 1

[CR09]

Emmanuel J. Candès and Benjamin Recht, Exact matrix completion via convex optimization, Foundations of Computational Mathematics 9 (2009), no. 6, 717–772. 2

[CRTV05]

Emmanuel J. Candès, Mark Rudelson, Terence Tao, and Roman Vershynin, Error
correction via linear programming, FOCS, IEEE Computer Society, 2005, pp. 295–308. 2

[CSV17]

Moses Charikar, Jacob Steinhardt, and Gregory Valiant, Learning from untrusted data,
STOC, ACM, 2017, pp. 47–60. 1

[DKK+ 16]

Ilias Diakonikolas, Gautam Kamath, Daniel M. Kane, Jerry Li, Ankur Moitra, and Alistair Stewart, Robust estimators in high dimensions without the computational intractability,
FOCS, IEEE Computer Society, 2016, pp. 655–664. 1, 5, 7, 10

[DKK+ 17a]
[DKK+ 17b]

, Being robust (in high dimensions) can be practical, CoRR abs/1703.00893 (2017). 5
, Robustly learning a gaussian:
abs/1704.03866 (2017). 6, 7
33

Getting optimal error, efficiently, CoRR

[DLCC07]

Lieven De Lathauwer, Joséphine Castaing, and Jean-François Cardoso, Fourth-order
cumulant-based blind identification of underdetermined mixtures, IEEE Trans. Signal Process. 55 (2007), no. 6, part 2, 2965–2973. MR 2473563 8, 26

[GLS81]

M. Grötschel, L. Lovász, and A. Schrĳver, The ellipsoid method and its consequences in
combinatorial optimization, Combinatorica 1 (1981), no. 2, 169–197. MR 625550 14, 15

[GM15]

Rong Ge and Tengyu Ma, Decomposing overcomplete 3rd order tensors using sum-of-squares
algorithms, APPROX-RANDOM, LIPIcs, vol. 40, Schloss Dagstuhl - Leibniz-Zentrum
fuer Informatik, 2015, pp. 829–849. 1

[Gro11]

David Gross, Recovering low-rank matrices from few coefficients in any basis, IEEE Trans.
Information Theory 57 (2011), no. 3, 1548–1566. 2

[GVX13]

Navin Goyal, Santosh Vempala, and Ying Xiao, Fourier PCA, CoRR abs/1306.5825
(2013). 8

[HK13a]

Daniel Hsu and Sham M. Kakade, Learning mixtures of spherical Gaussians: moment
methods and spectral decompositions, ITCS’13—Proceedings of the 2013 ACM Conference
on Innovations in Theoretical Computer Science, ACM, New York, 2013, pp. 11–19. MR
3385380 9

[HK13b]

Daniel J. Hsu and Sham M. Kakade, Learning mixtures of spherical gaussians: moment
methods and spectral decompositions, ITCS, ACM, 2013, pp. 11–20. 1, 10

[HL17]

Sam B. Hopkins and Jerry Li, Mixture models, robustness, and sum of squares proofs, 2017.
5

[HRRS11]

Frank R Hampel, Elvezio M Ronchetti, Peter J Rousseeuw, and Werner A Stahel, Robust
statistics: the approach based on influence functions, vol. 114, John Wiley & Sons, 2011. 1

[HSS15]

Samuel B. Hopkins, Jonathan Shi, and David Steurer, Tensor principal component analysis
via sum-of-square proofs, COLT, JMLR Workshop and Conference Proceedings, vol. 40,
JMLR.org, 2015, pp. 956–1006. 2

[Hub11]

Peter J Huber, Robust statistics, International Encyclopedia of Statistical Science,
Springer, 2011, pp. 1248–1251. 1

[KMV10]

Adam Tauman Kalai, Ankur Moitra, and Gregory Valiant, Efficiently learning mixtures
of two gaussians, STOC, ACM, 2010, pp. 553–562. 1

[KS17]

Pravesh K. Kothari and Jacob Steinhardt, Better agnostic clustering via relaxed tensor
norms, 2017. 4, 5

[Las01]

Jean B. Lasserre, New positive semidefinite relaxations for nonconvex quadratic programs,
Advances in convex analysis and global optimization (Pythagorion, 2000), Nonconvex
Optim. Appl., vol. 54, Kluwer Acad. Publ., Dordrecht, 2001, pp. 319–331. MR 1846160
14
34

[LRV16]

Kevin A. Lai, Anup B. Rao, and Santosh Vempala, Agnostic estimation of mean and
covariance, FOCS, IEEE Computer Society, 2016, pp. 665–674. 1, 6, 9

[MMY06]

RARD Maronna, R Douglas Martin, and Victor Yohai, Robust statistics, John Wiley &
Sons, Chichester. ISBN, 2006. 1

[MR05]

Elchanan Mossel and Sébastien Roch, Learning nonsingular phylogenies and hidden
markov models, STOC, ACM, 2005, pp. 366–375. 1

[MSS16a]

Tengyu Ma, Jonathan Shi, and David Steurer, Polynomial-time tensor decompositions with
sum-of-squares, FOCS, IEEE Computer Society, 2016, pp. 438–446. 1, 2, 8, 26, 27, 29

[MSS16b]

, Polynomial-time tensor decompositions with sum-of-squares, CoRR abs/1610.01980
(2016). 13, 29

[MV10]

Ankur Moitra and Gregory Valiant, Settling the polynomial learnability of mixtures of
gaussians, FOCS, IEEE Computer Society, 2010, pp. 93–102. 1

[Nes00]

Yurii Nesterov, Squared functional systems and optimization problems, High performance
optimization, Appl. Optim., vol. 33, Kluwer Acad. Publ., Dordrecht, 2000, pp. 405–440.
MR 1748764 14

[Par00]

Pablo A Parrilo, Structured semidefinite programs and semialgebraic geometry methods in
robustness and optimization, Ph.D. thesis, California Institute of Technology, 2000. 14

[Pea94]

Karl Pearson, Contributions to the mathematical theory of evolution, Philosophical Transactions of the Royal Society of London. A 185 (1894), 71–110. 1

[PS17]

Aaron Potechin and David Steurer, Exact tensor completion with sum-of-squares, COLT,
Proceedings of Machine Learning Research, vol. 65, PMLR, 2017, pp. 1619–1673. 2

[Rec11]

Benjamin Recht, A simpler approach to matrix completion, Journal of Machine Learning
Research 12 (2011), 3413–3430. 2

[SCV17]

Jacob Steinhardt, Moses Charikar, and Gregory Valiant, Resilience: A criterion for learning in the presence of arbitrary outliers, CoRR abs/1703.04940 (2017). 5

[Sho87]

N. Z. Shor, Quadratic optimization problems, Izv. Akad. Nauk SSSR Tekhn. Kibernet.
(1987), no. 1, 128–139, 222. MR 939596 14

[Tuk75]

John W. Tukey, Mathematics and the picturing of data, 523–531. MR 0426989 1

[VX15]

Santosh Vempala and Ying Xiao, Max vs min: Tensor decomposition and ICA with nearly
linear sample complexity, COLT, JMLR Workshop and Conference Proceedings, vol. 40,
JMLR.org, 2015, pp. 1710–1723. 8

35

A

Sum-of-squares toolkit

We will rely on the following sum-of-squares version of the inequality between the arithmetic and
geometric mean.
Lemma A.1 ([BKS15, Lemma A.1]).
w
k

{w 1 > 0, . . . , w n > 0}
Lemma A.2. For all even k,

a,b
k



(

Ö

wi 6

i

Í

i

w ik

k

)

.

(a + b)k 6 2k−1 (a k + b k ) .

(A.1)

Proof. The lemma holds for k  2, because
a,b
2



2a 2 + 2b 2 − (a + b)2  a 2 + b 2 − 2ab  (a − b)2 > 0 .

(A.2)

By iterating this inequality k/2 times, we obtain



a,b
k

(a + b)k  ((a + b)2 )k/2 6 2k/2 (a 2 + b 2 )k/2 .

(A.3)

By the binomial identity and Lemma A.1,
a,b
k

2

2 k/2

(a + b )



k/2  
Õ
k
i0

6

a 2i b k−2i

k/2   
Õ
k
i0



i

i

k/2   
Õ
k

i

i0
k/2−1

2

2i
k

1
2

(A.4)

k

k−2i
k

·a +
· ak +

1
2

· (a k + b k ) .

· bk

·b



k







k
.
Here, we used for the second equality that ki  k−i
Combining the previous two inequalities yields the desired inequality,
a,b
k

(a + b)k 6 2k/2 · (a 2 + b 2 )k/2 6 2k−1 · (a k + b k ) .

(A.5)


Lemma A.3. For any even k,



fk 6 1

f
k

36



f 61 .

(A.6)

Proof. We apply Lemma A.1 for k/2 and substitute w 1  f 2 and w 2 , w 3 , . . . , w k/2  1 to obtain
f
k

f2 6

f k +(k/2−1)
.
k/2

Thus,





f

fk 6 1

k

f2 6 1 .

(A.7)

At the same time, since 2 − 2 f  1 − f 2 + (1 − f )2 , we have





f

f2 6 1

2

1− f > 0 .

(A.8)

Combining (A.7) and (A.8) yields the desired conclusion.
Lemma A.4. For every even k ∈  and all δ < 0.1,



where δ′  100δ.

( f − 1)k 6 δ k ( f + 1)k



f
k



1 − δ′ 6 f 6 1 + δ′ ,

(A.9)

Proof. By Lemma A.2,
f

Therefore, for δ′ 




4δ k
1−2δ ,



k

( f + 1)k  (( f − 1) + 2)k 6 2k ( f − 1)k + 4k .

f

( f − 1)k 6 δ k ( f + 1)k

By Lemma A.3, ( f − 1)k 6 (δ′)k
equality gives the desired inequality.



k



( f − 1)k 6 (δ′)k .

1 − δ′ 6 f 6 1 + δ , which together with the previous in-



Proposition A.5. Let A be a set of polynomial equality axioms in variable x such that:
A

2t

p(x, u) > 0,

for a polynomial p with total degree at most 2t.
Then, for any pseudo-distribution D of degree 2t on x satisfying A,

u
2t

˜ D p(x, u) > 0.

Proof. Let A  {q 1  0, q 2  0, . . . , q r  0}. Since A 2t p(x, u) > 0, there are degree 6 2t
(in x) polynomials 11 (x, u), 12(x, u), . . . 1r ′ (x, u) and polynomials h1 (x, u), . . . , h r (x, u) such that
Í
Í
p(x, u)  i 1i (x, u)2 + j h j (x, u)q j (x) where the degree of any h j (x, u) · q j (x) in x is at most 2t.
Í
Taking pseudo-expectations with respect to D, we have: ˜ D p(x, u)  i ˜ D 1i (x, u)2 +
Í ˜
Í ˜
2
j D h j (x, u)q j (x) 
i D 1 i (x, u) , where we used that the second term must vanish as D
satisﬁes A.
Now, observe that by matching degrees on both sides as polynomial in u, no ˜ D 1i (x, u)2 has
degree more than 2t.
Further, each 1i (x, u)2 can be written as hG(x)G(x)⊤ , (1, u)⊗2t i for some matrix valued polynomial G(x).

37

Í ˜
⊤
⊗2t i.
i h D [G(x)G(x) ], (1, u)
Taking

pseudo-expectations

with

respect

to

D

yields

that:

˜ D p(x, u)



Now, we claim that ˜ D G(x)G(x)⊤ is psd. To see this, let w be any vector and consider the
quadratic form of this matrix in w. Then, the quadratic form equals ˜ D hG(x), wi 2 which is nonnegative by the positivity of pseudo-expectations.
u

Thus, ˜ D p(x, u) is a sum of squares polynomial in u implying that 2t ˜ D p(x, u) > 0.

38

