Learning Mixtures of Gaussians in High Dimensions

Introduction

Learning mixtures of Gaussians is a fundamental problem in statistics and learning theory, whose study dates back to Pearson (1894). Gaussian mixture models arise in numerous areas including physics, biology and the social sciences (McLachlan and Peel (2004); Titterington et al. (1985)), as well as in image processing (Reynolds and Rose (1995)) and speech (Permuter et al. (2003)).

In a Gaussian mixture model, there are $k$ unknown $n$-dimensional multivariate Gaussian distributions. Samples are generated by first picking one of the $k$ Gaussians, then drawing a sample from that Gaussian distribution. Given samples from the mixture distribution, our goal is to estimate the means and covariance matrices of these underlying Gaussian distributions This is different from the problem of density estimation considered in Feldman et al. (2006); Chan et al. (2014).

This problem has a long history in theoretical computer science. The seminal work of Dasgupta (1999) gave an algorithm for learning spherical Gaussian mixtures when the means are well separated. Subsequent works (Dasgupta and Schulman (2000); Sanjeev and Kannan (2001); Vempala and Wang (2004); Brubaker and Vempala (2008)) developed better algorithms in the well-separated case, relaxing the spherical assumption and the amount of separation required.

When the means of the Gaussians are not separated, after several works (Belkin and Sinha (2009); Kalai et al. (2010)), Belkin and Sinha (2010) and Moitra and Valiant (2010) independently gave algorithms that run in polynomial time and with polynomial number of samples for a fixed number of Gaussians. However, both running time and sample complexity depend super exponentially on the number of components $k$ In fact, it is in the order of $O({e^{O(k)}}^{k})$ as shown in Theorem 11.3 in Valiant (2012). . Their algorithm is based on the method of moments introduced by Pearson (1894): first estimate the $O(k)$-order moments of the distribution, then try to find the parameters that agree with these moments. Moitra and Valiant (2010) also show that the exponential dependency of the sample complexity on the number of components is necessary, by constructing an example of two mixtures of Gaussians with very different parameters, yet with exponentially small statistical distance.

Recently, Hsu and Kakade (2013) applied spectral methods to learning mixture of spherical Gaussians. When $n\geq k+1$ and the means of the Gaussians are linearly independent, their algorithm can learn the model in polynomial time and with polynomial number of samples. This result suggests that the lower bound example in Moitra and Valiant (2010) is only a degenerate case in high dimensional space. In fact, most (in general position) mixture of spherical Gaussians are easy to learn. This result is also based on the method of moments, and only uses second and third moments. Several follow-up works (Bhaskara et al. (2014); Anderson et al. (2013)) use higher order moments to get better dependencies on $n$ and $k$.

However, the algorithm in Hsu and Kakade (2013) as well as in the follow-ups all make strong requirements on the covariance matrices. In particular, most of them only apply to learning mixture of spherical Gaussians. For mixture of Gaussians with general covariance matrices, the best known result is still Belkin and Sinha (2010) and Moitra and Valiant (2010), which algorithms are not polynomial in the number of components $k$. This leads to the following natural question:

Question: Is it possible to learn most mixture of Gaussians in polynomial time using a polynomial number of samples?

In this paper, we give an algorithm that learns most mixture of Gaussians in high dimensional space (when $n\geq\Omega(k^{2})$), and the argument is formalized under the smoothed analysis framework first proposed in Spielman and Teng (2004).

In the smoothed analysis framework, the adversary first choose an arbitrary mixture of Gaussians. Then the mean vectors and covariance matrices of this Gaussian mixture are randomly perturbed by a small amount $\rho$ See Definition 3.2 in Section 3.1 for the details.. The samples are then generated from the Gaussian mixture model with the perturbed parameters. The goal of the algorithm is to learn the perturbed parameters from the samples.

The smoothed analysis framework is a natural bridge between worst-case and average-case analysis. On one hand, it is similar to worst-case analysis, as the adversary chooses the initial instance, and the perturbation allowed is small. On the other hand, even with small perturbation, we may hope that the instance be different enough from degenerate cases. A successful algorithm in the smoothed analysis setting suggests that the bad instances must be very “sparse” in the parameter space: they are highly unlikely in any small neighborhood of any instance. Recently, the smoothed analysis framework has also motivated several research work (Kalai et al. (2009) Bhaskara et al. (2014)) in analyzing learning algorithms.

In the smoothed analysis setting, we show that it is easy to learn most Gaussian mixtures:

(informal statement of Theorem 3.4) In the smoothed analysis setting, when $n\geq\Omega(k^{2})$, given samples from the perturbed $n$-dimensional Gaussian mixture model with $k$ components, there is an algorithm that learns the correct parameters up to accuracy $\epsilon$ with high probability, using polynomial time and number of samples.

An important step in our algorithm is to learn Gaussian mixture models whose components all have mean zero, which is also a problem of independent interest (Zoran and Weiss (2012)). Intuitively this is also a “hard” case, as there is no separation in the means. Yet algebraically, this case gives rise to a novel tensor decomposition algorithm. The ideas for solving this decomposition problem are then generalized to tackle the most general case.

(informal statement of Theorem 3.5) In the smoothed analysis setting, when $n\geq\Omega(k^{2})$, given samples from the perturbed mixture of zero-mean $n$-dimensional Gaussian mixture model with $k$ components, there is an algorithm that learns the parameters up to accuracy $\epsilon$ with high probability, using polynomial running time and number of samples.

Organization

The main part of the paper will focus on learning mixtures of zero-mean Gaussians. The proposed algorithm for this special case contains most of the new ideas and techniques. In Section 2 we introduce the notations for matrices and tensors which are used to handle higher order moments throughout the discussion. Then in Section 3 we introduce the smoothed analysis model for learning mixture of Gaussians and discuss the moment structure of mixture of Gaussians, then we formally state our main theorems. Section 4 outlines our algorithm for learning zero-mean mixture of Gaussians. The details of the steps are presented in Section 5. The detailed proofs for the correctness and the robustness are deferred to Appendix (Sections B to D). In Section 6 we briefly discuss how the ideas for zero-mean case can be generalized to learning mixture of nonzero Gaussians, for which the detailed algorithm and the proofs are deferred to Appendix F.

Notations

We use $[n]$ to denote the set $\{1,2,...,n\}$ and $[n]\times[n]$ to denote the set $\{(i,j):i,j\in[n]\}$. These are often used as indices of matrices.

Symmetric matrices

Linear subspaces

Tensors

In particular, the vector given by $X(\mathbf{e}_{i},\mathbf{e}_{j},I)$ is the one-dimensional slice of the 3-way array, with the index for the first dimension to be $i$ and the second dimension to be $j$.

Matrix Products

Main results

In this section, we first formally introduce the smoothed analysis framework for our problem and state our main theorems. Then we will discuss the structure of the moments of Gaussian mixtures, which is crucial for understanding our method of moments based algorithm.

As an interesting special case of the general model, we also consider the mixture of “zero-mean” Gaussians, which has $\mu^{(i)}=0$ for all components $i\in[k]$.

A sample $x$ from a mixture of Gaussians is generated in two steps:

Sample $h\in[k]$ from a multinomial distribution, with probability $\Pr[h=i]=\omega_{i}$ for $i\in[k]$.

The learning problem asks to estimate the parameters of the underlying mixture of Gaussians:

Given $N$ samples $x_{1},x_{2},...,x_{N}$ drawn i.i.d. from a mixture of Gaussians $\mathcal{G}=\{(\omega_{i},\mu^{(i)},\Sigma^{(i)})\}_{i\in[k]}$, an algorithm learns the mixture of Gaussians with accuracy $\epsilon$, if it outputs an estimation $\widehat{\mathcal{G}}=\{(\widehat{\omega}_{i},\widehat{\mu}^{(i)},\widehat{\Sigma}^{(i)})\}_{i\in[k]}$ such that there exists a permutation $\pi$ on $[k]$, and for all $i\in[k]$, we have $|\widehat{\omega}_{i}-\omega_{\pi(i)}|\leq\epsilon$, $\|\widehat{\mu}^{(i)}-\mu^{(\pi(i))}\|\leq\epsilon$ and $\|\widehat{\Sigma}^{(i)}-\Sigma^{(\pi(i))}\|\leq\epsilon$.

In the worst case, learning mixture of Gaussians is a information theoretically hard problem (Moitra and Valiant (2010)). There exists worst-case examples where the number of samples required for learning the instance is at least exponential in the number of components $k$ (McLachlan and Peel (2004)). The non-convexity arises from the hidden variable $h$: without knowing $h$ we cannot determine which Gaussian component each sample comes from.

The smoothed analysis framework provides a way to circumvent the worst case instances, yet still studying this problem in its most general form. The basic idea is that, with high probability over the small random perturbation to any instance, the instance will not be a “worst-case” instance, and actually has reasonably good condition for the algorithm.

Next, we show how the parameters of the mixture of Gaussians are perturbed in our setup.

For $\rho<1/n$, a $\rho$-smooth $n$-dimensional $k$-component mixture of Gaussians $\widetilde{\mathcal{G}}=\{(\widetilde{\omega}_{i},\widetilde{\mu}^{(i)},\widetilde{\Sigma}^{(i)})\}_{i\in[k]}\in\mathcal{G}_{n,k}$ is generated as follows:

Choose an arbitrary (could be adversarial) instance $\mathcal{G}=\{(\omega_{i},\mu^{(i)},\Sigma^{(i)})\}_{i\in[k]}\in\mathcal{G}_{n,k}$. Scale the distribution such that $0\preceq\Sigma^{(i)}\preceq{1\over 2}I_{n}$ and $\|\mu^{(i)}\|\leq{1\over 2}$ for all $i\in[k]$.

Set $\widetilde{\omega}_{i}=\omega_{i}$, $\widetilde{\mu}^{(i)}=\mu^{(i)}+\delta_{i}$, $\widetilde{\Sigma}^{(i)}=\Sigma^{(i)}+\Delta_{i}$.

Choose the diagonal entries of $\widetilde{\Sigma}^{(i)}$ arbitrarily, while ensuring the positive semi-definiteness of the covariance matrix $\widetilde{\Sigma}^{(i)}$, and the diagonal entries are upper bounded by $1$. The perturbation procedure fails if this step is infeasible Note that by standard random matrix theory, with high probability the 4-th step is feasible and the perturbation procedure in Definition 3.2 succeeds. Also, with high probability we have $\|\widetilde{\mu}^{(i)}\|\leq 1$ and $0\preceq\widetilde{\Sigma}^{(i)}\preceq I_{n}$ for all $i\in[k]$. .

A $\rho$-smooth zero-mean mixture of Gaussians is generated using the same procedure, except that we set $\widetilde{\mu}^{(i)}=\mu^{(i)}=0$, for all $i\in[k]$.

When the original matrix is of low rank, a simple random perturbation may not lead to a positive semidefinite matrix, which is why our procedure of perturbation is more restricted in order to guarantee that the perturbed matrix is still a valid covariance matrix.

There could be other ways of locally perturbing the covariance matrix. Our procedure actually gives more power to the adversary as it can change the diagonals after observing the perturbations for other entries. Note that with high probability if we just let the new diagonal to be $5\sqrt{n}\rho$ larger than the original ones, the resulting matrix is still a valid covariance matrix. In other words, the adversary can always keep the perturbation small if it wants to.

Instead of the worst-case problem in Definition 3.1, our algorithms work on the smoothed instance. Here the model first gets perturbed to $\widetilde{\mathcal{G}}=\{(\widetilde{\omega}_{i},\widetilde{\mu}^{(i)},\widetilde{\Sigma}^{(i)})\}_{i\in[k]}$, the samples are drawn according to the perturbed model, and the algorithm tries to learn the perturbed parameters. We give a polynomial time algorithm in this case:

Consider a $\rho$-smooth mixture of Gaussians $\widetilde{\mathcal{G}}=\{(\widetilde{\omega}_{i},\widetilde{\mu}^{(i)},\widetilde{\Sigma}^{(i)})\}_{i\in[k]}\in\mathcal{G}_{n,k}$ for which the number of components is at least Note that the algorithms of Belkin and Sinha (2010) and Moitra and Valiant (2010) run in polynomial time for fixed $k$. $k\geq C_{0}$ and the dimension $n\geq C_{1}k^{2}$, for some fixed constants $C_{0}$ and $C_{1}$. Suppose that the mixing weights $\widetilde{\omega}_{i}\geq\omega_{o}$ for all $i\in[k]$. Given $N$ samples drawn i.i.d. from $\widetilde{\mathcal{G}}$, there is an algorithm that learns the parameters of $\widetilde{\mathcal{G}}$ up to accuracy $\epsilon$, with high probability over the randomness in both the perturbation and the samples. Furthermore, the running time and number of samples $N$ required are both upper bounded by $\text{poly}(n,k,1/\omega_{o},1/\epsilon,1/\rho)$.

To better illustrate the algorithmic ideas for the general case, we first present an algorithm for learning mixtures of zero-mean Gaussians. Note that this is not just a special case of the general case, as with the smoothed analysis, the zero mean vectors are not perturbed.

Consider a $\rho$-smooth mixture of zero-mean Gaussians $\widetilde{\mathcal{G}}=\{(\widetilde{\omega}_{i},0,\widetilde{\Sigma}^{(i)})\}_{i\in[k]}\in\mathcal{G}_{n,k}$ for which the number of components is at least $k\geq C_{0}$ and the dimension $n\geq C_{1}k^{2}$, for some fixed constants $C_{0}$ and $C_{1}$. Suppose that the mixing weights $\widetilde{\omega}_{i}\geq\omega_{o}$ for all $i\in[k]$. Given $N$ samples drawn i.i.d. from $\widetilde{\mathcal{G}}$, there is an algorithm that learns the parameters of $\widetilde{\mathcal{G}}$ up to accuracy $\epsilon$, with high probability over the randomness in both the perturbation and the samples. Furthermore, the running time and number of samples $N$ are both upper bounded by $\text{poly}(n,k,1/\omega_{o},1/\epsilon,1/\rho)$.

Throughout the paper we always assume that $n\geq C_{1}k^{2}$ and $\widetilde{\omega}_{i}\geq\omega_{o}$.

2 Moment Structure of Mixture of Gaussians

Our algorithm is also based on the method of moments, and we only need to estimate the $3$-rd, the $4$-th and the $6$-th order moments. In this part we briefly discuss the structure of $4$-th and $6$-th moments in the zero-mean case ($3$-rd moment is always 0 in the zero-mean case). These structures are essential to the proposed algorithm. For more details, and discussions on the general case see Appendix A.

where $y^{(i)}$ corresponds to the $n$-dimensional zero-mean Gaussian distribution $\mathcal{N}(0,\Sigma^{(i)})$. The moments for each Gaussian component are characterized by Isserlis’s theorem as below:

Let $(y_{1},\dots,y_{2t})$ be a multivariate zero-mean Gaussian random vector $\mathcal{N}(0,\Sigma)$, then

where the summation is taken over all distinct ways of partitioning $y_{1},\dots,y_{2t}$ into $t$ pairs, which correspond to all the perfect matchings in a complete graph.

Ideally, we would like to obtain the following quantities (recall $n_{2}={n+1\choose 2}$):

Note that the entries in $X_{4}$ and $X_{6}$ are quadratic and cubic monomials of the covariance matrices, respectively. If we have $X_{4}$ and $X_{6}$, the tensor decomposition algorithm in Anandkumar et al. (2014) can be immediately applied to recover $\omega_{i}$’s and $\Sigma^{(i)}$’s under mild conditions. It is easy to verify that those conditions are indeed satisfied with high probability in the smoothed analysis setting.

By Isserlis’s theorem, the entries of the moments $M_{4}$ and $M_{6}$ are indeed quadratic and cubic functions of the covariance matrices, respectively. However, the structure of the true moments $M_{4}$ and $M_{6}$ have more symmetries, consider for example,

Note that due to symmetry, the number of distinct entries in $M_{4}$ ( ${n+3\choose 4}\approx n^{4}/24$) is much smaller than the number of distinct entries in $X_{4}$ (${n_{2}+1\choose 2}\approx n^{4}/8$). Similar observation can be made about $M_{6}$ and $X_{6}$.

Therefore, it is not immediate how to find the desired $X_{4}$ and $X_{6}$ based on $M_{4}$ and $M_{6}$. We call the moments $M_{4},M_{6}$ the folded moments as they have more symmetry, and the corresponding $X_{4},X_{6}$ the unfolded moments. One of the key steps in our algorithm is to unfold the true moments $M_{4},M_{6}$ to get $X_{4},X_{6}$ by exploiting special structure of $M_{4},M_{6}$.

In some cases, it is easier to restrict our attention to the entries in $M_{4}$ with indices corresponding to distinct variables. In particular, we define

Moreover $\mathcal{F}_{4}$ is a projection from a ${n_{2}+1\choose 2}$-dimensional subspace to a $n_{4}$-dimensional subspace, and $\mathcal{F}_{6}$ is a projection from a ${n_{2}+2\choose 3}$-dimensional subspace to a ${n_{6}}$-dimensional subspace.

Algorithm Outline for Learning Mixture of Zero-Mean Gaussians

In this section, we present our algorithm for learning zero-mean Gaussian mixture model. The algorithmic ideas and the analysis are at the core of this paper. Later we show that it is relatively easy to generalize the basic ideas and the techniques to handle the general case.

For simplicity we state our algorithm using the exact moments $\widetilde{M}_{4}$ and $\widetilde{M}_{6}$, while in implementation the empirical moments $\widehat{M}_{4}$ and $\widehat{M}_{6}$ obtained with the samples are used. In later sections, we verify the correctness of the algorithm and show that it is robust: the algorithm learns the parameters up to arbitrary accuracy using polynomial number of samples.

Span Finding: Find the span of covariance matrices .

For a set of indices $\mathcal{H}\subset[n]$ of size $|\mathcal{H}|=\sqrt{n}$, find the span:

Find the span of the covariance matrices with the columns projected onto $\mathcal{S}^{\perp}$, namely,

For two disjoint sets of indices $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, repeat Step 1 (a) and Step 1 (b) to obtain $\mathcal{U}_{1}$ and $\mathcal{U}_{2}$, namely the span of covariance matrices projected onto two subspaces $\mathcal{S}_{1}^{\perp}$ and $\mathcal{S}_{2}^{\perp}$. Merge $\mathcal{U}_{1}$ and $\mathcal{U}_{2}$ to obtain the span of covariance matrices $\mathcal{U}$:

The unfolded moments $\widetilde{X}_{4},\widetilde{X}_{6}$ are then given by $\widetilde{X}_{4}=U\widetilde{Y}_{4}U^{\top},\widetilde{X}_{6}=\widetilde{Y}_{6}(U^{\top},U^{\top},U^{\top}).$

Tensor Decomposition: learn $\widetilde{\omega}_{i}$ and $\widetilde{\Sigma}^{(i)}$ from $\widetilde{Y}_{4}$ and $\widetilde{Y}_{6}$. Given $U$, and given $\widetilde{Y}_{4}$ and $\widetilde{Y}_{6}$ which are relate to the parameters as follows:

we apply tensor decomposition techniques to recover $\widetilde{\Sigma}^{(i)}$’s and $\widetilde{\omega}_{i}$’s.

Implementing the Steps for Mixture of Zero-Mean Gaussians

In this part we show how to accomplish each step of the algorithm outlined in Section 4 and sketch the proof ideas.

For each step, we first explain the detailed algorithm, and list the deterministic conditions on the underlying parameters as well as on the exact moments for the step to work correctly. Then we show that these deterministic conditions are satisfied with high probability over the $\rho$-perturbation of the parameters in the smoothed analysis setting. In order to analyze the sample complexity, we further show that when we are given the empirical moments which are close to the exact moments, the output of the step is also close to that in the exact case.

In particular we show the correctness and the stability of each step in the algorithm with two main lemmas: the first lemma shows that with high probability over the random perturbation of the covariance matrices, the exact moments satisfy the deterministic conditions that ensure the correctness of each step; the second lemma shows that when the algorithm for each step works correctly, it is actually stable even when the moments are estimated from finite samples and have only inverse polynomial accuracy to the exact moments.

The detailed proofs are deferred to Section B to D in the appendix.

Given the 4-th order moments $\widetilde{M}_{4}$, Step 1 finds the span of covariance matrices $\mathcal{U}$ as defined in (6). Note that by definition of the unfolded moments $\widetilde{X}_{4}$ in (1), the subspace $\mathcal{U}$ coincides with the column span of the matrix $\widetilde{X}_{4}$.

By Lemma 3.7, we know that the entries in $\widetilde{M}_{4}$ are linear mappings of entries in $\widetilde{X}_{4}$. Since the matrix $\widetilde{X}_{4}$ is of low rank ($k\ll n_{2}$), this corresponds to the matrix sensing problem first studied in Recht et al. (2010). In general, matrix sensing problems can be hard even when we have many linear observations (Hardt et al. (2014b)). Previous works (Recht et al. (2010); Hardt et al. (2014a); Jain et al. (2013)) showed that if the linear mapping satisfy matrix RIP property, one can uniquely recover $\widetilde{X}_{4}$ from $\widetilde{M}_{4}$.

However, properties like RIP do not hold in our setting where the linear mapping is determined by Isserlis’ Theorem. We can construct two different mixtures of Gaussians with different unfolded moments $\widetilde{X}_{4}$, but the same folded moment $\widetilde{M}_{4}$ (see Section A.3). Therefore the existing matrix recovery algorithm cannot be applied, and we need to develop new tools by exploiting the special moment structure of Gaussian mixtures.

Step 1 (a). Find the Span of a Subset of Columns of the Covariance Matrices.

The key observation for this step is that if we hit $\widetilde{M}_{4}$ with three basis vectors, we get a vector that lies in the span of the columns of the covariance matrices:

For a mixture of zero-mean Gaussians $\mathcal{G}=\{(\omega_{i},0,\Sigma^{(i)})\}_{i\in[k]}\in\mathcal{G}_{n,k}$, the one-dimensional slices of the 4-th order moments $M_{4}$ are given by:

In particular, if we pick the indices $j_{1},j_{2},j_{3}$ in the index set $\mathcal{H}$, the vector $M_{4}(\mathbf{e}_{j_{1}},\mathbf{e}_{j_{2}},\mathbf{e}_{j_{3}},I)$ lies in the desired span $\mathcal{S}=\left\{\Sigma^{(i)}_{[:,j]}:i\in[k],j\in\mathcal{H}\right\}$.

We shall partition the set $\mathcal{H}$ into three disjoint subsets $\mathcal{H}^{(i)}$ of equal size $\sqrt{n}/3$, and pick $j_{i}\in H^{(i)}$ for $i=1,2,3$. In this way, we have $(|\mathcal{H}|/3)^{3}=\Omega(n^{1.5})$ such one-dimensional slices of $M_{4}$, which all lie in the desired subspace $\mathcal{S}$. Moreover, the dimension of the subspace $\mathcal{S}$ is at most $k|\mathcal{H}|\ll n^{1.5}$. Therefore, with the $\rho$-perturbed parameters $\widetilde{\Sigma}^{(i)}$’s, we can expect that with high probability the slices of $\widetilde{M}_{4}$ span the entire subspace $\mathcal{S}$.

We first show that this deterministic condition is satisfied with high probability by bounding the $k|\mathcal{H}|$-th singular value of $\widetilde{Q}_{S}$ with smoothed analysis.

Given the exact 4-th order moments $\widetilde{M}_{4}$, for any index set $\mathcal{H}$ of size $|\mathcal{H}|=\sqrt{n}$, With high probability, the $k|\mathcal{H}|$-th singular value of $\widetilde{Q}_{S}$ is at least $\Omega(\omega_{o}\rho^{2}n)$.

The proof idea involves writing the matrix $\widetilde{Q}_{S}$ as a product of three matrices, and using the results on spectral properties of random matrices Rudelson and Vershynin (2009) to show that with high probability the smallest singular value of each factor is lower bounded.

Since this step only involves the singular value decomposition of the matrix $\widetilde{Q}_{S}$, we then use the standard matrix perturbation theory to show that this step is stable:

Note that we need the high dimension assumption ($n\gg k$) to guarantee the correctness of this step: in order to span the subspace $\mathcal{S}$, the number of distinct vectors should be equal or larger than the dimension of the subspace, namely $|\mathcal{H}|^{3}\geq k|\mathcal{H}|$; and the subspace should be non-trivial, namely $k|\mathcal{H}|<n$. These two inequalities suggest that we need $n\geq\Omega(k^{1.5})$. However, we used the stronger assumption $n\geq\Omega(k^{2})$ to obtain the lower bound of the smallest singular value in the proof.

Step 1 (b). Find the Span of Projected Covariance Matrices.

In this step, we continue to use the structural properties of the 4-th order moments. In particular, we look at the two-dimensional slices of $M_{4}$ obtained by hitting it with two basis vectors:

For a mixture of zero-mean Gaussians $\mathcal{G}=\{(\omega_{i},0,\Sigma^{(i)})\}_{i\in[k]}\in\mathcal{G}_{n,k}$, the two-dimensional slices of the 4-th order moments $M_{4}$ are given by:

Note that if we take the indices $j_{1}$ and $j_{2}$ in the index set $\mathcal{H}$, the slice $M_{4}(\mathbf{e}_{j_{1}},\mathbf{e}_{j_{2}},I,I)$ is almost in the span of the covariance matrices, except $2k$ additive rank-one terms in the form of $\Sigma^{(i)}_{[:,j_{1}]}(\Sigma^{(i)}_{[:,j_{2}]})^{\top}$. These rank-one terms can be eliminated by projecting the slice to the subspace $\mathcal{S}^{\perp}$ obtained in Step 1 (a), namely,

and this projected two-dimensional slice lies in the desired span $\mathcal{U}_{S}$ as defined in (5). Moreover, there are ${|\mathcal{H}|+1\choose 2}=\Omega(n)$ such projected two-dimensional slices, while the dimension of the desired span $\mathcal{U}_{S}$ is at most $k$.

We show that this deterministic condition is satisfied by bounding the $k$-th singular value of $\widetilde{Q}_{U_{S}}$ in the smoothed analysis setting:

Given the exact 4-th order moments $\widetilde{M}_{4}$, with high probability, the $k$-th singular value of $\widetilde{Q}_{U_{S}}$ is at least $\Omega(\omega_{o}\rho^{2}n^{1.5})$.

Similar to Lemma 5.3, the proof is based on writing the matrix $Q_{U_{S}}$ as a product of three matrices, then bound their $k$-th singular values using random matrix theory. The stability analysis also relies on the matrix perturbation theory.

Given the empirical 4-th order moments $\widehat{M}_{4}=\widetilde{M}_{4}+E_{4}$, assume that the absolute value of entries of $E_{4}$ are at most $\delta_{2}$. Also, given the output $\text{Proj}_{\widehat{S}^{\perp}}$ from Step 1 (a), and assume that $\|\text{Proj}_{\widehat{S}^{\perp}}-\text{Proj}_{\widetilde{S}^{\perp}}\|\leq\delta_{1}$. When $\delta_{1}$ and $\delta_{2}$ are inverse polynomially small, we have $\|\text{Proj}_{\widehat{U}_{S}}-\text{Proj}_{\widetilde{U}_{S}}\|\leq O\left({n^{2.5}\left(\delta_{2}+2\delta_{1}\right)/\sigma_{k}(\widetilde{Q}_{U_{S}})}\right)$.

Note that for a given index set $\mathcal{H}$, the span $\mathcal{U}_{S}$ obtained in Step 1 (b) only gives partial information about the span of the covariance matrices. The idea of getting the span of the full covariance matrices is to obtain two sets of such partial information and then merge them.

In order to achieve that, we repeat Step 1 (a) and Step 1 (b) for two disjoint sets $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, each of size $\sqrt{n}$. The two subspace $S_{1}$ and $S_{2}$ thus correspond to the span of two disjoint sets of covariance matrix columns. Therefore, we can hope that $U_{1}$ and $U_{2}$, the span of covariance matrices projected to $S_{1}^{\perp}$ and $S_{2}^{\perp}$ contain enough information to recover the full span $U$.

In particular, we prove the following claim:

The main idea in the proof is that since $s$ is not too large, the two subspaces $S_{1}^{\perp}$ and $S_{2}^{\perp}$ have a large intersection. Using this intersection we can “align” the two basis $V_{1}$ and $V_{2}$ and obtain $V_{1}^{\dagger}V_{2}$, and then it is easy to merge the two projections of the same matrix (instead of a subspace).

Moreover, we show that when applying this result to the projected span of covariance matrices, we have $s=k|\mathcal{H}|\leq n/3$, and the two deterministic conditions $\sigma_{2s}([S_{1},S_{2}])>0$ and $\sigma_{k}(\text{Proj}_{S_{3}}V_{1})>0$ are indeed satisfied with high probability over the parameter perturbation. The detailed smoothed analysis (Lemma B.13 and B.14) and the stability analysis (Lemma B.11) are provided in Section B.3 in the appendix.

We show that given the span of covariance matrices $\mathcal{U}$ obtained from Step 1, finding the unfolded moments $\widetilde{X}_{4}$, $\widetilde{X}_{6}$ is reduced to solving two systems of linear equations.

Note that once we know $\widetilde{Y}_{4}$ and $\widetilde{Y}_{6}$, the unfolded moments $\widetilde{X}_{4}$ and $\widetilde{X}_{6}$ are given by $\widetilde{X}_{4}=U\widetilde{Y}_{4}U^{\top}$ and $\widetilde{X}_{6}=\widetilde{Y}_{6}(U^{\top},U^{\top},U^{\top})$. Therefore, after changing the variable, we need to solve the two linear equation systems given in (7) with the variables $\widetilde{Y}_{4}$ and $\widetilde{Y}_{6}$.

This change of variable significantly reduces the number of unknown variables. Note that the number of distinct entries in $\widetilde{Y}_{4}$ and $\widetilde{Y}_{6}$ are $k_{2}={k+1\choose 2}$ and $k_{3}={k+2\choose 3}$, respectively. Since $k_{2}\leq n_{4}$ and $k_{3}\leq n_{6}$, we can expect that the linear mapping from $\widetilde{Y}_{4}$ to $\overline{\widetilde{M}}_{4}$ and the one from $\widetilde{Y}_{6}$ to $\overline{\widetilde{M}}_{6}$ are linearly invertible. This argument is formalized below.

We show with smoothed analysis that the smallest singular value of the two coefficient matrices are lower bounded with high probability:

With high probability over the parameter random perturbation, the $k_{2}$-th singular value of the coefficient matrix $\widetilde{H}_{4}$ is at least $\Omega(\rho^{2}n/k)$, and the $k_{3}$-th singular value of the coefficient matrix $\widetilde{H}_{6}$ is at least $\Omega(\rho^{3}(n/k)^{1.5})$.

To prove this lemma we rewrite the coefficient matrix as product of two matrices and bound their smallest singular values separately. One of the two matrices corresponds to a projection of the Kronecker product $\widetilde{\Sigma}\otimes_{kr}\widetilde{\Sigma}$. In the smoothed analysis setting, this matrix is not necessarily incoherent. In order to provide a lower bound to its smallest singular value, we further apply a carefully designed projection to it, and then we use the concentration bounds for Gaussian chaoses to show that after the projection its columns are incoherent, finally we apply Gershgorin’s Theorem to bound the smallest singular value Note that the idea of unfolding using system of linear equations also appeared in the work of Jain and Oh (2014). However, in order to show the system of linear equations in their setup is robust, i.e., the coefficient matrix has full rank, they heavily rely on the incoherence assumption, which we do not impose in the smoothed analysis setting. .

When implementing this step with the empirical moments, we solve two least squares problems instead of solving the system of linear equations. Again using results in matrix perturbation theory and using the lower bound of the smallest singular values of the two coefficient matrices, we show the stability of the solution to the least squares problems:

Given the empirical moments $\widehat{M}_{4}=\widetilde{M}_{4}+E_{4}$, $\widehat{M}_{6}=\widetilde{M}_{6}+E_{6}$, and suppose that the absolute value of entries of $E_{4}$ and $E_{6}$ are at most $\delta_{1}$. Let $\widehat{U}$, the output of Step 1, be the estimation for the span of the covariance matrices, and suppose that $\|\widehat{U}-\widetilde{U}\|\leq\delta_{2}$. Let $\widehat{Y}_{4}$ and $\widehat{Y}_{6}$ be the least squares solution respectively. When $\delta_{1}$ and $\delta_{2}$ are inverse polynomially small, we have $\|\widetilde{Y}_{4}-\widehat{Y}_{4}\|_{F}\leq O(\sqrt{n_{4}}(\delta_{1}+\delta_{2}/\sigma_{min}(\widetilde{H}_{4}))$ and $\|\widetilde{Y}_{6}-\widehat{Y}_{6}\|_{F}\leq O(\sqrt{n_{6}}(\delta_{1}+\delta_{2}/\sigma_{min}(\widetilde{H}_{6}))$.

Step 3. Tensor Decomposition.

Given $\widetilde{Y}_{4}$, $\widetilde{Y}_{6}$ and $\widetilde{U}$, the symmetric tensor decomposition algorithm can correctly and robustly find the mixing weights $\widetilde{\omega}_{i}$’s and the vectors $\widetilde{\sigma}_{i}$’s, up to some unknown permutation over $[k]$, with high probability over both the randomized algorithm and the parameter perturbation.

The algorithm and its analysis mostly follow the algorithm of symmetric tensor decomposition in Anandkumar et al. (2014), and the details are provided in Section D in the appendix.

Proof Sketch for the Main Theorem of Zero-mean Case.

Theorem 3.5 follows from the previous smoothed analysis and stability analysis lemmas for each step.

First, exploiting the randomness of parameter perturbation, the smoothed analysis lemmas show that the deterministic conditions, which guarantee the correctness of each step, are satisfied with high probability. Then using concentration bounds of Gaussian variables, we show that with high probability over the random samples, the empirical moments $\widehat{M}_{4}$ and $\widehat{M}_{6}$ are entrywise $\delta$-close to the exact moments $\widetilde{M}_{4}$ and $\widetilde{M}_{6}$. In order to achieve $\epsilon$ accuracy in the parameter estimation, we choose $\delta$ to be inverse polynomially small, and therefore the number of samples required will be polynomial in the relevant parameters. The stability lemmas show how the errors propagate only “polynomially” through the steps of the algorithm, which is visualized in Figure 1.

A more detailed illustration is provided in Section E in the appendix.

Algorithm Outline for Learning Mixture of General Gaussians

In this section, we briefly discuss the algorithm for learning mixture of general Gaussians. Figure 2 shows the inputs and outputs of each step in this algorithm. Many steps share similar ideas to those of the algorithm for the zero-mean case in previous sections. We only highlight the basic ideas and defer the details to Section F in the appendix.

Similar to Step 1 in the zero-mean case, this step makes use of the structure of the 4-th order moments $\widetilde{M}_{4}$, and is achieved in three small steps:

For a subset $\mathcal{H}\subset[n]$ of size $|\mathcal{H}|=\sqrt{n}$, find the span:

Find the span of the covariance matrices with the columns projected onto $\mathcal{S}^{\perp}$, namely,

For disjoint subsets $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, repeat Step 1 (a) and Step 1 (b) to obtain $\mathcal{U}_{1}$ and $\mathcal{U}_{2}$, the span of the covariance matrices projected onto the subspaces $\mathcal{S}_{1}^{\perp}$ and $\mathcal{S}_{2}^{\perp}$. The intersection of the two subspaces $\mathcal{U}_{1}$ and $\mathcal{U}_{2}$ gives the span of the mean vectors $\widetilde{Z}=\text{span}\left\{\widetilde{\mu}^{(i)},i\in[k]\right\}$. Merge the two subspaces $\mathcal{U}_{1}$ and $\mathcal{U}_{2}$ to obtain the span of the covariance matrices projected to the subspace orthogonal to $\widetilde{Z}$, namely $\widetilde{\Sigma}_{o}=\text{span}\left\{\text{Proj}_{\widetilde{Z}^{\perp}}\widetilde{\Sigma}^{(i)}\text{Proj}_{\widetilde{Z}^{\perp}}:i\in[k]\right\}$.

The key observation of this step is that when the samples are projected to the subspace orthogonal to all the mean vectors, they are equivalent to samples from a mixture of zero-mean Gaussians with covariance matrices $\widetilde{\Sigma}^{(i)}_{o}=\text{Proj}_{\widetilde{Z}^{\perp}}\widetilde{\Sigma}^{(i)}\text{Proj}_{\widetilde{Z}^{\perp}}$ and with the same mixing weights $\widetilde{\omega}_{i}$’s. Therefore, projecting the samples to $\widetilde{Z}^{\perp}$, the subspace orthogonal to the mean vectors, and use the algorithm for the zero-mean case, we can obtain $\widetilde{\Sigma}^{(i)}_{o}$’s, the covariance matrices projected to this subspace, as well as the mixing weights $\widetilde{\omega}_{i}$’s.

Step 3. Find the means

With simple algebra, this step extracts the projected covariance matrices $\widetilde{\Sigma}_{o}^{(i)}$’s from the $3$-rd order moments $\widetilde{M}_{3}$, the mixing weights $\widetilde{\omega}_{i}$ and the projected covariance matrices $\widetilde{\Sigma}_{o}^{(i)}$’s obtained in Step 2.

Step 4. Find the full covariance matrices

In Step 2, we obtained $\widetilde{\Sigma}^{(i)}_{o}$, the covariance matrices projected to the subspace orthogonal to all the means. Note that they are equal to matrices $(\widetilde{\Sigma}^{(i)}+\widetilde{\mu}^{(i)}(\widetilde{\mu}^{(i)})^{\top})$ projected to the same subspace. We claim that if we can find the span of these matrices ($(\widetilde{\Sigma}^{(i)}+\widetilde{\mu}^{(i)}(\widetilde{\mu}^{(i)})^{\top})$’s), we can get each matrix $(\widetilde{\Sigma}^{(i)}+\widetilde{\mu}^{(i)}(\widetilde{\mu}^{(i)})^{\top})$, and then subtracting the known rank-one component to find the covariance matrix $\widetilde{\Sigma}^{(i)}$. This is similar to the idea of merging two projections of the same subspace in Step 1 (c) for the zero-mean case.

The idea of finding the desired span is to construct a $4$-th order tensor:

which corresponds to the 4-th order moments of a mixture of zero-mean Gaussians with covariance matrices $\widetilde{\Sigma}^{(i)}+\widetilde{\mu}^{(i)}(\widetilde{\mu}^{(i)})^{\top}$ and the same mixing weights $\widetilde{\omega}_{i}$’s. Then we can then use Step 1 of the algorithm for the zero-mean case to obtain the span of the new covariance matrices, i.e. $span\{\widetilde{\Sigma}^{(i)}+\widetilde{\mu}^{(i)}(\widetilde{\mu}^{(i)})^{\top}:i\in[k]\}$.

Conclusion

In this paper we give the first efficient algorithm for learning mixture of general Gaussians in the smoothed analysis setting. In the algorithm we developed new ways of extracting information from lower-order moment structure. This suggests that although the method of moments often involves solving systems of polynomial equations that are intractable in general, for natural models there is still hope of utilizing their special structure to obtain algebraic solution.

Smoothed analysis is a very useful way of avoiding degenerate examples in analyzing algorithms. In the analysis, we proved several new results for bounding the smallest singular values of structured random matrices. We believe the lemmas and techniques can be useful in more general settings.

Our algorithm uses only up to $6$-th order moments. We conjecture that using higher order moments can reduce the number of dimension required to $n\geq\Omega(k^{1+\epsilon})$, or maybe even $n\geq\Omega(k^{\epsilon})$.

We thank Santosh Vempala for many insights and for help in earlier attempts at solving this problem.

References

Appendix A Moment Structures

In this section we characterize the structure of the 3-rd, 4-th and 6-th moments of Gaussians mixtures.

where $y^{(i)}$ corresponds to the $n$-dimensional Gaussian distribution $\mathcal{N}(\mu^{(i)},\Sigma^{(i)})$.

Gaussian distribution is a highly symmetric distribution, and in the zero-mean case the higher moments are well-understood by Isserlis’ Theorem:

Let $\mathbf{y}=(y_{1},\dots,y_{2t})$ be a multivariate Gaussian random vector with mean zero and covariance $\Sigma$, then

where the summation is taken over all distinct ways of partitioning $y_{1},\dots,y_{2t}$ into $t$ pairs, which correspond to all the perfect matchings in a complete graph. Thus there are $(2t-1)!!$ terms in the sum, and each summand is a product of $t$ terms.

The non-zero mean case is a direct corollary using Isserlis’ Theorem and linearity of expectation.

Let $\mathbf{y}=(y_{1},\dots,y_{t})$ be a multivariate Gaussian random vector with mean $\mu$ and covariance $\Sigma$, then

where the summation is taken over all distinct ways of partitioning $y_{1},\dots,y_{t}$ into $p$ pairs of $(u,v)$ and $s$ singletons of $(w)$, where $p\geq 0$, $s\geq 0$ and $2p+s=t$.

We shall first prove Lemma 3.7 in Section 3.2. Recall that this lemma shows that for mixture of zero-mean Gaussians, the 4-th moments $\overline{M}_{4}$ and the 6-th moments $\overline{M}_{6}$ with distinct indices can be viewed as a linear projection of the unfolded moment $X_{4}$ and $X_{6}$ defined in (1).

By Isserlis Theorem A.1, the mapping $\sqrt{3}\mathcal{F}_{4}$ is characterized by: ($\forall 1\leq j_{1}<j_{2}<j_{3}<j_{4}\leq n$)

Therefore, with the normalization constant $\sqrt{3}$, the $(j_{1},j_{2},j_{3},j_{4})$-th mapping of $\mathcal{F}_{4}$ is a projection of the three elements in $X_{4}$. Similarly, we have for $\sqrt{15}\mathcal{F}_{6}$: ($\forall 1\leq j_{1}<j_{2}<\dots<j_{6}\leq n$)

Thus with the normalization constant $\sqrt{15}$, the mapping $\mathcal{F}_{6}$ is a linear projection. ∎

A.2 Slices of Moments

Next we shall characterize the slices of the moments of mixture of Gaussians.

For mixture of zero-mean Gaussians, a one-dimensional slice of the 4th moment tensor is a vector in the span of corresponding columns of the covariance matrices:

For a mixture of zero-mean Gaussians, the one-dimensional slices of the 4-th moments $M_{4}$ are given by:

By the definition of multilinear map, $M_{4}(e_{j_{1}},e_{j_{2}},e_{j_{3}},I)$ is a vector whose $p$-th entry is equal to $M_{4}(e_{j_{1}},e_{j_{2}},e_{j_{3}},e_{p})$. We can compute this entry by Isserlis’ Theorem:

For mixture of zero-mean Gaussians, a two-dimensional slice of the 4th moment $M_{4}$ is a matrix, and it is a linear combination of the covariance matrices with some additive rank one matrices:

For a mixture of zero-mean Gaussians, the two-dimensional slices of the 4-th moment $M_{4}$ are given by:

Again this follows from Isserlis’ theorem. By definition of multilinear map this is a matrix whose $(p,q)$-th entry is equal to

Similarly, for mixture of general Gaussians, we prove the following claims:

For a mixture of general Gaussians, the $(j_{1},j_{2},j_{3})$-th one-dimensional slice of $M_{4}$ is given by:

This is very similar to Claim 5.1 and follows from the corollary of Isserlis’s theorem (Corollary A.2). There are 10 ways to partition the indices $\{j_{1},j_{2},j_{3},j_{4}\}$ into pairs and singletons: $((j_{1}),(j_{2}),(j_{3}),(j_{4}))$, $((j_{1},j_{2}),(j_{3}),(j_{4}))$, $((j_{1},j_{3}),(j_{2}),(j_{4}))$, $((j_{1},j_{4}),(j_{2}),(j_{3}))$, $((j_{2},j_{3}),(j_{1}),(j_{4}))$, $((j_{2},j_{4}),(j_{1}),(j_{3}))$, $((j_{3},j_{4}),(j_{1}),(j_{2}))$, $((j_{1},j_{2}),(j_{3},j_{4}))$, $((j_{1},j_{3}),(j_{2},j_{4}))$, $((j_{1},j_{4}),(j_{2},j_{3}))$. From this enumeration, we can specify each element in the vector of the one-dimensional slice. ∎

Since $M_{4}$ gives linear observations on the symmetric low rank matrix $X_{4}$, it is natural to wonder whether we can use matrix completion techniques to recover $X_{4}$ from $M_{4}$. Here we show this is impossible by giving a counter example: there are two mixture of Gaussians that generates the same 4th moment $M_{4}$, but has different $X_{4}$ (even the span of $\Sigma^{(i)}$’s are different).

By $((a,b),(c,d))$ we denote a $5\times 5$ matrix $A$ which has $2$’s on diagonals, and the only nonzero off-diagonal entries are $A_{a,b}=A_{b,a}=A_{c,d}=A_{d,c}=1$. For example, $((1,2),(4,5))$ will be the following matrix:

where all the missing entries are 0’s. Now we construct two mixtures of 3 Gaussians, all with mean 0 and weight $1/3$. The covariance matrices are $((1,2),(4,5)),((1,3),(2,5)),((1,4),(3,5))$ for the first mixture and $((1,2),(3,5)),((1,3),(4,5)),((1,4),(2,5))$ for the second mixture. These are clearly different mixtures with different span of $\Sigma^{(i)}$’s: in the first mixture, $\Sigma^{(i)}_{1,2}=\Sigma^{(i)}_{4,5}$ for all matrices, but this is not true for the second mixture.

These two mixture of Gaussians have the same 4th moment $M_{4}$. This can be checked by using Isserlis’ theorem to compute the moments. Intuitively, this is true because all the pairs $(1,i)$ and $(i,5)$ appeared exactly twice in the covariance matrices for both mixtures; also, every 4-tuple $(1,i,j,5)$ appeared exactly once in the covariance matrices for both mixtures.

Appendix B Step 1: Span Finding

Recall that in Step 1 of the algorithm for learning mixture of zero-mean Gaussians, we find the span of the covariance matrices in three small steps. In this section, we prove the correctness and the robustness of each step with smoothed analysis.

For completeness we restate each substep and highlight the key properties we need, followed by the detailed proofs.

In Step 1 (a), for any set $\mathcal{H}$ of size $\sqrt{n}$, we want to show that the one-dimensional slices of $M_{4}$ span the entire subspace $\mathcal{S}=\text{span}\left\{\widetilde{\Sigma}^{(i)}_{[:,j]}:i\in[k],j\in\mathcal{H}\right\}$, which is the span of a subset of the columns in the covariance matrices.

This in particular means when $j_{1},j_{2},j_{3}\in\mathcal{H}$, the vector $\widetilde{M}_{4}(e_{j_{1}},e_{j_{2}},e_{j_{3}},I)$ is in $\mathcal{S}$. We need to show that the columns of the matrix $Q$ indeed span the entire subspace $\mathcal{S}$.

It is sufficient to show that a subset of the column span the entire subspace. Form a three-way even partition of the set $\mathcal{H}$, i.e., $|\mathcal{H}^{(1)}|=|\mathcal{H}^{(2)}|=|\mathcal{H}^{(3)}|=|\mathcal{H}|/3=\sqrt{n}/3$, and only consider the one-dimensional slices of $\widetilde{M}_{4}$ corresponding to the indices $j_{i}\in\mathcal{H}^{(i)}$ for $i=1,2,3$. In particular, we define matrix $\widetilde{Q}_{S}$ with these one-dimensional slices of $\widetilde{M}_{4}$:

Define matrix $\widetilde{P}_{S}$ with the corresponding columns of the covariance matrices, forming a basis (although not orthogonal) of the desired subspace $\mathcal{S}$:

In the following two lemmas, we show that with high probability over the random perturbation, the column span of $\widetilde{Q}_{S}$ is exactly equal to the column span of $\widetilde{P}_{S}$, and robustly so.

Given $\widetilde{M}_{4}$, the exact 4-th order moment of the $\rho$-smooth mixture of zero-mean Gaussians, for any subset $\mathcal{H}\in[n]$ with cardinality $|\mathcal{H}|=\sqrt{n}$, let $\widetilde{Q}_{S}$ be the matrix defined as in (12) with the one-dimensional slices of $\widetilde{M}_{4}$. For any $\epsilon>0$, and for some absolute constant $C_{1},C_{2},C_{3}>0$, with probability at least $1-(C_{1}\epsilon)^{C_{2}n}$, the $k|\mathcal{H}|$-th singular value of $\widetilde{Q}_{S}$ is bounded below by:

In order to prove this lemma, we first write $\widetilde{Q}_{S}$ as the product of three matrices.

Under the same assumptions of Lemma B.1, the matrix $\widetilde{Q}_{S}$ can be written as

With some fixed and known row permutation $\pi^{(2)}$ and $\pi^{(3)}$, the matrix $\widetilde{B}^{(2)}$ and $\widetilde{B}^{(3)}$ can be made block diagonal with identical blocks equal to $\widetilde{\Sigma}_{\mathcal{H}^{(1)},\mathcal{H}^{(3)}}$ and $\widetilde{\Sigma}_{\mathcal{H}^{(1)},\mathcal{H}^{(2)}}$, respectively. Note that the three parts $\widetilde{B}^{(1)},\widetilde{B}^{(2)},\widetilde{B}^{(3)}$ do not have any common entry, nor do they involve any diagonal entry of the covariance matrices, therefore the three parts are independent when the covariances are randomly perturbed in the smoothed analysis.

It is easier to understand the structure by picture, see Figure 3. The rows of the matrix should be indexed by $(j_{1},j_{2},j_{3})\in\mathcal{H}^{(1)}\times\mathcal{H}^{(2)}\times\mathcal{H}^{(3)}$, which can also be interpreted as a cube (in the right). The block structure in the first part $\widetilde{B}^{(1)}$ correspond to a slice in $\mathcal{H}^{(2)}\times\mathcal{H}^{(3)}$ direction (for each block, the element in $\mathcal{H}^{(1)}$ is fixed, the elements in $\mathcal{H}^{(2)}$ and $\mathcal{H}^{(3)}$ take all possible values). Similarly for $\widetilde{B}^{(2)}$ and $\widetilde{B}^{(3)}$ (as shown in figure).

(of Claim B.2 ) The proof of this claim is using Claim 5.1, the definition of matrices and the rule of matrix multiplication. Consider the column in $\widetilde{Q}_{S}$ corresponding to the index $(j_{1},j_{2},j_{3})$ for $j_{1}\in\mathcal{H}^{(1)},j_{2}\in\mathcal{H}^{(2)},j_{3}\in\mathcal{H}^{(3)}$, and the row of $\widetilde{B}_{S}$ together with the mixing wights specifies how this column is formed as a linear combination of $3k$ columns of $\widetilde{P}_{S}$. By the structure of $M_{4}$ in Claim 5.1, the $(j_{1},j_{2},j_{3})$-th row of $\widetilde{B}^{(1)}$ has exactly $k$ entries corresponding to $\widetilde{\Sigma}^{(i)}_{j_{2},j_{3}}$ for $i\in[k]$, these entries are multiplied by $\widetilde{\omega}_{i}$ in the middle (diagonal) matrix. Therefore, these directly correspond to the $k$ terms in Claim 5.1. Similarly the entries in $\widetilde{B}^{(2)}$ and $\widetilde{B}^{(3)}$ correspond to the other $2k$ terms. ∎

Using Claim B.2, we need to bound the smallest singular value for each of the matrices in order to bound the $k|\mathcal{H}|$-th singular value of $\widetilde{Q}_{S}$, this is deferred to the end of this part. The most important tool is a corollary (Lemma G.16) of the random matrix result proved in Rudelson and Vershynin (2009), which gives a lowerbound on the smallest singular value of perturbed rectangular matrices.

By Lemma B.1, we know $\widetilde{Q}_{S}$ has exactly rank $k|\mathcal{H}|$, and is robust in the sense that its $k|\mathcal{H}|$-th singular value is large (polynomial in the amount of perturbation $\rho$). By standard matrix perturbation theory, if we get $\widehat{Q}_{S}$ close to $\widetilde{Q}_{S}$ up to a high accuracy (inverse polynomial in the relevant parameters), the top $k|\mathcal{H}|$ singular vectors will span a subspace that is very close to the span of $\widetilde{Q}_{S}$. We formalize this in the following lemma.

Note that the columns of $S$ are the leading left singular vectors of $Q_{S}$. We apply the standard matrix perturbation bound of singular vectors. Recall that $S$ is defined to be the first $k|\mathcal{H}|$ left singular vector of $Q_{S}$, and we have

Therefore by Wedin’s Theorem (in particular the corollary Lemma G.5), we can conclude (16). ∎

We first use Claim B.2 to write $\widetilde{Q}_{S}=\widetilde{P}_{S}\left(D_{\widetilde{\omega}}\otimes_{kr}I_{|\mathcal{H}|}\right)(\widetilde{B}_{S})^{\top}$, note that the matrix $(D_{\widetilde{\omega}}\otimes_{kr}I_{|\mathcal{H}|})$ has dimension $k|\mathcal{H}|\times k|\mathcal{H}|$, therefore we just need to show with high probability each of the three factor matrix has large $k|\mathcal{H}|$-th singular value, and that implies a bound on the $k|\mathcal{H}|$-th singular value of $\widetilde{Q}_{S}$ by union bound. The smallest singular value of $\widetilde{P}_{S}$ and $\widetilde{B}_{S}$ are bounded below by the following two Claims.

With high probability $\sigma_{k|\mathcal{H}|}(\widetilde{P}_{S})\geq\Omega(\rho\sqrt{n})$.

Now $\widetilde{P}^{\prime}_{S}$ is a randomly perturbed rectangular matrix, whose smallest singular value can be lower bounded using Lemma G.16, and we conclude that with probability at least $1-(C\epsilon)^{0.25n}$,

Next, we bound the smallest singular value of $\widetilde{B}_{S}$.

With high probability $\sigma_{k|\mathcal{H}|}(\widetilde{B}_{S})\geq\Omega(\rho\sqrt{n})$.

We make use of the special structure of the three blocks of $\widetilde{B}_{S}$ to lower bound its smallest singular value.

First, we prove that the block diagonal matrix $\widetilde{B}^{(1)}$ has large singular values, even after projecting to the orthogonal subspace of the column span of $\widetilde{B}^{(2)}$ and $\widetilde{B}^{(3)}$. This idea appeared several times in our proof and is abstracted in Lemma G.12. Apply the lemma and we have:

where the $j$-th block of $[\widetilde{B}^{(2)},\widetilde{B}^{(3)}]$ has dimension $(|\mathcal{H}|/3)^{2}\times 2k|\mathcal{H}|/3$. Since

this means for each block, even after projection it has more than $3k$ rows. Note that by definition the three blocks $\widetilde{B}^{(1)}$, $\widetilde{B}^{(2)}$ and $\widetilde{B}^{(3)}$ are independent and do not involve any diagonal elements of the covariance matrices, so each block after the two projections is again a rectangular random matrix. We can apply Lemma G.15, for any $j$, for some absolute constant $C_{1},C_{2},C_{3}$ (not fixed throughout the discussion), with probability at least $1-(C_{1}\epsilon)^{C_{2}n}$ over the randomness of $\widetilde{\Sigma}_{\mathcal{H}^{(2)},\mathcal{H}^{(3)}}$, we have:

Now we can take a union bound over the blocks and conclude that with high probability, the smallest singular value of each block is large.

In order to bound $\sigma_{k(2|\mathcal{H}|/3)}([\widetilde{B}^{(2)},\widetilde{B}^{(3)}])$, we use the same strategy. Note that $\widetilde{B}^{(2)}$ also has a block structure that corresponds to the $\mathcal{H}^{(1)}\times\mathcal{H}^{(3)}$ faces (see Figure 3). Again check the condition on dimension $(|\mathcal{H}|/3)^{2}-k-k|\mathcal{H}|/3\geq\Omega(n)>3k$, we can apply Lemma G.12 again to show that for any $j$, with probability at least $1-(C_{1}\epsilon)^{C_{2}n}$ over the randomness of $\widetilde{\Sigma}_{\mathcal{H}^{(1)},\mathcal{H}^{(3)}}$, we have: