The Identification of Discrete Mixture Models Citation Gordon, Spencer Lane (2023) The Identification of Discrete Mixture Models. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/ebf5-0b48. https://resolver.caltech.edu/CaltechTHESIS:02072023-112938936 Abstract In this thesis we discuss a variety of results on learning and identifying discrete mixture models, i.e., distributions that are a convex combination of k from a known class C of distributions. We first consider the case where C is the class of binomial distributions, before generalizing to the case of product distributions. We provide a necessary condition for identifiability of mixture of products distributions as well as a generalization to structured mixtures over multiple latent variables. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Mixture Models; Identifiability; Algorithms; Complexity; Statistical Learning Theory Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Computer Science Thesis Availability: Public (worldwide access) Research Advisor(s): Vidick, Thomas G. (co-advisor) Schulman, Leonard J. (co-advisor) Thesis Committee: Yue, Yisong (chair) Vidick, Thomas G. Schulman, Leonard J. Wierman, Adam C. Defense Date: 19 January 2023 Record Number: CaltechTHESIS:02072023-112938936 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:02072023-112938936 DOI: 10.7907/ebf5-0b48 Related URLs: URL URL Type Description https://doi.org/10.1109/TIT.2022.3146630 DOI Article adapted for Chapter 1 ORCID: Author ORCID Gordon, Spencer Lane 0000-0002-7101-2370 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 15101 Collection: CaltechTHESIS Deposited By: Spencer Gordon Deposited On: 17 Feb 2023 17:37 Last Modified: 28 May 2025 22:17 Thesis Files PDF - Final Version See Usage Policy. 842kB Repository Staff Only: item control page

The Identification of Discrete Mixture Models Thesis by Spencer Lane Gordon In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in Computer Science CALIFORNIA INSTITUTE OF TECHNOLOGY Pasadena, California 2023 Defended January 19, 2023 ii © 2023 Spencer Lane Gordon ORCID: 0000-0002-7101-2370 All rights reserved iii ACKNOWLEDGEMENTS [Intentionally left blank] iv ABSTRACT In this thesis we discuss a variety of results on learning and identifying discrete mixture models, i.e., distributions that are a convex combination of k from a known class C of distributions. We first consider the case where C is the class of binomial distributions, before generalizing to the case of product distributions. We provide a necessary condition for identifiability of mixture of products distributions as well as a generalization to structured mixtures over multiple latent variables. v PUBLISHED CONTENT AND CONTRIBUTIONS [1] S. L. Gordon and L. J. Schulman, “Hadamard extensions and the identification of mixtures of product distributions,” IEEE Transactions on Information Theory, vol. 68, no. 6, pp. 4085–4089, 2022. doi: 10.1109/ TIT.2022.3146630, S.G. participated in all parts of the project. [2] S. Gordon, B. H. Mazaheri, Y. Rabani, and L. Schulman, “Source identification for mixtures of product distributions,” in Proceedings of Thirty Fourth Conference on Learning Theory, M. Belkin and S. Kpotufe, Eds., ser. Proceedings of Machine Learning Research, vol. 134, PMLR, Aug. 2021, pp. 2193–2216. [Online]. Available: https://proceedings.mlr. press/v134/gordon21a.html, S.G. participated in all parts of the project. [3] S. Gordon, B. Mazaheri, L. J. Schulman, and Y. Rabani, “The sparse hausdorff moment problem, with application to topic models,” CoRR, vol. abs/2007.08101, 2020. arXiv: 2007 . 08101. [Online]. Available: https://arxiv.org/abs/2007.08101, S.G. participated in all parts of the project. vi TABLE OF CONTENTS Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Published Content and Contributions . . . . . . . . . . . . . . . . . . . v Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Chapter I: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Chapter II: The k-Mix IID problem . . . . . . . . . . . . . . . . . . . . . 5 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Mixture Models and Other Definitions . . . . . . . . . . . . . . 11 2.3 Properties of Hankel Matrices . . . . . . . . . . . . . . . . . . . 13 2.4 The Empirical Moments . . . . . . . . . . . . . . . . . . . . . . . 14 2.5 Learning the Source . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.6 Implications for Topic Models . . . . . . . . . . . . . . . . . . . 18 2.7 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.8 Computing the Weights . . . . . . . . . . . . . . . . . . . . . . . 22 2.9 Deferred Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.10 Useful Theorems . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Chapter III: Sufficient Conditions for the Identifiability of Mixtures of Products . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Some Theory for Hadamard Products, and a Proof of Theorem 45 35 3.4 Combinatorics of the NAE Condition: Proof of Theorem 47(a) 37 3.5 From NAE to Rank: Proof of Theorem 47(b) . . . . . . . . . . . 39 Chapter IV: Source Identification for Mixtures of Products . . . . . . . 41 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.3 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.4 The Condition Number Bound . . . . . . . . . . . . . . . . . . . 56 4.5 Analysis of the Algorithm . . . . . . . . . . . . . . . . . . . . . . 57 Chapter V: The Identifiability of Uniform Mixtures of Binomial Distributions with Log-Linear Influences . . . . . . . . . . . . . . . . . . . 65 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 vii LIST OF ILLUSTRATIONS Number Page 3.1 A Bayesian network diagram relating H and X1 , . . . , Xn . . . . . 34 3.2 Argument for Theorem 47(a). Upper-left region is white. Entries (t, f (t)) (indicated with black dots) are not white. . . . . . 38 4.1 Graphical depiction of a k-MixProd . . . . . . . . . . . . . . . . 41 1 Chapter 1 INTRODUCTION In this thesis, we will investigate problems of the following form: We have a distribution over vectors X = (X1 , X2 , . . . , Xn ) ∈ [B]n for B ∈ N with P (X1 , . . . , Xn ) a mixture of k distributions from a class D of distributions with a given parameterization. We are given samples of X and are asked to determine the mixture weights, that is, the weight of each mixture constituent in the distribution, along with the parameters determining the mixture constituents. Our goal is to recover these parameter up to an input accuracy ε with probability 1 − δ, where δ is also part of our input. We’d like to use as few samples as possible. We will write our target distribution as a convex combination of k distributions P1 (θ1 ), . . . , Pk (θk ) ∈ D, with non-negative weights π1 , . . . , πk summing to 1. Definition 1. The (k − 1)-simplex, ∆k−1 ⊆ Rk , consists of all vectors π ∈ Rk≥0 P such that i πi = 1. Definition 2 (Mixture of discrete distributions). A k-mixture of distributions is parameterized by distributions P1 , . . . , Pk and mixing weights π ∈ ∆k−1 and is the marginal distribution over X = (X1 , X2 , . . . , Xn ) of the joint distribution on (X, U ) given by P (X = x | U = u) = Pu (X = x) and P (U = u) = πu for all u ∈ [k]. Notation 3. When D is a family of discrete distributions parameterized by θ ∈ Rn we can write a mixture of k distributions from D as P (X) = π1 P (X; θ1 ) + · · · + πk P (X; θk ). When we don’t care about the parameterization, we will write P (X) = π1 P1 (X) + · · · + πk Pk (X). The k-Mix D problem • Input: The input to this problem is a distribution P (X) that is a mixture of k distributions from D with unknown mixing weights π1 , . . . , πk where the ith mixture constituent is parameterized by unknown parameter θi , along with parameters ε, δ > 0. 2 • Goal: The goal is to design an algorithm that uses S samples from P and with probability at least 1 − δ outputs the following: a vector π̃ ∈ ∆k−1 along with distributions P (·; θ̃1 ), . . . , P (·; θ̃k ) ∈ D such that ∥π̃ − π∥∞ ≤ ε and θ̃i − θi ≤ ε for i = 1, . . . , k. ∞ Identifiability There is one way in which the k-Mix D problem is too hard, no matter which class of distributions D we are considering; permuting the mixture components by applying some permutation σ : [k] → [k] will not change any of the observed statistics. In particular, let P ′ (X) = π1′ P (X; θ1′ ) + π2′ P (X; θ2′ ) + · · · + πk′ P (X; θk′ ) where πu′ = πσ−1 (u) and θu′ = θσ′ −1 (u) for all u. Then P = P ′ and there is no way to distinguish between the distribution on X induced by these two distinct mixture models. Thus, we will always be interested in approximately recovering the parameters only up to the equivalence class of mixture models induced by permuting mixture components. • Revised Goal: Our revised goal is to design an algorithm that uses S samples from P and with probability at least 1−δ outputs the following: a vector π̃ ∈ ∆k−1 along with distributions P (·; θ̃1 ), . . . , P (·; θ̃k ) ∈ D such that for some permutation σ : [k] → [k] the vector πσ obtained by permuting coordinates according to σ satisfies ∥π̃ − πσ ∥∞ ≤ ε and θ̃i − θσ(i) ≤ ε for i = 1, . . . , k. ∞ Having weakened our requirement for a solution to the problem, there is another way in which the k-Mix D problem may still be too hard for certain classes of distributions D and certain concrete choices of distributions from D. There may be multiple equivalence classes of parameters that generate the same distribution on X. To illustrate this phenomenon, consider the following class of distributions, where n ∈ N is fixed: IID := {IID(θ) : θ ∈ [0, 1]} , where X ∈ {0, 1}n and P (X = x) = Qn i=1 θ xi (1 − θ)1−xi for P = IID(θ). Example 4. In both of the following, let D = IID and let k = 3. 3 Degenerate mixture Consider π = (0, 1/2, 1/2) and θ = (θ1 , θ2 , θ3 ) = (1/2, 1/4, 3/4). Then 1 1 × P (X = x; 1/4) + × P (X = x; 3/4) 2 2 1 1 = × P (X = x; 1/4) + × P (X = x; 3/4), 2 2 P (X = x) = 0 × P (X = x; θ1 ) + so any choice of θ1 gives the same distribution on X. Identical mixture components Consider π = (1/3, 1/3, 1/3) and θ = (1/2, 1/2, 1/4). Then P (X = x) = π1 × P (X = x; 1/2) + π2 × P (X = x; 1/2) + π3 × P (X = x; 1/4) = (π1 + π2 ) × P (X = x; 1/2) + π3 × P (X = x; 1/4), so any choice of π1 and π2 that retains π1 + π2 = 2/3 will give the same distribution on X. What we need for the k-Mix D problem to be solvable is identifiability of the mixture model. Definition 5 (Identifiability). We will say that a k-Mix D instance π, θ is identifiable when P (X; π, θ) ̸= P (X; π ′ , θ′ ) for any π ′ , θ′ not in the equivalence class of parameters containing π, θ. Note that identifiability is a property of instances and not of problems. In this work we will be very interested in understanding when a given instance is identifiable. 1.1 Organization The k-Mix IID Problem The first case, considered in Chapter 2, will be mixtures of IID distributions, i.e., D := IID. For the k-Mix IID problem, there is a complete characterization of identifiability in terms of conditions on the parameters and the number of copies of X1 available. Moreover, there are known algorithms for solving the problem. In this chapter, we present a new analysis of Prony’s method, a classical method for solving this problem, and obtain the best known upper bounds on the required sample complexity and runtime for solving the k-Mix IID problem. Finally, we present an algorithm for the non-binary case, where X ∈ [B]n for some B ∈ N. 4 The k-Mix Prod Problem The next case, considered in Chapters 3 and 4, is that of product distributions. That is, n Y Pu (X = x; θ) = θixi (1 − θi )1−xi i=1 n for θ ∈ [0, 1] . In Chapter 3 we present a new sufficient condition for identifiability that generalizes the identifiability condition for the k-Mix IID problem and in Chapter 4 we give an algorithm to solve the k-Mix Product problem under a stronger identifiability assumption. Generalizations In Chapter 5, we consider a slightly different model and prove an identifiability result for this model. 5 Chapter 2 THE k-MIX IID PROBLEM [1] S. Gordon, B. Mazaheri, L. J. Schulman, and Y. Rabani, “The sparse hausdorff moment problem, with application to topic models,” CoRR, vol. abs/2007.08101, 2020. arXiv: 2007 . 08101. [Online]. Available: https://arxiv.org/abs/2007.08101, 2.1 Introduction We consider the problem of learning a mixture of k IID distributions, k-Mix IID. This is equivalent to learning a mixture of k binomial distributions; motivated by the following analogy, we call this the k-coin problem. (In the literature this is also called the k-spike problem.) Consider a set of k visually indistinguishable coins. When tossed, coin j has (unknown) probability θj of coming up heads. There is an (unknown) probability distribution π on the coins. Our sampling regimen is this: a coin is picked according to π. We do not discover which coin we picked. We then toss this coin m times, yielding a sequence X = (X1 , . . . , Xm ) ∈ {0, 1}m . Repeat. The samples X are distributed according to Pr(X = x) = Pr(X1 = x1 , . . . , Xm = xm ) = k X j=1 πj m Y θjxi (1 − θj )1−xi . (2.1) i=1 We write M = (θ, π) for the parameters of this model. Our goal is to start from empirical statistics on X which are close to the probabilities in (2.1), and from that to produce reconstructed parameters π̃1 , . . . , π̃k and θ̃1 , . . . , θ̃k such that θ − θ̃ and ∥π − π̃∥∞ are small. ∞ Equivalent formulation. The sample statistics are of course invariant unP der permutation of [m], so only Xi matters. Therefore we can effectively think of the distribution in (2.1) as a mixture of binomial distributions: Binomial(m, θj ) with mixture weights (πj ). 6 We will show below that estimates of the probabilities in (2.1) give estimates of the following quantities, µi := X πj θji , i = 0, . . . , 2k, (2.2) j which we will refer to as the moments of the distribution P, where we define P as the probability measure on [0, 1] supported on k discrete points θj , each with weight πj . That is, we can write Z 1 µi := θi dP(θ) (i ≥ 0). (2.3) 0 Obviously as m → ∞ we are able to learn from a single sample the bias of the chosen coin; and thus can determine the distribution on coins with enough samples. On the other hand, even with exact statistics, we need m ≥ 2k in order to verify that P is supported on at most k points (this will be explained below). In this chapter we will solve the problem for m = 2k. Understanding the relationship between the sequence (µi ) and measures on [0, 1] is a classical problem in probability, known as the Hausdorff moment problem. The Sparse Hausdorff Moment Problem. The Hausdorff moment problem is that of determining what moment sequences (µi ) are possible for a probability measure P supported on [0, 1]. See [1], [2] for background on this classical problem. Our problem is the computational version of the Hausdorff moment problem, “sparse” because the measure P is supported on an unknown but discrete set {θ1 , . . . , θk }. Motivation. The problem of reconstructing the parameters (θ, π) is quite useful: (i) By a known reduction, any algorithm for this problem lifts to an algorithm for learning Topic Models. This will be discussed in Sec. 2.6. (ii) The moments we analyze correspond to the cumulative distribution function of an exponential distribution. Thus, the problem of learning mixtures of binomial distributions is equivalent to learning mixtures of exponential 7 distributions. [3] show that learning population histories from coalescence times reduces to learning mixtures of exponential distributions. The use of [3] of the Matrix Pencil Method is interchangeable with our version of Prony’s method, which improves time complexity. Other applications of learning mixtures of binomial or (equivalently) exponential distributions abound in various areas, including ecology, geology, social sciences, and computer systems. For example, see [4] for an application in network evaluation. (iii) This problem is a special case of the problem of identifying mixture models of k product distributions on binary variables, a problem on which there has been an impressive sequence of contributions in the last two decades, as we will discuss below. In chapter 4, our results for the iid case are used to improve the complexity of the k-Mix Prod problem. In particular, we use the algorithm for k-Mix IID as a subroutine in the solution to k-Mix Prod. (iv) Algorithms for identifying mixtures of product distributions are the simplest case of the yet-more-general problem of identifying distributions on directed graphical models [5]. There has been little work in this direction, [6] being a notable exception. However, even that work has to make strong assumptions about the distributions of the variables Xi , and in particular they cannot be binary. (Except for the case k = 2, but we are concerned here with complexity of the problem as a function of k.) Source identification in causal graphical models is an important direction for future research, and the dependence of the sample size and runtime complexity on k matters a great deal. (One may think of k as a bound on how much “confounding” there is among the observables in the model.) Prony’s method. The algorithm we analyze is essentially that of Prony, 1795 [7]. The idea is to (a) characterize the coin biases (the support of P) as the roots of a polynomial whose coefficient vector is the kernel of the Hankel matrix; (b) use polynomial root-finding to determine the empirical coin biases; (c) reconstruct the mixture weights by polynomial interpolation. Prior work on the sparse Hausdorff, i.e., k-coin mixture, problem. It has long been acknowledged in the numerical analysis literature (e.g., [1] §9.4, [8]) that the Prony method is highly sensitive to sample error (i.e., to errors in the moments). This instability is also inherent to our problem (the source identification of k-coin mixture models); a lower bound ([9] Thm 6.1) 8 shows that even for any constant c ≥ 2, if ck (rather than just the minimum 2k) noisy moments are available, accurate source identification is possible only if those moments are available to accuracy exp(−Ω(k)). To be specific here, the input to our problem is an empirical moment sequence µ̃0 , . . . , µ̃m ; its accuracy is maxi |tildeµi − µi | where µi is as in (2.3). In order to obtain this short list of moments to within accuracy ε one needs sample size roughly 1/ε2 . Throughout the paper, therefore, accuracies of exp(−Ω(k)) in the moments, translate to sample size bounds of exp(Ω(k)), and vice versa; we shall therefore have no need to comment separately on sample size.1 The prior upper bounds for the problem were: (i) [9] re-invented the Prony 2 method and solved the problem using moment accuracy min{ζ O(k) , k −O(k ) } and runtime poly(k). (ii) A different algorithm in [10] improved require2 ments in the moment accuracy to k −O(k) but required runtime k O(k ) . (iii) Motivated by a problem in population genetics that reduces to the k-coin mixture problem, [3] analyzed a solution using the Matrix Pencil Method 4 (MPM), which required moment accuracy πmin ζ O(k) . We note that [9], [10] do not depend on πmin and [10] does not depend on ζ. This is the result of analyzing error in terms of transportation norm, which allows coins with equivalent biases (ζ = 0) to be merged and improbable coins (πmin ≈ 0) to be ignored without severe consequences. The Matrix Pencil Method used by [3] gives the desired parameters of the model (θ) directly, yielding a straightforward stability analysis. While runtime is not discussed in [3], the method requires solving a generalized eigenvalue problem, which in practice requires time O(k 3 ).2 We propose improving runtime by instead using Prony’s Method, a close “relative” of the Matrix Pencil Method. Prony’s Method gives coefficients of a polynomial whose 1 Of course, collecting samples also takes time, so one might ask whether the sample complexity should be included within the runtime term. However, the process of sampling and computing the moments or frequencies is computationally trivial. It can often be done under a very restrictive computational model, such as streaming. Or samples might be collected in parallel. Or the frequencies might be otherwise available from an external source. Thus, we make a clear distinction between two resources which a source identification algorithm requires: the sample size, which as noted, translates directly into the accuracy of the moments; and the runtime, given those moments. 2 It is possible that the runtime could be improved to the time it takes to multiply two k × k matrices. This is still much worse than O(k 2 ), and the best guarantees hide impractical constants. 9 roots lie at the desired parameters (θ). This significantly complicates the stability analysis, which becomes the main undertaking of this paper. Our result. Our main result is to provide source identification of k-coin mixtures using the Prony method, simultaneously requiring moment accuracy 2 only πmin ζ O(k) , and achieving runtime O(k 2+o(1) ). The main technical contributions needed for the result are Theorem 16 and Lemma 28, which give quantitative characterizations of the error propagation occurring in Prony’s method. Recall that Prony’s method interprets the kernel of a Hankel matrix as a polynomial with roots corresponding to the support of P. We are able to show that the approximate Hankel matrix has a pseudo-kernel that is close to the kernel of the exact Hankel matrix; this pseudo-kernel has roots close to the roots of the exact kernel when interpreted as a polynomial. Our result also implies an improvement in identifying pure topic models, via the reductions in [9], [10]. These reductions require solving k binary instances, and the required accuracy of those solutions necessitates a postreduction moment accuracy of exp(−Ω(k 2 log k)) in both papers (hence, required sample size exp(O(k 2 log k))). Our result improves the required accuracy the moments to exp(−Ω(k log k)) (sample size exp(O(k log k))) and the runtime to O(k 3+o(1) ). A more detailed comparison with previous work on topic models is given in Section 2.6. We have posted a working implementation of the algorithm on the following public Jupyter Notebook: Online notebook implementation.3 (Tested in Chrome and Safari.) Related work. The k-coin problem becomes easier when m is superlinear in k, and trivial when m is Ω(k 2 log k). Therefore, we focus on the smallest m for which the problem is solvable, which is m = 2k − 1 if k is assumed, or m = 2k if k needs to be verified. As noted previously, three prior papers gave algorithms with worse performance than ours. Roughly stating the results (ignoring dependence on ζ and on πmin ), they are as follows. The 2 paper [9] solved the problem with moment accuracy k −O(k ) (sample size 3 https://colab.research.google.com/drive/1qR6VOYSjq08LPxqHhyY0ap_ VL1apt9yS?usp=sharing 10 2 k O(k ) ) and runtime poly(k). That paper also proves a lower bound of exp(k) on the sample size needed to solve the problem. Subsequently, a different solution using near optimal sample size s = k O(k) , but much worse runtime 2 of k O(k ) , was given in [10]. More recently, an algorithm achieving sample size s = k O(k) and runtime of poly(k) was analyzed in [3]. (Incidentally all of these papers use m = 2k − 1, and hence do not deal with verifying that the source is a k-coin distribution.) In [9], [10], the k-coin problem arises as the output of a reduction from the problem of identifying topic models, introduced in [11], [12]. A (pure) k-topic model is simply analogous to the k-coin problem with highly multisided coins. There has been ample work on learning pure and mixed topic models, under various restrictive assumptions on the model, and also without restrictions [9], [10], [13], [14]. The reductions of [9], [10] can be used in conjunction with our algorithm to reduce the sample size and runtime required to solve the topic model problem. This is discussed in Section 2.6. In [3], the k-coin problem arises as output of a reduction from the problem of inferring population histories (see the references therein). Our results improve both the sample size and the runtime complexity of the solution. We do note that the k-coin algorithm in [3] could have been used in conjunction with the reductions in [9], [10] to solve the topic model problem; the bounds derived this way would be worse than the bounds we prove in this paper. We also mention some generalizations of the k-coin problem that were considered in the literature. Most obvious is mixtures of k product distributions on {0, 1}m . That is, the formulation is the same as ours except that X1 , . . . , Xm are merely required to be independent, but not necessarily iid, conditional on the hidden variable H. This problem has been the focus of considerable research in the past two decades [15]–[20]. Clearly, in this case a larger m is no longer purely helpful, since the number of degrees of freedom of the problem also goes up with m. It should be noted, though, that the strongest results in this sequence, [19] and [20], do not address the problem of identifying the source model; rather, they learn a model which generates similar statistics. On the positive side, this task can sometimes be performed even under conditions where there is not enough information in the statistics for identification (i.e., when there are models with near-enough statistics that are far apart in, say, transportation distance); but on the negative side, since 11 these algorithms (as well as the algorithm in [10]) are forced to perform an exhaustive enumeration over a large grid of potential models, their computational efficiency does not much improve even when the statistics are known to sufficiently-good accuracy that only a very small-diameter (in transportation distance) set of models could generate them. The distinction between the “identification” and “learning” goals was made already in [16], who solved the identification problem for mixtures of k = 2 product distributions on {0, 1}m . Similar results for somewhat more general models were achieved at a similar time in [17]. The best result to date [20] 3 2 3 learns in time k k · mO(k ) , improving upon a previous result [19] of mO(k ) . √ The same paper [20] shows a lower bound of mΩ( k) on the sample size of the task. Beyond mixtures of product distributions, an even more complex but important class of source identification problems arises when the hidden variable (our “H”) may be just one of several such variables, and when a known directed causal structure exists among the observed variables (the “Xi ”). This is a very broad field of investigation and we point only to [5], [21] for background, and to [6] for an example of how (with some additional assumptions on the distributions of the Xi ) certain models can be handled. 2.2 Mixture Models and Other Definitions We specialize the definitions appearing in Chapter 1 to this particular problem to reduce notational clutter. Definition 6 (The k-coin model). A k-coin model M = (θ, π) is a mixture, with non-negative mixing weights π1 , . . . , πk , of k Bernoulli variables with success probabilities θ1 , . . . , θk , respectively. Definition 7 (m-snapshots of a k-coin model). Given a k-coin model M = (θ, π), an m-snapshot is a sample from the mixture of binomial distributions π1 Binomial(m, θ1 ) + . . . + πk Binomial(m, θk ). (The binomial is a sufficient statistic for m rv’s X1 , . . . , Xm because they are iid given the selected coin.) For a k-coin model, the moments defined in equation (2.3) can be written as follows where δθ is the Dirac measure at θ, P = π1 δθ1 + · · · + πk δθk , µi = k X j=1 θji πj . 12 Definition 8 (Separation for polynomials and mixtures). For a k-coin probability model M = (θ, π), define the separation by ζ(M) = mini̸=j |θi − θj |. For a degree k polynomial with roots β1 , . . . , βk ∈ C, define the root separation by mini̸=j |βi − βj |. (m) Definition 9. The rectangular Vandermonde matrix Vθ ated with a vector θ ∈ Ck is given by   1 1 1 ··· 1    θ1 θ2 θ3 · · · θk      (m) Vθ =  θ12 θ22 θ32 · · · θk2  .  . .. .. . . ..   .. . . . .    m m m m θ1 θ2 θ3 · · · θk ∈ R(m+1)×k associ- (k−1) We will denote the square Vandermonde matrix by Vθ := Vθ . Definition 10 (Hankel Matrix). The (k + 1) × (k + 1) Hankel matrix Hk+1 = Hk+1 (P) is defined by Hk+1 = [µi+j ]ki,j=0 . Note that if P is supported on a set of cardinality k (a.k.a. a k-coin distribution), then Hk+1 = k X (k) πj θj θjT = Vθ (k) T diag(π1 , . . . , πk )Vθ (2.4) j=1 where θjT = (1, θj , θj2 , θj3 , . . . , θjk ). This also shows that the Hankel matrix is positive semi-definite. Definition 11 (Polynomial associated with a vector). We associate to each P j vector q ∈ Rk a degree k − 1 polynomial q̂(x) = k−1 j=0 qj x . (For this reason we use zero indexing for the vector.) Definition 12. For a matrix M , let ∥M ∥2 denote the 2 → 2 operator norm of M . Thus, ∥M ∥2 = σmax (M ), the largest singular value of M . Definition 13. For a Hermitian matrix M , let λi (M ) denote the ith smallest eigenvalue of M . In particular λ1 (M ) is the smallest eigenvalue of M . Definition 14 (Euclidean projection onto a closed convex set). For a closed convex set S ⊆ Rk and any point x ∈ / S, the Euclidean projection of x onto S is ProjS (x) := arg miny∈S ∥y − x∥2 . This projection is unique. 13 2.3 Properties of Hankel Matrices We begin with some properties of Hankel matrices corresponding to finitely supported distributions that follow from results in Chihara [22]. (See Schmüdgen [23, Ch. 10] for a complete characterization.) For completeness, a proof is provided in Section 2.9. Lemma 15. Let P be a probability measure on [0, 1]. Then, 1. P is supported on a set of cardinality at most k iff Hk+1 is singular. 2. If the support of P is a set {θ1 , . . . , θk } ⊂ [0, 1] then the kernel of Hk+1 is Q spanned by the vector q ∈ Rk+1 where q̂(z) = ki=1 (z − θi ) is the unique monic polynomial with roots at the support of P. The above lemma is classic in the theory of Orthogonal Polynomials, but for completeness we provide a proof in Section 2.9. An essential contribution of this chapter is to strengthen the above lemma with a quantitative version: Theorem 16. Let P = (θ, π) be a k-coin distribution with separation ζ, and let Hk := Hk (P). Then, with I representing the k × k identity matrix and ⪰ being the PSD order,  2k−2 πmin ζ ·I . Hk ⪰ · k 8 Proof. We must show that for every monic degree k ′ ≤ k − 1 polynomial q̂, represented by q ∈ Rk , πmin q Hk q ≥ · k T  2k−2 ζ · ∥q∥22 . 8 Let β1 , β2 , . . . , βk′ be the roots (possibly complex) of the polynomial q̂, ordered so that |β1 | ≥ |β2 | ≥ · · · ≥ |βk′ |. Since q̂ is monic, we can write q̂(x) = Qk ′ j=1 (x − βj ). As the balls B(θi , ζ/2), i = 1, 2, . . . , k, are disjoint, by the pigeonhole principle, there exists an i ∈ {1, 2, . . . , k} such that B(θi , ζ/2) ∩ {β1 , β2 , . . . , βk′ } = ∅. The value of q̂ at θi is ′ k Y q̂(θi ) = (θi − βj ). j=1 14 ∥q∥22 k′ +1 There must be some ℓ ∈ {0, 1, 2, . . . , k ′ } such that |qℓ |2 ≥ . Notice that |qℓ | = |ek′ −ℓ (β1 , β2 , . . . , βk′ )|, where er is the r-th elementary symmetric polyP P nomial over k ′ variables. (e0 = 1, e1 = βi , e2 = i<j βi βj , etc. . . ) So, ek′ −ℓ ′
′ ∥q∥2 is the sum over k′k−ℓ ≤ 2k monomials, hence |β1 β2 · · · βk′ −ℓ | ≥ (√k′ +1)2 k′ . Eliminating from the product all the factors whose absolute value is below 2, ∥q∥2 we get that for some r ≤ k ′ − ℓ, |β1 β2 · · · βr | ≥ (√k′ +1)4 k′ . For j ∈ {1, 2, . . . , r}, |β | since |βj | ≥ 2 and θi ∈ [0, 1], it follows that |θi − βj | ≥ 2j . Also, by the definition of i we have that |θi − βj | > ζ/2 for all j ∈ {1, 2, . . . , k ′ }. Thus, we have that ! !  k′ −r r k′ Y Y |β1 β2 · · · βr | ζ |q̂(θi )| = |θi − βj | |θi − βj | ≥ 2r 2 j=1 j=r+1  k−1 ∥q∥2 1 ζ k′ √ · ≥ √ ′ ∥q∥2 . ζ ≥ ′ 8 ( k + 1)8k k Therefore, k X 1 q Hk q = πj · (q̂(θj )) ≥ πmin · (q̂(θi )) > πmin · · k j=1 T 2 2 Corollary 17. For a k-coin model (θ, π), λ2 (Hk+1 ) > πmin ·  2k−2 ζ · ∥q∥22 . 8
ζ 2k−2 . 16 Proof. By the Courant-Fischer-Weyl min-max principle, the smallest eigenvalue of Hk is given by minimizing the Rayleigh-Ritz quotient. Let q ̸= 0 T ′ kq be a minimizer of q qH T q . Let k be greatest such that qk ′ ̸= 0, and w.l.o.g. set qk′ = 1. Then by Theorem 16,  2k−2  2k−2 πmin ζ ζ q T Hk q λ1 (Hk ) = min T ≥ · ≥ πmin · , q̸=0 q q k 8 16 where the last inequality follows from observing that 1/k ≥ 1/22k−1 for k ≥ 2. Notice that Hk is a principal submatrix of Hk+1 . Therefore, by the Cauchy interlacing theorem (Theorem 40), λ2 (Hk+1 ) ≥ λ1 (Hk ). 2.4 The Empirical Moments We bound the sampling error as follows. Sample s coins and let each of the random variables hj , 0 ≤ j ≤ 2k, be the fraction of coins which came up “heads” exactly j times. Then by the additive deviation bound known as Hoeffding’s inequality [24], Pr(|hj − E(hj )| ≥ t) ≤ 2 exp(−2t2 s). Thus 15 Lemma 18. If we use s > 2t12 log(4k/δ) samples then with probability at least 1 − δ: ∀j, |hj − E(hj )| < t. We can convert between the normalized histogram h and the standard moments of the distribution by using the observation (Lemma 1 in [25]) that for any t ∈ R,
  n j X n j i
ti = t (1 − t)n−j . n × j j=i i This gives us a linear transformation for converting from h to the vector µ̃ = (µ̃0 , . . . , µ̃2k ). Define Pas ∈ R(2k+1)×(2k+1) (using zero-indexing) by  j  (2ki ) if j ≥ i Pasij = ( i )  0 otherwise; then µ̃ = Pas h. Lemma 19. ∥Pas∥2 ≤ 6k . Proof in Appendix 2.9. Now let µ = (µ0 , . . . , µ2k ) be the actual vector of moments of the distribution P. Lemma 20. For every ε > 0, using s = 2O(k) · ε12 · log(1/δ) samples gives us estimated moments µ̃ = (µ̃0 , . . . , µ̃2k ) satisfying ∥µ̃ − µ∥∞ ≤ ε with probability at least 1 − δ. Proof. Follows directly from Lemma 19 and Lemma 18. Given an s-sample as above with empirical moments µ̃0 , µ̃1 , . . . , µ̃2k , denote by H̃k+1 the empirical Hankel matrix,   µ̃0 µ̃1 µ̃2 · · · µ̃k   µ̃1 µ̃2 µ̃3 · · · µ̃k+1   H̃k+1 =  . (2.5) .. .. ..  .. . . . . . .  . µ̃k µ̃k+1 µ̃k+2 · · · µ̃2k Corollary 21. For every ε > 0, using s = 2O(k) · ε12 · log(1/δ) samples, we can obtain an empirical Hankel matrix satisfying H̃k+1 − Hk+1 ≤ ε with probability 2 at least 1 − δ. Proof. We have H̃k+1 − Hk+1 ≤ 2 ε Now use Lemma 20 with k+1 . H̃k+1 − Hk+1 F ≤ (k + 1) · ∥µ̃ − µ∥∞ . 16 2.5 Learning the Source In this section, we define our learning algorithm, and we state and prove our main result and applications. The auxiliary lemmas are stated and proved in Section 2.7. The algorithm is specified given k, lower bounds on the source parameters ζ and πmin , the empirical histogram h, and a parameter γ controlling the output accuracy. The algorithm is a straightforward implementation of Prony’s method with numerical approximations where needed. See Algorithm 1 for the full description of the algorithm (where the parameter for probability of success, 1 − δ, has been suppressed in favor of a constant “0.99”). On line 5, we take the output of Prony’s method and project the roots back into [0, 1] and on line 6 in RectifyWeights we do a similar correction to ensure that the weights are non-negative and still sum to one. While simple, we defer the implementation of RectifyWeights to Algorithm 2 in Appendix 2.8. Procedure LearnPowerDistribution(k, ζ, πmin , H̃k+1 , γ): ε1 −approx 1 v ←−−−−−− arg min{v T H̃k+1 v : v T v = 1} // ε1 = πmin · 2−γ · (ζ/16)2k 2 3 ε2 −approx β̃1 , β̃2 , . . . , β̃k ←−−−−−− roots(v̂) 1 · 2−γ · (ζ/2)k // ε2 = 6k 4 5 6 θ̃1 , θ̃2 , . . . , θ̃k ← Proj[0,1] (β̃1 ), . . . , Proj[0,1] (β̃k ) π̃ ← RectifyWeights(Vθ̃−1 µ̃) // see Algorithm 2 on page 24 7 return M̃ ← (θ̃, π̃) Procedure LearnCoinMixture(k, ζ, πmin , γ, h): 9 µ̃ ← Pas h 10 H̃k+1 ← Hankel(µ̃) 11 return M̃ ← LearnPowerDistribution(k, ζ, πmin , H̃k+1 , γ) Algorithm 1: The main learning algorithm. 8 Theorem 22. Let M = (θ, π) be a k-coin model with separation ζ = ζ(M). Given an empirical Hankel matrix H̃k+1 satisfying H̃k+1 − Hk+1 ≤ πmin ·2−γ ·(ζ/16)4k , 2 the procedure LearnPowerDistribution in Algorithm 1 outputs a model M̃ = (θ̃, π̃) satisfying θ − θ̃ ∞ , ∥π − π̃∥∞ ≤ 2−γ 17 −1 using O(k 2 log k + k log2 k · log(log ζ −1 + log πmin + γ)) arithmetic operations. Proof. Throughout the proof, we make no attempt to optimize the absolute constants that are used. Let u1 denote the unit vector spanning the kernel of Hk+1 , and let v1 denote the eigenvector corresponding to the smallest eigenvalue of H̃k+1 . Also, let ε0 > 0 be a sufficiently small constant, to be determined later. The analysis of LearnPowerDistribution can be broken down into steps, each of which degrades the initial accuracy obtained for the Hankel matrix. The outline is as follows. The auxiliary claims and proofs appear mostly in Section 2.7. 1. As H̃k+1 − Hk+1 2 ∥u1 − v1 ∥2 < ≤ πmin · 2−γ · (ζ/16)4k , by Lemma 25, p 1 2(k + 1) · 2−γ · (ζ/16)2k < · 2−γ · (ζ/8)2k . 2 2. We use Lemma 26 with ε = πmin · 2−γ · (ζ/16)2k , which satisfies the conditions of the lemma. We compute v ∈ Rk+1 such that ∥v − v1 ∥2 ≤ ε < 1 −γ · 2 · (ζ/16)2k , 2 −1 using O(k 2 log k + k log2 k · log(log ζ −1 + log πmin + γ)) arithmetic operations. 3. As ∥u1 − v1 ∥2 , ∥v − v1 ∥2 < 12 · 2−γ · (ζ/16)2k , we have that ∥u1 − v∥2 < 2−γ · (ζ/16)2k . So, by Lemma 27, √ ∥q − r∥∞ < 2k · k + 1 · 2−γ · (ζ/16)2k < 2−γ · (ζ/8)2k , where q := u1 / |(u1 )k | , r := v/ |(u1 )k |. 1 · 2−γ · (ζ/2)2k−1 , by Lemma 28 we 4. As ∥q − r∥∞ < 2−γ · (ζ/8)2k ≤ 24k(k+1) have that 1 d(θ, β) ≤ · 2−γ · (ζ/2)k 6(k + 1) (where θ is the vector of roots of q̂ and β is the vector of roots of r̂ and d is the matching distance, defined in Lemma 28). 1 5. We use Corollary 43 with ρ = ε = 6(k+1) · 2−γ · (ζ/2)k , which satisfy the conditions of the corollary. Thus, we can compute biases θ̃1 , . . . , θ̃k satisfying θ̃ − θ ∞ ≤ρ+ε≤ 1 · 2−γ · (ζ/2)k , 3(k + 1) 18 using O(k log2 k · (log(log ζ −1 + γ) + log2 k)) arithmetic operations. 6. Finally, line 6 can be executed in the time it takes to invert the Vandermonde matrix Vθ̃ (i.e., O(k 2 ) arithmetic operations, for instance using Parker’s Algorithm [26]; by Lemma 34, the procedure RectifyWeights takes O(k) operations). By Corollary 35, as ∥θ̃ − θ∥∞ , ∥µ̃ − µ∥∞ ≤ 1 · 2−γ · (ζ/2)k (the guarantee for µ̃ is implied with plenty of room 3(k+1) to spare by our assumption on the sample), we have ∥π̃ − π∥∞ ≤ 2−γ . Corollary 23. Let M = (θ, π) be a k-coin model with separation ζ = ζ(M). For any γ ≥ 1, the procedure LearnCoinMixture in Algorithm 1 uses a histogram h −2 for a sample of 2k-snapshots of size s = πmin · 2O(k+γ) · ζ −O(k) · log δ −1 , and outputs a model M̃ = (θ̃, π̃) satisfying θ − θ̃ ∞ , ∥π − π̃∥∞ ≤ 2−γ with probability at least 1 − δ. After sampling, LearnCoinMixture computes the −1 approximate model M̃ using O(k 2 log k + k log2 k · log(log ζ −1 + log πmin + γ)) arithmetic operations. Proof. By Lemma 20 and Corollary 21 we can achieve H̃k+1 − Hk+1 ≤ 2 −2 πmin · 2−γ · (ζ/16)4k for a sample of size s = πmin · 2O(k+γ) · ζ −O(k) · log δ −1 , with probability at least 1 − δ. The theorem immediately follows from Theorem 22. Notice that the proof actually gives a stronger guarantee for θ̃ − θ ∞ , which is relative to (ζ/2) . We can get a relative guarantee ∥π̃ − π∥∞ ≤ πmin · 2−γ by −2 increasing the sample size by a factor of πmin . k Corollary 24. Let W (M, M̃) denote the Wasserstein distance between models M and M̃ (viewed as metric measure spaces on [0, 1]). Then, W (M, M̃) ≤ (k+1)·2−γ with probability at least 0.99. Proof in Appendix 2.9. 2.6 Implications for Topic Models Corollary 23 improves upon the upper bound of Theorem 5.1 in [9], which n o 2 O(k) uses a sample of (2k − 1)-snapshots of size max (2/ζ) , (2γ k)O(k ) to achieve accuracy 2−γ with high probability, using runtime of O(k c ) arithmetic 19 operations, for a relatively large constant c (in particular, the algorithm solves a convex quadratic program whose representation uses k 3 bits). Corollary 23 also improves upon the upper bound of [10].4 That algorithm uses a sample 2 size comparable to ours, but requires runtime (2γ k)O(k ) to achieve accuracy 2−γ with high probability. These improvements imply immediately a similar improvement for learning pure k-topic models, using known reductions from k-topic models to k-coin models. The reduction in Theorem 4.1 of [9] uses a sample of 1- and 2

−1 snapshots of size O n · poly log n, k, πmin , ζ −1 , 2γ , and runtime polynomial in the sample size, to reduce the of the k-coin n problem to solving k instances o −O(1) −O(k) γ γ problem with accuracy min (2 k/πmin ζ) , (2 k) . The reduction in [10] 5 uses a sample of 1- and 2-snapshots of size poly (n, k, 2γ ), and runtime polynomial in the sample size, to reduce the problem to solving at most k instances of the k-coin problem with accuracy (2γ k)−O(k) . Notice that solving the k-coin outcome of either one of the two reductions using either one of 2 the two previous algorithms requires a sample size of at least k O(k ) (on account of the required accuracy). Our algorithm enables a solution to the outcome of these reductions using a sample size of k O(k) (and total runtime of O(k 3+o(1) )). We note that the accuracy in [9], [10] is stated in terms of Wasserstein distance, which is a weaker guarantee than the one we use here (see Corollary 24). 2.7 Analysis In this section we prove the lemmas that are needed in the proof of Theorem 22 and Corollary 23. We have to cope with the fact that roots of polynomials (and even, generally, of polynomials with well-separated roots), are notoriously ill-conditioned in terms of the polynomial coefficients [27]. For this reason we will be developing bounds specifically adapated to our situation. We begin with an estimate on the accuracy of the recovered kernel of the Hankel matrix. Lemma 25. Let P be any k-coin distribution with separation ζ. Then, for every
ζ 2k ε < πmin · 16 the following holds. Suppose that H̃k+1 − Hk+1 ≤ ε. Let u1 be 2 the unit vector in the kernel of Hk+1 and let v1 be the unit eigenvector corresponding to 4 5 See Theorem 3.9 in the ArXiv version See Theorem 6.1 in the ArXiv version. 20 λ1 (H̃k+1 ) (chosen so that uT1 v1 ≥ 0). Then ∥u1 − v1 ∥2 < Proof in Appendix 2.9.  p 2(k + 1) ·  2k 16 ζ ε · πmin .  Recall that λ1 (H̃k+1 ), v1 is an eigenpair of H̃k+1 . We need to compute a good approximation of v1 . This can be done using the following lemma. The result is implied by the algorithm of Pan and Chen (Theorem 1.2 of [28]). Extracting our lemma from the result in that paper is somewhat involved and we provide in Appendix 2.9 a brief outline of the argument (in particular, the parts that are not spelled out in that paper). Lemma 26. For every ε such that 0 < ε < 12 min{λ2 (H̃k+1 )−λ1 (H̃k+1 ), 1}, we can
compute a unit vector v satisfying ∥v − v1 ∥2 < ε using O k 2 log k + k log2 k log log(1/ε) arithmetic operations. We need to show that our computed first eigenvector of the empirical Hankel matrix is close to the kernel eigenvector of the true Hankel matrix. Lemma 27. Let P be any k-coin distribution with separation ζ. Let u1 be a unit vector in the kernel of Hk+1 . Let v be a unit vector satisfying ∥u1 − v∥2 < ε for some √ ε > 0. Let q = u1 / |(u1 )k | and let r = v/ |(u1 )k |. Then ∥q − r∥∞ < 2k k + 1 · ε. Proof in Appendix 2.9. We are going to use the roots of the polynomial r̂ as our guessed coin biases (after projecting the roots back to [0, 1]). We first need to show that the roots of q̂ are well-behaved with respect to perturbations of q so that when q and r are close the roots of q̂ are close to the roots of r̂. Lemma 28. Let q ∈ Rk+1 be the vector representing a degree-k monic polynomial with roots θ1 , θ2 , . . . , θk contained in [0, 1]. Let ζ be the root separation  for (ζ/2)k k+1 q̂. Let r ∈ R represent another degree-k polynomial. Let ε ∈ 0, 4k . If r satisfies ∥q − r∥∞ ≤ ε, then the (possibly complex) roots β1 , β2 , . . . , βk of r̂ sat4kε isfy d(θ, β) ≤ (ζ/2) k−1 where d(θ, β) is the optimal matching distance defined by d(θ, β) := minσ∈Sk maxi θi − βσ(i) . Proof. Fix any root θi of q̂, and consider the ball   4kε Bi = B θi , (ζ/2)k−1 21 ζ 4kε in the complex plane. By assumption, (ζ/2) k−1 < 2 , so there are no other roots of q̂, aside from θi , in Bi . Moreover, for any x ∈ Bi , and for any j ̸= i, we have that |x − θj | ≥ ζ2 . Thus for every x ∈ ∂Bi , we have Y |q̂(x)| = (x − θi ) (x − θj ) > j̸=i 4kε (ζ/2)k−1  k−1 ζ = 4kε. 2 On the other hand, we also have that Bi ⊂ B(0, (2k − 1)/(2k − 2)), as 1 . Therefore, |x| ≤ 2k−1 , and thus θ1 , . . . , θk ∈ [0, 1] and ζ ≤ k−1 2k−2 |q̂(x) − r̂(x)| = k X j=0 j (qj − rj )x ≤ k X j=0 j |qj − rj | · |x| ≤ (k + 1) ·  2k − 1 2k − 2 k · ∥q − r∥∞ ≤ 4kε. By Rouché’s theorem (Theorem 41), we conclude that there is exactly one zero of r̂ in Bi and the matching distance bound follows immediately. Our reconstructed coin biases will be denoted θ̃1 , θ̃2 , . . . , θ̃k . We compute these biases by finding the roots of v̂ (approximately), and then by projecting these roots onto the unit interval. To find the approximate roots we can use Corollary 43 in Appendix 2.10. Once we have recovered the parameters θ̃1 , . . . , θ̃k , we need to use those to recover mixture weights. This sequence of steps—first solving (approximately) for the roots, then for the mixture weights—is the essence of Prony’s method [7], [1] §9.4, [8]. In Appendix 2.8, we will show that this recovery can be done by solving a linear system without paying too great a price in terms of accuracy. 22 2.8 Computing the Weights Recovering the Weights We will begin by stating results characterizing the condition number of a Vandermonde matrix under perturbations that preserve the Vandermonde structure. Lemma 29 (Operator norm bound for a Vandermonde inverse; equation 3.2 k Q in [29]). Let θ ∈ Rk be entry-wise non-negative, and let q(z) = (z − θi ). Then i=1 Vθ−1 ∞ = |q(−1)| . mini {(1 + θi ) |q ′ (θi )|} Claim 30. For roots θ1 , . . . , θj satisfying |θi − θj | ≥ ζ, we have Vθ−1 ∞ ≤ 2k /ζ k−1 . Proof. We apply Lemma 29 and observe that |q(−1)| ≤ 2k and q ′ (θi ) ≥ ζ k−1 . We define the derivative matrix of the Vandermonde matrix by interpreting each entry as the evaluation of a polynomial at a point, [Va ]ij = pi (aj ), where pi (t) = ti . Then [Va′ ]ij = p′i (aj ) = iai−1 j . We will now define the condition number of the system,  cond∞ (a, b) := lim sup ε→0 ∥∆a∥ ≤ε ∞ ∥∆b∥∞ ≤ε ∥∆x∥∞ ε  V (a + ∆a)(x + ∆x) = b + ∆b . (2.6) We will utilize a bound from [30]. After instantiating the theorem with the parameters relevant to our problem, the bound is the following: Theorem 31 (Theorem 2.2 of [30]). cond∞ (a, b) ≤ ∥Va−1 ∥∞ +∥Va−1 Va′ diag(x)∥∞ . Lemma 32. Let θ ∈ [0, 1]k , and let w ∈ Rk be a probability distribution over [k]. Let µ = Vθ w. If ζ ≤ mini̸=j |θi − θj |, cond∞ (θ, µ) ≤ (k + 1)2k /ζ k−1 . 23 Proof. We observe that Vθ−1 Vθ′ diag(π) ∞ ≤ Vθ−1 ∞ ∥Vθ′ diag(π)∥∞ ≤ 2k /ζ k−1 ∥Vθ′ diag(π)∥∞ k = 2 /ζ k−1 max (i + 1) i∈[k−2] k X θji πj j=1 ≤ k2k /ζ k−1 . Applying the bound of Theorem 31 gives the conclusion. Lemma 33. Let θ ∈ [0, 1]k and let w ∈ Rk be a probability distribution over [k]. Let µ = Vθ w, and ζ ≤ mini̸=j |θi − θj |. Then π ′ := Vθ̃−1 µ̃ satisfies ∥π ′ − π∥∞ ≤ o n (k + 1)2k , ∥µ̃ − µ∥ max θ̃ − θ ∞ . ζ k−1 ∞ Proof. This follows from Lemma 32 and the definition of the condition number. P Lemma 34. Given any weights π ′ ∈ Rk satisfying ki=1 πi′ = 1, the procedure RectifyWeights(π ′ ) outputs in time O(k) a weight vector π̃ ∈ [0, 1]k satisfying the following conditions (i) Pk i=1 π̃i = 1. (ii) ∥π̃ − π∥∞ ≤ (k + 1) ∥π ′ − π∥∞ . Corollary 35. Letting π̃ ∈ [0, 1]k be the output of RectifyWeights(π ′ ) where π ′ is as in Lemma 33, n o (k + 1)2 2k ∥π̃ − π∥∞ ≤ max θ̃ − θ , ∥µ̃ − µ∥∞ . ζ k−1 ∞ P Proof. Notice that the first equation in the linear system defining π ′ is ki=1 πi′ = 1T π ′ = µ̃0 = 1. Thus, π ′ satisfies the hypothesis of Lemma 34 and the conclusion follows. 24 Procedure RectifyWeights(π ′ ): {i | πi′ ≥ 0} 1 I − ← {iP| πi′ < 0}, I + ←P 2 W − ← i∈I − πi′ , W + ← i∈I + πi′ for i = 1, . . . , k do ( 0  if i ∈ I −  3 π̃i ← − if i ∈ I + πi′ 1 + W W+ end 4 return π̃ Algorithm 2: Correct computed weights to be non-negative and sum to one. Fixing the weights The weights produced prior to calling RectifyWeights sum to 1 (since they P satisfy the equation 1 = µ0 = j π̃j ) but they may lie outside [0, 1]. To correct this, we simply make all negative weight zero and scale all non-negative weights so that the sum of the weights doesn’t change. Proof of Lemma 34. Note that in Algorithm 2, I − denotes the indices of the negative weights, and I + the positive weights. W − and W + denote the sums of the weights in the corresponding set of indices. We will now analyze π̃. First, note that we maintain property (i): k X   W− π̃i = 1+ + W i=1 i∈I +   W− + =W 1+ + W X πi′ = W + + W − = 1. Now we show that the weights are non-negative. Trivially, π̃i ≥ 0 for i ∈ I − . For i ∈ I + , W+ = 1 − W− = 1 + |W − | ≥ |W − |. − So πi′ (1 + W ) ≥ 0 if i ∈ I + as well. W+ We now prove (ii). We know that the true weights π lie in [0, 1], so increasing the negative weights to 0 only moves them closer to their true values. Thus, 25 we have |π̃i − πi | ≤ |πi′ − πi | for all i ∈ I − . We observe that W − ≤ ∥π ′ − π∥1 ≤ k ∥π ′ − π∥∞ and then that |π̃i − πi′ | =
πi′ /W + W − ≤ k ∥π ′ − π∥∞ . | {z } ≤1 It follows that ∥π̃ − π ′ ∥∞ ≤ k ∥π ′ − π∥∞ . Now we can apply the triangle inequality to get that ∥π̃ − π∥∞ ≤ ∥π̃ − π ′ ∥∞ + ∥π ′ − π∥∞ ≤ (k + 1) ∥π ′ − π∥∞ . To see that the runtime is O(k) we observe we can compute I − and I + in linear time and likewise for W − and W + . Each subsequent computation of π̃i takes constant time. 2.9 Deferred Proofs Proof of Lemma 15. (Part 1.) By Equation (2.4), the rank of Hk+1 for a t-coin distribution is at most t, and that implies that if t ≤ k, then Hk+1 is singular. So consider a distribution P on [0, 1] that has positive mass at k + 1 points or more. Let q ∈ Rk+1 be a non-zero vector. We have !2 Z 1 X Z 1 k T j q Hk+1 q = qj θ dP(θ) = q̂ 2 (θ) dP(θ). 0 0 j=0 There are at most k points in [0, 1] where the polynomial q̂ evaluates to 0, and the total P measure of those points is less than 1. Thus, q T Hk+1 q > 0, so Hk+1 is positive definite. (Part 2.) Since Hk+1 is symmetric, its kernel is spanned by q s.t. q T Hk+1 q = 0. In order for the above integral to evaluate to zero over P, we need that q̂ 2 (θ) = 0 for each point θ ∈ supp(P). As q̂ is of degree ≤ k, it is necessarily a Q scalar multiple of ki=1 (z − θi ). Proof of Lemma 19. We first observe that 

−1

1 2k −1 0 2k ··· 0 0
2k
−1  0 0 1  0 ··· 1 1   0 0 ··· Pas =   . . ..  .. .. .  0 0 ···

2k−1 2k −1 0
0
2k−1 2k −1 1
1
2k−1 2k −1 2 2 .. . 0

2k 2k −1 0
0
 2k 2k −1  1
1
 2k 2k −1   2 2  .. . 2k 2k
2k
−1 2k    26 which can be factored to obtain 
−1 2k 0 0 ···
−1  0 2k  0 0 ··· 1 
−1  2k 0 ··· Pas =  0 2  . . . ..  .. .. .. .  0 0 0 ··· Now 0 
0 1 0
1 1
 0  0   0  0  . ..   . .   .
2k −1 0 2k 0      −1 2k 2k 2k diag , ,..., 0 1 2k 0 .. . ··· ··· ··· .. . 0 ··· 2k−1 0
2k−1 1
2k−1 2
.. . 0
 2k 0
 2k  1 
2k  2 . ..  . 
 2k 2k ≤ 1. 2 The Frobenius norm of the latter matrix is ! !   1/2 j  2 1/2 2k X 2k   1/2 X X j 4k 2j = ≤ 2k ≤ (2k 42k )1/2 ≤ 6k i 2k j j=0 i=0 j=0 for k ≥ 2. Using the sub-multiplicativity of the operator norm and the fact that the Frobenius norm upper bounds the operator norm, we get that ∥Pas∥ ≤ 6k , as desired. Proof of Corollary 24:. Each θi can be matched to its corresponding θ̃i up to weight min{πi , π̃i }. The additional |πi − π̃i | must move an additional distance of at most 1. This gives W (M, M̃) ≤ ≤ k X i=1 k X |θi − θ̃i | · min{πi , π̃i } + k X |πi − π̃i | i=1 2−γ min{πi , π̃i } + i=1 k X 2−γ i=1 −γ ≤ (k + 1) · 2 , using Theorem 22 and the fact that Pk i=1 min{πi , π̃i } ≤ Pk i=1 πi = 1. Proof of Lemma 25. By Weyl’s inequality, we have that λ1 (H̃k+1 ) ≤ ε. By Corollary 17, the eigengap λ2 (Hk+1 ) − λ1 (H̃k+1 ) is at least  πmin ζ 16 2k−2  − πmin ζ 16 2k  > πmin ζ 16 2k . 27 Now we can use Corollary 38 to obtain 2  Hk+1 − H̃k+1  uT1 v1 = uT1 v1 ≥ 1 − F λ2 (Hk+1 ) − λ1 (H̃k+1 )  4k  2 16 ε = 1 − (k + 1) · · . ζ πmin 1/2  2 >1− (k + 1) · ε2
ζ 4k 2 wmin 16 Since ∥u1 − v1 ∥22 = 2 − 2uT1 v1 we get that ∥u1 − v1 ∥22 < 2(k + 1) ·  16 ζ 4k  2 ε · . πmin Proof sketch of Lemma 26. We follow the outline in the papers by Pan, Chen, and Zheng [28], [31]. As H̃k+1 is a Hankel matrix, a similarity transformation A = T H̃k+1 T −1 , where A is tridiagonal, can be computed in time O(k 2 log k). The characteristic polynomial cA (x) of A can then be computed in time O(k). Then, a root λ̃ that satisfies |λ̃ − λ1 (H̃k+1 )| < ε2 can be com
puted in time O (k log2 k)(log log(1/ε) + log2 k) (see Theorem 42; note that ∥A∥2 = ∥H̃k+1 ∥2 , thus it is trivially upper bounded by (k + 1)2 ). Next, proceed to compute v as follows. Pick an initial guess v (0) uniformly at random on the unit sphere (i.e., from the unit Haar measure on the sphere). We need v1T v (0) > √1k , which happens with constant probability. To boost the success probability to 1 − δ, we can repeat the entire process O(log(1/δ)) times. For constant δ, this does not affect the asymptotic bound. We compute v (1) , v (2) , . . . using the inverse power iteration (see, for instance,  Chapter 4 (t) in [32]): Solve for ṽ the system of linear equations λ̃I − H̃k+1 ṽ (t) = v (t−1) , (t) then set v (t) = ∥ṽṽ(t) ∥2 . As H̃k+1 is a Hankel matrix, this can be done using O(k 2 ) arithmetic operations. How many iterations are needed?—It is known that if λ1 (H̃k+1 ) is the unique eigenvalue of H̃k+1 that is closest to λ̃, and if v1T v (0) > 0, then tan θ(t) ≤ ρ · tan θ(t−1) , where θ(t) is the angle between 1 (H̃k+1 )| v1 and v (t) , and ρ = |λ̃−λ|λ̃−λ , where λ2 is an eigenvalue of H̃k+1 that is | 2 1 (H̃k+1 )| ε second-closest to λ̃. Notice that in our case ρ = min|λ̃−λ|λ̃−λ < ε−ε 2 < 2ε. i>1 i (H̃k+1 )| √ As tan θ(0) ≤ k, after t = O(log1/2ε k) iterations, we have tan θ(t) < ε. This implies that v (t) − v1 2 < ε. 2 28 Proof of Lemma 27. Notice that q and r are well-defined, as (u1 )k ̸= 0 by the second part of Lemma 15. Now each of the coefficients of q can be bounded by   k |qi | = |ek−i (θ1 , . . . , θk )| ≤ i where er is the r-th elementary symmetric polynomial over k variables. Now √ √ ∥q∥2 ≤ k + 1 · ∥q∥1 ≤ 2k k + 1. Since |(u1 )k | · ∥q∥2 = ∥u1 ∥2 = 1, we have |(u1 )k | ≤ 2k √1k+1 , and √ √ ∥q − r∥∞ ≤ ∥q − r∥2 ≤ 2k k + 1 · ∥u1 − v∥2 < 2k k + 1 · ε, as stipulated. Proof of Lemma 32. We observe that Vθ−1 Vθ′ diag(π) ∞ ≤ Vθ−1 ∞ ∥Vθ′ diag(π)∥∞ ≤ 2k /ζ k−1 ∥Vθ′ diag(π)∥∞ = 2k /ζ k−1 max (i + 1) i∈[k−2] k ≤ k2 /ζ k−1 k X θji πj j=1 . Applying the bound of Theorem 31 gives the conclusion. 2.10 Useful Theorems Consider two n × n Hermitian matrices A, B, with spectral decomposiPn Pn T T tions A = i=1 κi ui ui and B = i=1 λi vi vi , where the eigenvalues of both matrices are sorted in increasing order (i.e., κ1 ≤ κ2 ≤ · · · ≤ κn and λ1 ≤ λ2 ≤ · · · ≤ λn ). Also, let P = B − A and let ρ1 ≤ ρ2 ≤ · · · ≤ ρn be the eigenvalues of P in increasing order. Theorem 36 (Weyl’s inequality). For every i ∈ {1, 2, . . . , n}, κi + ρ1 ≤ λi ≤ κi + ρn . Theorem 37 (Davis-Kahan sin Θ theorem). Using the above definitions, let i0 , i1 be integers such that 1 ≤ i0 ≤ i1 ≤ n, and let g = inf{|κ − λ| : κ ∈ [κi0 , κi1 ] ∧ λ ∈ (−∞, λi0 −1 ] ∪ [λi1 +1 , +∞)}, 29 where we define λ0 = −∞ and λn+1 = ∞. Then, ∥sin Θ(U, V )∥F ≤ ∥P ∥F , g where U (V , respectively) is the n×i1 −i0 +1 matrix whose columns are ui0 , . . . , ui1 (vi0 , . . . , vi1 , respectively), Θ(U, V ) is the i1 − i0 + 1 × i1 − i0 + 1 diagonal matrix whose i-th diagonal entry is the i-th principal angle between the column spaces of U and V , and sin Θ(U, V ) is the diagonal matrix derived by applying the function sin entrywise to Θ(U, V ). The same inequality holds if the Frobenius norm is replaced by any orthogonally invariant norm, e.g., an operator norm ∥ · ∥op . Corollary 38. Using the same definitions, s |uT1 v1 | ≥ 1− ∥P ∥2 . |κ1 − λ2 |2 ∥ Proof. Take i0 = i1 = 1. By Theorem 37, | sin θ(u1 , v1 )| ≤ |κ∥P . The corollary 1 −λ2 | p 2 T follows as |u1 v1 | = | cos θ(u1 , v1 )| = 1 − sin θ(u1 , v1 ). Theorem 39 (Courant-Fischer-Weyl min-max principle). For every i = 1, 2, . . . , n,     T x Bx λi = minn max : x ̸= 0 : dim(U ) = i U ⪯R x∈U xT x     T x Bx = maxn min : x ̸= 0 : dim(U ) = n − i + 1 . U ⪯R x∈U xT x Let C be an m × m Hermitian matrix with eigenvalues ν1 ≤ ν2 ≤ · · · ≤ νm , where m ≤ n. Theorem 40 (Cauchy’s interlacing theorem). If C = Π∗ BΠ for an orthogonal projection Π, then for all i = 1, 2, . . . , m it holds that λi ≤ νi ≤ λn−m+i . Theorem 41 (Rouché’s theorem). Let f and g be two complex-valued functions that are holomorphic inside a region R with a closed simple contour ∂R. If for every x ∈ ∂R we have that |g(x)| < |f (x)|, then f and f + g have the same number of zeros inside R, counting multiplicities. Theorem 42 (Pan’s Algorithm: Theorem 1.1 of [33]). Given a monic degree k polynomial p̂ with roots ρ1 , . . . , ρk ∈ B(0, 1) and an accuracy parameter γ > 1, we can compute approximate roots ρ̃1 , . . . , ρ̃k satisfying ∥ρ − ρ̃∥∞ ≤ 2−γ in time O(k log2 k · (log γ + log2 k)). 30 Qk Corollary 43. Let q ∈ Rk+1 represent the polynomial q̂(z) = i=1 (z − θi ) where θ1 , . . . , θk ∈ [0, 1] are ζ-separated, and let r ∈ Rk+1 represent a polynomial of degree k with roots β1 , . . . , βk satisfying d(θ, β) ≤ ρ < ζ/2. For every ε ∈ (0, ζ/2 − ρ), we can reconstruct biases θ̃1 , . . . , θ̃k satisfying θ̃ − θ ≤ ρ + ε using O(k log2 k · (log log(1/ε) + log2 k)) arithmetic operations. ∞ Proof. We will first approximate the roots β̃1 , . . . , β̃k of r̂ using Theorem 42.
Since the roots of r̂ are in B 0, 2k−1 instead of B(0, 1), we will actually find 2k−2
2k−2 the roots of r̂ 2k−1 z and then multiply by 2k−1 to get the roots of r̂ up to 2k−2 2 2 accuracy ε in time O(k log k · (log log(1/ε) + log k)). (Notice that in order to get the desired accuracy we need to run Pan’s algorithm to get the rescaled roots to within distance 2k−2 · ε; this doesn’t matter for the purposes of 2k−1 runtime.) Our output is θ̃i := Proj[0,1] (β̃i ) for i = 1, . . . , k, where we label the roots β̃1 , . . . , β̃k by the permutation achieving the matching distance, i.e., the ordering of coordinates so that ∥θ − β∥∞ = d(θ, β). Now θi − θ̃i ≤ θi − ℜ(β̃i ) ≤ |θi − ℜ(βi )|+ ℜ(βi ) − ℜ(β̃i ) ≤ |θi − βi |+ βi − β̃i ≤ ρ+ε. 31 Chapter 3 SUFFICIENT CONDITIONS FOR THE IDENTIFIABILITY OF MIXTURES OF PRODUCTS [1] S. L. Gordon and L. J. Schulman, “Hadamard extensions and the identification of mixtures of product distributions,” IEEE Transactions on Information Theory, vol. 68, no. 6, pp. 4085–4089, 2022. doi: 10.1109/ TIT.2022.3146630, 3.1 Introduction The Hadamard product for row vectors u = (u1 , . . . , uk ), v = (v1 , . . . , vk ) is the mapping ⊙ : Rk × Rk → Rk given by u ⊙ v := (u1 v1 , . . . , uk vk ). The identity for this product is the all-ones vector 1. We associate with vector v the linear operator v⊙ = diag(v), a k × k diagonal matrix, so that u · v⊙ = v ⊙ u. Throughout this chapter m is a real matrix with row set [n] := {1, . . . , n} and column set [k]; write mi for a row and mj for a column. As a matter of notation, for a matrix Q and nonempty sets R of rows and C of columns, let Q|C R be the restriction of Q to those rows and columns (with either index omitted if all rows or columns are retained). Definition 44. The Hadamard Extension of m, written H(m), is the 2n × k matrix with rows mS for all S ⊆ [n], where, for S = {i1 , . . . , iℓ }, mS = Q mi1 ⊙ · · · ⊙ miℓ ; equivalently mjS = i∈S mji . (In particular m∅ = 1.) This construction originated recently in learning theory [20], [34] where it is arises naturally and unavoidably when we wish to perform source identification (i.e., parameter estimation) given data from a mixture (convex combination) of k product distributions on n binary random variables. We explain the connection further in Section 3.2. Motivated by this application, we are interested in the following two questions: 32 (1) If H(m) has full column rank, must there exist a subset R of the rows, of bounded size, such that H(m|R ) has full column rank? (2) In each row of m, assign distinct colors to the distinct real values. Is there a condition on the coloring that ensures H(m) has full column rank? In answer to the first question we show: Theorem 45. If H(m) has full column rank then there is a set R of no more than k − 1 of the rows of m, such that H(m|R ) has full column rank. Considering the more combinatorial second question, observe that if m possesses two identical columns then the same is true of H(m), and so the latter cannot have full column rank. Extending this further, suppose there are three columns C in which only one row r has more than one color. Then rowspace H(m|C ) is spanned by 1|C and r|C , so again H(m) cannot have full column rank. Motivated by these necessary conditions, set: Definition 46. For a matrix Q let NAE(Q) be the set of nonconstant rows of Q (NAE=“not all equal”); let ε(Q|C ) = | NAE(Q|C )| − |C|; and let ε(Q) = minC̸=∅ ε(Q|C ). If ε(Q) ≥ −1 we say Q satisfies the NAE condition. In answer to the second question we have the following: Theorem 47. If m satisfies the NAE condition then (a) There is a restriction of m to some k − 1 rows R such that ε(m|R ) = −1. (b) H(m) is full column rank. (As a consequence also H(m|R ) is full column rank.) Apparently the only well-known example of the NAE condition is when m contains k−1 rows which are identical and whose entries are all distinct. Then the vectors m∅ , m{1} , m{1,2} , . . . , m{1,...,k−1} form a nonsingular Vandermonde matrix. This example shows that the bound of k − 1 in (a) is best possible. For another example in which the NAE condition ensures that rank H(m) = k, take the (k − 1)-row matrix with mji = 1 for i ≤ j and mji = 1/2 for i > j. Here the NAE condition is only minimally satisfied, in that for every ℓ ≤ k there are ℓ columns C s.t. ε(m|C ) = −1. 33 For k > 3 the NAE condition is no longer necessary in order that H(m) have full column rank. E.g., for k = 2ℓ , the ℓ × k “Hamming matrix” mji = (−1)ji where j is an ℓ-bit string j = (j1 , . . . , jℓ ), forms H(m) = the Fourier transform for the group (Z/2)ℓ (often called a Walsh or Hadamard transform), which is invertible. Furthermore, for k ≤ 2ℓ , almost all (in the sense of Lebesgue measure) ℓ × k matrices m form a full-column-rank H(m). (For k = 2ℓ this is because det H(m) is a polynomial in the entries of m, and the Walsh example shows that this polynomial is nonzero. For k < 2ℓ , consideration of the same 2ℓ × 2ℓ Walsh transform implies that there are some k rows of H(m) such that the determinant of the minor they form is a nonzero polynomial.) Despite this observation, the Vandermonde case, in which k − 1 rows are required, is very typical, as it is what arises in H(m) for a mixture model of observables Xi that are iid conditional on a hidden variable. Another class of examples that is far from Lebesgue-typical, and furthermore also far from being “separated” (see next section), is this. There are two possible coins, with biases p1 ̸= p2 . A hidden variable H is sampled in {0, . . . , k − 1}, and then the process is that you observe the result of H independent tosses of coin 1, followed by k − 1 − H independent tosses of coin 2. The NAE condition implies that here H(m) has full column rank. As a consequence (applying [34]) the following model is identifiable: a hidden H is sampled (from unknown prior) in {0, . . . , k − 1}, and then you observe the result of 2H independent tosses of coin 1 followed by 2k − 1 − 2H independent tosses of coin 2. A similar class of examples (sometimes identifiable but in general not) are the “subcube mixtures” studied in [20], where all coin biases must be one of {0, 1/2, 1}. 3.2 Motivation Consider observable random variables X1 , . . . , Xn that are statistically independent conditional on H, a hidden or latent random variable H supported on {1, . . . , k}. (See Figure 3.1 for an illustration.) The most fundamental case is that the Xi are binary. Then we denote mji = Pr(Xi = 1|H = j). The model parameters are m along with a probability distribution (the mixture distribution) π = (π1 , . . . , πk ) on H. The study of finite mixture models was pioneered in the late 1800s in [35], 34 H X1 X2 ... Xi ... Xn Figure 3.1: A Bayesian network diagram relating H and X1 , . . . , Xn . [36]. The problem of learning such distributions has drawn a great deal of attention. For surveys see, e.g., [37]–[40]. For some algorithmic papers on discrete-valued Xi , see [6], [9], [10], [15]–[20], [34]. The source identification (or parameter estimation) problem is that of computing (m, π) from the joint statistics of the Xi . Put another way, the problem is to invert the multilinear moment map [n] µ : (m, π) → R2 µ(m, π)S = Pr(XS = 1) where S ⊆ [n], XS = Y i∈S Xi = mS · π ⊤ . Since mjS = Pr(XS = 1|H = j), this shows the essential role of H(m) in the mixture model. Connection to rank H(m) In general µ is not injective (even allowing for permutation among the values of π and columns of m). For instance it is clearly not injective if m has two identical columns (unless π places no weight on those). More generally, and assuming all πj > 0, it cannot be injective unless H(m) has full column rank. (Suppose α ∈ Rk is nonzero s.t. H(m)α = 0. Since (H(m)α)|{∅} = 0, P j αj = 0. So for sufficiently small δ > 0, π + δα is a mixture distribution, distinct from π, with identical statistics.) One sufficient condition for injectivity, due to [41], is that there be 2k − 1 “separated” observables Xi . Xi is separated if all mji are distinct, or in our terminology, if no color recurs in mi . (Further it is shown in [34], Theorem 1, that one can lower bound the distance between µ(m, π) and any µ(m′ , π ′ ) ′ in terms of ζ = mini minj̸=j ′ |mji − mji | and the distance between (m, π) and (m′ , π ′ ).) There are examples with X1 , . . . , X2k−1 where the mapping is injective but is no longer so if any single Xi is omitted [9]. 35 A weaker and still sufficient condition for injectivity of µ, due to [34], is that for every i ∈ [n] there exist two disjoint sets A, B ⊆ [n] − {i} such that H(m|A ) and H(m|B ) have full column rank. (It is not known whether two disjoint such A, B are strictly necessary.) Observable Xi with larger finite range If each Xi can take on one of say L values, m can be considered as a nonnegP ative n × k × L real array, with mji,ℓ = Pr(Xi = ℓ|H = j), Lℓ=1 mji,ℓ = 1; the multivariate moments are indexed not by sets S but by mappings S : [n] → [L], with mS = mS(1) ⊙ · · · ⊙ mS(n) and µ : (m, π) → R[L] [n] µ(m, π)S = Pr(XS = 1) where XS = Yn i=1 δXi ,S(i) (Kronecker delta) = mS · π ⊤ For any given k, if L is sufficiently large and m satisfies a certain nonsingularity condition, the mixture learning problem becomes easier; this insight is due to [6]. It will be interesting to explore what conditions exactly m must satisfy for identifiability (for positive π), for arbitrary L. But in this paper we study only the most extreme, and hardest for identification, case L = 2. 3.3 Some Theory for Hadamard Products, and a Proof of Theorem 45 For v ∈ Rk and U a subspace, extend the definition of v⊙ to v⊙ (U ) = {u · v⊙ : u ∈ U } and introduce the notation v⊙¯ (U ) = span(U ∪ v⊙ (U )). We want to understand which subspaces U are invariant under v⊙¯ . Let v have distinct values λ1 > . . . > λℓ for ℓ ≤ k. Let the polynomials pv,i (i = 1, . . . , ℓ) of degree ℓ − 1 be the Lagrange interpolation polynomials for these values, so pv,i (λj ) = δij (Kronecker delta). Let B(v) denote the partition of [k] into blocks B(v)(i) = {j : vj = λi }. Let V(i) be the space spanned by the elementary basis vectors in B(v)(i) , and P(i) the projection onto V(i) w.r.t. the standard 36 inner product. Since v⊙ is diagonal with entries λi in B(v)(i) , we have the matrix equation pv,i (v⊙ ) = P(i) , (3.1) where pv,i is interpreted as a matrix polynomial. The collection of all linear combinations of the matrices P(i) is a commutative algebra, the B(v) projection P algebra, which we denote AB(v) . The identity of the algebra is I = P(i) . Definition 48. A subspace of Rk respects B(v) if it has a basis in which each vector lies in some V(i) . For a subspace U we let U ⊥ be its orthogonal complement w.r.t. the standard inner product. S L For U respecting B(v) write U = span( U(i) ) for U(i) ⊆ V(i) . (Thus U = U(i) L ⊥ ⊥ and U(i) = P(i) U .) Let D(i) = (U(i) ) ∩ V(i) . Then (U(i) ) = D(i) ⊕ j̸=i V(j) . Lemma 49. A subspace U ⊥ respects B(v) if U does. Proof. The subspaces of an inner product space form an orthocomplemented lattice in which the meet operation is intersection, and the negation operation is orthogonal complement. So for any subspaces W, W ′ we have De Morgan’s T L law (span(W ∪ W ′ ))⊥ = W ⊥ ∩ W ′⊥ . Thus U ⊥ = (U(i) )⊥ = D(i) . L Lemma 50. A subspace U respects B(v) iff U = (P(i) U ). Proof. (⇐): Because this gives an explicit representation of U as a direct sum of subspaces each restricted to some V(i) . ′ (⇒): By definition U is spanned by some collection of subspaces V(i) ⊆ V(i) ; L ′ since these subspaces are necessarily orthogonal, U = V(i) . Moreover, since P(i) annihilates V(j) , j ̸= i, and is the identity on V(i) , it follows that each ′ V(i) = P(i) U . Theorem 51. A subspace U is invariant under v⊙¯ iff U respects B(v). P Proof. (⇐): Let w ∈ U and write w = wi for wi ∈ U(i) . Then v ⊙ wi = P L λi wi ∈ U(i) . So v ⊙ w = v ⊙ wi ∈ U(i) = U . (⇒): If U = v⊙¯ (U ) then these also equal v⊙¯ (v⊙¯ (U )), etc., so U is an invariant space of AB(v) , meaning, aU ⊆ U for any a ∈ AB(v) . In particular, applyL ing (3.1), this holds for a = P(i) . So U ⊇ (P(i) U ). On the other hand, 37 P L L since P(i) = I, U = ( P(i) )U ⊆ (P(i) U ). So U = (P(i) U ). Now apply Lemma 50. P The symbol ⊂ is reserved for strict inclusion. Lemma 52. If R, T ⊆ [n] and rowspace H(m|R ) ⊂ rowspace H(m|R∪T ), then there is a row t ∈ T such that rowspace H(m|R ) ⊂ rowspace H(m|R∪{t} ). Proof. Without loss of generality R, T are disjoint. Let T ′ ⊆ T be a smallest set s.t. ∃R′ ⊆ R s.t. mR′ ⊙ mT ′ ∈ / rowspace H(m|R ). Select any t ∈ T ′ and write mR′ ⊙ mT ′ = mR′ ⊙ mT ′ −{t} ⊙ mt . By minimality of T ′ , mR′ ⊙ mT ′ −{t} ∈ rowspace H(m|R ). But then mR′ ⊙ mT ′ ∈ rowspace H(m|R∪{t} ), so rowspace H(m|R ) ⊂ rowspace H(m|R∪{t} ). Proof of Theorem 45. This is now a consequence of Lemma 52. Start with R = ∅, and repeatedly use the Lemma to adjoin to R a row from [n] \ R which will increase the rank of H(m|R ) by at least 1. Remark rank H(m), along with a basis (using only rows of H(m)) for rowspace H(m), can be computed in time O(nk 3 ) using Chen and Moitra’s “GrowByOne” procedure [20]. For completeness here is a version of that procedure: For ℓ ≥ 0 let Wℓ = span(m|[ℓ] ), and let rℓ = rank Wℓ . Wℓ is spanned by some vectors mSℓ,1 , . . . , mSℓ,rℓ , with all Sℓ,i ⊆ [ℓ], which we compute as follows. For ℓ = 0 we have r0 = 1, S0,1 = ∅. For ℓ > 1 form the matrix with rows mSℓ−1,1 , . . . , mSℓ−1,rℓ−1 followed by rows mSℓ−1,1 ∪{ℓ} , . . . , mSℓ−1,rℓ−1 ∪{ℓ} . Perform Gaussian elimination to zero-out all but rℓ − rℓ−1 of the second batch of rows. The first batch, together with the non-eliminated rows of the second batch, become mSℓ,1 , . . . , mSℓ,rℓ . 3.4 Combinatorics of the NAE Condition: Proof of Theorem 47(a) Recall we are to show: If ε(m) ≥ −1 then m has a restriction to some k − 1 rows on which ε = −1. Proof of Theorem 47(a). We induct on k. The (vacuous) base-case is k = 1. For k > 1, we proceed by way of contradiction. Suppose the theorem fails for k, and let m be a k-column counterexample with the least possible number of rows, n. So n > k − 1 ≥ 1. Necessarily every row of m is in NAE(m). Our 38 strategy is to show m has a restriction m′ to n − 1 rows, for which ε(m′ ) ≥ −1; this will imply a contradiction because, by minimality of the number of rows of m, m′ has a restriction to k − 1 rows on which ε = −1. If ε(m) ≥ 0 then we can remove any single row of m and still satisfy ε ≥ −1. Otherwise, ε(m) = −1, so there is a nonempty S such that | NAE(m|S )| = |S| − 1; choose a largest such S. It cannot be that S = [k] (as then n = k − 1). Arrange the rows NAE(m|S ) as the bottom |S| − 1 rows of the matrix. As discussed earlier, for the NAE condition one may regard the distinct real values in each row of m simply as distinct colors; relabel the colors in each row above NAE(m|S ) so the color above S is called “white.” (There need be no consistency among the real numbers called white in different rows.) See Fig. 3.2. f R' f(T1) R1 m R'' R1 NAE(m|S) m' R2 T2 S T1 [k] - S Figure 3.2: Argument for Theorem 47(a). Upper-left region is white. Entries (t, f (t)) (indicated with black dots) are not white. Due to the maximality of |S| and the fact that ε(m) ≥ −1, there is no set of columns S ′ with S ⊂ S ′ such that for some set of rows A ⊆ [n] − NAE(m|S ), ′ with |A| = n − |S ′ | + 1, m|SA is all white. That is to say, if we form a bipartite graph on “right” vertices corresponding to the columns [k] − S, and “left” 39 vertices corresponding to the rows [n]−NAE(m|S ), with non-white cells being edges, then any subset of the right vertices of size ℓ ≥ 1 has at least ℓ + 1 neighbors within the left vertices. By the induction on k (since S ̸= ∅), for the set of columns [k] − S there is [k]−S a set R′′ of k − |S| − 1 rows such that ε(m|R′′ ) = −1. Together with the rows of NAE(m|S ) this amounts to at most k − 2 rows, so since n ≥ k, we can find two rows outside this union; delete either one of them, leaving a matrix m′ with n − 1 rows. This matrix has the rows NAE(m|S ) at the bottom, and n − |S| remaining rows which we call R′ . The lemma will follow by showing that ε(m′ ) ≥ −1. In m′ , the induced bipartite graph on right vertices [k] − S and left vertices R′ has the property that any right subset of size ℓ ≥ 1 has a neighborhood of size at least ℓ in R′ . Applying Hall’s Marriage Theorem, there is an injective f : [k] − S → R′ employing only edges of the graph. Now consider any nonempty set of columns T , and write it as T = T1 ∪ T2 for T1 ⊆ [k] − S and T2 ⊆ S. We need to show that ε(m′ |T ) ≥ −1. Let R1 = NAE(m′ |T1 ) ∩ R′′ and R2 = NAE(m′ |T2 ). We have that |R1 | ≥ |T1 | − 1 [k]−S because ε(m|R′′ ) = −1. We further have that |R2 | ≥ |T2 | − 1 because ε(m) ≥ −1 and because NAE(m|T2 ) ⊆ NAE(m|S ) = NAE(m′ |S ), so no row of NAE(m|T2 ) has been removed in m′ . If T2 = ∅, the rows R1 witness that ε(m′ |T ) ≥ −1. Likewise if T1 = ∅, the rows R2 witness the same conclusion. Lastly suppose both T1 and T2 are nonempty. Nonemptyness of T2 gives | NAE(m|T2 )| ≥ |T2 | − 1. Now use the matching f . The set of rows f (T1 ) lies in R′ and is therefore disjoint from NAE(m|T2 ), which as noted is a subset of NAE(m′ |S ). Moreover since T2 ̸= ∅, every entry (t, j) for t ∈ T2 , j ∈ R′ is white. On the other hand due to the construction of f , for every t ∈ T1 the entry (t, f (t)) is non-white. Therefore every row in f (T1 ) is in NAE(m′ |T1 ∪T2 ). So | NAE(m′ |T1 ∪T2 )| ≥ |T2 | − 1 + |T1 |, which is to say ε(m′ |T ) ≥ −1. Thus ε(m′ ) ≥ −1. 3.5 From NAE to Rank: Proof of Theorem 47(b) Recall we are to show: H(m) has full column rank if ε(m) ≥ −1. 40 Proof of Theorem 47(b). The case k = 1 is trivial. Now suppose k ≥ 2 and that Theorem 47(b) holds for all k ′ < k. Any constant rows of m affect neither the hypothesis nor the conclusion, so remove them, leaving m with at least k − 1 rows. Now pick any set, C, of k − 1 columns of m. By Theorem 47(a) there are some k − 2 rows of m, call them R′ , on which ε(m|C R′ ) = −1. Let v be a row of m outside R′ . Let R′′ denote the set of rows of m other than v. Since R′′ contains R′ , by induction dim rowspace H(m|C R′′ ) = k − 1. Therefore U := rowspace H(m|R′′ ) ⊆ Rk is of dimension at least k − 1. We claim now that dim v⊙¯ (U ) = k. (Note that v⊙¯ (U ) = rowspace H(m).) Suppose to the contrary that dim v⊙¯ (U ) = k − 1. It must then be that dim U = k − 1 and v⊙¯ (U ) = U . So as proven in Theorem 51, U respects B(v). Since v is nonconstant, B(v) is a partition of [k] into ℓ ≥ 2 nonempty blocks L B(v)(i) , and U = ℓi=1 U(i) with U(i) = P(i) U(i) . So there is some i0 for which U(i0 ) ⊂ V(i0 ) ; specifically, U(i) = V(i) for all i ̸= i0 , and dim U(i0 ) = dim V(i0 ) − 1. Since |B(v)(i0 ) | < k, we know by induction that P(i0 ) rowspace H(m) = V(i0 ) . But since rowspace H(m) = v⊙¯ (U ) = U , this means that P(i0 ) U = V(i0 ) . Contradiction. 41 Chapter 4 SOURCE IDENTIFICATION FOR MIXTURES OF PRODUCTS [1] S. Gordon, B. H. Mazaheri, Y. Rabani, and L. Schulman, “Source identification for mixtures of product distributions,” in Proceedings of Thirty Fourth Conference on Learning Theory, M. Belkin and S. Kpotufe, Eds., ser. Proceedings of Machine Learning Research, vol. 134, PMLR, Aug. 2021, pp. 2193–2216. [Online]. Available: https://proceedings.mlr. press/v134/gordon21a.html, 4.1 Introduction A k-MixProd is a mixture (that is, a convex combination) of k product distributions on a set of n random variables X . In the notation of Bayesian networks, this situation is represented graphically by a single unobservable random variable U with edges to each of the variables Xi ∈ X . U is referred to as a “confounding” variable with range 1, . . . , k and the variables in X are referred to as the “observables.” The main complexity parameter of the problem is k, the number of mixture constituents or “sources.” See Fig. 4.1. U X1 X2 X3 ... Xn Figure 4.1: Graphical depiction of a k-MixProd In this chapter we study the identification problem for k-MixProds. Specifically, given a joint distribution P on the variables (vertices), recover up to small statistical error (a) the mixture weights, up to a permutation of the constituents, and (b) the conditional distribution of each vertex X within each mixture source. Problem Statement and Notation. The hidden variable U ranges in [k] = {1, . . . , k}, with unknown distribution πj := Pr(U = j). There are n binary observable variables X1 , X2 , . . . , Xn ; their distributions are given by the rows mi of an unknown matrix n × k matrix 42 m, with the meaning that mij := Pr[Xi = 1 | U = j]. Thus, the identification problem is to recover (πj )j∈[k] and (mij )i∈[n],j∈[k] (up to permuting the set [k]) to within high accuracy. Identification is not always possible, as there can be disparate models sharing the same distribution on the observables. To give sufficient conditions for identifiability, we need the following definition, which allows us to quantify the degree to which an observable variable has distinct behavior in each of the different mixture constituents. Definition 53. An observable Xi is ζ-separated if minj̸=j ′ |mij − mij ′ | > ζ. It is known that a k-MixProd is identifiable (from perfect statistics) if it has at least 2k − 1 0-separated observables [42] and πmin > 0. In general (although not for every m), 2k − 1 observables are also necessary [9]. Moreover some separation assumption is certainly necessary for identification: e.g., a mixture of several identical sources generates the same statistics as a single source. Clearly, a sample size bound requires a further, quantitative assumption on separation; this is the role of the parameter ζ in our work. For our nearoptimal running time and sample complexity guarantee, we also require a slightly larger-than-minimal number of observables: 3k − 3 rather than 2k − 1. (A similar algorithm works with just 2k − 1 ζ-separated observables, and its post-sampling runtime is about the same, but its sample size requirement is larger.) It should be noted that the post-sampling runtime, although exponential in k, is hardly affected by ζ. The bottleneck resource is sample size. Main Result Theorem 54. The algorithm given later w.h.p. identifies a k-MixProd on n ≥ 3k − 3 binary observables, each of which is ζ-separated, with runtime (and, a fortiori, sample complexity) (1/ζ)O(k log k) (n log n)(1/ min πi )O(log k) (1/ε)2 . and w.h.p. (over the distribution of samples) computes π̃ and m̃ such that maxi |π̃i − πi | ≤ ε and maxij |m̃ji − mji | ≤ ε. 43 It is actually sufficient that there merely exist a subset of 3k − 3 ζ-separated random variables. Under this weaker assumption, the runtime above is multiplied by nO(k) (the sample size is unaffected). Discussion The k-MixProd problem for finite-range observables has been studied for nearly 30 years [15]–[19], [42]–[44] There are two versions of the problem: (1) Learning the model, namely, producing any model consistent with (or close to) the observed statistics; (2) Identifying the model, namely, producing the true model (or one close to it). Any algorithm for “learning” will also of course achieve identification if it happens that the model is uniquely specified by the statistics, but the algorithm may not, and existing algorithms do not, provide a certificate of uniqueness under such conditions. Our result is concerned entirely with identification. Since our algorithm computes, as part of its operation, the condition number of the matrices that it uses for inverting the model→statistics mapping, it will only return an output under conditions which ensure the output is indeed unique (to within the allowed ε error). Hence, as compared with a “learning” algorithm, we are effectively providing a stronger output guarantee under stronger conditions on the input. One of the chief motivations for our work is the characterization of interventional distributions in Bayesian networks (causal DAGs). Causal inference is a sizable and growing field; in that setting, we want to understand the distribution on some “downstream” variable(s) that will result if we disconnect a specified “intervention variable” from its parents in the DAG, and instead assign it to a chosen value. This task is easy when all variables in the network are observable, but challenging when there are latent or “hidden” variables whose statistics we cannot observe. In the most extreme case there may be a latent variable U which can directly affect all observables. It turns out, as we show in separate work, that it is possible to solve that problem when U takes on a bounded number k of values. In that work, a key subroutine has to solve k-MixProd problems. The role of the present work is to provide and analyze the needed algorithm for that problem. We should emphasize that “learning” the model, without “identification,” is (so far as we know) useless for causal inference applications. Fortunately, 44 we show that when identification is possible, it can also be performed much faster than the best current results for “learning” (which suffer a k 3 in the exponent, as compared with our k log k). Of course it remains of great interest to ask whether a similar runtime can be achieved for learning. Our algorithm, in any event, does not achieve this. Comments on the assumptions in Theorem 54. 1. The observable variables are binary. The binary case is actually the most difficult. (The problem also becomes significantly easier when each variable Xi takes on at least k values and the k distributions on Xi span a k-dimensional space. This kind of assumption was first to our knowledge used in, in a somewhat different setting, in [45].) 2. The mixture is ζ-separated. A separation condition of some kind, on at least some of the variables, is necessary, since we certainly cannot identify the distribution of the latent variable if it does not have sufficient effect upon the observable variables. The 0-separation condition fails only on a set of Lebesgue measure 0 in the parameter space. Likewise for small ζ, most distributions are ζ-separated. The value of ζ-separation will be, as we show below, that it ensure that matrices representing the parameters for each source are well-conditioned. ζ-separation is not always a necessary condition; characterizing necessary conditions is a difficult open question, tackled in part in [46], where it is shown, essentially, that some “batching” of not-fully-separated variables is also sufficient to imply invertibility of the relevant matrices (however, no quantitative bounds on condition number are known). 3. The required value of n. Depending on the actual distributions (as expressed by m), the number of ζ-separated variables n required to identify the source may be anywhere from lg n to 2k − 1. Less than lg n is not possible: in that case a k-MixProd can be any probability distribution on X , and it is not hard to see that parameters are not unique. Prior work. In [19] a seminal learning algorithm for k-MixProd was given. 3 Its running time, for mixtures on n binary variables (n suff. large), is nO(k ) . 3 2 This was improved in [43] to k O(k ) nO(k ) . The most recent algorithm [44] identifies a mixture of k product distributions on at least 3k − 3 ζ-separated 2 variables in time 2O(k ) n. In this work we improve the latter results, giving 45 upper bounds with exponent O(k log k) rather than O(k 2 ). This is nearly optimal, as it is known that even in the special case that the observables are known to be iid conditional on U —call this the k-MixIID problem—exponent Ω(k) is unavoidable [9]. Moreover n ≥ 2k − 1 observable variables are generally necessary. (There are interesting special cases in which better efficiency is possible. [43] study “mixtures of sub-cubes,” i.e., k-MixProd with all bit probabilities in {0, 12 , 1}. For this problem they achieve complexity nO(log k) .) It turns out that the k-MixIID problem plays a special role at the base of a tower of reductions. On the one hand it can be solved using a relatively simple, two-century-old method of Prony [7], which connects it to the classical Hausdorff moment problem; it can also be solved by the Matrix Pencil Method [3], [9], [10], [47]. Both analyses are quite recent. In turn, the k-MixProd problem is solved by a nontrivial reduction to k-MixIID, first with exp(O(k 2 )) complexity [44] and here with exp(O(k log k)) complexity. In forthcoming work the authors have shown in turn how to reduce to k-MixProd from the yet more general k-MixBND problem, in which the input distribution is a k-component mixture over distributions consistent with a (known) Bayesian Network. The dominant term in the complexity of that algorithm is that of the k-MixProd instances to which the problem reduces. Hence the key role played by the improved sample- and run-time complexities of the present work. Comparison with the parametric case. The literature on mixture models for parametric families (exponential distributions, gaussians in R or Rd , etc.) is even more extensive and older than for the discrete case. It is necessary however to point out a fundamental difference between the types of problems. In general, in a k-MixIID problem, a mixture source is chosen from distribution π), and then we sample some n independent samples from that source. In almost every parametric scenario (think e.g., of a mixture of k unit-variance gaussians on the line), n = 1 is sufficient in order to (in the limit of many repetitions) identify the model. This is fundamentally untrue in the non-parametric case. To see this, consider an instance of 2-MixIID for binary variables with two equiprobable sources (i.e., π1 = π2 = 1/2), one of which has Pr(X = 1 | U = 0) = 43 and the other Pr(X = 1 | U = 1) = 14 . With only a single sample from each source, it is impossible to distinguish between a 46 mixture in which Pr(X = 1 | U = 0) = 1 and the other Pr(X = 1 | U = 1) = 0. More specifically, it is only possible to identify the mean of these two parameters. With access to multiple samples from the same source, however, we can begin to infer the parameters associated with each mixture component. For the k-MixIID problem for binary rv’s, it is known that n = 2k − 1 such samples [9] are needed. (This value is called there the threshold “aperture” of the problem.) Unlike k-MixIID, the k-MixProd problem does not allow us access to multiple samples of iid variables from a source. However it does give access to multiple variables which are independent conditional on the source. An approach introduced in [44] is to create synthetic copies of a single variable via linear combinations of the other variables; since we need to modify the method, we describe below how this is done. This approach enables a reduction from the k-MixProd problem to the k-MixIID problem. Organization. Section 4.2 sets up the key mathematical objects needed for our work. Section 4.3 describes the algorithm. The algorithm is similar to that in [44], differing in one crucial way which will be described. Along with the algorithm pointers are provided to the main steps of the analysis in Appendix 4.5. The part of this work which is a full departure from the previous literature is the condition number bound in Section 4.4, in which the condition number of a key linear mapping (namely, the matrix of multilinear moments of the distribution on observables), is bounded through a novel argument characterizing the images of rank-1 tensors under mapping by Hadamard extensions. 4.2 Preliminaries Hadamard Extensions As in earlier works [43], [44], the algorithm and its analysis make extensive use of the Hadamard extension of a matrix m. Definition 55. Given a matrix m of any dimensions, let mi∗ denote the ith row of m and m∗j the jth column of m. Where clear from context we write mi instead of mi∗ . Definition 56 (Hadamard product). The Hadamard product is the mapping ⊙ : Rk × Rk → Rk which, for row vectors u = (u1 , . . . , uk ) and v = (v1 , . . . , vk ), 47 is given by u ⊙ v := (u1 v1 , . . . , uk vk ). Equivalently, using the linear operator v⊙ = diag(v), the Hadamard product is u ⊙ v = u · v⊙ . The identity element for the Hadamard product is the all-ones vector 1. The Hadamard product of u with itself k times is written u⊙k . Definition 57 (Hadamard extension). For n ∈ R[n]×[p] , the Hadamard extension of n, written H(n), is the 2n × p matrix with rows H(n)S for all S ⊆ [n], where, Q for S = {i1 , . . . , iℓ }, H(n)S = ni1 ⊙ · · · ⊙ niℓ ; equivalently H(n)jS = i∈S nji . In particular H(n)∅ = 1, and for all i ∈ [n], H(n){i} = ni . Notation for subsets and collections of subsets We will reserve calligraphic fonts for collections of subsets of [n], i.e., S ⊆ 2[n] . In contrast, the same variable in a standard math font will represent a subset of [n], i.e., S ⊆ [n]. Definition 58 (Sum of two collections). The sum of two collections of subsets X , Y ⊆ 2[n] is defined as X + Y := {X ∪ Y : X ∈ X , Y ∈ Y}. Definition 59 (Kronecker product of vectors). Let S, T ⊂ [n] be two disjoint sets. Let x and y, respectively, be vectors indexed by the subsets in X = 2S and Y = 2T , respectively. Then the Kronecker product x ⊗ y is the vector indexed by the subsets in X + Y given by (x ⊗ y)X∪Y := xX yY , for all X ∈ X and Y ∈ Y. Note that X + Y = 2S∪T . (Notice that every set Z ∈ X + Y can be written uniquely as Z = X ∪ Y for X ∈ X and Y ∈ Y. We don’t define the Kronecker product when X ∩ Y = ̸ ∅.) Definition 60 (Kronecker product of vectors — alternate definition). Let S T S, T ⊂ [n] be two disjoint sets. Let x ∈ R2 and y ∈ R2 , which is to say, x has a coordinate xR for every R ⊆ S, and similarly for y. Then the Kronecker product x ⊗ y is the vector with coordinates xR for every R ⊆ S ∪ T , given by (x ⊗ y)R := xR∩S yR∩T . In the algorithm and analysis we will make heavy use of the singular values of a matrix M, which we will denote σmax (M) = σ1 (M) ≥ σ2 (M) ≥ · · · .  k When v (i) i=1 is a set of row vectors, we will write (v (1) ; v (2) ; . . . ; v (k) ) for the matrix obtained by vertically concatenating each row vector in order. 48 Finally, for a matrix M, M+ will denote the Moore-Penrose inverse (i.e., the pseudoinverse), given by (MT M)−1 MT in the case that M has full column rank. Definition 61 (Column-wise Khatri-Rao product of matrices). The columnwise Khatri-Rao product of matrices A ∈ Rn×k , B ∈ Rm×k is denoted A ∗ B; A ∗ B has dimensions nm × k with row ij given by Ai∗ ⊙ Bj∗ . Definition 62 (Vandermonde matrix). The Vandermonde matrix Vdm(m) ∈ Rk×k associated with a vector m ∈ Rk has entries Vdm(m)ji = mji for j ∈ {0, 1, . . . , k − 1} and i ∈ {1, 2, . . . , k}. Multilinear Moments Using the notation from Sec. 4.1, observe that the expectation of each Xi is given by mi π T . The model parameters are (πj )j∈[k] and (mij )i∈[n],j∈[k] (up to a permutation of the range of U ). When k > 1, observable variables Xi and Xi′ will not in general be independent. We will make use of information about the dependencies between observables through measurements of E[Xi Xi′ ]. Generalizing this, for any Q S ⊆ [n], we will define the random variable XS := i∈S Xi . The data we obtain from our samples will be estimates of E[XS ] for different subsets S ⊆ [n]. These are known as the multilinear moments of the distribution, since they are multilinear polynomials in the parameters mij . Our goal will be to eventually compute power moments for the single variable (2) (3) (ℓ) (j) X1 , i.e. the values E[X1 X1 X1 · · · X1 ] for ℓ ≤ 2k where each X1 is a copy of X1 that is iid conditioned on U . Such copies are not guaranteed to exist among the Xi , but we can still define (and our algorithm will compute) the power moments. The ℓth power moment is defined as E[X1⊙ℓ ], the expectation of the product of ℓ copies of X1 that are iid conditioned on U . P Q Now E[XS ] = j πj i∈S mij , or equivalently, E[XS ] = (mi1 ⊙ mi2 ⊙ · · · ⊙ mis )π T where S = {i1 , i2 , . . . , is }. If we replace m with the Hadamard extension M := H(m), we can simplify the above equation to E[XS ] = MS π T . We can also write power moments in terms of Hadamard products. In particP T ular, we have E[Xi⊙ℓ ] = j πj mℓij = m⊙ℓ i π . 49 Observe that source identification is not possible if M has less than full column rank, i.e., rank M < k, as then the mixing weights cannot be unique. For a collection of subsets S ⊆ 2[n] , let M[S] denote the restriction of M to the rows MS , S ∈ S. E.g., M = M[2[n] ]. If we have non-overlapping collections of subsets A and B, we can build a matrix of multilinear moments as follows: Definition 63 (Observable matrices). For collections A, B ⊆ 2[n] with A ⊆ 2S and B ⊆ 2T for disjoint S, T ⊆ [n], we define the matrix CB,A := M[B]π⊙ M[A]T . For any A ∈ A, B ∈ B, we have (CB,A )B,A = MB π⊙ MTA which is the multilinear moment E[XA∪B ] and therefore observable. The Empirical Multi-linear Moments For a finite sample drawn from the model, we let g̃(S) be the empirical estimate of E[XS ], i.e., the fraction of Q samples for which i∈S Xi = 1. These g̃(S) for S ⊆ [n] are the complete list of “observables” of the model. Each converges, in the infinite-sample limit, to the value g(S) := E[XS ], g(S) = MS π T = MS π⊙ 1T . Our algorithm will heavily utilize (the empirical estimates of) matrices of observable moments. 4.3 The Algorithm The algorithm appears below as Algorithm 3. It uses disjoint sets S, T, T ′ ⊆ [n] of ζ-separated variables, with |S| = |T | = |T ′ | = k − 1 and 1 ∈ T , and takes as input the desired accuracy ε. We will use tildes in the algorithm to denote the empirical version of quantities we define with respect to exact moments. To avoid clutter however we will describe most steps of the algorithm as though we had access to the exact moments we need (i.e. g(S) instead of g̃(S)). In the actual algorithm we use our empirical estimates of the moments. 50 Definition 64. Define the following collections from S, T, T ′ : A := 2S , ′ B ′ := 2T . B := 2T , We introduce shorthand for the corresponding submatrices of M: A := M[A], B := M[B], B′ := M[B]. An essential difference between this algorithm and the one given in [44] is ′ that our A, B, B ′ comprise of the entirety of 2S , 2T , and 2T . Consequently, we use the entire 2k × 2k matrix for CAB rather than performing an expensive 2 2O(k ) search within it for a well-conditioned k × k submatrix. In the earlier paper this search cost was dominated by other costs, but that is no longer the case here, due to the stronger bounds we provide (coming from Section 4.4). Our goal will be to compute the power moments E[X1⊙ℓ ]. To do this, we will alternate between computing vectors of moments vi , and vectors of coefficients ui , u′i , defined below: Definition 65. The sequence of vectors v1 , . . . , v2k ∈ RA is given by (vi )A := E[XA X1⊙i ] for each A ∈ A. In particular, (vi )∅ = E[X1⊙i ]. Proposition 66. vi = m⊙i 1 π⊙ A. That is, each entry of the vector vi will be a mixed moment, with a multilinear part given by XA and a power moment part given by X1⊙ℓ . Our algorithms begins by computing ṽ0 and ṽ1 , the empirical counterparts to v0 and v1 , respectively. These empirical counterparts are available directly as empirical multilinear moments. Definition 67. The sequence of vectors u1 , . . . , u2k ∈ RB is given by the equations ui M[B] := m⊙i 1 . ′ The sequence of vectors u′1 , . . . , u′2k ∈ RB is defined analogously with B ′ replacing B everywhere. 51 + ′ Proposition 68. ui = vi C+ BA and ui = vi CB′ A . That is, ui is the vector of coefficients of linear combinations of parameter vectors from M[B] (or M[B ′ ]) needed to create the parameter vector for the variable X1⊙i . Once we have computed v1 , v2 , . . . , v2k , we can take the entries (vi )∅ = E[X1⊙i ] to obtain the first 2k power moments for X1 , which we can then pass as input to an algorithm for learning k-MixIID. At this point, we will have recovered m1 and π and we can use these to recover the remaining parameters with some more straightforward linear algebra. Empirical counterparts In the pseudocode for Algorithm 3, we will exclusively be working with empirically computed (hence, approximate) versions of the quantities above. The empirical counterpart to any quantity is distinguished by a tilde, e.g., Ã is the empirical counterpart to A. The Steps of the Algorithm We now outline what each non-trivial line of the algorithm accomplishes. • Line 1: As stated in Theorem 54, we write the algorithm under the assumption that all observable bits are ζ-separated (and consequently in this line, S, T, T ′ are arbitrary disjoint sets of size k−1). This assumption is not strictly necessary; all we actually need is that the disjoint sets S, T, T ′ create ĈBA and ĈB′ A with σk above the threshold in Line 6, something that is guaranteed if S, T, T ′ are of size k − 1 and consist of ζ-separated variables, but may also hold without either of these conditions. In any case, provided one assumes that there exist some 3k − 3 ζ-separated variables, an exhaustive search running in time nO(k) will suffice to find S, T, T ′ passing the threshold in Line 6. • Line 3: For the algorithm outlined here to produce m and π within additive error ε (in each coordinate), we shall need the g̃(S) that we use to be accurate to within ε(min πi )O(log k) ζ O(k log k) . This can be achieved with high probability using a sample of size (1/ε)2 (1/ min πi )O(log k) (1/ζ)O(k log k) . 52 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 input: disjoint sets of ζ−separated S, T, T ′ of size k − 1 each, with 1 ∈ T , as well as the desired accuracy ε. ′ A ← 2S , B ← 2T , B ′ ← 2T . Use (1/ε)2 (1/ min πi )O(log k) (1/ζ)O(k log k) samples to estimate {g̃(A ∪ B ∪ B ′ )}A∈A,B∈B,B ′ ∈B′ . Construct C̃BA , C̃B′ A , and C̃B+B′ ,A from the empirical moments. Compute ĈBA and ĈB′ A by taking the SVD of C̃BA and C̃B′ A and truncating to the first k singular values. if min{σk (ĈBA ), σk (ĈB′ A )} < πmin ζ O(k) then terminate ṽ0 ← (g̃(A))A∈A ṽ1 ← (g̃(A ∪ {1}))A∈A ũ1 ← ṽ1 Ĉ+ BA + ′ ũ1 ← ṽ1 ĈB′ A for i = 1, . . . , 1 + lg k for j = 1, . . . , 2i−1 ṽ2i−1 +j ← (ũj ⊗ ũ′2i−1 )ĈB+B′ ,A ũ2i−1 +j ← ṽ2i−1 +j Ĉ+ BA + ′ i ũ2i ← ṽ2 ĈB′ A . Let Hk+1 be the (k + 1) × (k + 1) Hankel matrix with entries given by [Hk+1 ]i,j = (ṽi+j ){1} if the second−smallest eigenvalue of Hk+1 is below πmin ζ O(k) then terminate m̃1 , π̃ ← LearnPowerDistribution(Hk+1 ). Ṽ ← (ṽ0 ; . . . ; ṽk−1 ) ⊙(k−1) Vdm(m̃1 ) ← (m̃⊙0 ) 1 ; . . . ; m̃1 −1 T −1 Ã ← π̃⊙ (Vdm(m̃1 )) Ṽ −1 B̃ ← ĈBA (ÃT )+ π⊙ .  + T −1 . for every i ∈ [n] \ S, m̃i ← ((g̃(R ∪ {i}))R∈A ) ÃT π̃⊙  + −1 for every i ∈ S, m̃i ← ((g̃(R ∪ {i}))R∈B )T B̃T π̃⊙ . Algorithm 3: Identifies a mixture of product distributions given 3k − 3 ζ-separated observable bits. 53 • Line 5: This line is a critical place where it is necessary to acknowledge the empirical nature of the quantities we are computing with, and account for this explicitly. The matrix C̃BA will w.h.p. be close in Frobenius norm to CBA when the moments are computed from a large sample, but C̃BA will almost certainly have rank(C̃BA ) = 2k−1 > rank(CBA ) = k, which will introduce instability in the pseudoinverse. To avoid this we replace C̃BA (resp. C̃B′ A ) with ĈBA (resp. ĈB′ A ), the rank k approximation obtained by truncating to the first k singular values. The SVD of an m × m matrix can be computed in time O(m3 ) [48] which in our setting is O(23k ). Note that computing the pseudoinverse of ĈBA takes no more time than computing ĈBA since we use the same SVD of C̃BA to compute both. (There are also methods which are faster for a low-rank approximation.) • Line 6: If σk (CBA ) is too small, we have a contradiction to our assumption that S and T consisted of ζ-separated observables. Likewise, σk (CBA ) being too small contradicts the ζ-separation of the variables in S and T ′ . • Lines 7-8: Here we create v0 and v1 , both of which are directly constructed from observable quantities. • Lines 9-15: In lines 9-10, we do the first step of the iteration that we continue in lines 11-15. We initialize u1 and u′1 by multiplying v1 by the pseudoinverses of CBA and CB′ A , respectively. Note here a difference from [44], where the CBA matrices were k × k and full-rank. We then compute vi , ui , and u′i for larger values of i. We use the fact that vi+j = (ui ⊗uj )CB+B′ ,A , proven in Lemma 69 to do this computation. We analyze the error introduced by this step in Lemma 80. This crucially relies on the condition number bound in Theorem 71, which is proved in Section 4.4. Each time we multiply by Ĉ+ BA we introduce some additional error and we are only able to bound the error by showing that + the operator norm of Ĉ+ BA − CBA is small, which requires this condition number bound. There are lg k steps in the iteration, each of which requires 2O(k) time, since we multiply by matrices of size 2O(k) constantly many times each iteration. (The Kronecker product computation also requires 2O(k) time.) 54 • Lines 16-18: We assemble the power moments of X1 into a Hankel matrix and invoke an algorithm for identifying a k-MixIID. (If the Hankel matrix has insufficiently large second smallest singular value, we conclude that one of our assumptions has not been satisfied—in particular, our estimated moments could be outside the required accuracy; or that there was a rare statistical error in the sampling.) The recovery of m1 and π can be done in time k 2+o(1) + O(k(log2 k) log log(ε−1 )) (see Corollary 83 for the runtime and the error bound). • Lines 19-21: We compute A by multiplying the matrix of the first k vi vectors by the inverse of the matrix Vdm(m̃1 ), the Vandermonde matrix generated from the empirical version of row m1 . (See Lemma 86 for the resulting error bound.) • Line 22: We compute B by multiplying CBA by the pseudoinverse of A. Again, the pseudo-inverse here replaces the regular inverse used in [44] because A is no longer square. (See Lemma 88 for the resulting error bound.) • Lines 23-24: We compute the remaining parameters in a similar fashion, using A+ or B+ as needed, depending on if the observable in question could overlap with an observable used by A or B. Lemma 69. For all ℓ, (vℓ )∅ = E[X1⊙ℓ ]. T Proof. We show by induction on ℓ that vℓ = m⊙ℓ 1 π⊙ A , hence (vℓ ){1} = T ℓ m⊙ℓ 1 π⊙ 1 = E[X1 ]. The base case of ℓ = 1 follows trivially from the definition of v1 . So suppose that this is true of vℓ′ for ℓ′ < ℓ. Let ℓ = 2i + j for j ∈ {1, 2, . . . , 2i }. Now, A, B, and B′ have full column rank, so + vℓ = (vj C+ BA ⊗ v2i CB′ A )CB+B′ ,A i T T + ⊙2 T ′ T + ′ T = (m⊙j 1 π⊙ A (Bπ⊙ A ) ⊗ m1 π⊙ A (B π⊙ A ) )(B ∗ B )π⊙ A (using Defn. 61) ⊙2i + ′ + ′ T = (m⊙j 1 (B) ⊗ m1 (B ) )(B ∗ B )π⊙ A = mℓ1 π⊙ AT , When actually running the algorithm, we will only have access to the approximations C̃BA , C̃B′ A , and C̃B+B′ ,A , and, starting with ṽ1 , we will compute 55 + ′ approximate vectors ũj = ṽj Ĉ+ BA , and ũ2i = ṽ2i ĈB′ A , and use those to compute an approximate vector ṽℓ = (ũj ⊗ ũ′2i )C̃B+B′ ,A . Notice the advantage of the recurrence in Line 12; we are able to get away with performing at most 1 + lg k iterations to compute any of ṽ1 , . . . , ṽ2k . Each iteration consists of at + most two matrix multiplications by Ĉ+ BA and ĈB′ A , followed by a convolution, followed by multiplication by CB+B′ ,A . Each such step can increase the −O(1) initial sampling error by at most a factor of ζ −O(k) πmin . By starting with Ω(log k) empirical moments accurate to within εζ Ω(k log k) πmin , we can ensure that the resulting vectors ṽi are sufficiently close to the vectors vi to start solving for m1 and π. The conclusion of our analysis can be stated as follows: Theorem 70. 1. If all multilinear moments used in the computation are within ±ε of their true values, the output π̃ and m̃ will satisfy −O(lg k) −O(k lg k) ∥m̃ − m∥∞ , ∥π̃ − π∥∞ ≤ πmin ζ ε. 2. With high probability, all empirical multilinear moments appearing in Algorithm 3 will be within ±ε of their true value. Proof. 1. This is an immediate consequence of Theorem 89 and Corollary 83. 2. With N samples, Hoeffding’s inequality upper bounds the probability that each moment g(S) differs from its empirical estimate g̃(S) by at least ε as h i ˜ − g(S)| > ε ≤ 2e−2ε2 N . Pr |g(S) A union bound over 3(2k ) such moments gives us a bound for the failure event, h i ˜ − g(S)| > ε ≤ 6(2k )e−2ε2 N . Pr ∃S : |g(S) This becomes constant when the number of samples N is O(k/ε2 ). We note that a simpler version of the above procedure uses only 2k − 1 ζ-separated variables, but needs 2k − 1 iterations to compute the ṽi -s, hence the required initial accuracy would be exponential in k 2 rather than in k log k. 56 4.4 The Condition Number Bound The key to the sample-complexity and run-time bounds for our algorithm lies in the following condition number bound for the Hadamard Extension. Theorem 71. Let m have k − 1 rows, every row ζ-separated, and let M = H(m). List the singular values of Mπ⊙ in decreasing order: σ1 ≥ σ2 ≥ . . ., and note σk+1 = . . . = σ2k−1 = 0. Then √ (min πi )(ζ/2 5)k−1 σk ≥ p . k(k 2 + 1)k−1 Definition 72. Hp = {H(n) : n ∈ R[k−1]×[p] }. (So H1 consists of rank-1 tensors of order k − 1.) The proof of Theorem 71 relies on the following insight. Since M has dimensions 2k−1 × k, σk characterizes the least norm of Mv ranging over any unit vector v, but it does not characterize the least norm of vectors of the form hM; the left-kernel of M is of course very large. The insight is that it does become possible to bound σk in terms of such vectors h, provided h is restricted to rank 1 tensors. Proof. Let τ (m, π) = min0̸=h∈H1 hT Mπ⊙ / ∥h∥. The proof rests on two key lemmas: q Q Lemma 73. σk (Mπ⊙ ) ≥ τ (m, π)/ k 1k−1 (1 + ∥mi∗ ∥2 ). Proof. Consider v ∈ Rk , ∥v∥ = 1, achieving σk , i.e.,r := H(m)π⊙ v satisfies  √ Pk−1 1 ∥r∥ = σk . W.l.o.g. ∥vk ∥ ≥ 1/ k. Now H(m)∗k = πk vk r − j=1 πj vj H(m)∗j . [k−1] Now we carefully choose by ni := Q h ∈ H1 based on m. Define n ∈ R k−1 −1/mii ; and let h := mii H(n). For j ̸= k we have 1 ! ! ! ! k−1 k−1 Y X Y X Y Y hT H(m)∗j = mii nS mS,j = mii (−1)|S| 1/mii mij 1 1 S ! = X Y (−1)|S| mii S T i∈S / S ! Y i∈S mij = i∈S i∈S k−1 Y (mii − mij ) = 0. i=1 So, h H(m)∗j πj = 0 for j = 1, . . . , k − 1. For j = k we can also show that hT H(m)∗j πj is small: ! k−1 k−1 X

1 1 T X 1 T T h H(m)π⊙ k = h r − πj vj H(m)∗j = h r− πj vj hT H(m)∗j = hT r. vk vk vk j=1 j=1 57 The norm of hT H(m)π⊙ is then upper bounded by hT H(m)π⊙ √ 1 (hT H(m)π⊙ )k = hT r ≤ k ∥h∥ ∥r∥ v v k u k−1 u Y √ = k ∥h∥ σk = σk tk (1 + m2ii ) due to h being a rank 1 tensor = 1 √ Lemma 74. τ (m, π) ≥ (min πi )(ζ/2 5)k−1 . T Proof. Consider any G ∈ H1 , say G = H(g), g ∈ R[k−1] . Then (G H(m)π⊙ )j = q Qk−1 P Qk−1 2 πj S gS mS,j = πj 1 (1+gi mi,j ). We also note that ∥G∥ = 1 (1 + gi ). We now show that there is some j such that Qk−1 1+gi mij √ 2 is large. First, for i=1 1+gi any i for which gi ≥ 21 , there is at most one j s.t. mij < ζ; exclude these j’s. Next, for each i for which gi < 12 , there is at most one j s.t. g1i + mij < ζ/2; exclude these j’s. For the remainder of the argument fix any j which has not been excluded. Since m has k columns while g ∈ Rk−1 , such a j exists. We 1+g m consider three cases for i. In each case we lower bound √ i ij2 . 1+gi 1+g m 1. gi ≥ 1/2, mij ≥ ζ. Then √ i ij2 ≥ mij ≥ ζ. 1+gi q √ (mij −2)2 1+gi mij √ 2. −1/2 ≤ gi < 1/2. Then ≥ ≥ 1/ 5. (This does not 2 5 1+gi depend on j.) 1+g m 3. gi < −1/2, | g1i + mij | ≥ ζ/2. Then √ i ij2 = 1+gi gi ( g1 +mij ) √i gi 1+1/gi2 ≥ 2√ζ 5 . √ We therefore have τ (m, π) ≥ (min πi )(ζ/2 5)k−1 . The proof of Theorem 71 now follows from Lemmas 73, 74, and the fact that ∥mi∗ ∥ ≤ k for any i. 4.5 Analysis of the Algorithm Here we will bound the errors introduced in each step of Algorithm 3, referring often to both the empirical and idealized versions of each quantity in question. First, we will introduce some notation to simplify some intermediate calculations: Definition 75. ζ1 := ζ/9k 3/2 . 58 We will work under the following assumptions when doing the analysis. Assumption 1. All empirical multilinear moments appearing in computations are within an additive ε of their true values. Assumption 2. The sets S, T, T ′ contain ζ-separated observables. Finally, for a matrix M , ∥M ∥ will denote the spectral norm. We will frequently use (implicitly) the upper-bound ∥M ∥ ≤ ∥M ∥F , where ∥·∥F is the Frobenius norm. Bounding the Operator Norm and Condition Number for the Ideal Matrices and Vectors First, we argue that the ideal versions of the matrices under consideration are well-behaved, crucially relying on Theorem 71. Lemma 76. The matrices A, B, and B′ satisfy 1. A∅∗ = B∅∗ = B′∅∗ = 1, the all ones row vector. 2. σk (A), σk (B), σk (B′ ) ≥ ζ1k . 3. σmax (A), σmax (B), σmax (B′ ) ≤ k2k−1 . Moreover, the matrices CBA and CB′ A satisfy 1. σmax (CBA ), σmax (CB′ A ) ≤ 22k−2 . 2. σk (CBA ), σk (CB′ A ) ≥ πmin ζ12k . Proof. This follows immediately from Theorem 71, the definitions of CBA and CB′ A , and the min-max characterization of the first and last singular values. −1 −2k Corollary 77. ∥(CBA )+ ∥ , ∥(CB′ A )+ ∥ ≤ πmin ζ1 . −1 −2k ζ1 , and ∥vi ∥ ≤ ζ1−k . Lemma 78. ∥ui ∥ , ∥u′i ∥ ≤ πmin Proof. Clearly, ∥vi ∥ ≤ 2k−1 ≤ ζ1−k , as vi is a vector of 2k−1 moments of prod−k ucts of Bernoulli random variables. Now ∥ui ∥ = vi C+ C+ BA ≤ ζ1 BA ≤ −1 −2k ′ πmin ζ1 . A similar argument bounds ∥ui ∥. 59 Bounding Error in Derived Quantitites Lemma 79. If ε is sufficiently small, ĈBA − CBA 2 , ĈB′ A − CB′ A 2 , C̃B+B′ ,A − CB+B′ ,A 2 ≤ 23k ε < ζ1−k ε, and + Ĉ+ BA − CBA 2 + , Ĉ+ B′ A − CB′ A 2 −2 −5k ≤ 2πmin ζ1 ε. Proof. To prove the first inequalities we use the fact that ĈBA − CBA ≤ 2 2 C̃BA − CBA since ĈBA − C̃BA ≤ C̃BA − CBA (as the nearest rank k matrix to C̃BA can be no further than the distance from C̃BA to CBA ). We then observe that ∥·∥2 ≤ ∥·∥F and the fact that every entry in the matrix of differences is at most ε in magnitude. For the final inequality we use Lemmas 90 and 91 to get + + Ĉ+ BA − CBA ≤ 2 CBA 2 −2 −5k ĈBA − CBA ≤ 2πmin ζ1 ε. + The same argument bounds Ĉ+ B ′ A − CB ′ A . 2 i Lemma 80. For any i ∈ lg(2k) and j ≤ 2 , −2i −11ik ∥ṽj − vj ∥ , ∥ũj − uj ∥ , ũ′j − u′j ≤ πmin ζ1 ε. Proof. Recall that we initialize the algorithm with ṽ1 ← (g̃(R ∪ {1}))R∈A , ũ1 ← ṽ1 Ĉ+ BA , ũ′1 ← ṽ1 Ĉ+ BA . First, we observe that ∥ṽ1 − v1 ∥ ≤ ε by assumption. Since ũ1 , ũ′1 are computed in the same manner here as in the loop, we will bound that error in the induction. Now assume that the claim holds up to i − 1, and let εi−1 be the bound obtained in the i − 1st step. Recall that in each iteration of the outer loop we compute ṽ2i ← (ũ2i−1 ⊗ ũ′2i−1 )C̃B+B′ ,A , ũ2i ← ṽ2i Ĉ+ BA , ũ′2i ← ṽ2i Ĉ+ B′ A . We also observe that we choose ε so that ∥ũj − uj ∥ ≪ ∥uj ∥ and ∥ṽi − vi ∥ ≪ −1 −2k ∥vj ∥ for all i, j so that we can upper bound ∥ũj ∥ ≤ 2 ∥uj ∥ ≤ 2πmin ζ1 and −k ∥ṽi ∥ ≤ 2 ∥vi ∥ ≤ 2ζ1 . We will first focus on bounding ∥ṽ2i − v2i ∥. To do this we write ṽ2i − v2i = (ũ2i−1 ⊗ ũ′2i−1 )C̃B+B′ ,A − (u2i−1 ⊗ u′2i−1 )CB+B′ ,A . 60 We now let w = ũ2i−1 − u2i−1 , w′ = ũ′2i−1 − u2i−1 , and E = C̃B+B′ ,A − CB+B′ ,A , and use the triangle inequality to obtain ∥ṽ2i − v2i ∥ ≤ (w ⊗ ũ′2i−1 )C̃B+B′ ,A + (ũ2i−1 ⊗ w′ )C̃B+B′ ,A +∥(ũ2i−1 ⊗ ũ′2i−1 )E∥ . The first two terms can each be bounded by −1 −2k 3k ζ1 2 εi−1 , ∥w∥ 2 ∥u1 ∥ C̃B+B′ ,A ≤ 2πmin and the final term is bounded by −2 −4k −k 4 ∥u1 ∥2 ∥E∥ ≤ 4πmin ζ1 ζ1 ε. As a result we get that −2 −5k ζ1 (εi−1 + ε). ∥ṽ2i − v2i ∥ ≤ 4πmin Now we can bound ∥ũ2i − u2i ∥ by observing that + ũ2i − u2i = ṽ2i Ĉ+ BA − v2i CBA . + Let z = ṽ2i − v2i and let D = Ĉ+ BA − CBA . The above equation becomes + + ũ2i − u2i = (v2i + z)(C+ BA + D) − v2i CBA = v2i D + zCBA + zD, and after taking norms and using the triangle inequality we obtain ∥ũ2i − u2i ∥ ≤ ∥v2i D∥ + zC+ BA + ∥zD∥ . By Corollary 77, Lemma 78 and the induction hypothesis, we get −2 −k −1 −2k −2 −5k ∥ũ2i − u2i ∥ ≤ πmin ζ1 2ζ1−5k ε + ∥ṽ2i − v2i ∥ πmin ζ1 + 2πmin ζ1 ε
−2 −6k −2 −5k −1 −2k −2 −5k ≤ 2πmin ζ1 ε + 4πmin ζ1 (εi−1 + ε)(πmin ζ1 + 2πmin ζ1 ε) −2 −10k ≤ 32πmin ζ1 εi−1 −2 −11k ≤ πmin ζ1 εi−1 , where we use the fact that εi−1 ≥ ε. For j not a power of 2, we can do the same analysis, and since the error bound is increasing in j, the result will follow. Corollary 81. Algorithm 3 will produce vectors ṽi for i ≤ 2k satisfying −O(lg k) −O(k lg k) ∥ṽi − vi ∥ ≤ πmin ζ ε. 61 Applying the k-MixIID Algorithm. We use the following theorem. Theorem 82 (A slight restatement of Theorem 22 from Chapter 2). Given a mixture M = (m, π) of k Bernoulli random variables with probabilities m1 , . . . , mk and mixing probabilities π1 , . . . , πk , respectively, let [Hk+1 ]ki,j=0 = µi+j be the matrix of moments of the distribution. If m is ζ-separated, then there is an algorithm, LearnPowerDistribution, that takes as input a Hankel matrix [H̃k+1 ]ki,j=0 = µ̃i+j of approximate moments of M satisfying H̃k+1 − Hk+1 2 ≤ πmin 2−γ ζ 16k (for some γ ≥ 1), and outputs a model M̃ = (m̃, π̃) satisfying ∥m̃ − m∥∞ , ∥π̃ − π∥∞ ≤ 2−γ −1 using O(k 2 log k + k log2 k · log(log ζ −1 + log πmin + γ)) arithmetic operations. Corollary 83. The output (m̃1 , π̃) of LearnPowerDistribution in line 14 of Algorithm 3 satisfies −O(lg k) −O(k lg k) ∥m̃1 − m1 ∥ , ∥π̃ − π∥ ≤ πmin ζ ε. −O(lg k) −O(k lg k) Proof. Every entry (ṽi ){1} satisfies ∥(ṽi )1 − (vi )1 ∥ ≤ πmin −O(lg k) −O(k lg k) H̃k+1 − Hk+1 ≤ πmin ζ ζ ε so ε which implies that −O(lg k) −O(k lg k) ∥m̃1 − m1 ∥∞ , ∥π̃ − π∥∞ ≤ πmin ζ ε. Finally, we need to multiply by k ≤ ζ −1 to account for the conversion to the Euclidean norm to get the stated bound. Solving for the Rest of the Model. Once we have computed m̃1 and π̃, we will use them to compute the remaining model parameters. In this section we bound the additional error introduced by these computations. Proposition 84. ∥Vdm(m̃1 ) − Vdm(m1 )∥ ≤ k ∥m̃1 − m1 ∥ ≤ ζ −1 ∥m̃1 − m1 ∥. Claim 85 (Claim 26 in [47]). ∥Vdm(m1 )−1 ∥ ≤ 2k /ζ k−1 ≤ ζ −2k when m1 is ζ-separated. 62 Lemma 86. The computed à produced by Algorithm 3 will satisfy −O(lg k) −O(k lg k) à − A ≤ πmin ζ ε. Proof. Recall Ṽ = (ṽ0 ; . . . ; ṽk−1 ) from Algorithm 3 and V = (v0 ; . . . ; vk−1 ) is its real-value analog. First, we observe that Ṽ ≤ ζ −3k . Also, by Lemma 90 −1 −1 and Claim 85, ∥Vdm(m̃1 )−1 ∥ ≤ ζ −3k . Now, by Corollary 83, π̃⊙ − π⊙ ≤ −O(lg k) −O(k lg k) πmin ζ ε. Thus, −1 −1 π̃⊙ − π⊙ (Vdm(m̃1 ))−1 −O(lg k) −O(k lg k) ≤ πmin Ṽ ζ ε. −O(lg k) −O(k lg k) Now, by Lemma 91, ∥(Vdm(m̃1 ))−1 − (Vdm(m1 ))−1 ∥ ≤ πmin so, −1 π̃⊙ (Vdm(m̃1 ))−1 − (Vdm(m1 ))−1 Finally, Ṽ − V −1 π⊙ −O(lg k) −O(k lg k) ≤ πmin ζ (Vdm(m̃1 ))−1 Ṽ ζ −O(lg k) −O(k lg k) ≤ πmin ζ ε, ε. ε, so that Ṽ − V −O(lg k) −O(k lg k) ≤ πmin ζ ε. Putting these together, we easily obtain −1 −1 (Vdm(m1 ))−1 V à − A = π̃⊙ (Vdm(m̃1 ))−1 Ṽ − π⊙ −1 −1 ≤ π̃⊙ − π⊙ −1 + π⊙ (Vdm(m̃1 ))−1 (Vdm(m̃1 ))−1 + ∥w∥ Vdm(m̃1 )−1 −O(lg k) −O(k lg k) ≤ πmin ζ −1 Ṽ + π̃⊙ (Vdm(m̃1 ))−1 − (Vdm(m1 ))−1 Ṽ − V ( Vdm(m1 )−1 E2 + ∥E1 V ∥∞ ) Ṽ ε. Lemma 87. (ÃT )+ ≤ ζ −O(k) and (B̃T )+ ≤ ζ −O(k) . Proof. We prove this for (ÃT )+ , the argument for (B̃T )+ is the same. Notice that Lemma 86 together with the lower bound on σk (A) from Theorem 71 imply that ÃT has full row rank. Now, (ÃT )+ ≤ (AT )+ + (ÃT )+ − (AT )+ . Ṽ 63 We bound each term separately. We have (AT )+ = A(AAT )−1 ≤ ∥A∥ (AAT )−1 ≤ k2k−1 /(σk (A))2 ≤ ζ −O(k) . As for the second term, notice that ÃÃT − AAT = Ã(ÃT − AT ) + (à − A)AT ≤ à AT ÃT − AT + à − A −O(lg k) −O(k lg k) ≤ πmin ζ ε. −O(lg k) −O(k lg k) Therefore, by Lemma 91 also (ÃÃT )−1 − (AAT )−1 Thus, we get for the second term ≤ πmin ζ ε. (ÃT )+ − (AT )+ = Ã(ÃÃT )−1 − A(AAT )−1 = (à − A)(ÃÃT )−1 + A((ÃÃT )−1 − (AAT )−1 ) (ÃÃT )−1 + ∥A∥ (ÃÃT )−1 − (AAT )−1 ≤ à − A −O(lg k) −O(k lg k) ≤ πmin ζ ε. −O(lg k) −O(k lg k) Lemma 88. B̃ − B ≤ πmin ζ ε. Proof. We can bound B̃ − B using the same tools as in the previous bounds. First, we bound C̃BA − CBA −1 −O(k) −1 π̃⊙ ≤ πmin ζ ε. (ÃT )+ Next, C̃BA (ÃT )+ − (AT )+ C̃BA (ÃT )+ −O(lg k) −O(k lg k) −1 π̃⊙ ≤ πmin ζ ε. Finally, −O(lg k) −O(k lg k) −1 −1 π̃⊙ − π⊙ ≤ πmin ζ ε. The resulting bound is −1 B̃ − B = C̃BA (ÃT )+ π̃⊙ − CBA (AT )+ π⊙ ≤ C̃BA − CBA (ÃT )+ −1 π̃⊙ + + C̃BA (ÃT )+ − (AT )+ + C̃BA (ÃT )+ −O(lg k) −O(k lg k) ≤ πmin ζ ε. −1 π̃⊙ + −1 −1 π̃⊙ − π⊙ 64 Theorem 89. Algorithm 3 will compute m̃i satisfying, for all i ∈ [n], −O(lg k) −O(k lg k) ∥m̃i − mi ∥∞ ≤ πmin ζ ε. Proof. We will compute the bound using the inversion of B̃T since this will give us the worst case. Let ỹ = (g̃(R ∪ {i}))R∈B , and let y = (g(R ∪ {i}))R∈B . We note that ∥ỹ − y∥ ≤ ζ −O(k) ε, by assumption on the sample size, and ∥ỹ∥ ≤ 2k−1 . Then, ∥ỹ − y∥ (B̃T )+ ∥ỹ∥ (B̃T )+ − (BT )+ −O(1) −1 π̃⊙ ≤ πmin ζ −O(k) ε, −O(lg k) −O(k lg k) −1 π̃⊙ ≤ πmin ζ ε, and ∥ỹ∥ (B̃T )+ −O(lg k) −O(k lg k) −1 −1 − π⊙ ≤ πmin π̃⊙ ζ ε. So that we get −1 −1 ∥m̃i − mi ∥ = ỹ(B̃T )+ π̃⊙ − y(BT )+ π⊙ −1 π̃⊙ + ∥ỹ∥ (B̃T )+ − (BT )+ ≤ ∥ỹ − y∥ (B̃T )+ −O(lg k) −O(k lg k) ≤ πmin ζ −1 π̃⊙ + ∥ỹ∥ (B̃T )+ ε. Miscellaneous Proofs Lemma 90. Let M, M̃ ∈ Rn×n have rank k ≤ n. Let M̃ − M = ε ≤ σk (M )/2. Then M̃ + ≤ 2 ∥M + ∥. Proof. We observe that M̃ + = 1 1 ≤ ≤ 2 M+ . σk (M ) − σk (M )/2 σk (M̃ ) Lemma 91 (Theorem 3.4 in [49]). Let M, M̃ ∈ Rn×n have rank k ≤ n. Then M + − M̃ + ≤ 2 ∥M + ∥ M̃ + M − M̃ . −1 −1 − π⊙ π̃⊙ 65 Chapter 5 THE IDENTIFIABILITY OF UNIFORM MIXTURES OF BINOMIAL DISTRIBUTIONS WITH LOG-LINEAR INFLUENCES In the work thus far we have been exclusively concerned with the setting where there is a single latent variable U which induces independence among X1 , . . . , Xn when it is conditioned upon. In general Bayesian networks we are unlikely to have a single predecessor for most nodes and so this setting is insufficiently general to handle many cases of interest. However, this setting does allow us to identify the parameters of the joint distribution on U and the Xi using only O(k) moments, when U is supported on k states. Towards the goal of eventually understanding identifiability in more general graphs, we will return the k-Mix IID setting and start generalizing from there. We begin by introducing multiple latent predecessors U1 , . . . , Uℓ ; we will have each visible variable X be conditionally independent from the others given the values of each of the Ui . With ℓ hidden variables having k states each, we then obtain a setting in which there k ℓ possible assignments to the latent variables, and each such assignment induces a (potentially different) probability that X = 1. One method for identifying parameters in this case (excluding the measure zero set of parameters for which the probability that X = 1 is the same across two or more assignments to U1 , . . . , Uℓ ) is to run a k ℓ -Mix IID algorithm. Unfortunately, doing this requires observing O(k ℓ ) copies of X and requires learning k ℓ parameters. Here, there is no advantage to be gained by knowing that there are multiple independent latent variables, as opposed to a single latent variable with k ℓ states. We will call this the naïve setting. We would like to understand when we can hope to take advantage of the independence structure among the latent variables, which requires that there be some observable trace of this indpendence structure in the Bernoulli parameters induced by the assignments to the latent variables. Towards that end, we will consider the case where the contributions to the Bernoulli parameters induced by an assignment can be decoupled into a contribution 66 for each latent variable. The simplest way to do this is to have Pr(X = 1 | Q u1 , . . . , uℓ ) = ℓi=1 αi,ui for parameters {αi,j }i,j∈[ℓ]×[k] so that the contributions from each latent variable are multiplied together; changing one latent variable only changes that term in the product. This model is already interestingly different from the naïve setting when k = 2. Then we have 2ℓ + ℓ parameters describing the joint distribution in this setting, compared to 2ℓ + ℓ parameters in the naïve setting. We can hope to do better than just running Prony’s algorithm with a mixture of 2ℓ distributions in this case. This assumption of “independent influence” of the latent variables on the observed variables makes the resulting model into a log-linear model; we have ℓ X log Pr(X = 1 | u1 , . . . , uℓ ) = log αi,ui . i=1 While we have a log-linear dependence of the probability on the parameters, these are unnormalized, in that the probability of X = 1 | u1 , . . . , uℓ does not depend on the sum of the probabilities of X = 1 across all assignments to U1 , . . . , Uℓ . As a result, these are in fact distinct from many of the models that go under the heading of “log-linear models”, such as restricted Boltzmann machines or Markov random fields. Formalizing our model We are concerned with the following mixture model: We have n binary random variables X1 , . . . , Xn that are iid conditioned on a collection of ℓ independent latent binary variables U1 , . . . , Uℓ . Each Ui is uniform over {0, 1}.1 Finally, we have that the the probability that X = 1 given U = u is given by the product of terms corresponding to each Q of the hidden variables, namely P (X = x | U = u) = ℓj=1 αj,uj . We have 2ℓ parameters αib for i = 1, 2, . . . , ℓ and b = 0, 1. The probability that X = 1 is given by 2−ℓ times ℓ Y µ1 := (αi0 + αi1 ). i=1 With enough copies of X we have access to ℓ Y j j µj := (αi0 + αi1 ) i=1 1 Relaxing this assumption makes the analysis quite a bit more difficult, and we do not know if the same result will hold in that setting. 67 for arbitrary j.2 We are going to throw out the restrictions on the statistical model, and consider the problem of solving the polynomial system. We first observe that µ1 is invariant to the transformation (αi0 , αi1 ) 7→ (λαi0 , λαi1 ), (αj0 , αj1 ) 7→ (λ−1 αi0 , λ−1 αi1 ) for any λ > 0 and i ̸= j (where all other parameters are unchanged). As a result, we can without loss of generality solve for parameters satisfying 1 1 αi0 + αi1 = µ1ℓ for all i ∈ [ℓ]. Let γ := µ1ℓ . To solve for αi0 , αi1 it now suffices to solve for ai := αi0 αi1 , the product of our two parameters. In particular, we get that αi0 (γ − αi0 ) = ai which we can easily solve for αi0 . We will obtain the following theorem: Theorem 92. The mapping from {ai } to µ2 , . . . , µℓ+1 is identifiable. Our immediate goal will now be to understand the polynomial system of equations in the variables {ai } that we obtain after this transformation. First, we observe that 2 2 αi0 + αi1 = (αi0 + αi1 )2 − 2αi0 αi1 = γ 2 − 2ai so that ℓ Y µ2 = (γ 2 − 2ai ). i=1 Moreover, 2 2 3 3 = (αi0 + αi1 )3 − 3(αi0 αi1 + αi0 αi1 ) αi0 + αi1 = γ 3 − 3αi0 αi1 (αi0 + αi1 ) = γ 3 − 3γai . In subsequent computations, I’ll omit the subscript ‘i’ everywhere since it remain constant throughout. 2 The actual observation we can make is of µj /2ℓ , which we can multiply by 2ℓ to get µj . 68 Definition 93. Let pm (a) denote the univariate polynomial in a obtained by simplifying α0n + α1n using α0 + α1 = γ and a = α0 α1 . Proposition 94. p0 = 2, (5.1) p1 = γ, (5.2) pm = γpm−1 − apm−2 , m ≥ 2. (5.3) Proof. First, p0 = α00 + α10 = 2 and p1 = α01 + α11 = γ. Now pm = α0m + α1m = (α0m−1 + α1m−1 )(α0 + α1 ) − α0 α1 (α0m−2 + α1m−2 ) = γpm−1 − apm−2 . Observation 95. The polynomials defined above satisfy the following: 1. The degree of pm is ⌊m/2⌋. 2. pm (0) = γ m for m > 0. 3. The leading term of pm has a negative sign if ⌊m/2⌋ is odd, and a positive sign if ⌊m/2⌋ is even. Proof. By induction over m. Proposition 96. For any m > 2, 1. the roots r1 , . . . , r⌊m/2⌋ of pm are all real, distinct, and contained in (0, ∞). 2. Let s1 , . . . , s⌊(m−1)/2⌋ be the roots of pm−1 . Then 0 < r1 < s1 and si−1 < ri < si for i = 2, . . . , ⌊m/2⌋ and s⌊(m−1)/2⌋ < r⌊m/2⌋ when pm has degree greater than pm−1 . Proof. We will induct over m. The claim clearly holds for m = 0, 1. Now p2 = −2a + γ 2 which has a root at γ 2 /2 and p3 = γ(−2a + γ 2 ) − aγ = −2γa + γ 3 − γa = −3γa + γ 3 which has a root at γ 2 /3. 69 Now fix an arbitrary m ∈ N. Let d := deg(pm−2 ), so d + 1 = deg(pm ). Let s1 , . . . , sd be the roots of pm−2 and let r1 , . . . , rd be the first d roots of pm−1 . By the inductive hypothesis, 0 < r1 < s1 < r2 < · · · < rd < sd . If deg(pm−1 ) = d + 1 > deg(pm ), then there is an additional root rd+1 of pm−1 with sd < rd+1 . Since we have accounted for every root of pm−1 and pm−2 , it must be the case that the value of pm−1 alternates between strictly positive and strictly negative on the sequence of open intervals (−∞, r1 ), (r1 , r2 ), (r2 , r3 ), . . . , (r⌊(m−1)/2⌋ , ∞) and pm−2 alternates in sign on the intervals (−∞, s1 ), (s1 , s2 ), . . . , (sd , ∞). Now we compute pm (r1 ) = γpm−1 (r1 ) − r1 pm−2 (r1 ) = −r1 pm−2 (r1 ) < 0; pm (s1 ) = γpm−1 (s1 ) − s1 pm−2 (s1 ) = γpm−1 (s1 ) < 0. By the preceding observation, pm (0) = γ m > 0 so there must be a root of pm in (0, r1 ). For 1 < i < d, ri ∈ (si−1 , si ) and si ∈ (ri , ri+1 ). Moreover, pm (ri ) = −ri pm−2 (ri ) and pm (si ) = γpm−1 (si ), so sign(pm (ri )) = sign(pm (si )) = − sign(pm (ri+1 )). We conclude that there is a root of pm in the interval (si , ri+1 ) for 1 < i < d. we have thus shown that there are roots of pm in each of the intervals (0, r1 ), (s1 , r2 ), (s2 , r3 ), . . . , (sd−1 , rd ). If deg(pm−1 ) = d+1, then we also get that there is a root in (sd , rd+1 ) and we are done. If deg(pm−1 ) = d, then the leading term of pm has a different sign than the leading terms of pm−1 and pm−2 . Since sign(pm (sd )) = sign(pm−1 (sd )) and sd is greater than all the roots of pm−1 it must be the case that sign(pm−1 (x)) = sign(pm−1 (sd )) for all x ∈ [sd , ∞). But limx→∞ pm (x) = − limx→∞ pm−1 (x), so there must be a root of pm in (sd , ∞). we have thus accounted for all d + 1 roots of pm , proving the claim. The identifiability of symmetric systems of polynomial equations Consider a sequence of ℓ univariate polynomials p1 , . . . , pℓ . Let x1 , . . . , xℓ be indeterminates. We can construct symmetric polynomials in the xi by taking products as follows: ℓ Y qi (x) := pi (xj ). j=1 70 Proposition 97. Suppose that for all i ∈ [ℓ], pi has a smallest positive real root αi with multiplicity 1 that is not a root of p1 , p2 , . . . , pi−1 . Then the mapping (y1 , . . . , yℓ ) 7→ (q1 (y), . . . , qℓ (y)) is invertible, except on a set of measure zero. Proof. In what follows, let #(αj , pi ) denote the multiplicity of root αj in pi . By assumption, #(αj , pi ) = 0 for j > i and #(αj , pj ) = 1 for all j. We will now construct a sequence of rational functions r1 , . . . , rℓ satisfying #(αj , ri ) = 1i=j , where #(αj , ri ) < 0 if αj is a pole of ri . First, we set r1 := p1 , since #(αj , p1 ) = 0 for all j > 1. Now inductively we construct ri as follows: ri := pi i−1 Y −#(α ,p ) ri′ i′ i ! . i′ =1 By construction, #(αj , ri ) = 0 for j < i and #(αi , ri ) = 1. Moreover, #(αj , ri ) = 0 for j > i since #(αj , pi ) = #(αj , r1 ) = · · · = #(αj , ri−1 ) = 0. Q Now define si := ℓj=1 ri (yj ) for i = 1, . . . , ℓ so that si is the symmetric product of ri evaluated at each indeterminate, just as qi is the symmetric product of pi evaluated at each indeterminate. In fact, we have si = qi ! i−1 Y −#(αi′ ,pi ) si′ . i′ =1 Let F and G be the following mappings: F G (y1 , . . . , yℓ ) 7− → (q1 , . . . , qℓ ) 7− → (s1 , . . . , sℓ ). Finally, let H = G ◦ F . We now consider the Jacobian of H, evaluated at the point α = (α1 , . . . , αℓ ). 71 By construction we have that ∂si (α) = ∂yj ! Y ri (αj ′ ) ri′ (αj ). j ′ ̸=j Now Y ri (αj ′ ) ̸= 0 ⇐⇒ i = j j ′ ̸=j since ri (αi ) = 0 and #(αj , ri ) = 0 for any j ̸= i. Moreover, ri′ (αi ) ̸= 0 since ∂si αi is a simple root of ri , so we conclude that ∂y (α) ̸= 0 ⇐⇒ i = j. Thus, j the Jacobian is a diagonal matrix with non-zero diagonal entries and is thus invertible. We conclude that there is an open neighborhood around α in which H is invertible, which implies the same for F . Since F is a polynomial map, we conclude it is generically identifiable. Proof of Theorem 92. Finally, Theorem 92 follows immediately from Propositions 94 and 97. Q Conjecture 98. Let p1 , . . . , pℓ be univariate polynomials and let qi := ℓj=1 pi (xj ) for i = 1, . . . , ℓ. Let α1 , . . . , αL be the set of all roots appearing in any pi . Consider the matrix M ∈ Rℓ×L with entry Mij ∈ Z≥0 being the multiplicity with which αj appears as a root in pi . Then the mapping (x1 , . . . , xℓ ) 7→ (q1 (x), . . . , qℓ (x)) is invertible if and only if M has rank ℓ over Q. 72 BIBLIOGRAPHY [1] F. B. Hildebrand, Introduction to Numerical Analysis, 2nd. McGraw-Hill, 1974. [2] B. Simon, A Comprehensive Course in Analysis. American Mathematical Society, 2015. [3] Y. Kim, F. Koehler, A. Moitra, E. Mossel, and G. Ramnarayan, “How many subpopulations is too many? Exponential lower bounds for inferring population histories,” in Int’l Conference on Research in Computational Molecular Biology, ser. Lecture Notes in Computer Science, L. Cowen, Ed., vol. 11457, Springer, 2019, pp. 136–157. doi: 10.1007/9783-030-17083-7_9. [4] A. Feldmann and W. Whitt, “Fitting mixtures of exponentials to longtail distributions to analyze network performance models,” Performance Evaluation, vol. 31, no. 3, pp. 245–279, 1998. [5] J. Pearl, Causality, 2nd. Cambridge, 2009. [6] A. Anandkumar, D. J. Hsu, and S. M. Kakade, “A method of moments for mixture models and hidden Markov models,” in Proceedings of the 25th Annual Conference on Learning Theory - COLT, ser. JMLR Proceedings, vol. 23, 2012, pp. 33.1–33.34. [Online]. Available: http: //proceedings.mlr.press/v23/anandkumar12/anandkumar12.pdf. [7] R. de Prony, “Essai expérimentale et analytique,” J. Écol. Polytech., vol. 1, no. 2, pp. 24–76, 1795. [8] R. Kumaresan, D. W. Tufts, and L. L. Scharf, “A Prony method for noisy data: Choosing the signal components and selecting the order in exponential signal models,” Proceedings of the IEEE, vol. 72, no. 2, pp. 230–233, 1984. [9] Y. Rabani, L. J. Schulman, and C. Swamy, “Learning mixtures of arbitrary distributions over large discrete domains,” in Proceedings of the 5th Conference on Innovations in Theoretical Computer Science, 2014, pp. 207–224. doi: 10.1145/2554797.2554818. [10] J. Li, Y. Rabani, L. J. Schulman, and C. Swamy, “Learning arbitrary statistical mixtures of discrete distributions,” in Proceedings of the 47th Annual ACM Symposium on Theory of Computing, 2015, pp. 743–752. doi: 10.1145/2746539.2746584. [11] T. Hofmann, “Probabilistic latent semantic analysis,” in Proceedings of the 15th Conference on Uncertainty in Artificial Intelligence, 1999, pp. 289– 296. 73 [12] C. Papadimitriou, P. Raghavan, H. Tamaki, and S. Vempala, “Latent semantic indexing: A probabilistic analysis,” Journal of Computer and System Sciences, vol. 61, no. 2, pp. 217–235, 2000. [13] S. Arora, R. Ge, and A. Moitra, “Learning topic models — going beyond SVD,” in Proceedings of the 53rd Annual IEEE Symposium on Foundations of Computer Science, 2012. [14] A. Anandkumar, D. P. Foster, D. J. Hsu, S. M. Kakade, and Y.-K. Liu, “A spectral algorithm for latent Dirichlet allocation,” in Advances in Neural Information Processing Systems 25, F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, Eds., Curran Associates, Inc., 2012, pp. 917–925. [15] M. Kearns, Y. Mansour, D. Ron, R. Rubinfeld, R. Schapire, and L. Sellie, “On the learnability of discrete distributions,” in Proceedings of the 26th Annual ACM Symposium on Theory of Computing, 1994, pp. 273–282. doi: 10.1145/195058.195155. [16] Y. Freund and Y. Mansour, “Estimating a mixture of two product distributions,” in Proceedings of the 12th Annual Conference on Computational Learning Theory, Jul. 1999, pp. 53–62. doi: 10.1145/307400.307412. [17] M. Cryan, L. Goldberg, and P. Goldberg, “Evolutionary trees can be learned in polynomial time in the two state general Markov model,” SIAM J. Comput., vol. 31, no. 2, pp. 375–397, 2001. doi: 10 . 1137 / S0097539798342496. [18] K. Chaudhuri and S. Rao, “Learning mixtures of product distributions using correlations and independence,” in Proceedings of the 21st Annual Conference on Learning Theory - COLT, Omnipress, 2008, pp. 9–20. [Online]. Available: http://colt2008.cs.helsinki.fi/papers/7Chaudhuri.pdf. [19] J. Feldman, R. O’Donnell, and R. A. Servedio, “Learning mixtures of product distributions over discrete domains,” SIAM J. Comput., vol. 37, no. 5, pp. 1536–1564, 2008. doi: 10.1137/060670705. [20] S. Chen and A. Moitra, “Beyond the low-degree algorithm: Mixtures of subcubes and their applications,” in Proceedings of the 51st Annual ACM Symposium on Theory of Computing, 2019, pp. 869–880. [21] J. Peters, D. Janzing, and B. Schölkopf, Elements of Causal Inference. MIT Press, 2017. [22] T. S. Chihara, An Introduction to Orthogonal Polynomials. Gordon and Breach, 1978. [23] K. Schmüdgen, The Moment Problem (Graduate Texts in Mathematics). Springer International Publishing, 2017, vol. 277. 74 [24] Wikipedia contributors, Hoeffding’s inequality — Wikipedia, the free encyclopedia, [Online; accessed 6-April-2020], 2020. [Online]. Available: https://en.wikipedia.org/w/index.php?title=Hoeffding%27s_ inequality&oldid=946643344. [25] T. J. Rivlin, “Bounds on a polynomial,” Journal of Research of the National Bureau of Standards - B. Math. Sci., vol. 74B, no. 1, pp. 47–54, Jan. 1970. [26] F. D. Parker, “Inverses of Vandermonde matrices,” The American Mathematical Monthly, vol. 71, no. 4, pp. 410–411, 1964. [27] J. H. Wilkinson, “The perfidious polynomial,” in Studies in Numerical Analysis, ser. Studies in Mathematics, G. H. Golub, Ed., vol. 24, Mathematical Association of America, 1984, pp. 1–28. [28] V. Y. Pan and Z. Q. Chen, “The complexity of the matrix eigenproblem,” in Proceedings of the 31st Annual ACM Symposium on Theory of Computing, 1999, pp. 507–516. [29] W. Gautschi, “How (un)stable are Vandermonde systems,” Asymptotic and computational analysis, vol. 124, pp. 193–210, 1990. [30] S. G. Bartels and D. J. Higham, “The structured sensitivity of Vandermondelike systems,” Numerische Mathematik, vol. 62, pp. 17–33, 1992. [31] V. Y. Pan, Z. Q. Chen, and A. Zheng, “The complexity of the algebraic eigenproblem,” Math. Science Research Institute, Berkeley California, Tech. Rep. 1998-071, 1998. [32] B. N. Parlett, The Symmetric Eigenvalue Problem. Prentice-Hall, Inc., 1998. [33] V. Y. Pan, “Optimal and nearly optimal algorithms for approximating polynomial zeros,” Computers & Mathematics with Applications, vol. 31, no. 12, pp. 97–138, 1996. [34] S. L. Gordon, B. Mazaheri, Y. Rabani, and L. J. Schulman, “Source identification for mixtures of product distributions,” in Proceedings of the 34th Annual Conference on Learning Theory - COLT, ser. Proceedings of Machine Learning Research, vol. 134, PMLR, 2021, pp. 2193–2216. [Online]. Available: http://proceedings.mlr.press/v134/gordon21a. html. [35] S. Newcomb, “A generalized theory of the combination of observations so as to obtain the best result,” American Journal of Mathematics, vol. 8, no. 4, pp. 343–366, 1886. [36] K. Pearson, “Contributions to the mathematical theory of evolution III,” Philosophical Transactions of the Royal Society of London (A.), vol. 185, pp. 71–110, 1894. 75 [37] B. S. Everitt and D. J. Hand, “Mixtures of discrete distributions,” in Finite Mixture Distributions, Dordrecht: Springer Netherlands, 1981, pp. 89–105. [38] D. M. Titterington, A. F. M. Smith, and U. E. Makov, Statistical Analysis of Finite Mixture Distributions. John Wiley and Sons, Inc., 1985. [39] B. G. Lindsay, Mixture Models: Theory, Geometry and Applications. 1995, pp. i–163. [40] G. J. McLachlan, S. X. Lee, and S. I. Rathnayake, “Finite mixture models,” Annual Review of Statistics and Its Application, vol. 6, no. 1, pp. 355–378, 2019. doi: 10.1146/annurev-statistics-031017-100325. [41] B. Tahmasebi, S. A. Motahari, and M. A. Maddah-Ali, “On the identifiability of finite mixtures of finite product measures,” (Also in “On the identifiability of parameters in the population stratification problem: A worst-case analysis,” Proceedings of the ISIT 2018 pp. 1051-1055), 2018. [Online]. Available: https://arxiv.org/abs/1807.05444. [42] B. Tahmasebi, S. A. Motahari, and M. A. Maddah-Ali, “On the identifiability of finite mixtures of finite product measures,” (Also in “On the identifiability of parameters in the population stratification problem: A worst-case analysis,” Proceedings of the ISIT’18 pp. 1051-1055.), 2018. [Online]. Available: https://arxiv.org/abs/1807.05444. [43] S. Chen and A. Moitra, “Beyond the low-degree algorithm: Mixtures of subcubes and their applications,” in Proceedings of the 51st Annual ACM Symposium on Theory of Computing, 2019, pp. 869–880. doi: 10. 1145/3313276.3316375. [44] S. L. Gordon, B. Mazaheri, Y. Rabani, and L. J. Schulman, “Source identification for mixtures of product distributions,” in Proceedings of the 34th Annual Conference on Learning Theory - COLT, ser. Proceedings of the Machine Learning Research, vol. 134, PMLR, 2021, pp. 2193– 2216. [Online]. Available: http://proceedings.mlr.press/v134/ gordon21a.html. [45] D. Hsu, S. M. Kakade, and T. Zhang, “A spectral algorithm for learning hidden Markov models,” J. Comput. Syst. Sci., vol. 78, no. 5, pp. 1460– 1480, Sep. 2012. doi: 10.1016/j.jcss.2011.12.025. [46] S. L. Gordon and L. J. Schulman, “Hadamard extensions and the identification of mixtures of product distributions,” IEEE Transactions on Information Theory, vol. 68, no. 6, pp. 4085–4089, 2022. doi: 10.1109/ TIT.2022.3146630. [47] S. L. Gordon, B. Mazaheri, Y. Rabani, and L. J. Schulman, “The sparse Hausdorff moment problem, with application to topic models,” 2020. [Online]. Available: https://arxiv.org/abs/2007.08101. 76 [48] G. H. Golub and C. F. V. Loan, Matrix Computations, 4th. The Johns Hopkins University Press, 2013. [49] G. W. Stewart, “On the perturbation of pseudo-inverses, projections and linear least squares problems,” SIAM Review, vol. 19, no. 4, pp. 634– 662, 1977. doi: 10.1137/1019104.