Outlier-robust moment-estimation via sum-of-squares

Introduction

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 effort are computationally efficient only for low-dimensional distributions (because the running time is say exponential in the dimension $d$) [Ber06].

A recent line of research developed the first robust estimators for basic parameter estimation problems (e.g., estimating the mean and covariance matrix of a Gaussian distribution) that are computationally efficient 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 efficient algorithms to estimate low-degree moments of distributions. Our estimators succeed under significantly 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 (affine transformations of) product distributions.

Besides these qualitative differences, our algorithms also offer 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 $\varepsilon$ of corruptions—match, for the first time in this generality, information-theoretic lower bounds.

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 effect on the empirical moments. Instead, we apply method-of-moment-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 (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-efficient provable parameter estimation that has been developed in the context of the sum-of-squares method. We say that a parameter estimation problem satisfies identifiability if it is information-theoretically possible to recover the desired parameter from the available data (disregarding computational efficiency). The key idea of this paradigm is that a proof of identifiability can be turned into an efficient 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].

Next we describe the form of our identifiability proofs for outlier-robust estimation.

Identifiability in the presence of outliers

$x_{i}^{\prime}=x_{i}$ for all but an $\varepsilon$ fraction of the indices $i\in[n]$,

the uniform distribution over $X^{\prime}$ is in $\mathcal{C}$.

Note that given $X$ and $X^{\prime}$, we can efficiently check the above conditions, assuming that the conditions on the low-degree moments that define $\mathcal{C}$ are efficiently checkable (which they will be).

Also note that the above notion of proof is complete in the following sense: if $X$ is indeed an $\varepsilon$-corruption of a typicalThe 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 satisfied with high probability. sample $X^{0}$ from a distribution $D\in\mathcal{C}$, then there exists a set $X^{\prime}$ that satisfies 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 $\varepsilon$-corruption of a (typical) sample $X^{0}$ from a distribution $D\in\mathcal{C}$ and $X^{\prime}$ satisfies the above conditions, then the empirical mean $\mu^{\prime}\mathrel{\mathop{:}}=\tfrac{1}{n}\sum_{i=1}^{n}x^{\prime}_{i}$ of the uniform distribution over $X^{\prime}$ is close to $\mu$. We can rephrase this soundness as the following concise mathematical statement, which we prove in Section 2: if $D$ and $D^{\prime}$ are two distributions in $\mathcal{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 efficient algorithm for outlier-robust mean estimation (see Section 4).

Outlier-robustness and (certifiable) subgaussianity

Here, $\mu_{D}$ is the mean of the distribution $D$.

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 polynomial-time algorithm to estimate the covariance matrix $\Sigma_{D}$ of $D$ up to constant factors (in the Löwner order sense) assuming again that $Ck\cdot\varepsilon^{1-2/k}\ll 1$.

In Section 5, we show that a wide range of distributions are certifiably subgaussian and that many operations on distributions preserve this property. In particular, (affine transformations of) products of scalar-valued subgaussian distributions and mixtures thereof satisfy this property. In a subsequent work, Kothari and Steinhardt [KS17] show that certifiable subgaussianity holds for distributions that satisfy a Poincaré inequality (which includes all strongly log-concave distributions).

Sum-of-squares and quantifier alternation

We remark that previous work on computationally-efficiently outlier-robust estimation also used convex optimization techniques albeit in different ways. For example, Diakonikolas et al. solve an implicitly-defined convex optimization problem using a customized separation oracle [DKK+16]. (Their optimization problem is implicit in the sense that it is defined in terms of the uncorrupted sample which we do not observe.)

An unusual feature of the aforementioned system of equations $\mathcal{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^{\prime}$ such that the uniform distribution $X^{\prime}$ is certifiably subgaussian. It is interesting to note that in this way we can use sum-of-squares as an approach to solve $\exists\,\forall$-problems as opposed to just the usual $\exists$-problems. (The current problem is an $\exists\,\forall$-problem in the sense that we want to find $X^{\prime}$ such that for all vectors $u$ the inequality Eq. 1.2 holds for the uniform distribution over $X^{\prime}$.)

We remark that the idea of using sum-of-squares to solve problems with quantifier alternation also plays a role in control theory (where the goal is find 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 first that uses this idea for the design of computationally-efficient 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 Results

Without any assumptions about the underlying distribution, the best known efficient algorithms for robust mean estimation incur an estimation error that depends on the spectral norm of the covariance matrix $\Sigma$ of the underlying distribution and is proportional to $\sqrt{\varepsilon}$ (where $\varepsilon>0$ is the fraction of outliers) [DKK+17a, SCV17]. Concretely, given an $\varepsilon$-corrupted sample of sufficiently larger polynomial size from a distribution $D$, they compute an estimate $\hat{\mu}$ for the mean $\mu$ of $D$ such that with high probability $\lVert\hat{\mu}-\mu\rVert\leqslant O(\sqrt{\varepsilon})\cdot\lVert\Sigma\rVert^{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 $\varepsilon$ and the spectral norm of $\Sigma$.

In the following theorem, we show that better bounds for the mean estimation error are possible for large classes of distributions. Concretely, we assume (certifiable) 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 $\varepsilon$ (instead of a $\sqrt{\varepsilon}$ bound as for the unconditional mean estimation before we obtain an $\varepsilon^{1-1/k}$ if we assume a bound on the degree-$k$ 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 certifiable subgaussianity does not restrict the covariance matrix in any way.)

For the last two bounds, we assume in addition $Ck\cdot\varepsilon^{1-2/k}\leqslant\Omega(1)$.This notation means that we require $Ck\cdot\varepsilon^{1-2/k}\leqslant c_{0}$ for some absolute constant $c_{0}>0$ (that could in principle be extracted from the proof).

Note that the second guarantee for the mean estimation error $\mu-\hat{\mu}$ is stronger because $\lVert\mu-\hat{\mu}\rVert\leqslant\lVert\Sigma\rVert^{1/2}\cdot\lVert\Sigma^{-1/2}(\mu-\hat{\mu})\rVert$. We remark that $\lVert\Sigma^{-1/2}(\mu-\hat{\mu})\rVert$ is the Mahalanobis distance between $\hat{\mu}$ and $D$.

Previous work for robust covariance estimation [LRV16, DKK+17b] work with Frobenius norms for measuring the estimation error $\Sigma-\hat{\Sigma}$ and obtain in this way bounds that can be stronger than ours. However, it turns out that assuming only $k$-certifiable 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.

Multiplicative vs. additive estimation error

Another benefit 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\pm\delta)$ 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.

Robust Estimation of Higher Moments

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, it turns out that our algorithm not only guarantees an injective-norm bound for the error but also a good certificate for this bound, in form of a low-degree sum-of-squares proof. It turns out that this kind of certificate 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 certified by a low-degree sum-of-squares proof.

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 $\varepsilon$. For example, we show that there are two $k$-certifiably $O(1)$-subgaussian distributions with statistical distance $\varepsilon$ but means that are $\Omega(\sqrt{k}\varepsilon^{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 efficient algorithms for ICA [GVX13, VX15].

Outlier robust version of ICA was considered as an application of the outlier-robust mean and covariance estimation problems in [LRV16]. They sketched an algorithm with the guarantee that the relative error in the columns of $A$ is at most $\varepsilon\sqrt{\log{(d)}}\operatorname{poly}(\kappa)$.

In particular, this guarantee is meaningful only if the fraction of outliers $\varepsilon\ll\frac{1}{\sqrt{\log{(d)}}\operatorname{poly}(\kappa)}$. Here, 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 $\kappa$ 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.

The quantity $\langle A^{-1}\hat{a}_{i},A^{-1}a_{\pi(i)}\rangle$ is closely related to the Mahalanobis distance between $\hat{a}_{i}$ and $a_{\pi(i)}$ with respect to the distribution $\{A\mathbf{x}\}$

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

Let $D$ be mixtures of $\mathcal{N}(\mu_{i},I)$ for $i\leqslant q$ with uniformWhile our algorithm generalizes naturally to arbitrary mixture weights, we restrict to this situation for simplicity mixture weights. Assume that $\mu_{i}$s are linearly independent and, further, assume that $\kappa$, the smallest non-zero eigenvalue of $\frac{1}{q}\sum_{i}\mu_{i}\mu_{i}^{\top}$ is $\Omega(1)$.

Given an $\varepsilon$-corrupted sample of size $n\geqslant n_{0}=\Omega((d\log{(d)})^{k/2}/\varepsilon^{2})$, for every $k\geqslant 4$, there’s a $\operatorname{poly}(n)d^{O(k)}$ time algorithm that recovers $\hat{\mu}_{1},\hat{\mu}_{2},\ldots\hat{\mu}_{q}$ so that there’s a permutation $\pi:[q]\rightarrow[q]$ satisfying

Diakonikolas et. al. [DKK+16] gave an outlier-robust algorithm that learns mixtures of $q$ gaussians with error $\approx q\sqrt{\varepsilon}$ 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 difference 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 fixed polynomial time [HK13b].

Robust identifiability of low-degree moments

In this section, we show that low-degree moments of certifiably subgaussian distributions are identifiable in the presence of outliers. As we will formalize in Section 4, these proofs of identifiability 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).

find a certifiably subgaussian distribution $D^{\prime}$ that has statistical distance at most $\varepsilon$ from the uniform distribution over $Y$,

output the low-degree moments of $D^{\prime}$.

We can typically find a distribution $D^{\prime}$ as above, because with high-probability $D$ satisfies the conditions. What remains to prove is that no matter what distribution $D^{\prime}$ 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^{\prime}$, whose low-degree moments are (certifiably) subgaussian, have small statistical distance, then their low-degree moments are close.

Our first bound of this kind, controls the distance between the means of $D$ and $D^{\prime}$ in terms of their covariances $\Sigma$ and $\Sigma^{\prime}$. (Later bounds will also relate $\Sigma$ and $\Sigma^{\prime}$.)

Since we assume that $\varepsilon<0.9$, we can conclude that $\lvert\langle u,\mu-\mu^{\prime}\rangle\rvert\leqslant\delta\cdot\langle u,(\Sigma+\Sigma^{\prime})u\rangle^{1/2}$ for $\delta\leqslant O\left(\varepsilon^{1-1/k}\cdot\sqrt{Ck\,}\right)$ as desired. ∎

Next we show that if two certifiably subgaussian distributions have small statistical distance, then their (raw) second moments are close in the Löwner order sense.

In particular, \bigl{(}1-O(\delta)\bigr{)}M^{\prime}\preceq M\preceq\bigl{(}1+O(\delta)\bigr{)}M^{\prime} if $\delta<0.9$.

which implies the desired bound $-\delta\cdot(M+M^{\prime})\preceq M-M^{\prime}\preceq\delta\cdot(M+M^{\prime})$ for $\delta\leqslant O\left(Ck\cdot\varepsilon^{1-1/k}\right)$. For the third inequality, we use Lemma 5.2 (moment bounds for shifts of subgaussian distributions).

For $\delta\leqslant 1$, we can now rearrange this bound to get the following relationship between $M$ and $M^{\prime}$,

which implies the desired bound \bigl{(}1-O(\delta)\bigr{)}M^{\prime}\preceq M\preceq\bigl{(}1+O(\delta)\bigr{)}M^{\prime} for $\delta<0.9$. ∎

Next, we combine the previous bounds for the first and second moments in order to establish the identifiability of the covariance matrix.

Under the same conditions as the previous Lemma 2.2, the covariance matrices $\Sigma$ and $\Sigma^{\prime}$ of the distributions $D$ and $D^{\prime}$ satisfy,

As before, in the case that $\delta\leqslant 0.9$, the above bounds imply a strong multiplicative approximation of the covariance so that $(1-\delta^{\prime})\cdot\Sigma^{\prime}\preceq\Sigma\preceq(1+\delta^{\prime})\Sigma^{\prime}$ for some $\delta^{\prime}\leqslant O(\delta)$.

Here, $\delta_{1}\leqslant O(Ck\cdot\varepsilon^{1-1/k})$. By Lemma 2.1 (mean identifiability), $\langle u,(\mu-\mu^{\prime})(\mu-\mu^{\prime}){}^{\mkern-1.5mu\mathsf{T}}u\rangle\leqslant\delta_{2}\cdot\langle u,(\Sigma+\Sigma^{\prime})u\rangle$, where $\delta_{2}\leqslant O(Ck\cdot\varepsilon^{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^{\prime}$).

Under the same conditions as the previous Lemma 2.2 and the additional condition $Ck\cdot\varepsilon^{1-2k}\leqslant\Omega(1)$, the means $\mu$ and $\mu^{\prime}$ of the distributions $D$ and $D^{\prime}$ satisfy,

Next we bound the distances of higher order moments for certifiably subgaussian distributions that are close in statistical distance. Compared to first 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^{\prime}$ as two polynomials $p$ and $p^{\prime}$ of degree $2r$ and we bound their difference $p(u)-p^{\prime}(u)$ relative to the variance of the distribution in direction $u$.

which implies the desired bound for $\delta\leqslant\delta\leqslant O(Ck)^{r/2}\cdot\varepsilon^{1-\frac{r}{k}}$. For the third inequality, we rely on Lemma 5.2 (moment bounds for shifts of subgaussian distributions). ∎

Preliminaries

In this section, we define 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.

This fact, together with the equivalence of weak separation and optimization [GLS81] allows us to efficiently optimize over pseudo-distributions (approximately)—this algorithm is referred to as the sum-of-squares algorithm.

We remark that if $D$ is an actual (discrete) probability distribution, then we have $D\sdtstile{}{}\mathcal{A}$ if and only if $D$ is supported on solutions to the constraints $\mathcal{A}$.

We say that a system $\mathcal{A}$ of polynomial constraints is explicitly bounded if it contains a constraint of the form $\{\|x\|^{2}\leqslant M\}$. The following fact is a consequence of Fact 3.1 and [GLS81],

2 Sum-of-squares proofs

Let $f_{1},f_{2},\ldots,f_{r}$ and $g$ be multivariate polynomials in $x$. A sum-of-squares proof that the constraints $\{f_{1}\geqslant 0,\ldots,f_{m}\geqslant 0\}$ imply the constraint $\{g\geqslant 0\}$ consists of polynomials $(p_{S})_{S\subseteq[m]}$ such that

Low-degree sum-of-squares proofs are sound and complete if we take low-level pseudo-distributions as models.

Concretely, sum-of-squares proofs allow us to deduce properties of pseudo-distributions that satisfy some constraints.

If the pseudo-distribution $D$ satisfies $\mathcal{A}$ only approximately, soundness continues to hold if we require an upper bound on the bit-complexity of the sum-of-squares $\mathcal{A}\sststile{r^{\prime}}{}B$ (number of bits required to write down the proof).

The following fact shows that every property of low-level pseudo-distributions can be derived by low-degree sum-of-squares proofs.

Suppose $d\geqslant r^{\prime}\geqslant r$ and $\mathcal{A}$ is a collection of polynomial constraints with degree at most $r$, and $\mathcal{A}\vdash\{\sum_{i=1}^{n}x_{i}^{2}\leqslant B\}$ for some finite $B$.

Let $\{g\geqslant 0\}$ be a polynomial constraint. If every degree-$d$ pseudo-distribution that satisfies $D\sdtstile{r}{}\mathcal{A}$ also satisfies $D\sdtstile{r^{\prime}}{}\{g\geqslant 0\}$, then for every $\varepsilon>0$, there is a sum-of-squares proof $\mathcal{A}\sststile{d}{}\{g\geqslant-\varepsilon\}$.

Moment estimation algorithm

Our robust moment estimation algorithm is a direct translation of the robust identifiability results from Section Section 2.

The following fact shows that the solutions to $\mathcal{A}_{Y,\varepsilon}$ correspond to all $\varepsilon$-corruptions of the (already corrupted) sample $Y$. In particular, if $Y=\{y_{1},\ldots,y_{n}\}$ is an $\varepsilon$-corruption of a sample $X=\{x_{1},\ldots,x_{n}\}$ of $D$, then one solution to $\mathcal{A}$ consists of the uncorrupted sample $X$ so that $x^{\prime}_{1}=x_{1},\ldots,x^{\prime}_{n}=x_{n}$.

An assignment to the variables $x^{\prime}_{1},\ldots,x^{\prime}_{n}$ can be extended to a solution for $\mathcal{A}_{Y,\varepsilon}$ if and only if $x^{\prime}_{i}=y_{i}$ holds for all but at most an $\varepsilon$-fraction of the indices $i\in[n]$.

Note that in order to extend such an assignment, we can choose $w_{i}=1$ for $(1-\varepsilon)n$ indices $i\in[n]$ such that $x_{i}=y_{i}$ and $w_{i}=0$ for the remaining indices.

The following fact shows that the solutions to the system $\mathcal{B}_{C,k}$ correspond to certifiably subgaussian distributions.

The following algorithm is the key ingredient of our efficient outlier-robust estimators.

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

where $\delta=O\left(\sqrt{Ck\,}\cdot\varepsilon^{1-1/k}\right)$. Furthermore, the bit complexity of the sum-of-squares proof is polynomial.

Choose $r_{1},\ldots,r_{n}\in\{0,1\}$ such that $\sum_{i}r_{i}=(1-\varepsilon)n$ and $r_{i}(x_{i}-y_{i})=0$ for all $i\in[n]$. (Such a choice exists because $Y$ is an $\varepsilon$-corruption of $X$.)

Letting $z_{i}=r_{i}\cdot w_{i}$, we claim that

Indeed, since $\{w_{i}^{2}=w_{i}\}\sststile{2}{w_{i}}\{(1-z_{i})\leqslant(1-w_{i})+(1-r_{i})\}$, we have $\mathcal{A}_{Y,\varepsilon}\sststile{2}{}\{\sum_{i=1}^{n}(1-z_{i})\leqslant 2\varepsilon n\}$

Since we assume that $\varepsilon<0.01$, we can conclude that

for $\delta_{1}\leqslant O\left(\varepsilon^{1-1/k}\cdot\sqrt{Ck\,}\right)$ as desired.

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

Furthermore, if $\delta\leqslant 0.0001$,

Furthermore, the bit complexity of the sum-of-squares proof is polynomial.

By an argument analogous to the one in the proof of Lemma 4.4, we can obtain:

for $\delta\leqslant O\left(\varepsilon^{1-2/k}\cdot Ck\right)$. Lemma A.4 together with (4.8) implies for some $\delta^{\prime}\leqslant 100\delta$,

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

Furthermore, the bit complexity of the sum-of-squares proof is polynomial.

Here, $\delta_{2}\leqslant O(Ck\cdot\varepsilon^{1-2/k})$. Further,

By Lemma 2.1 (mean estimation), $\langle u,(\mu-\mu^{\prime})(\mu-\mu^{\prime}){}^{\mkern-1.5mu\mathsf{T}}u\rangle^{k/2}\leqslant\delta_{1}^{2}\cdot\langle u,(\Sigma+\Sigma^{\prime})u\rangle^{k/2}$, where $\delta_{1}\leqslant O(Ck\cdot\varepsilon^{1-1/k})$. Combining the above inequalities, we obtain:

Using this in conjunction with (4.13), we obtain:

Now, using Lemma A.3 with $f=\frac{\langle u,(\Sigma-\Sigma^{\prime})u\rangle}{O(\delta_{2})^{k/2}\langle u,\Sigma u\rangle}$ completes the proof. ∎

The following is the analog of the Corollary 2.4.

Combining Lemma 4.4 and Corollary 4.6 we have for $\delta_{1}=O(\sqrt{Ck}\varepsilon^{1-1/k})$ and $\delta_{2}=O(Ck)\varepsilon^{1-2/k}$,

Observe that $\langle u,\Sigma u\rangle^{k/2}$ is a constant, i.e., does not depend on any variable in the SoS proof. Using Lemma A.3 with $f=\frac{\langle u,\mu-\mu^{\prime}(X^{\prime})\rangle}{(6\delta_{1})^{2}\langle u,\Sigma u\rangle^{1/2}}$ now completes the proof. ∎

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.

where $\delta=O\left(C^{r}k^{r}\right)\cdot\varepsilon^{1-r/k}$. Consequently,

Furthermore, the bit complexity of the sum-of-squares proof is polynomial.

We point out a subtle technical difference 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 $\leqslant 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.

We can now complete the proof of Theorem 1.2.

By Cauchy-Schwarz inequality for pseudo-distributions, we have:

Using $u\rightarrow\Sigma^{-1/2}u$ thus immediately yields:

Taking square roots thus implies $\langle u,\Sigma^{-1/2}(\mu-\hat{\mu})\rangle\leqslant\delta\|u\|$ or equivalently, $\|\Sigma^{-1/2}(\mu-\hat{\mu})\|\leqslant\delta$.

On the other hand, taking square roots in (4.17) gives: $|\langle u,\mu-\hat{\mu}\rangle|=\|\mu-\hat{\mu}\|\leqslant\delta\|\Sigma\|^{1/2}$ 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.

Certifiably subgaussian distributions

In this section, we give sufficient conditions for distributions to be certifably subgaussian. As concrete consequences, we derive that (affine transformations of) products of scalar-valued subgaussian distributions are certifiably subgaussian and that certain mixture models are certifiably subgaussian.

For the benefit of the reader, we restate here the definition of certifiable subgaussianity.

We first note some basic invariance properties of certifiable subgaussianity.

Next, we construct examples of various classes of distributions that are certifiably subgaussian. We start by noting that the standard gaussian distribution is $k$-certifiably $O(1)$-subgaussian for every $k$.

Next, we prove a similar fact about uniform mixtures of arbitrary mean gaussians.

Let $\mu_{1},\mu_{2},\ldots\mu_{q}$ be arbitrary. Let $D$ be the mixture of $\mathcal{N}(\mu_{i},I)$ with mixture weights $\frac{1}{q}$ for each component. Then, $D$ is $(k,k)$-certifiably $C$-subgaussian with $C=O(q)$.

We have using Lemma 5.2 and certifiable $C$-subgaussianity of the gaussian distribution for $C=O(1)$,

We will use the following matrix concentration inequality in the proof.

Now, consider the affine transformation $x\rightarrow M_{2}^{-1/2}x$ which makes second moment is now $I$. By Lemma 5.1, it is enough to show certifiable subgaussianity of the distribution after this transformation. Thus, without loss of generality, we assume that the $D$ is isotropic, that is, $M_{2}=I.$

Now, by sum-of-squares version of the Cauchy-Schwarz inequality

Let $D$ be a mixture of $D_{1}$ and $D_{2}$ with mixture weights $1-\lambda$ and $\lambda$ respectively such that $\lambda<\Delta^{-k/2}$.

Assume that $D_{1}$ and $D_{2}$ are mean zero. We have for any $u$,

The proof of this lemma follows immediately from the following proposition.

Call a multi-set $S$ of $[d]$ even if every individual element in $S$ appears even number of times.

Applications

In this section, we apply our robust moment estimation algorithm to give an efficient algorithm for outlier-robust independent component analysis. We recall the definition of subgaussian random variables first.

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

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

Our algorithm is simple. It first estimates the 2nd and 4th moments of the observed distribution from an $\varepsilon$-corrupted sample using Algorithm 4.3 and then applies a blackbox tensor decomposition algorithm to a “whitened” version of 4th moments.

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 define the following relaxation of the injective tensor norm in describing the guarantees of this algorithm.

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

Let $T=\hat{M}_{4}^{\prime}-3I\otimes I$ for the whitenened 4th moments $\hat{M}_{4}^{\prime}$. Let $a_{1},a_{2},\ldots,a_{d}$ be the columns of $A$ and let $b_{i}=(AA^{\top})^{-1/2}a_{i}$ be the whitened versions. Then,

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.

The whitened and shifted 4th moment $T$ then, by definition, satisfies

We will use $\Sigma,M_{4},M_{4}^{\prime},T,b_{i}=\Sigma^{-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$, $\langle u,\hat{\Sigma}u\rangle=(1\pm\delta_{2})\langle u,\Sigma u\rangle$ for some $\delta_{2}<O(Ck)\varepsilon^{1-2/k}$. Thus:

By Theorem 1.3, for $\delta_{4}=O(C^{2}k^{2})\varepsilon^{1-4/k}$,

Since $M_{4}^{\prime}-\hat{M}_{4}^{\prime}=T-\hat{T}$, this gives:

The proof of Theorem 1.4 then follows easily.

We are now ready to apply Fact 6.6 to obtain that decomposing the tensor $T$ recovers components $c_{1},c_{2},\ldots,c_{d}$ such that there’s a permutation $\pi:[d]\rightarrow[d]$ such that $\|c_{i}-b_{\pi(i)}\|_{2}^{2}\leqslant O(1)\cdot\frac{1}{|\gamma|}\cdot\delta_{4}(2+3\delta_{2})$.

Applying the reverse whitening transform $\hat{a}_{i}=\hat{\Sigma}^{1/2}c_{i}$ then recovers the set of vectors $a_{1},a_{2},\ldots,a_{d}$ satisfying the requirements of Theorem 1.4. ∎

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.

Let $D$ be mixtures of $\mathcal{N}(\mu_{i},I)$ for $i\leqslant q$ with uniformWhile our algorithm generalizes naturally to arbitrary mixture weights, we restrict to this situation for simplicity mixture weights. Assume that $\mu_{i}$s are linearly independent and, further, let $\kappa$ be the smallest non-zero eigenvalue of $\frac{1}{q}\sum_{i}\mu_{i}\mu_{i}^{\top}$.

Given an $\varepsilon$-corrupted sample of size $n\geqslant n_{0}=\Omega((\kappa d\log{(d)})^{k/2})$, for every $k\geqslant 4$, there’s a $\operatorname{poly}(n)d^{O(k)}$ time algorithm that recovers $\hat{\mu}_{1},\hat{\mu}_{2},\ldots\hat{\mu}_{q}$ so that there’s a permutation $\pi:[q]\rightarrow[q]$ satisfying

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

The following is the main technical lemma required in the proof of 6.9.

Let $Y$ be an $\varepsilon$-corrupted sample of size $n>n_{0}=\Omega((\kappa d\log d)^{k}/\varepsilon^{2})$ from $D_{0}=\sum_{i}\frac{1}{q}\mathcal{N}(\mu_{i},I)$. Let $\hat{M}_{i}$ be the estimated ith moment of $D_{0}$ given by Algorithm 4.3 when run on $Y$. Let $\hat{M}_{i}^{\prime}$ be the whitening of the estimated $i$th moment of $D_{0}$ as in Algorithm 6.11.

Then, with high probability over the draw of the $\varepsilon$-corrupted sample $Y$, for $\hat{T}=\hat{M}^{\prime}_{3}-3\hat{M}_{1}^{\prime}\otimes I$ and $\nu_{i}=(M_{2}-I)^{-1/2}\mu_{i}$ for $i\leqslant q$ and $\delta_{i}=O(qk)^{i/2}\varepsilon^{1-i/k}$ for $i=1,2,3$,

As before, we will first analyze the tensor $T$ assuming all estimates of the moments are exactly correct and then account for the estimation errors. We will write $M_{1},M_{2},M_{3}$ for the true first three moments of the input mixture of gaussian and $\Sigma=M_{2}-I$. We will let $\hat{M}_{1},\hat{M}_{2},\hat{M}_{3},\hat{\Sigma}$ be the corresponding quantities estimated by our algorithm. We will write $M_{1}^{\prime},M_{3}^{\prime}$ for the true whitened moments and $\hat{M}_{1}^{\prime},\hat{M}_{3}^{\prime}$, the estimated counterparts.

By direct computtion, $\langle M_{3},u^{\otimes 3}\rangle=\frac{1}{q}\sum_{i}(\langle\mu_{i},u\rangle^{3}+3\langle\mu_{i},u\rangle)$. Thus, $\langle M_{3}-3M_{1}\otimes I,u^{\otimes 3}\rangle=\frac{1}{q}\sum_{i}\langle\mu_{i},u\rangle^{3}$.

Now, by definition, $\langle M_{3}^{\prime}-3M_{1}^{\prime}\otimes I,u^{\otimes 3}\rangle=\langle M_{3}-3M_{1}\otimes I,(\Sigma^{-1/2}u)^{\otimes 3}\rangle$.

From Theorem 1.2, for $\delta_{1}=O(\sqrt{qk})\varepsilon^{1-1/k},\delta_{2}=O(qk)\varepsilon^{1-2/k}$, $\langle M_{1}-\hat{M}_{1},\Sigma^{-1/2}u\rangle\leqslant\delta_{1}\langle u,M_{2}\Sigma^{-1}u\rangle^{1/2}$.

Using that $M_{2}=I+\frac{1}{q}\sum_{i}\mu_{i}\mu_{i}^{\top}$, $M_{2}\Sigma^{-1}=(\frac{1}{q}\sum_{i}\mu_{i}\mu_{i}^{\top})^{-1}+I$. Thus, the first term in (6.6) can be upper bounded by $\delta_{1}\|u\|_{2}+\delta_{1}\langle u,\Sigma^{-1},u\rangle^{1/2}\leqslant\delta_{1}(1+1/\sqrt{\kappa})\|u\|_{2}$.

From Theorem 1.2, for any $u$, $\langle u,(M_{2}-\hat{M}_{2})u\rangle\leqslant\delta_{2}\langle u,M_{2}u\rangle$. Writing $M_{2}=\Sigma+I$ for $\Sigma=\frac{1}{q}\mu_{i}\mu_{i}^{\top}$, we obtain that $\langle u,\hat{\Sigma}u\rangle=(1\pm O(\delta_{2}/\kappa))\langle u,\Sigma u\rangle$.

Thus, $\langle(\Sigma^{-1/2}-\hat{\Sigma}^{-1/2})\hat{M}_{1},u\rangle\leqslant O(\delta_{2}^{1/2}/\sqrt{\kappa})\cdot\langle\Sigma^{-1/2}\hat{M}_{1},u\rangle$. Writing $\langle\Sigma^{-1/2}\hat{M}_{1},u\rangle=\langle\Sigma^{-1/2}(\hat{M}_{1}-M_{1})u\rangle+\langle\Sigma^{-1/2}M_{1}u\rangle$ and using that $\frac{1}{q}\sum_{i}\langle\Sigma^{-1/2}\mu_{i},u\rangle\leqslant\frac{1}{\sqrt{q}}\|u\|_{2}$, we finally obtain anupper bound on $\eqref{eq:err-mean}$ of: $\delta_{1}(1+1/\sqrt{\kappa})+O(\delta_{2}^{1/2}/\sqrt{\kappa})\delta_{1}\left(1+1/\sqrt{\kappa}+\frac{1}{\sqrt{q}}\right)\|u\|_{2}\leqslant O(\delta_{2}^{1/2}\delta_{1}/\kappa+\delta_{1}/\sqrt{\kappa})\|u\|_{2}$.

For the first term of (6.5), we apply Theorem 1.3 to get that for $\delta_{3}=O(qk)^{3/2}\varepsilon^{1-3/k}$,

As a result, $\sststile{6}{u}\langle\hat{M_{3}}^{\prime}-M_{3}^{\prime},u^{\otimes 3}\rangle^{2}\leqslant 2\delta_{3}^{2}\|u\|_{2}^{3}$ giving us an upper bound on the first term in the RHS of (6.5).

We can now complete the proof of Theorem 6.9.

Applying Fact 6.10 to $\hat{T}$ as in Lemma 6.12 yields that the resulting components satisfy $\min_{\pi:[q]\rightarrow[q]}\max_{i}\|\nu_{i}-\hat{\nu}_{\pi(i)}\|_{2}\leqslant O(1)\delta_{1}/\sqrt{\kappa}+\delta_{1}\delta_{2}^{1/2}/\kappa+2\delta_{3}^{1/3}$. The theorem now follows. ∎

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 (certifiably) 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)$-certifiable $C$-subgaussianity and 2) are at most $\varepsilon$ apart in total variation distance but have low-degree 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 certifiability 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.

$D_{1}$, $D_{2}$ are both $2$-subgaussian in all moments up to $k$.

$d_{TV}(D_{1},D_{2})\leqslant\varepsilon$.

Let $D_{1}$ be the standard gaussian distribution $\mathcal{N}(0,1)$. Let $D_{2}$ be the distribution that with probability $1-\varepsilon$ outputs a sample from $\mathcal{N}(0,1)$ and with probability $\varepsilon$, outputs $\sqrt{k}\varepsilon^{-1/k}$. Observe that by construction, $d_{TV}(D_{1},D_{2})\leqslant\varepsilon$.

The mean of $D_{2}$ is easily seen to be $\sqrt{k}\varepsilon^{1-1/k}$. The variance of $D_{2}$ can be computed to be between $1$ and $(1+k\varepsilon^{1-2/k})$.

$D_{1}$, $D_{2}$ are both $2$-subgaussian in all moments up to $k$.

$d_{TV}(D_{1},D_{2})\leqslant\varepsilon$.

Let $D_{1}$ be the standard gaussian distribution $\mathcal{N}(0,1)$. Let $D_{2}$ be the distribution that with probability $1-\varepsilon$ outputs a sample from $\mathcal{N}(0,1)$ and with probability $\varepsilon$, outputs a uniform element of $\pm\sqrt{k}\varepsilon^{-1/k}$. Observe that by construction, $d_{TV}(D_{1},D_{2})\leqslant\varepsilon$.

The mean of $D_{1},D_{2}$ is easily seen to be . The same computation as in the proof of Lemma 7.1 establishes that $D_{1}$, $D_{2}$ are both $2$-subgaussian.

The variance of $D_{1}$ is $1$ and for $D_{2}$ can be computed to be at least $(1-\varepsilon)(1+k\varepsilon^{1-2/k}-k\varepsilon^{2-2/k})\geqslant 1+\frac{k}{2}\varepsilon^{1-2/k}$.

The $2r$th moment of $D_{1}$ is $\sim k^{r}$ and for $D_{2}$ can be computed to as: $k^{r}(1-\varepsilon+\varepsilon^{1-2r/k})\geqslant k^{r}(1+0.5\varepsilon^{1-2r/k})$ for $\varepsilon<2^{-k/2r}$.

References

Appendix A Sum-of-squares toolkit

We will rely on the following sum-of-squares version of the inequality between the arithmetic and geometric mean.

By iterating this inequality $k/2$ times, we obtain

Here, we used for the second equality that $\binom{k}{i}=\binom{k}{k-i}$.

Combining the previous two inequalities yields the desired inequality,

We apply Lemma A.1 for $k/2$ and substitute $w_{1}=f^{2}$ and $w_{2},w_{3},\ldots,w_{k/2}=1$ to obtain $\sststile{k}{f}f^{2}\leqslant\frac{f^{k}+(k/2-1)}{k/2}$. Thus,

At the same time, since $2-2f=1-f^{2}+(1-f)^{2}$, we have

Combining (A.7) and (A.8) yields the desired conclusion. ∎

Therefore, for $\delta^{\prime}=\left(\frac{4\delta}{1-2\delta}\right)^{k}$,

By Lemma A.3, $\left\{(f-1)^{k}\leqslant(\delta^{\prime})^{k}\right\}\sststile{}{}\left\{1-\delta^{\prime}\leqslant f\leqslant 1+\delta\right\}$, which together with the previous inequality gives the desired inequality. ∎

Let $\mathcal{A}$ be a set of polynomial equality axioms in variable $x$ such that:

for a polynomial $p$ with total degree at most $2t$.

Let $\mathcal{A}=\{q_{1}=0,q_{2}=0,\ldots,q_{r}=0\}$. Since $\mathcal{A}\sststile{2t}{}p(x,u)\geqslant 0$, there are degree $\leqslant 2t$ (in $x$) polynomials $g_{1}(x,u),g_{2}(x,u),\ldots g_{r^{\prime}}(x,u)$ and polynomials $h_{1}(x,u),\ldots,h_{r}(x,u)$ such that $p(x,u)=\sum_{i}g_{i}(x,u)^{2}+\sum_{j}h_{j}(x,u)q_{j}(x)$ where the degree of any $h_{j}(x,u)\cdot q_{j}(x)$ in $x$ is at most $2t$.

Further, each $g_{i}(x,u)^{2}$ can be written as $\langle G(x)G(x)^{\top},(1,u)^{\otimes 2t}\rangle$ for some matrix valued polynomial $G(x)$.