Robustly Learning General Mixtures of Gaussians

1.1 Key Challenges

It is believed that the sum-of-squares hierarchy might actually be the key to solving the more general problem of robustly learning a mixture of any constant number of Gaussians. However there are some obstacles that need to be overcome:

Robust Identifiability. Behind every polynomial time algorithm for learning the parameters of a mixture model is an argument for why there cannot be two mixtures with noticeably different parameters that produce almost the same distribution. In fact, we need quantitative bounds that say any two mixtures that are \(\epsilon\)-different must be at least \(\mbox{poly}(\epsilon , 1/d)\) different according to some natural family of test functions, usually the set of low degree polynomials. Here d is the dimension. This is called polynomial identifiability [46, 51]. Because we allow a polynomial dependence on \(1/d\), it often does not matter too much how we measure the differences between two mixtures, either in terms of some natural parameter distance between their components or in terms of the total variation distance, again between their components.

However, we need much stronger bounds when it comes to robust learning problems where we want to be able to tolerate a constant fraction of corruptions that is dimension independent. In particular, we need a family of test functions with the stronger quantitative property that whenever we have two mixtures whose components are \(\epsilon\)-different in total variation distance there is some function in the family that has at least \(\mbox{poly}(\epsilon)\) discrepancy when we take the difference between its expectations with respect to the two distributions (and its variance must also be bounded). In particular, this relationship cannot involve dependence on the dimension d. We will call this robust identifiability. Recall that the non-robust learning algorithms for mixtures of Gaussians reduce to a series of one-dimensional problems. Unfortunately this strategy inherently introduces polynomial factors in d and it cannot give what we are after. For the special case of clusterable mixtures of Gaussians, Bakshi and Kothari [6] and Diakonikolas et al. [21] proved robust identifiability and their approach was based on classifying the ways in which two single Gaussians can have total variation distance close to one. When it comes to the more general problem of handling mixtures where the components can overlap non-trivially, it seems difficult to follow the same route because we can no longer match components from the two mixtures to each other and almost cancel them both out. To prove robust identifiability for general mixtures of Gaussians, we use a carefully chosen set of low-degree polynomials as test functions, namely the (a multivariate analog) of the Hermite polynomials. Our proof of robust identifiability is explained in more detail in the following subsection.

Reasoning About the Sum-of-Squares Relaxation. A proof of robust identifiability alone does not give an efficient algorithm. However, we leverage the now standard sum-of-squares framework [9] to convert our proof of robust identifiability into an actual learning algorithm. The sum-of-squares hierarchy is a general framework for coming up with large and powerful semidefinite programming relaxations that can be applied to many sorts of problems [10, 11, 12, 33]. The semidefinite programs automatically enforce a large family of constraints and proofs of statements that use only these types of sum-of-squares amenable constraints can be directly converted into algorithms. However, it is often quite challenging to understand whether or not it works and/or to identify, out of all the constraints that are enforced on the pseudo-expectation, which ones are actually useful in the analysis [30, 32]. What makes matters especially challenging in our setting is that it is clear the structure of the higher moments of a mixture of Gaussians must play a major role. But how exactly do we leverage them in our analysis?

1.2 Our Techniques and Main Result

Actually, we overcome both obstacles using the same approach. We store the relevant moments and variables we would like to solve for in certain generating functions. Then by manipulating the generating functions using differential operators, we are able to reason about an SOS relaxation of a natural polynomial system that allows us to solve for the parameters that we want to learn.

Let us describe the setup. For simplicity, assume that the mixture is in isotropic position. First, we have the unknown parameters of the mixture. Let \(\begin{equation*} \mathcal {M} = w_1 N(\mu _1, I + \Sigma _1) + \dots + w_k N(\mu _k, I + \Sigma _k) \end{equation*}\) where \(w_i\) are the mixing weights and \(\mu _i\) and \(I + \Sigma _i\) denote the mean and covariance of the ith component. Second, we have the indeterminates we would like to solve for. These will be denoted \(\widetilde{\mu _i}\) and \(\widetilde{\Sigma _i}\) and our intention is for these to be the means and covariances of a hypothetical mixture of Gaussians.¹ We will also guess the mixing weights of the hypothetical mixture \(\widetilde{w_i}\). Finally, we have a d-dimensional vector \(X = (X_1, \dots , X_d)\) of formal variables and one auxiliary formal variable y. These will mostly be used to help us organize everything in a convenient way. Roughly, we would like to solve for \(\widetilde{\mu _i}\) and \(\widetilde{\Sigma _i}\) so that the hypothetical mixture \(\begin{equation*} \widetilde{\mathcal {M}} = \widetilde{w_1} N(\widetilde{\mu _1}, I + \widetilde{\Sigma _1}) + \dots + \widetilde{w_k} N(\widetilde{\mu _k}, \widetilde{I + \Sigma _k}) \end{equation*}\) is close to \(\mathcal {M}\) on the family of test functions (which will be low-degree multivariate Hermite polynomials). It turns out that this amounts to solving a polynomial system for the indeterminates.

Now we explain in more detail how to actually reason about and solve the polynomial system. It will be useful to work with the following generating functions. First, let \(\begin{equation*} F(y) = \sum _{i = 1}^k w_i e^{\mu _i(X)y + \frac{1}{2}\Sigma _i(X)y^2} \end{equation*}\) Here we have used the notation that \(\mu _i(X)\) denotes the inner product of \(\mu _i\) with the d-dimensional vector X and that \(\Sigma _i(X)\) denotes the quadratic form of X on \(\Sigma _i\). Second, let \(\begin{equation*} \widetilde{F}(y) = \sum _{i = 1}^k \widetilde{w_i} e^{\widetilde{\mu _i}(X)y + \frac{1}{2}\widetilde{\Sigma _i}(X)y^2} \end{equation*}\) As is familiar from elementary combinatorics, we can tease out important properties of the generating function by applying carefully chosen operators that involve differentiation. This requires a lot more bookkeeping than usual because there are unknown parameters of the mixture, indeterminates, and formal variables. But it turns out that there are simple differential operators we can apply which can isolate components. To gain some intuition, consider the operator \(\begin{equation*} \mathcal {D}_i = \partial _y - (\mu _i(X) + \Sigma _i(X)y) \,. \end{equation*}\) Note that \(\begin{equation*} \mathcal {D}_i\left(e^{\mu _i(X)y + \frac{1}{2}\Sigma _i(X)y^2} \right) = 0 \,, \end{equation*}\) in other words, this operator annihilates the i^th component. Thus, by composing such operators, we can annihilate all but one of the components in F.² On the other hand, note that applying differential operators is just a rearrangement of the polynomials that show up in the infinite sum representation of the generating function but using differential operators and generating functions in exponential form gives a particularly convenient way to derive useful expressions that would otherwise be extremely complex to write down.

Ultimately we derive a symbolic identity (1) \(\begin{equation} \widetilde{w_k}\prod _{i=1}^k (\widetilde{\Sigma _k}(X) - \Sigma _i(X))^{2^{i-1}}\prod _{i=1}^{k-1}(\widetilde{\Sigma _k}(X) - \widetilde{\Sigma _i}(X))^{2^{k+i-1}} = \sum _{i = 1}^m P_i(X)(\widetilde{h_i}(X) - h_i(X)) \end{equation}\) where m is a function of k. In the above, the \(h_i(X)\)’s are the expectations of the multivariate Hermite polynomials for the true mixture \(\mathcal {M}\) and the \(\widetilde{h_i}\)’s are the expectations of the multivariate Hermite polynomials for the hypothetical mixture \(\widetilde{\mathcal {M}}\). A more detailed explanation of Hermite polynomials is given in Section 3, but for now we may simply think of them as modified moments. The reason that we use Hermite polynomials instead of standard moments is that they can be robustly estimated without losing dimension-dependent factors (see e.g., [39]).

The above identity (when combined with a few others of a similar form) allows us to deduce robust identifiability. Roughly, this is because if we have a mixture \(\widetilde{\mathcal {M}}\) with means \(\widetilde{\mu _i}\) and covariances \(I + \widetilde{\Sigma _i}\) that don’t match those of \(\mathcal {M}\), then the LHS of (1) will be bounded away from 0, implying that some term on the RHS must also be bounded away from 0. This means that there must be some \(i \le m = O_k(1)\) such that \(h_i(X)\) and \(\widetilde{h_i(X)}\) are substantially different — i.e., there is some test function that is a low-degree multivariate Hermite polynomial that distinguishes the two mixtures.

Robust identifiability alone does not give us a polynomial time learning algorithm. However, it turns out that we can use SOS to obtain a polynomial time learning algorithm from the argument for robust identifiability i.e., we essentially get the learning algorithm “for free”. The key point is that the \(h_i(X)\) can be estimated using our samples and the coefficients of the \(\widetilde{h_i}\) are explicit polynomials in the indeterminates that we can write down. Also the \(P_i\)’s are polynomials in everything: the unknown parameters, the indeterminates, and the formal variables except for y. To set up an SOS system, we obtain robust estimates \(\overline{h_i}\) for the expectations of the Hermite polynomials \(h_i(X)\) for the true mixture that we can compute from existing techniques in the literature. We then enforce that the expectations of the Hermite polynomials for the hypothetical mixture \(\widetilde{h_i}(X)\) are close to these robust estimates where closeness is defined in terms of the distance between their coefficient vectors.

It is not immediately clear why the expression in (1) ought to be useful for solving the SOS system that we set up. After all, we cannot explicitly compute it because it depends on things we do not have, like the true parameters, the \(\mu _i\)’s and \(\Sigma _i\)’s. However, the sum-of-squares relaxation enforces that the pseudo-expectation operator assigns values to polynomials in the indeterminates in a way that behaves like an actual distribution on solutions when we are evaluating certain types of low degree expressions that contain the one above. So even though we do not know the actual polynomials in the identity, they exist and the fact that they are enforced is enough to ensure that we can estimate the covariance of the kth component We stress that this is just the high-level picture and many more details are needed to fill it in.

Using these techniques, we come to the main result of our paper, which is a polynomial-time algorithm for robustly learning the parameters of a high-dimensional mixture of a constant number of Gaussians. Our main theorem is (informally) stated below. A formal statement can be found in Theorem 8.3.

1.3 Proof Overview

The proof of our main theorem can be broken down into several steps. We first present our main contribution, an algorithm for learning mixtures of Gaussians when no pair of components are too far apart. We introduce the necessary generating function machinery in Section 3 and then present our algorithm in Sections 4 and 5. Specifically, in Section 4 we show how to learn the parameters once we have estimates for the Hermite polynomials of the true mixture. And in Section 5, we show how to robustly estimate the Hermite polynomials, using similar techniques to [39].

To complete our full algorithm for learning general mixtures of Gaussians, we combine our aforementioned results with a clustering algorithm similar to [21]. Combining these algorithms, we prove that our algorithm outputs a mixture that is close to the true mixture in TV distance. This is done in Sections 6 and 7. We then prove identifiability in Section 8, implying that our algorithm actually learns the true parameters.

1.4 Concurrent and Subsequent Work

There are three main pieces of work that we discuss. The first by Bakshi et al. [4], was independent and concurrent to this one. The next is a subsequent work by Bakshi et al. [5] that improves their earlier result but also borrows techniques from our paper. We also discuss our follow-up work [45] that was after [4] but before [5] (and the contributions are essentially disjoint).

Bakshi et al. in [4] obtain a result that is similar to our main result Theorem 1.1, but using rather different techniques. There are a few key differences, which we discuss below. We learn the mixture to accuracy \(\epsilon ^{\Omega _k(1)}\) while their result only achieves accuracy \((1/ \log (1/\epsilon))^{\Omega _k(1)}\), an exponentially worse guarantee. Also, our result solves parameter learning — i.e., we estimate the parameters of the true mixture — while their algorithm solves proper density estimation — i.e., it outputs a mixture of k Gaussians that is close to the true density. Their algorithm does not need any lower bound on the mixing weights or on the pairwise separation of the components. In fact, lower bounds on these quantities are necessary for parameter learning. However, our original result makes an even stronger assumption about the bounded fractionality of the mixing weights and pairwise separation.

Subsequently in [5], Bakshi et al. improve their earlier result to achieve parameter learning and accuracy \(\epsilon ^{\Omega _k(1)}\). In fact, the analysis of their new parameter learning algorithm relies crucially on the robust identifiability result from our paper (see Section 9 in [5]). Thus, all known algorithms for robust parameter learning go through our machinery and robust identifiability. Compared to our result, the main improvement in [5] is an improved, polynomial dependence on the minimum mixing weight.

Our follow-up work [45] improves the results here in a different direction. We are able to obtain an algorithm that solves density estimation to accuracy \(\widetilde{O}(\epsilon)\) (instead of \(\epsilon ^{\Omega _k(1)}\)). Note that parameter learning to this accuracy is information-theoretically impossible. The work in [45] does need a stronger assumption that the components have variances in all directions lower and upper bounded by a constant and the learning algorithm is improper, outputting a mixture of polynomials times Gaussians, instead of just a mixture of Gaussians. The main relevance of [45] to this paper is that as a first step, [45] shows how to improve the separation assumption in Theorem 1.1 from some constant to \(\epsilon ^{\Omega _k(1)}\).

2.1 Problem Setup

We use \(N(\mu , \Sigma)\) to denote a Gaussian with mean \(\mu\) and covariance \(\Sigma\). We use \(d_{\textsf {TV}}(\mathcal {D}, \mathcal {D^{\prime }})\) to denote the total variation distance between two distributions \(\mathcal {D}, \mathcal {D^{\prime }}\). We begin by formally defining the problem that we will study. First we define the contamination model. This is a standard definition from robust learning (see e.g., [21]).

In this paper, we study the following problem. There is an unknown mixture of Gaussians \(\begin{equation*} \mathcal {M} = w_1G_1 + \dots + w_kG_k \end{equation*}\) where \(G_i = N(\mu _i, \Sigma _i)\). We receive an \(\epsilon\)-corrupted sample \(Y_1, \dots , Y_n\) from \(\mathcal {M}\) where \(n = {\rm poly}(d/\epsilon)\). The goal is to output a set of parameters \(\widetilde{w_1}, \dots , \widetilde{w_k}\) and \((\widetilde{\mu _1}, \widetilde{\Sigma _1}), \dots , (\widetilde{\mu _k}, \widetilde{\Sigma _k})\) that are \({\rm poly}(\epsilon)\) close to the true parameters in the sense that there exists a permutation \(\pi\) on \([k]\) such that for all i \(\begin{equation*} |w_i - \widetilde{w_{\pi (i)}}|, d_{\textsf {TV}}\left(N(\mu _i, \Sigma _i) , N(\widetilde{\mu _{\pi (i)}}, \widetilde{\Sigma _{\pi (i)}})\right) \le {\rm poly}(\epsilon). \end{equation*}\)

Throughout our paper, we will assume that all of the Gaussians that we consider have variance at least \({\rm poly}(\epsilon /d)\) and at most \({\rm poly}(d/\epsilon)\) in all directions i.e., they are not too flat. This implies that their covariance matrices are invertible so we may write expressions such as \(\Sigma _i^{-1}\). We will also make the following assumptions about the mixture:

for some positive constants \(A,b\). Note that a lower bound on the minimum mixing weight and a lower bound on the TV distance between components is necessary for parameter learning. Throughout this paper, we treat \(k,A,b\) as constants — i.e., A and b could be any function of k — and when we say polynomial, the exponent may depend on these parameters. We are primarily interested in dependence on \(\epsilon\) and d (the dimension of the space).

2.2 Sum of Squares Proofs

We will make repeated use of the Sum of Squares (SOS) proof system. We review a few basic facts here (see [9] for a more extensive treatment). Our exposition here closely mirrors [21].

The dual objects to SOS proofs are pseudoexpectations. We will repeatedly make use of pseudoexpectations later on.

Definition 2.4.

Let \(x_1, \dots , x_n\) be indeterminates. A degree-k pseudoexpectation \(\widetilde{\mathbb {E}}\) is a linear map

from degree-k polynomials to \(\mathbb {R}\) such that \(\widetilde{\mathbb {E}}[p(x)^2] \ge 0\) for any p of degree at most \(k/2\) and \(\widetilde{\mathbb {E}}[1] = 1\). For a set of polynomial constraints \(\mathcal {A} = \lbrace p_1(x) \ge 0, \dots , p_m(x) \ge 0 \rbrace\), we say that \(\widetilde{\mathbb {E}}\) satisfies \(\mathcal {A}\) if

for all polynomials \(s(x)\) and \(i \in [m]\) such that \(s(x)^2p_i(x)\) has degree at most k.

The key fact is that given a set of polynomial constraints, we can solve for a constant-degree pseudoexpectation that satisfies those constraints (or determine that none exist) in polynomial time as it reduces to solving a polynomially sized SDP.

Note that there are a few technical details with regards to only being able to compute a pseudoexpectation that nearly satisfies the constraints. These technicalities do not affect our proof (as \(2^{-n}\) errors will be negligible) so we will simply assume that we can compute a pseudoexpectation that exactly satisfies the constraints. See [9] for more details about these technicalities.

Finally, we state a few simple inequalities for pseudoexpectations that will be used repeatedly later on.

Claim 2.6

(Cauchy Schwarz for Pseudo-distributions).

Let \(f,g\) be polynomials of degree at most k in indeterminates \(x = (x_1, \dots , x_n)\). Then for any degree k pseudoexpectation,

Corollary 2.7.

Let \(f_1,g_1, \dots , f_m,g_m\) be polynomials of degree at most k in indeterminates \(x = (x_1, \dots , x_n)\). Then for any degree k pseudoexpectation,

Proof.

Note

where the first inequality follows from Cauchy Schwarz for pseudoexpectations and the second follows from standard Cauchy Schwarz.□

3.1 Polynomial Factorizations

The analysis of our SOS-based learning algorithm will rely on manipulations of Hermite polynomials. An important piece of our analysis is understanding how the coefficients of polynomials behave under addition and (polynomial) multiplication. Specifically, if we have two polynomials \(f(X), g(X)\) and we have bounds on the coefficients of f and g, we now want to give bounds on the coefficients of the polynomials \(f(X) + g(X)\) and \(f(X)g(X)\). Most of these bounds are easy to obtain. The one that is somewhat nontrivial is lower bounding the coefficients of \(f(X)g(X)\) i.e., if the coefficients of f and g are not all small, then the product \(f(X)g(X)\) cannot have all of its coefficients be too small.

First, we prove a simple result about the norm of the vectorization of a sum of polynomials.

Next, we upper bound the norm of the vectorization of a product of polynomials.

Claim 3.15.

Let \(f,g,h_1, \dots , h_k\) be polynomials in \(X_1, \dots , X_d\) of degree at most k with coefficients that are polynomials in formal variables \(u_1, \dots , u_n\) of degree \(O_k(1)\) Then for any pseudoexpectation \(\widetilde{\mathbb {E}}\) of degree \(C_k\) for some sufficiently large constant \(C_k\) depending only on k,

where the pseudoexpectation operates on polynomials in \(u_1, \dots , u_n\).

Before we prove the final result in this section, we introduce a few definitions.

Now we can prove a lower bound on the norm of the vectorization of the product of polynomials.

Claim 3.18.

Let \(f,g,h_1, \dots , h_k\) be polynomials in \(X_1, \dots , X_d\) of degree at most k with coefficients that are polynomials in formal variables \(u_1, \dots , u_n\) of degree \(O_k(1)\). Then for any pseudoexpectation \(\widetilde{\mathbb {E}}\) of degree \(C_k\) for some sufficiently large constant \(C_k\) depending only on k,

where the pseudoexpectation operates on polynomials in \(u_1, \dots , u_n\).

Proof.

We will first prove the statement for \(h_1 = \dots = h_k = 1\).

Let S be the set of all types that can be obtained by taking the sum of two degree vectors for monomials of degree at most k and let T be the set of all types that can be obtained by taking the difference of two degree vectors for monomials of degree at most k. Note that \(|S|,|T| = O_k(1)\). Now \(\begin{align*} \widetilde{\mathbb {E}}[\left\Vert v(fg)\right\Vert ^2] = \widetilde{\mathbb {E}}\left[\sum _{a} \left(\sum _{u+v = a}f_ug_v\right)^2\right] = \widetilde{\mathbb {E}}\left[ \sum _{u_1 + v_1 = u_2 +v_2} f_{u_1}g_{v_1}f_{u_2}g_{v_2}\right] \\ = \widetilde{\mathbb {E}}\left[ \sum _{u_1 - v_2 = u_2 - v_1} f_{u_1}g_{v_1}f_{u_2}g_{v_2}\right] = \widetilde{\mathbb {E}}\left[\sum _{b} \left(\sum _{u-v = b}f_ug_v \right)^2\right] \end{align*}\) where the sums in the above expression are over all a and all b that are vectors in \(\mathbb {Z}^d\) for which the inner summands are nonempty. Let \(T = \lbrace t_1, \dots , t_n \rbrace\) where the types \(t_1, \dots , t_n\) are sorted in non-increasing order of their \(L^2\) norm. Recall that T consists of all types that can be obtained by taking the difference of two degree vectors corresponding to monomials of degree at most k. Now first note

since \(t_1\) corresponds to the type \((k,-k)\) and each of the inner summands only contains one term. Now consider \(t_i\) for \(i \gt 1\). \(\begin{align*} \widetilde{\mathbb {E}}\left[\sum _{b,\tau (b) = t_i}\left(\sum _{u-v = b}f_ug_v \right)^2 \right] = \widetilde{\mathbb {E}}\left[\sum _{b,\tau (b) = t_i}\left(\sum _{u-v = b}(f_ug_v)^2 \right) \right] + 2\widetilde{\mathbb {E}}\left[\sum _{b,\tau (b) = t_i}\left(\sum _{\lbrace u_1,v_1\rbrace \ne \lbrace u_1,v_2 \rbrace \\ u_1 - v_1 = u_2 - v_2 = b}f_{u_1}f_{u_2}g_{v_1}g_{v_2}\right) \right] \end{align*}\) Note that in the second sum, either \(u_1 - v_2 \in t_j\) for \(j \lt i\) or \(u_2 - v_1 \in t_j\) for \(j \lt i\). To see this, let \(a = v_1 - v_2\). Then \(u_1 - v_2 = b + a\) and \(u_2 - v_1 = b - a\). Now \(\begin{equation*} \left\Vert b - a\right\Vert _2^2 + \left\Vert b + a\right\Vert _2^2 \gt \left\Vert b\right\Vert _2^2 \end{equation*}\) since \(a \ne 0\) so one of the differences must be of an earlier type.

Next, note that for a fixed \(u_1, v_2\), there are at most \(O_k(1)\) possible values for \(u_2^{\prime },v_1^{\prime }\) such that the term \(f_{u_1}f_{u_2^{\prime }}g_{v_1^{\prime }}g_{v_2}\) appears. This is because we must have \(u_1 + v_2 = u_2^{\prime } + v_1^{\prime }\) and there are only \(O_k(1)\) ways to achieve this. Thus, by Cauchy Schwarz

Now combining the above with the fact that \(|T| = n = O_k(1)\) and that

we can complete the proof. To see this, for each i, let

Also, normalize so that

Let \(\delta\) be some suitably chosen constant depending only on k. If \(Q_1 \ge \delta\) then we are done. Otherwise, we have an upper bound on the square root terms that are subtracted in the expression for \(R_2\). If \(Q_2 \ge \Omega _1(k)\sqrt {\delta }\) then we are again done (since we have now reduced to the case where \(Q_1 \le \delta\)). Iteratively repeating this procedure, we are done whenever one of the \(Q_i\) is sufficiently large compared to \(Q_1, \dots , Q_{i-1}\). However, not all of \(Q_1, \dots , Q_n\) can be small since their sum is 1. Choosing \(\delta\) to be a sufficiently small constant but depending only on k we conclude that

as desired.

For the general case when not all of the \(h_i\) are 1, we can multiply the insides of all of the pseudoexpectations above by \(\left\Vert v(h_1)\right\Vert ^2 \dots \left\Vert v(h_k)\right\Vert ^2\) and the same argument will work.□

4.1 Reducing to All Pairs of Parameters Equal or Separated

We claim that it suffices to work under the following assumption. All pairs of parameters are either separated of equal. More specifically, for each pair of parameters \(\mu _i, \mu _j\) (and same for \(\Sigma _i, \Sigma _j\)), either \(\mu _i = \mu _j\) or \(\begin{equation*} ||{\mu _i - \mu _j}||_2 \ge c \end{equation*}\)

We now prove that it suffices to work with the above simplification. For any function \(0 \lt f(k) \lt 1\) depending only on k, there is some \(C \ge (f(k))^{k^2}\) such that there is no pair of parameters \(\mu _i, \mu _j\) or \(\Sigma _i, \Sigma _j\) whose distance is in the interval \([ (\epsilon ^{\prime })^{C}, (\epsilon ^{\prime })^{ f(k)C} ]\). Now consider the graph on the k nodes where \(i,j\) are connected if and only if \(\begin{equation*} ||{\mu _i - \mu _j}|| \le (\epsilon ^{\prime })^{f(k)C} \end{equation*}\) We now construct a new mixture \(N(\mu _i^{\prime }, \Sigma _i^{\prime })\). For each connected component say \(\lbrace i_1, \dots , i_j \rbrace\) , pick a representative and set \(\mu _{i_1}^{\prime } = \mu _{i_2}^{\prime } = \dots = \mu _{i_j}^{\prime } = \mu _{i_1}\). Do this for all connected components and similar in the graph on covariance matrices. For all i, we have \(\begin{equation*} \left\Vert \mu _i^{\prime } - \mu _i\right\Vert , \left\Vert \Sigma _i^{\prime } - \Sigma _i\right\Vert \le O_k(1)(\epsilon ^{\prime })^C \end{equation*}\) because there is a path of length at most k connecting i to the representative in its component that it is rounded to, and all edges correspond to pairs within distance of \((\epsilon ^{\prime })^C\).

The Hermite polynomials of this new mixture satisfy \(\begin{equation*} \left\Vert v(h_m^{\prime } - h_m)\right\Vert ^2 \le O_k(1) \Delta ^{O_k(1)} (\epsilon ^{\prime })^{\Omega _k(1)C} \end{equation*}\) as long as m is bounded as a function of k by Corollary 3.7. If we pretend that the new mixture is the true mixture, we have estimates \(\overline{h}_m(X)\) such that \(\begin{equation*} \left\Vert v(h_m^{\prime } - \overline{h_m})\right\Vert ^2 \le O_k(1) \Delta ^{O_k(1)} (\epsilon ^{\prime })^{\Omega _k(1)C} \end{equation*}\) and all pairs of parameters in the new mixture are either equal or \((\epsilon ^{\prime })^{f(k)C}\) separated. If we prove Theorem 4.1 with the assumption that the pairs of parameters are separated or equal, then we can choose \(f(k)\) accordingly and then we deduce that the theorem holds in the general case (with worse, but still polynomial, bounds on \(\Delta ,c, w_{\min }\) and the accuracy of our output as a function of \(\epsilon ^{\prime }\)).

From now on we will work with the assumption that each pair of parameters is either equal or separated by c.

4.2 SOS Program Setup

Our algorithm for learning the parameters when given estimates of the Hermite polynomials involves solving an SOS program. Here we set up the SOS program that we will solve.

We will let \(D = \binom{d}{2} + d\). We think of mapping between symmetric \(d \times d\) matrices and \(\mathbb {R}^D\) as \(\begin{equation*} \begin{bmatrix}x_{11} & \dots & x_{1d} \\ \vdots & \ddots & \vdots \\ x_{d1} & \dots & x_{dd} \end{bmatrix} \leftrightarrow (x_{11}, 2x_{12} , 2x_{13}, \dots , x_{dd}) \end{equation*}\)

We will solve for an orthonormal basis for the span of the means \(\mu _i\) and the span of the (vectorized) covariances \(\Sigma _i\). We will have SOS variables for the elements of this basis. Note that regardless of the basis, since it is k-dimensional, we can brute-force search over the representations of the true means \(\mu _i\) and covariances \(\Sigma _i\) in terms of elements of the basis. To set up our SOS-program, we will first guess these coefficients (which are real numbers) and then set up the SOS-program to solve for the desired bases i.e., there is a distinct program for each guess of the coefficients.

Note that we can explicitly write down the last set of constraints because we have estimates \(\overline{h_i}\).

It is important to note that the \(\widetilde{w_i}, A_i, B_i\) are real numbers. We will attempt to solve the system for each of our guesses and show that for some set of guesses, we obtain a solution from which we can recover the parameters. We can brute-force search over an \(\epsilon ^{\prime }\)-net because there are only \(O_k(1)\) parameters to guess. We call the SOS program that we set up \(\mathcal {S}\).

4.3 Analysis

We now prove a set of properties that must be satisfied by any pseudoexpectation of degree \(C_k\) satisfying \(\mathcal {S}\) where \(C_k\) is a sufficiently large constant depending only on k. What we would ideally want to show is that

However, it appears to be difficult to prove a statement of the above form within an SOS framework. Instead, we will look at the pseudoexpectations of the matrices

Roughly Lemma 4.11 says that any singular vector that corresponds to a large singular value of \(M_i\) must be close to the span of the \(\lbrace \Sigma _i \rbrace\). Lemma 4.12 says that any vector v that has large projection onto the subspace spanned by the \(\lbrace \Sigma _i \rbrace\) must have the property that \(v^TM_iv\) is large for some i. Putting these together, we can take the the top-k principal components of each of \(M_1, \dots , M_k\) and show that the span of these essentially contains the span of the \(\lbrace \Sigma _i \rbrace\) (this last step is done outside the SOS framework). We can now brute-force over an \(\epsilon ^{\prime }\)-net and guess the \(\Sigma _i\) (since we have narrowed them down to an \(O_k(1)\)-dimensional subspace). We can then plug in real values for the covariances and solve for the means using a similar method.

4.3.2 Warm-up: All Pairs of Parameters are Separated.

As a warm-up, we first analyze the case where all pairs of true parameters \(\mu _i, \mu _j\) and \(\Sigma _i, \Sigma _j\) satisfy \(||{\mu _i - \mu _j}||_2 \ge c\) and \(||{\Sigma _i - \Sigma _j}||_2 \ge c\). We will show how to deal with the general case where parameters may be separated or equal in Section 4.3.4.

We can assume that our guesses satisfy \(||{A_i - A_j}||_2 \ge c/2\) and \(||{B_i - B_j}||_2 \ge c/2\) for all \(i,j\). The key expressions to consider are applying the following differential operators \(\begin{align*} \mathcal {D} = \widetilde{\mathcal {D}_{k}}^{2^{2k-1} - 1}\widetilde{\mathcal {D}_{k-1}}^{2^{2k-2}} \dots \widetilde{\mathcal {D}_{1}}^{2^k} \mathcal {D}_k^{2^{k-1}} \dots \mathcal {D}_1^{1} \\ \widetilde{\mathcal {D}} = \mathcal {D}_{k}^{2^{2k-1}- 1}\mathcal {D}_{k-1}^{2^{2k-2}} \dots \mathcal {D}_{1}^{2^k} \widetilde{\mathcal {D}_k}^{2^{k-1}} \dots \widetilde{\mathcal {D}_1}^{1} \end{align*}\) to F and \(\widetilde{F,}\) respectively. The reason these differential operators are so useful is that \(\mathcal {D}\) zeros out the generating function for the true mixture and also zeros out all but one component of the generating function for the hypothetical mixture with parameters \(\widetilde{w_i} , \widetilde{\mu _i}, I + \widetilde{\Sigma _i}\). For the one component that is not zeroed out, only the leading coefficient remains and we can use Claim 3.11 to explicitly compute the leading coefficient. Thus, we can compare the results of applying these operators on the generating functions for the true and hypothetical mixtures and, using the fact that the Hermite polynomials for these mixtures must be close, we obtain algebraic relations that allow us to extract information about individual components.

We begin by explicitly computing the relevant leading coefficients.

Claim 4.7.

Write \(\begin{equation*} \widetilde{\mathcal {D}_{k}}^{2^{2k-1} - 1}\widetilde{\mathcal {D}_{k-1}}^{2^{2k-2}} \dots \widetilde{\mathcal {D}_{1}}^{2^k} \mathcal {D}_k^{2^{k-1}} \dots \mathcal {D}_1^{1}(\widetilde{F}) \end{equation*}\) as a formal power series in y. Its evaluation at \(y = 0\) is \(\begin{equation*} C_k\widetilde{w_k}\prod _{i=1}^k (\widetilde{\Sigma _k}(X) - \Sigma _i(X))^{2^{i-1}}\prod _{i=1}^{k-1}(\widetilde{\Sigma _k}(X) - \widetilde{\Sigma _i}(X))^{2^{k+i-1}} \end{equation*}\) where \(C_k\) is a constant depending only on k.

Proof.

Write \(\begin{equation*} \widetilde{F}(y) = \sum _{i = 1}^k \widetilde{w_i} e^{\widetilde{\mu _i}(X)y + \frac{1}{2}\widetilde{\Sigma _i}(X)y^2} \end{equation*}\) When applying the differential operator, by Corollary 3.10, all of the terms become 0 except for \(\begin{equation*} \widetilde{w_k} e^{\widetilde{\mu _k}(X)y + \frac{1}{2}\widetilde{\Sigma _i}(X)y^2}. \end{equation*}\) We now use Claims 3.11 and 3.9 to analyze what happens when applying the differential operator to this term. We know that \(\begin{equation*} \widetilde{\mathcal {D}_{k-1}}^{2^{2k-2}} \dots \widetilde{\mathcal {D}_{1}}^{2^k} \mathcal {D}_k^{2^{k-1}} \dots \mathcal {D}_1^{1}(\widetilde{F}) = P(y)e^{\widetilde{\mu _k}(X)y + \frac{1}{2}\widetilde{\Sigma _i}(X)y^2} \end{equation*}\) where P has leading coefficient \(\begin{equation*} \widetilde{w_k}\prod _{i=1}^k (\widetilde{\Sigma _k}(X) - \Sigma _i(X))^{2^{i-1}}\prod _{i=1}^{k-1}(\widetilde{\Sigma _k}(X) - \widetilde{\Sigma _i}(X))^{2^{k+i-1}} \end{equation*}\) and degree \(2^{2k-1} - 1\). Thus, \(\begin{align*} &\widetilde{\mathcal {D}_{k}}^{2^{2k-1} - 1}\widetilde{\mathcal {D}_{k-1}}^{2^{2k-2}} \dots \widetilde{\mathcal {D}_{1}}^{2^k} \mathcal {D}_k^{2^{k-1}} \dots \mathcal {D}_1^{1}(\widetilde{F}) \\ &= (2^{2k-1} - 1)!\widetilde{w_k}\prod _{i=1}^k (\widetilde{\Sigma _k}(X) - \Sigma _i(X))^{2^{i-1}}\prod _{i=1}^{k-1}(\widetilde{\Sigma _k}(X) - \widetilde{\Sigma _i}(X))^{2^{k+i-1}}e^{\widetilde{\mu _k}(X)y + \frac{1}{2}\widetilde{\Sigma _i}(X)y^2} \end{align*}\) and plugging in \(y = 0\), we are done.□

Claim 4.8.

Write \(\begin{equation*} \mathcal {D}_{k}^{2^{2k-1}- 1}\mathcal {D}_{k-1}^{2^{2k-2}} \dots \mathcal {D}_{1}^{2^k} \widetilde{\mathcal {D}_k}^{2^{k-1}} \dots \widetilde{\mathcal {D}_1}^{1} (F) \end{equation*}\) as a formal power series in y. Its evaluation at \(y = 0\) is \(\begin{equation*} C_k w_k\prod _{i=1}^{k}(\Sigma _k(X) - \widetilde{\Sigma _i}(X))^{2^{i-1}} \prod _{i=1}^{k-1} (\Sigma _k(X) - \Sigma _i(X))^{2^{k + i-1}} \end{equation*}\) where \(C_k\) is a constant depending only on k.

Proof.

This can be proved using the same method as Claim 4.7.□

Combining the previous two claims with Claims 4.5 and 4.6, we can write the expressions for the leading coefficients as polynomial combinations of the Hermite polynomials.

Lemma 4.9.

Consider the polynomial \(\begin{equation*} \widetilde{w_k}\prod _{i=1}^k (\widetilde{\Sigma _k}(X) - \Sigma _i(X))^{2^{i-1}}\prod _{i=1}^{k-1}(\widetilde{\Sigma _k}(X) - \widetilde{\Sigma _i}(X))^{2^{k+i-1}} \end{equation*}\) It can be written in the form \(\begin{equation*} P_0(X) + P_1(X)\widetilde{h_1}(X) + \dots + P_{m}(X)\widetilde{h_{m}}(X) \end{equation*}\) where

•	m is a function of k
•	Each of the \(P_i\) is m-simple
•	We have \(\begin{equation*} P_0(X) + P_1(X)h_1(X) + \dots + P_{m}(X)h_{m}(X) = 0 \end{equation*}\) as an algebraic identity over formal variables \(X_1, \dots , X_d, \lbrace u_i \rbrace , \lbrace v_i \rbrace\).

Proof.

Consider the power series \(\begin{align*} \widetilde{\mathcal {D}_{k}}^{2^{2k-1} - 1}\widetilde{\mathcal {D}_{k-1}}^{2^{2k-2}} \dots \widetilde{\mathcal {D}_{1}}^{2^k} \mathcal {D}_k^{2^{k-1}} \dots \mathcal {D}_1^{1}(\widetilde{F}) \end{align*}\) Now using Claim 4.7 and repeating the proof of Claim 4.5, we get the desired identity.□

Similarly, we have:

Lemma 4.10.

Consider the polynomial \(\begin{equation*} w_k\prod _{i=1}^{k}(\Sigma _k(X) - \widetilde{\Sigma _i}(X))^{2^{i-1}} \prod _{i=1}^{k-1} (\Sigma _k(X) - \Sigma _i(X))^{2^{k + i-1}} \end{equation*}\) It can be written in the form \(\begin{equation*} P_0(X) + P_1(X)h_1(X) + \dots + P_{m}(X)h_{m}(X) \end{equation*}\) where

•	m is a function of k
•	Each of the \(P_i\) is m-simple
•	We have \(\begin{equation*} P_0(X) + P_1(X)\widetilde{h_1}(X) + \dots + P_{m}(X) \widetilde{h_{m}}(X) = 0 \end{equation*}\) as an algebraic identity over formal variables \(X_1, \dots , X_d, \lbrace u_i \rbrace , \lbrace v_i \rbrace\).

Everything we’ve done so far has been symbolic manipulations and the claims in this section are all true as algebraic identities. We are now ready to analyze the SOS program. Note the polynomials \(P_0, \dots , P_m\) in Lemma 4.9 are unknown because they depend on the true parameters. This is fine because we will simply use their existence to deduce properties of pseudoexpectations that solve the SOS-system \(\mathcal {S}\).

Let U be the subspace spanned by the true \(\mu _1, \dots , \mu _k\) and let V denote the subspace spanned by the true (flattened) \(\Sigma _1, \dots , \Sigma _k\). We will use \(\Gamma _{V}, \Gamma _{V^{\perp }}\) to denote projections onto V and the orthogonal complement of V (and similar for \(U, U^{\perp }\)). Note that these are linear maps.

Our goal now will be to show that V is essentially contained within the span of the union of the top k principal components of the matrices

This gives us a \(k^2\)-dimensional space that essentially contains V and then we can guess the true covariance matrices via brute force search. In the first key lemma, we prove that the matrix \(\widetilde{\mathbb {E}}[\widetilde{\Sigma _i} \widetilde{\Sigma _i}^T]\) lives almost entirely within the subspace V.

Lemma 4.11.

Let \(\widetilde{\mathbb {E}}\) be a pseudoexpectation of degree \(C_k\) for some sufficiently large constant \(C_k\) depending only on k that solves \(\mathcal {S}\). Consider the matrix

where by this we mean we construct the \(D \times D\) matrix \(\widetilde{\Sigma _k} \widetilde{\Sigma _k}^T\) whose entries are quadratic in the variables \(\lbrace u \rbrace , \lbrace v \rbrace\) and then take the entry-wise pseudoexpectation. Then \(\begin{equation*} \text{Tr}_{V^{\perp }}(M) \le (\epsilon ^{\prime })^{2^{-k}}O_k(1) \left(\frac{\Delta }{w_{\min }c} \right)^{O_k(1)} \end{equation*}\) where \(\text{Tr}_{V^{\perp }}(M)\) denotes the trace of M on the subspace \(V^{\perp }\).

Proof.

Using Lemma 4.9, we may write \(\begin{equation*} \widetilde{w_k}\prod _{i=1}^k (\widetilde{\Sigma _k}(X) - \Sigma _i(X))^{2^{i-1}}\prod _{i=1}^{k-1}(\widetilde{\Sigma _k}(X) - \widetilde{\Sigma _i}(X))^{2^{k+i-1}} = P_1(X)(\widetilde{h_1}(X) - h_1(X)) + \dots + P_{m}(X)(\widetilde{h_{m}}(X) - h_m(X)) \end{equation*}\) where \(m = O_k(1)\)

Now we bound

Using Claims 3.15 and 3.14,

Where the last step is true because Claim 3.14 allows us to write the difference between the two sides as a sum of squares.

Now \(||{v(\overline{h_i}(X) - h_i(X))}||^2\) is just a real number and is bounded above by \(\epsilon ^{\prime }\) by assumption. We also have the constraint that \(\begin{equation*} \left\Vert v(\widetilde{h_i}(X) - \overline{h_i}(X))\right\Vert ^2 \le 100\epsilon ^{\prime } \end{equation*}\) so

Now we use the properties from Lemma 4.9 that each of the \(P_i\) can be written as a linear combination of a constant number of terms that are a product of some of \(\lbrace \mu _i(X) \rbrace , \lbrace \Sigma _i(X) \rbrace , \lbrace \widetilde{\mu _i}(X) \rbrace , \lbrace \widetilde{\Sigma _i}(X) \rbrace\) where

•	The coefficients in the linear combination are bounded by a constant depending only on k
•	The number of terms in the sum depends only on k
•	The number of terms in each product depends only on k.

Note for each \(\widetilde{\mu _i}(X)\), since we ensured that our guesses for the coefficients that go with the orthonormal basis \(u_1, \dots , u_k\) are at most \(\Delta\) and we have the constraints \(\left\Vert u_i\right\Vert _2^2 = 1, u_i \cdot u_j = 0\), we have \(\begin{equation*} \left\Vert v(\widetilde{\mu _i}(X))\right\Vert ^2 \dashv _{O(1)} O_k(1)\Delta ^2 \end{equation*}\) where \(\dashv _{O(1)}\) means the difference can be written as a sum of squares of polynomials of constant degree. We can make similar arguments for \(\widetilde{\Sigma _i}(X), \mu _i(X), \Sigma _i(X)\). Now using Claims 3.14 and 3.15 we can deduce

Overall, we have shown

Now we examine the expression

By Claim 3.18 (recall \(\widetilde{w_k}\) is a constant that we guess),

Note that \(\begin{equation*} \left\Vert v(\widetilde{\Sigma _k}(X) - \Sigma _i(X))\right\Vert ^2 \vdash _{O(1)} \left\Vert \Gamma _{V^{\perp }}(\widetilde{\Sigma _k})\right\Vert ^2 \end{equation*}\) (recall that \(\Gamma _{V^{\perp }}\) is a projection map with unknown but constant coefficients). Next, since we ensure that the coefficients \(B_i\) that we guess for the orthonormal basis satisfy \(||{B_i - B_j}||_2 \ge \frac{c}{2}\), we have \(\begin{equation*} \left\Vert v(\widetilde{\Sigma _k}(X) - \widetilde{\Sigma _i}(X))\right\Vert ^2 \vdash _{O(1)} \frac{c^2}{4} \end{equation*}\) where we use the constraints in \(\mathcal {S}\) that \(\left\Vert v_i\right\Vert _2^2 = 1, v_i \cdot v_j = 0\). Overall, we conclude

Note

because the inner expressions are equal symbolically. Thus,

Thus,

It remains to note that

and we are done.□

In the next key lemma, we prove that any vector that has nontrivial projection onto V must also have nontrivial projection onto \(\widetilde{\mathbb {E}}[\widetilde{\Sigma _i} \widetilde{\Sigma _i}^T]\) for some i.

Lemma 4.12.

Let \(\widetilde{\mathbb {E}}\) be a pseudoexpectation of degree \(C_k\) for some sufficiently large constant \(C_k\) depending only on k that solves \(\mathcal {S}\). Consider the matrix

where by this we mean we construct the \(D \times D\) matrix whose entries are quadratic in the variables \(\lbrace u \rbrace , \lbrace v \rbrace\) and then take the entry-wise pseudoexpectation. Then for any unit vector \(z \in \mathbb {R}^D\), \(\begin{equation*} z^TNz \ge \left(\frac{w_{\min }(z \cdot \Sigma _k)^{O_k(1)} - O_k(1)\epsilon ^{\prime } \Delta ^{O_k(1)}}{O_k(1)\Delta ^{O_k(1)} }\right)^2 \end{equation*}\) as long as \(\begin{equation*} w_{\min }(z \cdot \Sigma _k)^{O_k(1)} \gt O_k(1)\epsilon ^{\prime } \Delta ^{O_k(1)}. \end{equation*}\)

Proof.

Using Lemma 4.10, we may write \(\begin{equation*} w_k\prod _{i=1}^{k}(\Sigma _k(X) - \widetilde{\Sigma _i}(X))^{2^{i-1}} \prod _{i=1}^{k-1} (\Sigma _k(X) - \Sigma _i(X))^{2^{k + i-1}} = P_1(X)(\widetilde{h_1}(X) - h_1(X)) + \dots + P_{m}(X)(\widetilde{h_{m}}(X) - h_m(X)) \end{equation*}\) where \(m = O_k(1)\).

Using the same method as the proof in Lemma 4.11, we have

Now by Claim 3.18,

where the second inequality is true because \(\begin{equation*} \left\Vert v\left(\Sigma _k(X) - \widetilde{\Sigma _i}(X)\right)\right\Vert ^2 \vdash _{O(1)} (z \cdot \Sigma _k - z \cdot \widetilde{\Sigma _i})^2 \end{equation*}\) Now we claim

To see this, first recall that \(z \cdot \Sigma _k\) is just a constant. Next, we can expand the LHS into a sum of monomials in the \(z \cdot \widetilde{\Sigma _i}\). In particular, we can write the expansion in the form \(\begin{equation*} \left(z \cdot \Sigma _k \right)^{O_k(1)} + \sum _{i} (z \cdot \widetilde{\Sigma _i})P_i(z \cdot \widetilde{\Sigma _1}, \dots , z \cdot \widetilde{\Sigma _k}) \end{equation*}\) where P is some polynomial in k variables. We can upper bound the coefficients of the polynomial in terms of \(\Delta , k\) and we also know that \(\begin{equation*} (z \cdot \widetilde{\Sigma _i})^2 \dashv _{O(1)} O_k(1)\Delta ^{O(1)} \end{equation*}\) due to the constraints in our system. Thus, we can bound the pseudoexpectation

via Cauchy Schwarz. Putting everything together the same way as in Lemma 4.11, we deduce

and now we are done.□

Putting Lemmas 4.11 and 4.12 together, we now prove that V is essentially contained within the span of the union of the top principal components of \(\widetilde{\mathbb {E}}[\widetilde{\Sigma _i} \widetilde{\Sigma _i}^T]\) over all i.

Lemma 4.13.

For each i, let \(M_i\) be the \(D \times D\) matrix given by

Assume that for a sufficiently small function f depending only on k, \(\begin{align*} \Delta \le (\epsilon ^{\prime })^{-f(k)} \\ w_{\min }, c \ge (\epsilon ^{\prime })^{f(k)} \end{align*}\) Let \(V_i\) be the subspace spanned by the top k singular vectors of \(M_i\). Then for all i, the projection of the true covariance matrix \(\Sigma _i\) onto the orthogonal complement of \(\textsf {spn}(V_1, \dots , V_k)\) has length at most \((\epsilon ^{\prime })^{\Omega _k(1)}\).

Proof.

Assume for the sake of contradiction that the desired statement is false for \(\Sigma _i\). Let z be the projection of \(\Sigma _i\) onto the orthogonal complement of \(\textsf {spn}(V_1, \dots , V_k)\). By Lemma 4.12, (2) \(\begin{equation} \sum _j z^TM_jz \ge \left(\frac{w_{\min }(z \cdot \Sigma _i)^{O_k(1)} - O_k(1)\epsilon ^{\prime } \Delta ^{O_k(1)}}{O_k(1)\Delta ^{O_k(1)} }\right)^2 \end{equation}\) so there is some j for which \(\begin{equation*} z^TM_j z \ge \frac{1}{k}\left(\frac{w_{\min }(z \cdot \Sigma _i)^{O_k(1)} - O_k(1)\epsilon ^{\prime } \Delta ^{O_k(1)}}{O_k(1)\Delta ^{O_k(1)} }\right)^2 \end{equation*}\) On the other hand, Lemma 4.11 implies that the sum of the singular values of \(M_j\) outside the top k is at most \(\begin{equation*} (\epsilon ^{\prime })^{2^{-k}}O_k(1) \left(\frac{\Delta }{w_{\min }c} \right)^{O_k(1)} \end{equation*}\) Since z is orthogonal to the span of the top-k singular vectors of \(M_j\), we get (3) \(\begin{equation} z^T M_j z \le (\epsilon ^{\prime })^{2^{-k}}O_k(1) \left(\frac{\Delta }{w_{\min }c} \right)^{O_k(1)} \left\Vert z\right\Vert _2^2 \end{equation}\) Note \(z \cdot \Sigma _i = \left\Vert z\right\Vert _2^2\) since z is a projection of \(\Sigma _i\) onto a subspace. Now combining (2) and (3) we get a contradiction unless \(\begin{equation*} \left\Vert z\right\Vert _2 \le (\epsilon ^{\prime })^{\Omega _k(1)}. \end{equation*}\)□

5.1 Distance between Gaussians

As mentioned earlier, we will first introduce a few tools for relating parameter distance and TV distance between Gaussians.

The following is a standard fact.

Next, we will prove a bound in the opposite direction, that when Gaussians are not too far apart in TV distance, then their parameters also cannot be too far apart.

As a corollary to the previous lemma, in a \(\delta\)-well conditioned mixture, all component means and covariances are close to the mean and covariance of the overall mixture.

5.2 Hermite Polynomial Estimation

Now we show how to estimate the Hermite polynomials of a \(\delta\)-well-conditioned mixture \(\mathcal {M}\) if we are given an \(\epsilon\)-corrupted sample (where \(\delta \ge \epsilon ^{f(k)}\) for some sufficiently small function \(f(k) \gt 0\) depending only on k). Our algorithm will closely mirror the algorithm in [39].

The first step will be to show that we can robustly estimate the mean and covariance of the mixture \(\mathcal {M}\) and then we will use these estimates to compute a linear transformation to place the mixture in isotropic position.

Using the above, we can place our mixture in isotropic position. This mirrors Proposition 4.2 in [39].

Once our mixture is placed in isotropic position, we will estimate the Hermite polynomials and then we will be able to use Theorem 4.1. The following lemma can be easily derived from the results in [39] (see Lemmas 2.7, 2.8, and 5.2 there).

Lemma 5.8.

Let \(\mathcal {M}\) be a mixture of Gaussians \(w_1 N(\mu _1, I + \Sigma _1) + \dots + w_k N(\mu _k, I + \Sigma _k)\). Then

where \(H_{m}(X,z)\) is defined as in Definition 3.4 and \(v_X(H_{m}(X, z))\) denotes vectorizing as a polynomial in X so that the entries of the vector are polynomials in z.

Kane [39] works with Hermite polynomial tensors, which are tensorized versions of the Hermite polynomials we are using. It is clear that these two notions are equivalent up to \(O_k(1)\) factors as long as m is \(O_k(1)\) (writing them as formal polynomials instead of tensors simply collapses symmetric entries of the tensor but this collapses at most \(O_m(1)\) entries together at once).

We can now combine everything in this section with Theorem 4.1 to complete the proof of Theorem 5.2.

6.1 SOS Program

To solve the clustering problem, we set up the same polynomial constraints as in Diakonikolas et al. [21]. Recall that Definition 2.4 gives a recipe for turning this into an SDP relaxation.

We will work with the same set of deterministic conditions on the samples as in Diakonikolas et al. [21]. These conditions hold with high probability for the uncorrupted points.

We will use the following key lemmas from [21]. The setup is exactly the same. Let \(X_1, \dots , X_n \in \mathbb {R}^d\) satisfy the \((\delta , \psi , t)\)-deterministic conditions (Definition 6.7) with respect to Gaussians \(G_1, \dots , G_k\). Let \(S_1, \dots , S_k\) be the partition guaranteed in the definition. Let \(Y_1, \dots , Y_n\) be an \(\epsilon\)-corruption of \(X_1, \dots , X_n\) and let \(\mathcal {A}\) be the clustering program (Definition 6.6) for \(Y_1, \dots , Y_n\). For indeterminates \(w_1, \dots , w_n\), define \(\begin{equation*} \alpha _i(w) = \sum _{j \in S_i}w_j. \end{equation*}\) Below we will assume \(\psi , \tau\) are smaller than some universal constants \(\psi _0, \tau _0 \gt 0\).

Recall in Claim 6.5 that there are essentially three different ways that two Gaussians can be separated in TV distance. We call these mean separation, variance separation, and covariance separation. The lemmas below roughly assert that if two Gaussians are separated in one of these ways, then a valid solution to the clustering program \(\mathcal {A}\) cannot assign significant weight to both of them.

Lemma 6.9 (Mean Separation, Restated from [21]).

For every \(\tau \gt 0\), there is \(s = \widetilde{O}(1/\tau ^2)\) such that if \(\epsilon , \delta \le s^{-O(s)}k^{-20}\) then for all \(a,b \in [k]\), all \(v \in \mathbb {R}^d\) and all sufficiently small \(\rho \gt 0\), if

then \(\begin{equation*} \mathcal {A} \vdash _{O(s)} \left(\frac{\alpha _a(w)\alpha _b(w)}{n^2} \right)^s \le (s\log 1/\psi)^{O(s)} \cdot \left(\frac{\langle v, \Sigma _av \rangle + \langle v, \Sigma _b v \rangle }{\langle \mu _a - \mu _b, v \rangle ^2} \right)^{\Omega (s)} + \rho ^{-O(s)}(\tau ^{\Omega (s)} + \epsilon ^{\Omega (s)}k^{O(s)}s^{O(s^2)} + \psi ^{\Omega (s)}). \end{equation*}\)

Lemma 6.10 (Variance Separation, Restated from [21]).

For every \(\tau \gt 0\), there is \(s = \widetilde{O}(1/\tau ^2)\) such that if \(\epsilon , \delta \le s^{-O(s)}k^{-20}\) then for all \(a,b \in [k]\), all \(v \in \mathbb {R}^d\) and all sufficiently small \(\rho \gt 0\), if

then \(\begin{equation*} \mathcal {A} \vdash _{O(s)} \left(\frac{\alpha _a(w)\alpha _b(w)}{n^2} \right)^s \le \psi ^{-O(s)} \cdot \left(s^{O(s)}\left(\frac{\langle v, \Sigma _a v \rangle }{\langle v, \Sigma _b v \rangle }\right)^{\Omega (s)} + \rho ^{-O(s)}(\tau ^{\Omega (s)} + \epsilon ^{\Omega (s)}k^{O(s)}s^{O(s^2)})\right) + \rho ^{-O(s)}\psi ^{\Omega (s)}. \end{equation*}\)

6.2 Clustering Algorithm

We use essentially the same clustering algorithm as [21].

Proof of Theorem 6.2

We can use Claim 6.8 to ensure that with \(1 - \gamma /2\) probability, the deterministic conditions in Definition 6.7 are satisfied for all submixtures and the various values of \(\delta , \psi , t\) that we will need.

First, if all pairs of components are \(D^{\prime }\) close, then returning the entire sample as one cluster suffices. Now, we may assume that there is some pair that is not \(D^{\prime }\)-close. We apply Claim 6.5 and let \(U,V\) be the partition of the components given by the claim. Let C be a sufficiently large function of \(k,D, \gamma\) that we will set later. We can do this as long as we ensure that \(D^{\prime }\) is a sufficiently large function of \(k,C\). We ensure that \(k^C \gt D\). Note that each of the pieces \(A_1, \dots , A_l\), must be entirely in U or in V because of our assumption about closeness between the components. We claim that

(4)

where we can make \(\gamma ^{\prime }\) sufficiently small in terms of \(\gamma ,D, k\) by choosing \(D^{\prime }\) and the functions \(f, F\) suitably.

Below we will let \(a,b\) be indices such that \(a \in U\) and \(b \in V\). If the first clause of Claim 6.5 is satisfied, then we can take \(\rho = poly(1/k)\) and for \(\tau\) sufficiently small in terms of \(\gamma , k, D, C\), we have

Summing over all \(a \in U, b \in V\), this gives (4).

If the second clause of Claim 6.5 is satisfied, then we can take \(\begin{equation*} \rho = \min _{a \in U}\frac{v^T\Sigma _a v}{v^T \Sigma v} \ge \frac{1}{k^{2Ck}} \end{equation*}\) We choose \(\tau\) sufficiently small in terms of \(\gamma , k, C, D\) and combining with the fact that \((v^t\Sigma _a v) \ge k^{C}(v^T\Sigma _bv)\) for all \(a \in U, b \in V\) we get

Finally, when the third clause of Claim 6.5 is satisfied it follows similarly after setting \(A = A_{ab}\). In all cases, we now have (4). The next step will be to analyze our random sampling to select the subset R. First note

Next we analyze the intersections with the two sides of the partition \(U,V\). We will slightly abuse notation and use \(i \in U\) when \(i \in \cup _{j \in U}S_j\) and it is clear from context that we are indexing the samples. Conditioned on the first index that is randomly chosen satisfying \(i \in U\) then

repeating the same argument for when \(i \in V\), we have \(\mathbb {E}[\min (|R \cap U| , |R \cap V|)] \le \gamma ^{\prime }kn\). Finally, we lower bound the expected number of new elements that R adds to the list L. This quantity is

where by \(j \notin L\) we mean j is not in the union of all previous subsets in the list L. Note that indicator functions of the components \(S_1, \dots , S_k\) are all valid pseudoexpectations and since we are picking the pseudoexpectation that maximizes \(\sum _{j \notin L}\widetilde{\mathbb {E}}[w_j]\), the expected number of new elements added to L is at least \(\begin{equation*} \frac{n - |\cup _{R \in L}R|}{k} \end{equation*}\) Now we analyze the recombination step once we finalize \(L = \lbrace R_1, \dots , R_m \rbrace\). For any sufficiently small function \(h(k, \gamma , D)\), we claim that by choosing \(D^{\prime }\) and the functions \(f, F\) appropriately, we can ensure with \(1 - h(k, D, \gamma)\) probability, there is some recombination that gives a \(1 - h(k, D, \gamma)\)-corrupted sample of the submixture corresponding to U. To see this, it suffices to set \(\eta \lt h(k, D, \gamma)\) and then look at the first \(m^{\prime } = 100 k \log 1/h(k,D, \gamma)\) subsets in L. Their union has expected size \((1 - h(k,D, \gamma)^{100})n\). Next, among \(R_1, \dots , R_{m^{\prime }}\),

If we ensure that \(\gamma ^{\prime }\) is sufficiently small in terms of \(\gamma ,D, k\), then using Markov’s inequality, with \(1 - h(k,D, y)\) probability, there is some recombination that gives a \(1 - h(k, D, \gamma)\)-corrupted sample of the submixture corresponding to U. We can make the same argument for V. Now we can recurse and repeat the argument because each of these submixtures only contains at most \(k-1\) true components.□

6.3 Improved Clustering Result from [4]

In [4], the authors obtain a rough clustering result similar to Theorem 6.2 but are able to remove the bounded fractionality assumption. Their result is restated using our notation below.

Observe that Theorem 6.12 is the same as Theorem 6.2 but with the bounded fractionality assumption removed. Theorem 6.2 is the only source of the bounded fractionality assumption in our paper. In the subsequent sections, replacing all uses of Theorem 6.2 with Theorem 6.12 allows us to remove the bounded fractionality assumption from our main result.

7.1 Distance Between Gaussians

We will need to prove a few preliminary results. The main lemma we prove in this section is the following, which gives a stronger bound than the triangle inequality for TV distance between Gaussians.

Note that this result is not true for arbitrary distributions \(A,B,C\). We actually need to exploit the fact that \(A,B,C\) are Gaussian.

Our proof will parallel results in Section 8 of [21]. First, a definition:

The key ingredient is the following result from [21]:

7.2 Full Algorithm

We are now ready to prove Theorem 7.1. We begin by describing the algorithm. Our full algorithm consists of several phases.

The algorithm Learn Parameters (well-conditioned) for learning the parameters of a well-conditioned mixture of Gaussians (see Theorem 5.2) is summarized in Algorithm 2.

Our full algorithm is described in the next block Algorithm 3.

7.3 Analysis of Full Algorithm

The first step will be to show that among the first set of candidate components that we output, there are some that are within constant distance (say \(\lt c(k)\) for some sufficiently small function c) of the true components.

Our next step is to show that if our algorithm starts with component estimates that are accurate within some constant and guesses a good set of clusters, then the resulting subsamples (after assigning according to maximum likelihood) are equivalent to \({\rm poly}(\epsilon)\)-corrupted samples from the corresponding submixtures. First, we prove a preliminary claim which implies that a good set of clusters exists.

We will also need the following results about the VC-dimension of hypotheses obtained by comparing the density functions of two mixtures of Gaussians. The reason we need these VC dimension bounds is that we will need to argue that given any constant-accuracy estimates, we can obtain a clustering that is \({\rm poly}(\epsilon)\) accurate. While naively this would require union bounding over infinitely many possibilities for the initial estimates, the VC dimension bound allows to get around this and obtain uniform convergence over all possible initial estimates.

Technically for our clustering result, we only need the VC dimension bound for single Gaussians (instead of mixtures of k Gaussians). However, we will need the VC dimension bound for mixtures of Gaussians later when we do hypothesis testing so we state the full result below. First we need a definition.

It is a standard result in learning theory that for a hypothesis class with bounded VC dimension, taking a polynomial number of samples suffices to get a good approximation for all hypotheses in the class.

Lemma 7.9 ([54]).

Let \(\mathcal {H}\) be a hypothesis class of functions from some domain \(\mathcal {X}\) to \(\lbrace 0,1 \rbrace\) with VC dimension V. Let \(\mathcal {D}\) be a distribution on \(\mathcal {X}\). Let \(\epsilon , \delta \gt 0\) be parameters. Let S be a set of \(n = {\rm poly}(V, 1/\epsilon , \log 1/\delta)\) i.i.d samples from \(\mathcal {D}\). Then with \(1 - \delta\) probability, for all \(f \in \mathcal {H}\)

Now we can prove our lemma about obtaining a \({\rm poly}(\epsilon)\)-accurate clustering into submixtures when given constant-accuracy estimates for the components.

Combining Lemma 7.5, Claim 7.6, Lemma 7.10, and Theorem 5.2, we can show that at least one of the sets of candidate parameters that our algorithm outputs is close to the true parameters.

We have shown that our algorithm recovers a list of candidate mixtures, at least one of which is close to the true mixture. The last result that we need is a hypothesis testing routine. This is similar to the hypothesis testing result in [39]. However, there is a subtle difference that the samples we use to hypothesis test may not be independent of the hypotheses. This is because the adversary sees all of the data points and may corrupt the data in a way to affect the list of hypotheses that we output. Thus, we must prove that given an \(\epsilon\)-corrupted sample and any list of hypotheses with the promise that at least one of the hypotheses is close to the true distribution, we must output a hypothesis that is close to the true distribution.

We can now complete the proof of our main theorem.

8.1 Improving the Separation Assumption

With simple modifications to the analysis, we obtain the following improvement of Theorem 7.1 in [45].

To go from Theorem 7.1 to 8.4, the main idea to remove the constant separation assumption is just that we can find a scale \(\delta\) such that all pairs of components either have \(d_{\textsf {TV}}(G_i, G_j) \le \delta\) or \(d_{\textsf {TV}}(G_i, G_j) \ge \delta ^{\prime }\) for \(\delta ^{\prime } \gg \delta\). This is possible because the number of components k is a constant. We can then merge components whose TV distance is less than \(\delta\), treating them as the same component and the remaining components will be sufficiently separated. See [45] for more details.

The results of [45] were obtained using the previous clustering subroutine of [21]. If we instead plug in the updated clustering results of [4] (see Theorem 6.12 vs Theorem 6.2), we can remove the bounded fractionality assumption on the mixing weights. Also, we can ensure that the algorithm outputs a unique mixture instead of a list by running the hypothesis testing routine in Lemma 7.12.

Finally, combining the above with identifiability (Theorem 8.1), we immediately get an improved version of our main theorem for parameter learning.