A Method of Moments for Mixture Models and Hidden Markov Models

Animashree Anandkumar, Daniel Hsu, Sham M. Kakade

Introduction

Mixture models are a fundamental tool in applied statistics and machine learning for treating data taken from multiple subpopulations (Titterington et al., 1985). In a mixture model, the data are generated from a number of possible sources, and it is of interest to identify the nature of the individual sources. As such, estimating the unknown parameters of the mixture model from sampled data—especially the parameters of the underlying constituent distributions—is an important statistical task. For most mixture models, including the widely used mixtures of Gaussians and hidden Markov models (HMMs), the current practice relies on the Expectation-Maximization (EM) algorithm, a local search heuristic for maximum likelihood estimation. However, EM has a number of well-documented drawbacks regularly faced by practitioners, including slow convergence and suboptimal local optima (Redner and Walker, 1984).

An alternative to maximum likelihood and EM, especially in the context of mixture models, is the method of moments approach. The method of moments dates back to the origins of mixture models with Pearson’s solution for identifying the parameters of a mixture of two univariate Gaussians (Pearson, 1894). In this approach, model parameters are chosen to specify a distribution whose pp-th order moments, for several values of pp, are equal to the corresponding empirical moments observed in the data. Since Pearson’s work, the method of moments has been studied and adapted for a variety of problems; their intuitive appeal is also complemented with a guarantee of statistical consistency under mild conditions. Unfortunately, the method often runs into trouble with large mixtures of high-dimensional distributions. This is because the equations determining the parameters are typically based on moments of order equal to the number of model parameters, and high-order moments are exceedingly difficult to estimate accurately due to their large variance.

This work develops a computationally efficient method of moments based on only low-order moments that can be used to estimate the parameters of a broad class of high-dimensional mixture models with many components. The resulting estimators can be implemented with standard numerical linear algebra routines (singular value and eigenvalue decompositions), and the estimates have low variance because they only involve low-order moments. The class of models covered by the method includes certain multivariate Gaussian mixture models and HMMs, as well as mixture models with no explicit likelihood equations. The method exploits the availability of multiple indirect “views” of a model’s underlying latent variable that determines the source distribution, although the notion of a “view” is rather general. For instance, in an HMM, the past, present, and future observations can be thought of as different noisy views of the present hidden state; in a mixture of product distributions (such as axis-aligned Gaussians), the coordinates in the output space can be partitioned (say, randomly) into multiple non-redundant “views”. The new method of moments leads to unsupervised learning guarantees for mixture models under mild rank conditions that were not achieved by previous works; in particular, the sample complexity of accurate parameter estimation is shown to be polynomial in the number of mixture components and other relevant quantities. Finally, due to its simplicity, the new method (or variants thereof) also offers a viable alternative to EM and maximum likelihood for practical deployment.

Gaussian mixture models. The statistical literature on mixture models is vast (a more thorough treatment can be found in the texts of Titterington et al. (1985) and Lindsay (1995)), and many advances have been made in computer science and machine learning over the past decade or so, in part due to their importance in modern applications. The use of mixture models for clustering data comprises a large part of this work, beginning with the work of Dasgupta (1999) on learning mixtures of kk well-separated dd-dimensional Gaussians. This and subsequent work (Arora and Kannan, 2001; Dasgupta and Schulman, 2007; Vempala and Wang, 2002; Kannan et al., 2005; Achlioptas and McSherry, 2005; Chaudhuri and Rao, 2008; Brubaker and Vempala, 2008; Chaudhuri et al., 2009) have focused on efficient algorithms that provably recover the parameters of the constituent Gaussians from data generated by such a mixture distribution, provided that the distance between each pair of means is sufficiently large (roughly either dcd^{c} or kck^{c} times the standard deviation of the Gaussians, for some c>0c>0). Such separation conditions are natural to expect in many clustering applications, and a number of spectral projection techniques have been shown to enhance the separation (Vempala and Wang, 2002; Kannan et al., 2005; Brubaker and Vempala, 2008; Chaudhuri et al., 2009). More recently, techniques have been developed for learning mixtures of Gaussians without any separation condition (Kalai et al., 2010; Belkin and Sinha, 2010; Moitra and Valiant, 2010), although the computational and sample complexities of these methods grow exponentially with the number of mixture components kk. This dependence has also been shown to be inevitable without further assumptions (Moitra and Valiant, 2010).

Method of moments. The latter works of Belkin and Sinha (2010), Kalai et al. (2010), and Moitra and Valiant (2010) (as well as the algorithms of Feldman et al. (2005, 2006) for a related but different learning objective) can be thought of as modern implementations of the method of moments, and their exponential dependence on kk is not surprising given the literature on other moment methods for mixture models. In particular, a number of moment methods for both discrete and continuous mixture models have been developed using techniques such as the Vandermonde decompositions of Hankel matrices (Lindsay, 1989; Lindsay and Basak, 1993; Boley et al., 1997; Gravin et al., 2012). In these methods, following the spirit of Pearson’s original solution, the model parameters are derived from the roots of polynomials whose coefficients are based on moments up to the Ω(k)\Omega(k)-th order. The accurate estimation of such moments generally has computational and sample complexity exponential in kk.

Spectral approach to parameter estimation with low-order moments. The present work is based on a notable exception to the above situation, namely Chang’s spectral decomposition technique for discrete Markov models of evolution (Chang, 1996) (see also Mossel and Roch (2006) and Hsu et al. (2009) for adaptations to other discrete mixture models such as discrete HMMs). This spectral technique depends only on moments up to the third-order; consequently, the resulting algorithms have computational and sample complexity that scales only polynomially in the number of mixture components kk. The success of the technique depends on a certain rank condition of the transition matrices; but this condition is much milder than separation conditions of clustering works, and it remains sufficient even when the dimension of the observation space is very large (Hsu et al., 2009). In this work, we extend Chang’s spectral technique to develop a general method of moments approach to parameter estimation, which is applicable to a large class of mixture models and HMMs with both discrete and continuous component distributions in high-dimensional spaces. Like the moment methods of Moitra and Valiant (2010) and Belkin and Sinha (2010), our algorithm does not require a separation condition; but unlike those previous methods, the algorithm has computational and sample complexity polynomial in kk.

Some previous spectral approaches for related learning problems only use second-order moments, but these approaches can only estimate a subspace containing the parameter vectors and not the parameters themselves (McSherry, 2001). Indeed, it is known that the parameters of even very simple discrete mixture models are not generally identifiable from only second-order moments (Chang, 1996)See Appendix G for an example of Chang (1996) demonstrating the non-identifiability of parameters from only second-order moments in a simple class of Markov models.. We note that moments beyond the second-order (specifically, fourth-order moments) have been exploited in the methods of Frieze et al. (1996) and Nguyen and Regev (2009) for the problem of learning a parallelepiped from random samples, and that these methods are very related to techniques used for independent component analysis (Hyvärinen and Oja, 2000). Adapting these techniques for other parameter estimation problems is an enticing possibility.

Multi-view learning. The spectral technique we employ depends on the availability of multiple views, and such a multi-view assumption has been exploited in previous works on learning mixtures of well-separated distributions (Chaudhuri and Rao, 2008; Chaudhuri et al., 2009). In these previous works, a projection based on a canonical correlation analysis (Hotelling, 1935) between two views is used to reinforce the separation between the mixture components, and to cancel out noise orthogonal to the separation directions. The present work, which uses similar correlation-based projections, shows that the availability of a third view of the data can remove the separation condition entirely. The multi-view assumption substantially generalizes the case where the component distributions are product distributions (such as axis-aligned Gaussians), which has been previously studied in the literature (Dasgupta, 1999; Vempala and Wang, 2002; Chaudhuri and Rao, 2008; Feldman et al., 2005, 2006); the combination of this and a non-degeneracy assumption is what allows us to avoid the sample complexity lower bound of Moitra and Valiant (2010) for Gaussian mixture models. The multi-view assumption also naturally arises in many applications, such as in multimedia data with (say) text, audio, and video components (Blaschko and Lampert, 2008; Chaudhuri et al., 2009); as well as in linguistic data, where the different words in a sentence or paragraph are considered noisy predictors of the underlying semantics (Gale et al., 1992). In the vein of this latter example, we consider estimation in a simple bag-of-words document topic model as a warm-up to our general method; even this simpler model illustrates the power of pair-wise and triple-wise (i.e., bigram and trigram) statistics that were not exploited by previous works on multi-view learning.

2 Outline

Section 2 first develops the method of moments in the context of a simple discrete mixture model motivated by document topic modeling; an explicit algorithm and convergence analysis are also provided. The general setting is considered in Section 3, where the main algorithm and its accompanying correctness and efficiency guarantee are presented. Applications to learning multi-view mixtures of Gaussians and HMMs are discussed in Section 4. All proofs are given in the appendix.

3 Notations

Warm-up: bag-of-words document topic modeling

We first describe our method of moments in the simpler context of bag-of-words models for documents.

The generative process for a document is given as follows:

The document’s topic is drawn according to the multinomial distribution specified by the probability vector w=(w1,w2,,wk)Δk1\vec{w}=(w_{1},w_{2},\dotsc,w_{k})\in\Delta^{k-1}. This is modeled as a discrete random variable hh such that

This probabilistic model has the conditional independence structure depicted in Figure 2(a) as a directed graphical model.

We assume the following condition on w\vec{w} and MM.

wj>0w_{j}\kern-2.0pt>\kern-2.0pt0 for all j[k]j\kern-2.0pt\in\kern-2.0pt[k], and MM has rank kk.

This condition requires that each topic has non-zero probability, and also prevents any topic’s word distribution from being a mixture of the other topics’ word distributions.

2 Pair-wise and triple-wise probabilities

Since x1\vec{x}_{1}, x2\vec{x}_{2}, and x3\vec{x}_{3} are conditionally independent given hh,

3 Observable operators and their spectral properties

The pair-wise and triple-wise probabilities can be related in a way that essentially reveals the conditional probability matrix MM. This is achieved through a matrix called an “observable operator”. Similar observable operators were previously used to characterize multiplicity automata (Schützenberger, 1961; Jaeger, 2000) and, more recently, for learning discrete HMMs (via an operator parameterization) (Hsu et al., 2009).

The matrix B(η)B(\vec{\eta}) is called “observable” because it is only a function of the observable variables’ joint probabilities (e.g., Pr[x1=ei,x2=ej]\Pr[\vec{x}_{1}=\vec{e}_{i},\vec{x}_{2}=\vec{e}_{j}]). In the case η=ex\vec{\eta}=\vec{e}_{x} for some x[d]x\in[d], the matrix B(ex)B(\vec{e}_{x}) is similar (in the linear algebraic sense) to the diagonal matrix diag(Mex)\operatorname{diag}(M^{\scriptscriptstyle\top}\vec{e}_{x}); the collection of matrices {diag(Mex):x[d]}\{\operatorname{diag}(M^{\scriptscriptstyle\top}\vec{e}_{x}):x\in[d]\} (together with w\vec{w}) can be used to compute joint probabilities under the model (see, e.g., Hsu et al. (2009)). Note that the columns of UMU^{\scriptscriptstyle\top}M are eigenvectors of B(ex)B(\vec{e}_{x}), with the jj-th column having an associated eigenvalue equal to Pr[xv=xh=j]\Pr[\vec{x}_{v}=x|h=j]. If the word xx has distinct probabilities under every topic, then B(ex)B(\vec{e}_{x}) has exactly kk distinct eigenvalues, each having geometric multiplicity one and corresponding to a column of UMU^{\scriptscriptstyle\top}M.

4 Topic-word distribution estimator and convergence guarantee

The spectral properties of the observable operators B(η)B(\vec{\eta}) implied by Lemma 2.2 suggest the estimation procedure (Algorithm A) in Figure 1. The procedure is essentially a plug-in approach based on the equations relating the second- and third-order moments in Lemma 2.2. We focus on estimating MM; estimating the mixing weights w\vec{w} is easily handled as a secondary step (see Appendix B.5 for the estimator in the context of the general model in Section 3.1).

The following theorem establishes the convergence rate of Algorithm A.

The proof of Theorem 2.1, as well as some illustrative empirical results on using Algorithm A, are presented in Appendix A. A few remarks about the theorem are in order.

On boosting the confidence. Although the convergence depends polynomially on 1/δ1/\delta, where δ\delta is the failure probability, it is possible to boost the confidence by repeating Step 3 of Algorithm A with different random η\vec{\eta} until the eigenvalues of B^(η)\hat{B}(\vec{\eta}) are sufficiently separated (as judged by confidence intervals).

On the scaling factors cjc_{j}. With a larger sample complexity that depends on dd, an error bound can be established for μ^jμτ(j)1\|\hat{\mu}_{j}-\vec{\mu}_{\tau(j)}\|_{1} directly (without the unknown scaling factors cjc_{j}). We also remark that the scaling factors can be estimated from the eigenvalues of B^(η)\hat{B}(\vec{\eta}), but we do not pursue this approach as it is subsumed by Algorithm B anyway.

A method of moments for multi-view mixture models

We now consider a much broader class of mixture models and present a general method of moments in this context.

We assume the following conditions on w\vec{w} and the MvM_{v}.

Because the conditional distribution of xv\vec{x}_{v} is not specified beyond its conditional means, it is not possible to develop a maximum likelihood approach to parameter estimation. Instead, as in the document topic model, we develop a method of moments based on solving polynomial equations arising from eigenvalue problems.

2 Observable moments and operators

Lemma 3.1 and Lemma 3.2 are straightforward generalizations of Lemma 2.1 and Lemma 2.2.

In particular, the kk roots of the polynomial λdet(B1,2,3(η)λI)\lambda\mapsto\det(B_{1,2,3}(\vec{\eta})-\lambda I) are {η,μ3,j:j[k]}\{\langle\vec{\eta},\vec{\mu}_{3,j}\rangle:j\in[k]\}.

Recall that Algorithm A relates the eigenvectors of B(η)B(\vec{\eta}) to the matrix of conditional means MM. However, eigenvectors are only defined up to a scaling of each vector; without prior knowledge of the correct scaling, the eigenvectors are not sufficient to recover the parameters MM. Nevertheless, the eigenvalues also carry information about the parameters, as shown in Lemma 3.2, and it is possible to reconstruct the parameters from different the observation operators applied to different vectors η\vec{\eta}. This idea is captured in the following lemma.

Observe that the unknown parameters M3M_{3} are expressed as the solution to a linear system in the above equation, where the elements of the right-hand side LL are the roots of kk-th degree polynomials derived from the second- and third-order observable moments (namely, the characteristic polynomials of the B1,2,3(U3θi)B_{1,2,3}(U_{3}\vec{\theta}_{i}), i[k]\forall i\in[k]). This template is also found in other moment methods based on decompositions of a Hankel matrix. A crucial distinction, however, is that the kk-th degree polynomials in Lemma 3.3 only involve low-order moments, whereas standard methods may involve up to Ω(k)\Omega(k)-th order moments which are difficult to estimate (Lindsay, 1989; Lindsay and Basak, 1993; Gravin et al., 2012).

3 Main result: general estimation procedure and sample complexity bound

The lemmas in the previous section suggest the estimation procedure (Algorithm B) presented in Figure 3.

As stated, the Algorithm B yields an estimator for M3M_{3}, but the method can easily be applied to estimate MvM_{v} for all other views vv. One caveat is that the estimators may not yield the same ordering of the columns, due to the unspecified order of the eigenvectors obtained in the third step of the method, and therefore some care is needed to obtain a consistent ordering. We outline one solution in Appendix B.4.

The sample complexity of Algorithm B depends on the specific concentration properties of x1,x2,x3\vec{x}_{1},\vec{x}_{2},\vec{x}_{3}. We abstract away this dependence in the following condition.

There exist positive scalars N0N_{0}, C1,2C_{1,2}, C1,3C_{1,3}, C1,2,3C_{1,2,3}, and a function f(N,δ)f(N,\delta) (decreasing in NN and δ\delta) such that for any NN0N\geq N_{0} and δ(0,1)\delta\in(0,1),

\Pr\Bigl{[}\|\hat{P}_{a,b}-P_{a,b}\|_{2}\leq C_{a,b}\cdot f(N,\delta)\Bigr{]}\geq 1-\delta for {a,b}{{1,2},{1,3}}\{a,b\}\in\{\{1,2\},\{1,3\}\},

Moreover (for technical convenience), P^1,3\hat{P}_{1,3} is independent of P^1,2,3\hat{P}_{1,2,3} (which may be achieved, say, by splitting a sample of size 2N2N).

For the discrete models such as the document topic model of Section 2.1 and discrete HMMs (Mossel and Roch, 2006; Hsu et al., 2009), Condition 3.2 holds with N0=C1,2=C1,3=C1,2,3=1N_{0}=C_{1,2}=C_{1,3}=C_{1,2,3}=1, and f(N,δ)=(1+ln(1/δ))/Nf(N,\delta)=(1+\sqrt{\ln(1/\delta)})/\sqrt{N}. Using standard techniques (e.g., Chaudhuri et al. (2009); Vershynin (2012)), the condition can also be shown to hold for mixtures of various continuous distributions such as multivariate Gaussians.

Now we are ready to present the main theorem of this section (proved in Appendix B.6).

where P1,2,32:=maxv0P1,2,3(v)2\|P_{1,2,3}\|_{2}:=\max_{\vec{v}\neq\vec{0}}\|P_{1,2,3}(\vec{v})\|_{2}, then with probability at least 15δ1-5\delta, Algorithm B returns M^3=[μ^3,1μ^3,2μ^3,k]\hat{M}_{3}=[\hat{\mu}_{3,1}|\hat{\mu}_{3,2}|\dotsb|\hat{\mu}_{3,k}] with the following guarantee: there exists a permutation τ\tau on [k][k] such that for each j[k]j\in[k],

Applications

In addition to the document clustering model from Section 2, a number of natural latent variable models fit into this multi-view framework. We describe two such cases in this section: Gaussian mixture models and HMMs, both of which have been (at least partially) studied in the literature. In both cases, the estimation technique of Algorithm B leads to new learnability results that were not achieved by previous works.

Condition 3.1 requires that the conditional mean matrix Mv=[μv,1μv,2μv,k]M_{v}=[\vec{\mu}_{v,1}|\vec{\mu}_{v,2}|\dotsb|\vec{\mu}_{v,k}] for each view vv have full column rank. This is similar to the non-degeneracy and spreading conditions used in previous studies of multi-view clustering (Chaudhuri and Rao, 2008; Chaudhuri et al., 2009). In these previous works, the multi-view and non-degeneracy assumptions are shown to reduce the minimum separation required for various efficient algorithms to learn the model parameters. In comparison, Algorithm B does not require a minimum separation condition at all. See Appendix D.3 for details.

While Algorithm B recovers just the means of the mixture components (see Appendix D.4 for details concerning Condition 3.2), we remark that a slight variation can be used to recover the covariances as well. Note that

Under the setting of Lemma 3.2, the matrix given by

satisfies F1,2,3(ϕ,ψ)=(U1M1)diag(ϕ,μ3,tψ,μ3,t+ϕ,Σ3,tψ:t[k])(U1M1)1F_{1,2,3}(\vec{\phi},\vec{\psi})=(U_{1}^{\scriptscriptstyle\top}M_{1})\operatorname{diag}(\langle\vec{\phi},\vec{\mu}_{3,t}\rangle\langle\vec{\psi},\vec{\mu}_{3,t}\rangle+\langle\vec{\phi},\varSigma_{3,t}\vec{\psi}\rangle:t\in[k])(U_{1}^{\scriptscriptstyle\top}M_{1})^{-1} and hence is diagonalizable (in fact, by the same matrices as B1,2,3(η)B_{1,2,3}(\vec{\eta})).

of conditional means and covariances have full rank. This requirement can be met even if the means μv,j\vec{\mu}_{v,j} of the mixture components are all the same.

2 Hidden Markov models

The vector πΔk1\vec{\pi}\in\Delta^{k-1} is the initial state distribution:

The restriction of the HMM to three time steps, say t{1,2,3}t\in\{1,2,3\}, is an instance of the three-view mixture model.

If the hidden variable hh (from the three-view mixture model of Section 3.1) is identified with the second hidden state h2h_{2}, then {x1,x2,x3}\{\vec{x}_{1},\vec{x}_{2},\vec{x}_{3}\} are conditionally independent given hh, and the parameters of the resulting three-view mixture model on (h,x1,x2,x3)(h,\vec{x}_{1},\vec{x}_{2},\vec{x}_{3}) are

From Proposition 4.2, it is easy to verify that B3,1,2(η)=(U3OT)diag(Oη)(U3OT)1B_{3,1,2}(\vec{\eta})=(U_{3}^{\scriptscriptstyle\top}OT)\operatorname{diag}(O^{\scriptscriptstyle\top}\vec{\eta})(U_{3}^{\scriptscriptstyle\top}OT)^{-1}. Therefore, after recovering the observation conditional mean matrix OO using Algorithm B, the Markov chain transition matrix can be recovered using the matrix of right eigenvectors RR of B3,1,2(η)B_{3,1,2}(\vec{\eta}) and the equation (U3O)1R=T(U_{3}^{\scriptscriptstyle\top}O)^{-1}R=T (up to scaling of the columns).

We thank Kamalika Chaudhuri and Tong Zhang for many useful discussions, Karl Stratos for comments on an early draft, David Sontag and an anonymous reviewer for some pointers to related work, and Adel Javanmard for pointing out a problem with Theorem D.1 in an earlier version of the paper.

References

Appendix A Analysis of Algorithm A

In this appendix, we give an analysis of Algorithm A (but defer most perturbation arguments to Appendix C), and also present some illustrative empirical results on text data using a modified implementation.

where the first inequality follows by Cauchy-Schwarz. ∎

A.2 Proof of Theorem 2.1

By Lemma C.6 and the fact U^M(eiej)22σk(U^M)\|\hat{U}^{\scriptscriptstyle\top}M(\vec{e}_{i}-\vec{e}_{j})\|_{2}\geq\sqrt{2}\sigma_{k}(\hat{U}^{\scriptscriptstyle\top}M), we have Pr[E2E1]1δ\Pr[E_{2}|E_{1}]\geq 1-\delta, and thus Pr[E1E2](12δ)(1δ)13δ\Pr[E_{1}\cap E_{2}]\geq(1-2\delta)(1-\delta)\geq 1-3\delta. So henceforth condition on this joint event E1E2E_{1}\cap E_{2}.

The conditions on NN and the bounds in (2), (3), (4), (5), and (7) imply that ε3<12\varepsilon_{3}<\frac{1}{2}. By Lemma C.3, there exists a permutation τ\tau on [k][k] such that, for all j[k]j\in[k],

where sj:=sign(ξ^j,U^μτ(j))s_{j}:=\operatorname{sign}(\langle\hat{\xi}_{j},\hat{U}^{\scriptscriptstyle\top}\vec{\mu}_{\tau(j)}\rangle) and cj:=U^μτ(j)2μτ(j)2c_{j}^{\prime}:=\|\hat{U}^{\scriptscriptstyle\top}\vec{\mu}_{\tau(j)}\|_{2}\leq\|\vec{\mu}_{\tau(j)}\|_{2} (the eigenvectors ξ^j\hat{\xi}_{j} are unique up to sign sjs_{j} because each eigenvalue has geometric multiplicity 11). Since μτ(j)range(U)\vec{\mu}_{\tau(j)}\in\operatorname{range}(U), Lemma C.1 and the bounds in (8) and (6) imply

Therefore, for cj:=sjcj1,U^ξ^jc_{j}:=s_{j}c_{j}^{\prime}\langle\vec{1},\hat{U}\hat{\xi}_{j}\rangle, we have

Making all of the substitutions into the above bound gives

A.3 Some illustrative empirical results

The first and fourth topics appear to be about computer graphics (comp.graphics), the fifth and sixth about baseball (rec.sports.baseball), the third about encryption (sci.crypt), and the second about Christianity (soc.religion.christian).

Appendix B Proofs and details from Section 3

In this section, we provide ommitted proofs and discussion from Section 3.

B.2 Proof of Lemma 3.2

We have U1P1,2U2=(U1M1)diag(w)(M2U2)U_{1}^{\scriptscriptstyle\top}P_{1,2}U_{2}=(U_{1}^{\scriptscriptstyle\top}M_{1})\operatorname{diag}(\vec{w})(M_{2}^{\scriptscriptstyle\top}U_{2}) by Lemma 3.1, which is invertible by the assumptions on UvU_{v} and Condition 3.1. Moreover, also by Lemma 3.1,

B.3 Proof of Lemma 3.3

B.4 Ordering issues

Although Algorithm B only explicitly yields estimates for M3M_{3}, it can easily be applied to estimate MvM_{v} for all other views vv. The main caveat is that the estimators may not yield the same ordering of the columns, due to the unspecified order of the eigenvectors obtained in the third step of the method, and therefore some care is needed to obtain a consistent ordering. However, this ordering issue can be handled by exploiting consistency across the multiple views.

The first step is to perform the estimation of M3M_{3} using Algorithm B as is. Then, to estimate M2M_{2}, one may re-use the eigenvectors in R^1\hat{R}_{1} to diagonalize B^1,3,2(η)\hat{B}_{1,3,2}(\vec{\eta}), as B1,2,3(η)B_{1,2,3}(\vec{\eta}) and B1,3,2(η)B_{1,3,2}(\vec{\eta}) share the same eigenvectors. The same goes for estimating MvM_{v} for other all other views vv except v=1v=1.

It remains to provide a way to estimate M1M_{1}. Observe that M2M_{2} can be estimated in at least two ways: via the operators B^1,3,2(η)\hat{B}_{1,3,2}(\vec{\eta}), or via the operators B^3,1,2(η)\hat{B}_{3,1,2}(\vec{\eta}). This is because the eigenvalues of B3,1,2(η)B_{3,1,2}(\vec{\eta}) and B1,3,2(η)B_{1,3,2}(\vec{\eta}) are the identical. Because the eigenvalues are also sufficiently separated from each other, the eigenvectors R^3\hat{R}_{3} of B^3,1,2(η)\hat{B}_{3,1,2}(\vec{\eta}) can be put in the same order as the eigenvectors R^1\hat{R}_{1} of B^1,3,2\hat{B}_{1,3,2} by (approximately) matching up their respective corresponding eigenvalues. Finally, the appropriately re-ordered eigenvectors R^3\hat{R}_{3} can then be used to diagonalize B^3,2,1(η)\hat{B}_{3,2,1}(\vec{\eta}) to estimate M1M_{1}.

B.5 Estimating the mixing weights

Given the estimate of M^3\hat{M}_{3}, one can obtain an estimate of w\vec{w} using

B.6 Proof of Theorem 3.1

The proof is similar to that of Theorem 2.1, so we just describe the essential differences. As before, most perturbation arguments are deferred to Appendix C.

for all i[k]i\in[k]. Therefore by Condition 3.2 and a union bound, we have Pr[E1]13δ\Pr[E_{1}]\geq 1-3\delta. Second, let E2E_{2} be the event in which

Since each θi\vec{\theta}_{i} is distributed uniformly over Sk1\mathcal{S}^{k-1}, it follows from Lemma C.6 and a union bound that Pr[E2E1]12δ\Pr[E_{2}|E_{1}]\geq 1-2\delta. Therefore Pr[E1E2](13δ)(12δ)15δ\Pr[E_{1}\cap E_{2}]\geq(1-3\delta)(1-2\delta)\geq 1-5\delta.

Also define R:=U^1M1diag(U^1M1e12,U^1M1e22,,U^1M1ek2)1R:=\hat{U}_{1}^{\scriptscriptstyle\top}M_{1}\operatorname{diag}(\|\hat{U}_{1}^{\scriptscriptstyle\top}M_{1}\vec{e}_{1}\|_{2},\|\hat{U}_{1}^{\scriptscriptstyle\top}M_{1}\vec{e}_{2}\|_{2},\dotsc,\|\hat{U}_{1}^{\scriptscriptstyle\top}M_{1}\vec{e}_{k}\|_{2})^{-1}. Using most of the same arguments in the proof of Theorem 2.1, we have

Finally, by Lemma C.1 (as applied to P^1,3\hat{P}_{1,3} and P1,3P_{1,3}),

Making all of the substitutions into the above bound gives

Appendix C Perturbation analysis for observable operators

The following lemma establishes the accuracy of approximating the fundamental subspaces (i.e., the row and column spaces) of a matrix XX by computing the singular value decomposition of a perturbation X^\hat{X} of XX.

σk(X^)=σk(U^X^V^)(1ε0)σk(X)>0\sigma_{k}(\hat{X})=\sigma_{k}(\hat{U}^{\scriptscriptstyle\top}\hat{X}\hat{V})\geq(1-\varepsilon_{0})\cdot\sigma_{k}(X)>0;

σk(U^U)1ε12\sigma_{k}(\hat{U}^{\scriptscriptstyle\top}U)\geq\sqrt{1-\varepsilon_{1}^{2}};

σk(V^V)1ε12\sigma_{k}(\hat{V}^{\scriptscriptstyle\top}V)\geq\sqrt{1-\varepsilon_{1}^{2}};

σk(U^XV^)(1ε12)σk(X)\sigma_{k}(\hat{U}^{\scriptscriptstyle\top}X\hat{V})\geq(1-\varepsilon_{1}^{2})\cdot\sigma_{k}(X);

The next lemma bounds the error of the observation operator in terms of the errors in estimating the second-order and third-order moments.

U^X^V^\hat{U}^{\scriptscriptstyle\top}\hat{X}\hat{V} and U^XV^\hat{U}^{\scriptscriptstyle\top}X\hat{V} are both invertible, and (U^X^V^)1(U^XV^)12ε2σk(X)\|(\hat{U}^{\scriptscriptstyle\top}\hat{X}\hat{V})^{-1}-(\hat{U}^{\scriptscriptstyle\top}X\hat{V})^{-1}\|_{2}\leq\frac{\varepsilon_{2}}{\sigma_{k}(X)};

(U^Y^V^)(U^X^V^)1(U^YV^)(U^XV^)12ϵY(1ε0)σk(X)+Y2ε2σk(X)\|(\hat{U}^{\scriptscriptstyle\top}\hat{Y}\hat{V})(\hat{U}^{\scriptscriptstyle\top}\hat{X}\hat{V})^{-1}-(\hat{U}^{\scriptscriptstyle\top}Y\hat{V})(\hat{U}^{\scriptscriptstyle\top}X\hat{V})^{-1}\|_{2}\leq\frac{\epsilon_{Y}}{(1-\varepsilon_{0})\cdot\sigma_{k}(X)}+\frac{\|Y\|_{2}\cdot\varepsilon_{2}}{\sigma_{k}(X)}.

where the first inequality follows from the triangle inequality, the second follows from the sub-multiplicative property of the spectral norm, and the last follows from Lemma C.1 and the first claim. ∎

The following lemma establishes standard eigenvalue and eigenvector perturbation bounds.

The Bauer-Fike theorem (Lemma E.3) implies that for every eigenvalue λ^i\hat{\lambda}_{i} of A^\hat{A}, there exists an eigenvalue λj\lambda_{j} of AA such that λ^iλjR1(A^A)R2ε3γA|\hat{\lambda}_{i}-\lambda_{j}|\leq\|R^{-1}(\hat{A}-A)R\|_{2}\leq\varepsilon_{3}\cdot\gamma_{A}. Therefore, the assumption on ε3\varepsilon_{3} implies that there exists a permutation τ\tau such that λ^τ(i)λiε3γA<γA2|\hat{\lambda}_{\tau(i)}-\lambda_{i}|\leq\varepsilon_{3}\cdot\gamma_{A}<\frac{\gamma_{A}}{2}. In particular,

Since A^\hat{A} is real, all non-real eigenvalues of A^\hat{A} must come in conjugate pairs; so the existence of a non-real eigenvalue of A^\hat{A} would contradict (12). This proves the first claim.

again by the triangle inequality. Therefore, it suffices to show ci,j2R12ε3|c_{i,j}|\leq 2\|R^{-1}\|_{2}\cdot\varepsilon_{3} for jij\neq i to prove the second claim.

Observe that Aξ^i=A(i=1kci,iξi)=i=1kci,iλiξiA\hat{\xi}_{i}=A(\sum_{i^{\prime}=1}^{k}c_{i,i^{\prime}}\vec{\xi}_{i^{\prime}})=\sum_{i^{\prime}=1}^{k}c_{i,i^{\prime}}\lambda_{i^{\prime}}\vec{\xi}_{i^{\prime}}, and therefore

Multiplying through the above equation by ζj\vec{\zeta}_{j}^{\scriptscriptstyle\top}, and using the fact that ζjξi=\mathds1{j=i}\vec{\zeta}_{j}^{\scriptscriptstyle\top}\vec{\xi}_{i^{\prime}}=\mathds{1}\{j=i^{\prime}\} gives

The above equation rearranges to (λjλi)ci,j=(λ^iλi)ζjξ^i+ζj(AA^)ξ^i(\lambda_{j}-\lambda_{i})c_{i,j}=(\hat{\lambda}_{i}-\lambda_{i})\vec{\zeta}_{j}^{\scriptscriptstyle\top}\hat{\xi}_{i}+\vec{\zeta}_{j}^{\scriptscriptstyle\top}(A-\hat{A})\hat{\xi}_{i} and therefore

by the Cauchy-Schwarz and triangle inequalities and the sub-multiplicative property of the spectral norm. The bound ci,j2R12ε3|c_{i,j}|\leq 2\|R^{-1}\|_{2}\cdot\varepsilon_{3} then follows from the first claim.

The third claim follows from standard comparisons of matrix norms. ∎

The next lemma gives perturbation bounds for estimating the eigenvalues of simultaneously diagonalizable matrices A1,A2,,AkA_{1},A_{2},\dotsc,A_{k}. The eigenvectors R^\hat{R} are taken from a perturbation of the first matrix A1A_{1}, and are then subsequently used to approximately diagonalize the perturbations of the remaining matrices A2,,AkA_{2},\dotsc,A_{k}. In practice, one may use Jacobi-like procedures to approximately solve the joint eigenvalue problem.

The matrix R^\hat{R} is invertible and its inverse satisfies R^1Rτ12R12ε41ε4\|\hat{R}^{-1}-R_{\tau}^{-1}\|_{2}\leq\|R^{-1}\|_{2}\cdot\frac{\varepsilon_{4}}{1-\varepsilon_{4}};

For all i{2,3,,k}i\in\{2,3,\dotsc,k\} and all j[k]j\in[k], the (j,j)(j,j)-th element of R^1A^iR^\hat{R}^{-1}\hat{A}_{i}\hat{R}, denoted by λ^i,j:=ejR^1A^iR^ej\hat{\lambda}_{i,j}:=\vec{e}_{j}^{\scriptscriptstyle\top}\hat{R}^{-1}\hat{A}_{i}\hat{R}\vec{e}_{j}, satisfies

If ε412\varepsilon_{4}\leq\frac{1}{2}, then λ^i,jλi,τ(j)3ε3γA+4κ(R)ε4λmax|\hat{\lambda}_{i,j}-\lambda_{i,\tau(j)}|\leq 3\varepsilon_{3}\cdot\gamma_{A}+4\kappa(R)\cdot\varepsilon_{4}\cdot\lambda_{\max}.

Observe that ζj2R12\|\vec{\zeta}_{j}\|_{2}\leq\|R^{-1}\|_{2}, ξj2R2\|\vec{\xi}_{j}\|_{2}\leq\|R\|_{2}, ζ^jζj2R^1Rτ12R12ε41ε4\|\hat{\zeta}_{j}-\vec{\zeta}_{j}\|_{2}\leq\|\hat{R}^{-1}-R_{\tau}^{-1}\|_{2}\leq\|R^{-1}\|_{2}\cdot\frac{\varepsilon_{4}}{1-\varepsilon_{4}}, ξ^jξj24kR12ε3\|\hat{\xi}_{j}-\vec{\xi}_{j}\|_{2}\leq 4k\cdot\|R^{-1}\|_{2}\cdot\varepsilon_{3} (by Lemma C.3), and Ai2R2(maxjλi,j)R12\|A_{i}\|_{2}\leq\|R\|_{2}\cdot(\max_{j}|\lambda_{i,j}|)\cdot\|R^{-1}\|_{2}. Therefore, continuing from (13), λ^i,jλi,τ(j)|\hat{\lambda}_{i,j}-\lambda_{i,\tau(j)}| is bounded as

Rearranging gives the claimed inequality. ∎

We have R=Vdiag(Ve12,Ve22,,Vek2)1R=V\operatorname{diag}(\|V\vec{e}_{1}\|_{2},\|V\vec{e}_{2}\|_{2},\dotsc,\|V\vec{e}_{k}\|_{2})^{-1}, so by the sub-multiplicative property of the spectral norm, R2V2/minjVej2V2/σk(V)=κ(V)\|R\|_{2}\leq\|V\|_{2}/\min_{j}\|V\vec{e}_{j}\|_{2}\leq\|V\|_{2}/\sigma_{k}(V)=\kappa(V). Similarly, R12V12maxjVej2V12V2=κ(V)\|R^{-1}\|_{2}\leq\|V^{-1}\|_{2}\cdot\max_{j}\|V\vec{e}_{j}\|_{2}\leq\|V^{-1}\|_{2}\cdot\|V\|_{2}=\kappa(V). ∎

\displaystyle\Pr\biggl{[}\min_{i\neq j}|\langle\vec{\theta},A(\vec{e}_{i}-\vec{e}_{j})\rangle|>\frac{\min_{i\neq j}\|A(\vec{e}_{i}-\vec{e}_{j})\|_{2}\cdot\delta}{\sqrt{em}{n\choose 2}}\biggr{]}\geq 1-\delta.

\displaystyle\Pr\biggl{[}\forall i\in[m],\ |\langle\vec{\theta},A\vec{e}_{i}\rangle|\leq\frac{\|A\vec{e}_{i}\|_{2}}{\sqrt{m}}\Bigl{(}1+\sqrt{2\ln(m/\delta)}\Bigr{)}\biggr{]}\geq 1-\delta.

For the first claim, let δ0:=δ/(n2)\delta_{0}:=\delta/{n\choose 2}. By Lemma F.2, for any fixed pair {i,j}[n]\{i,j\}\subseteq[n] and β:=δ0/e\beta:=\delta_{0}/\sqrt{e},

Therefore the first claim follows by a union bound over all (n2){n\choose 2} pairs {i,j}\{i,j\}.

For the second claim, apply Lemma F.2 with β:=1+t\beta:=1+t and t:=2ln(m/δ)t:=\sqrt{2\ln(m/\delta)} to obtain

Therefore the second claim follows by taking a union bound over all i[m]i\in[m]. ∎

Appendix D Proofs and details from Section 4

In this section, we provide ommitted proofs and details from Section 4.

As in the proof of Lemma 3.1, it is easy to show that

The claim then follows from the same arguments used in the proof of Lemma 3.2.

D.2 Proof of Proposition 4.2

The conditional independence properties follow from the HMM conditional independence assumptions. To check the parameters, observe first that

The rest of the parameters are similar to verify.

D.3 Learning mixtures of product distributions

Chaudhuri and Rao (2008) show that under a similar condition (which they call a spreading condition), a random partitioning of the coordinates into two “views” preserves the separation between the means of kk component distributions. They then follow this preprocessing with a projection based on the correlations across the two views (similar to CCA). However, their overall algorithm requires a minimum separation condition on the means of the component distributions. In contrast, Algorithm B does not require a minimum separation condition at all in this setting.

Follows from Lemma D.1 (below) together with a union bound. ∎

D.4 Empirical moments for multi-view mixtures of subgaussian distributions

The required concentration behavior of the empirical moments used by Algorithm B can be easily established for multi-view Gaussian mixture models using known techniques (Chaudhuri et al., 2009). This is clear for the second-order statistics P^a,b\hat{P}_{a,b} for {a,b}{{1,2},{1,3}}\{a,b\}\in\{\{1,2\},\{1,3\}\}, and remains true for the third-order statistics P^1,2,3\hat{P}_{1,2,3} because x3\vec{x}_{3} is conditionally independent of x1\vec{x}_{1} and x2\vec{x}_{2} given hh. The magnitude of U^3θi,x3\langle\hat{U}_{3}\vec{\theta}_{i},\vec{x}_{3}\rangle can be bounded for all samples (with a union bound; recall that we make the simplifying assumption that P^1,3\hat{P}_{1,3} is independent of P^1,2,3\hat{P}_{1,2,3}, and therefore so are U^3\hat{U}_{3} and P^1,2,3\hat{P}_{1,2,3}). Therefore, one effectively only needs spectral norm error bounds for second-order statistics, as provided by existing techniques.

Indeed, it is possible to establish Condition 3.2 in the case where the conditional distribution of xv\vec{x}_{v} given hh (for each view vv) is subgaussian. Specifically, we assume that there exists some α>0\alpha>0 such that for each view vv and each component j[k]j\in[k],

Appendix E General results from matrix perturbation theory

The lemmas in this section are standard results from matrix perturbation theory, taken from Stewart and Sun (1990).

See Theorem 4.11, p. 204 in Stewart and Sun (1990). ∎

See Theorem 4.4, p. 262 in Stewart and Sun (1990). ∎

See Theorem 3.3, p. 192 in Stewart and Sun (1990). ∎

See Theorem 2.5, p. 118 in Stewart and Sun (1990). ∎

Appendix F Probability inequalities

Fix μ=(μ1,μ2,,μn)Δm1\vec{\mu}=(\mu_{1},\mu_{2},\dotsc,\mu_{n})\in\Delta^{m-1}. Let x\vec{x} be a random vector for which Pr[x=ei]=μi\Pr[\vec{x}=\vec{e}_{i}]=\mu_{i} for all i[m]i\in[m], and let x1,x2,,xn\vec{x}_{1},\vec{x}_{2},\dotsc,\vec{x}_{n} be nn independent copies of x\vec{x}. Set μ^:=(1/n)i=1nxi\hat{\mu}:=(1/n)\sum_{i=1}^{n}\vec{x}_{i}. For all t>0t>0,

This is a special case of Lemma 2.2 from Dasgupta and Gupta (2003). ∎

This is a direct corollary of Theorem 19 from Ahlswede and Winter (2002). ∎

Appendix G Insufficiency of second-order moments

Chang (1996) shows that a simple class of Markov models used in mathematical phylogenetics cannot be identified from pair-wise probabilities alone. Below, we restate (a specialization of) this result in terms of the document topic model from Section 2.1.

1Q=1\vec{1}^{\scriptscriptstyle\top}Q=\vec{1}^{\scriptscriptstyle\top};

MQ1MQ^{-1}, Qdiag(w)Mdiag(Mw)1Q\operatorname{diag}(\vec{w})M^{\scriptscriptstyle\top}\operatorname{diag}(M\vec{w})^{-1}, and QwQ\vec{w} have non-negative entries;

Qdiag(w)QQ\operatorname{diag}(\vec{w})Q^{\scriptscriptstyle\top} is a diagonal matrix.

A simple example for d=k=2d=k=2 can be obtained from

for some p(0,1)p\in(0,1). We take p=0.25p=0.25, in which case QQ satisfies the conditions of Proposition G.1, and

However, the triple-wise probabilities, for η=(1,0)\eta=(1,0), differ: for (M,w)(M,\vec{w}), we have