Learning Latent Tree Graphical Models

We study the problem of learning a latent tree graphical model where samples are available only from a subset of variables. We propose two consistent and computationally efficient algorithms for learning minimal latent trees, that is, trees without any redundant hidden nodes. Unlike many existing methods, the observed nodes (or variables) are not constrained to be leaf nodes. Our first algorithm, recursive grouping, builds the latent tree recursively by identifying sibling groups using so-called information distances. One of the main contributions of this work is our second algorithm, which we refer to as CLGrouping. CLGrouping starts with a pre-processing procedure in which a tree over the observed variables is constructed. This global step groups the observed nodes that are likely to be close to each other in the true latent tree, thereby guiding subsequent recursive grouping (or equivalent procedures) on much smaller subsets of variables. This results in more accurate and efficient learning of latent trees. We also present regularized versions of our algorithms that learn latent tree approximations of arbitrary distributions. We compare the proposed algorithms to other methods by performing extensive numerical experiments on various latent tree graphical models such as hidden Markov models and star graphs. In addition, we demonstrate the applicability of our methods on real-world datasets by modeling the dependency structure of monthly stock returns in the S&P index and of the words in the 20 newsgroups dataset.

Learning Latent Tree Graphical Models Learning Latent Tree Graphical Models arXiv:1009.2722v1 [stat.ML] 14 Sep 2010 Myung Jin Choi† Vincent Y. F. Tan† Animashree Anandkumar‡ Alan S. Willsky† myungjin@mit.edu vtan@mit.edu a.anandkumar@uci.edu willsky@mit.edu † Stochastic Systems Group, Laboratory for Information and Decision Systems, Massachusetts Institute of Technology. ‡ Center for Pervasive Communications and Computing, Electrical Engineering and Computer Science, University of California, Irvine. Editor: Abstract We study the problem of learning a latent tree graphical model where samples are available only from a subset of variables. We propose two consistent and computationally efficient algorithms for learning minimal latent trees, that is, trees without any redundant hidden nodes. Unlike many existing methods, the observed nodes (or variables) are not constrained to be leaf nodes. Our first algorithm, recursive grouping, builds the latent tree recursively by identifying sibling groups using so-called information distances. One of the main contributions of this work is our second algorithm, which we refer to as CLGrouping. CLGrouping starts with a pre-processing procedure in which a tree over the observed variables is constructed. This global step groups the observed nodes that are likely to be close to each other in the true latent tree, thereby guiding subsequent recursive grouping (or equivalent procedures) on much smaller subsets of variables. This results in more accurate and efficient learning of latent trees. We also present regularized versions of our algorithms that learn latent tree approximations of arbitrary distributions. We compare the proposed algorithms to other methods by performing extensive numerical experiments on various latent tree graphical models such as hidden Markov models and star graphs. In addition, we demonstrate the applicability of our methods on real-world datasets by modeling the dependency structure of monthly stock returns in the S&P index and of the words in the 20 newsgroups dataset. Keywords: Graphical Models, Hidden Variables, Latent Tree Models, Structure Learning 1. Introduction The inclusion of latent variables in modeling complex phenomena and data is a wellrecognized and a valuable construct in a variety of applications, including bioinformatics and computer vision, and the investigation of machine-learning methods for models with latent variables is a substantial and continuing direction of research. There are three challenging problems in learning a model with latent variables: learning the number of latent variables; inferring the structure of how these latent variables relate to each other and to the observed variables; and estimating the parameters characterizing 1 Choi, Tan, Anandkumar, and Willsky those relationships. Issues that one must consider in developing a new learning algorithm include developing tractable methods; incorporating the tradeoff between the fidelity to the given data and generalizability; deriving theoretical results on the performance of such algorithms; and studying applications that provide clear motivation and contexts for the models so learned. One class of models that has received considerable attention in the literature is the class of latent tree models, i.e., graphical models Markov on trees, in which variables at some nodes represent the original (observed) variables of interest while others represent the latent variables. The appeal of such models for computational tractability is clear: with a tree-structured model describing the statistical relationships, inference - processing noisy observations of some or all of the original variables to compute the estimates of all variables - is straightforward and scalable. Although the class of tree-structured models, with or without latent variables, is a constrained one, there are interesting applications that provide strong motivation for the work presented here. In particular, a very active avenue of research in computer vision is the use of context - e.g., the nature of a scene to aid the reliable recognition of objects (and at the same time to allow the recognition of particular objects to assist in recognizing the scene). For example, if one knows that an image is that of an office, then one might expect to find a desk, a monitor on that desk, and perhaps a computer mouse. Hence if one builds a model with a latent variable representing that context (“office”) and uses simple, noisy detectors for different object types, one would expect that the detection of a desk would support the likelihood that one is looking at an office and through that enhance the reliability of detecting smaller objects (monitors, keyboards, mice, etc.). Work along these lines, including by some of the authors of this paper (Parikh and Chen, 2007; Choi et al., 2010), show the promise of using tree-based models of context. This paper considers the problem of learning tree-structured latent models. If all variables are observed in the tree under consideration, then the well-known algorithm of Chow and Liu (1968) provides a tractable algorithm for performing maximum likelihood (ML) estimation of the tree structure. However, ML estimation of latent tree models is NP-hard (Roch, 2006). This has motivated a number of investigations of other tractable methods for learning such trees as well as theoretical guarantees on performance. Our work represents a contribution to this area of investigation. There are three main contributions in our paper. Firstly, by adopting a statistical distance-based framework, we develop two new algorithms, recursive grouping and CLGrouping, for the learning of latent trees, which applies equally well to discrete and Gaussian models. Secondly, we provide consistency guarantees (both structural and parametric) as well as very favorable computational and sample complexity characterizations for both of our algorithms. Thirdly, through extensive numerical experiments on both synthetic and real-world data, we demonstrate the superiority of our approach for a wide variety of models ranging from ones with very large tree diameters (e.g., hidden Markov models (HMMs)) to star models and complete trees.1 Our first algorithm, which we refer to as recursive grouping, constructs a latent tree in a bottom-up fashion, grouping nodes into sibling groups that share the same parent 1. A complete k-ary tree (or k-complete tree) is one in which all leaf nodes are at same depth and all internal nodes have degree k. 2 Learning Latent Tree Graphical Models node, recursively at each level of the resulting hierarchy (and allowing for some of the observed variables to play roles at arbitrary levels in the resulting hierarchy). Our second algorithm, CLGrouping first implements a global construction step, namely producing the Chow-Liu tree for the observed variables without any hidden nodes. This global step then provides guidance for groups of observed nodes that are likely to be topologically close to each other in the latent tree, thereby guiding subsequent recursive grouping or neighborjoining (Saitou and Nei, 1987) computations. Each of these algorithms is consistent and has excellent sample and computational complexity.2 As Pearl (1988) points out, the identification of latent tree models has some built-in ambiguity, as there is an entire equivalence class of models in the sense that when all latent variables are marginalized out, each model in this class yields the same joint distribution over the observed variables. For example, we can take any such latent model and add another hidden variable as a leaf node connected to only one other (hidden or observed) node. Hence, much as one finds in fields such as state space dynamic systems (e.g., Luenberger (1979, Section 8)), there is a notion of minimality that is required here, and our results are stated in terms of consistent learning of such minimal latent models. 1.1 Related Work The relevant literature on learning latent models is vast and in this section, we summarize the main lines of research in this area. The classical latent cluster models (LCM) consider multivariate distributions in which there exists only one latent variable and each state of the variable corresponds to a cluster in the data (Lazarsfeld and Henry, 1968). Hierarchical latent class (HLC) models (Zhang and Kočka, 2004; Zhang, 2004; Chen et al., 2008) generalize these models by allowing multiple latent variables. HLC allows latent variables to have different number of states, but assume that all observed nodes are at the leaves of the tree. Their learning algorithm is based on a greedy approach of making one local move at a time (e.g., introducing one hidden node, or replacing an edge), which is computationally expensive and does not have consistency guarantees. Another greedy learning algorithm called BIN (Harmeling and Williams, 2010) is computationally more efficient, but enforces that each internal node is hidden and has three neighboring nodes. In contrast, we fix the number of states in each hidden node, but allow observed nodes to be internal nodes. Our algorithms are guaranteed to recover the correct structure when certain (mild) conditions are met. Many authors also propose reconstructing latent trees using the expectation maximization (EM) algorithm (Elidan and Friedman, 2005; Kemp and Tenenbaum, 2008). However, as with all other EM-based methods, these approaches depend on the initialization and suffer from the possibility of being trapped in local optima and thus no consistency guarantees can be provided. At each iteration, a large number of candidate structures need to be evaluated, these methods assume that all observed nodes are the leaves of the tree and to reduce the number of candidate structures. Algorithms have been proposed (Hsu et al., 2. As we will see, depending on the true latent tree model, one or the other of these may be more efficient. Roughly speaking, for smaller diameter graphs (such as the star), recursive grouping is faster, and for larger diameter graphs (such as an HMM), CLgrouping is more efficient. 3 Choi, Tan, Anandkumar, and Willsky 2009) with sample complexity guarantees for learning HMMs under the condition that the joint distribution of the observed variables generated by distinct hidden states are distinct. The reconstruction of latent trees has been studied extensively by the phylogenetic community where sequences of extant species are available and the unknown phylogenetic tree is to be inferred from these sequences. See Durbin et al. (1999) for a thorough overview. Efficient algorithms with provable performance guarantees are available (Erdős et al., 1999; Daskalakis et al., 2006). However, the works in this area mostly assume that only the leaves are observed and each internal node (which is hidden) has the same degree except for the root. The most popular algorithm for constructing phylogenetic trees is the neighbor-joining (NJ) method by Saitou and Nei (1987). Like recursive grouping, the input to the algorithm is a set of statistical distances between observed variables. The algorithm proceeds by recursively pairing two nodes that are the closest neighbors in the true latent tree and introducing a hidden node as the parent of the two nodes. For more details on NJ, the reader is referred to Durbin et al. (1999, Section 7.3). Another popular class of reconstruction methods used in the phylogenetic community is the family of quartet-based distance methods (Bandelth and Dress, 1986; Erdős et al., 1999; Jiang et al., 2001).3 Quartet-based methods first construct a set of quartets for all subsets of four observed nodes. Subsequently, these quartets are then combined to form a latent tree. However, when we only have access to the samples at the observed nodes, then it is not straightforward to construct a latent tree from a set of quartets since the quartets may be not be consistent.4 In fact, it is known that the problem of determining a latent tree that agrees with the maximum number of quartets is NP-hard (Steel, 1992) but many heuristics have been proposed (Farris, 1972; Sattath and Tversky, 1977). Also, when only samples are available, quartet-based methods are usually much less accurate than NJ (St. John et al., 2003) so we only compare our proposed algorithms to NJ. For further comparisons between (the sample complexity and other aspects of) quartet methods and NJ, the reader is referred to Csűrös (2000) and St. John et al. (2003). Another distance-based algorithm was proposed in Pearl (1988, Section 8.3.3). This algorithm is very similar in spirit to quartet-based methods but instead of finding quartets for all subsets of four observed nodes, it finds just enough quartets to determine the location of each observed node in the tree. Although the algorithm is consistent, it performs poorly when only the samples of observed nodes are available (Pearl, 1988, Section 8.3.5). The learning of phylogenetic trees is related to the emerging field of network tomography (Castro et al., 2004) in which one seeks to learn characteristics (such as structure) from data which are only available at the end points (e.g., sources and sinks) of the network. However, again observations are only available at the leaf nodes and usually the objective is to estimate the delay distributions corresponding to nodes linked by an edge (Tsang et al., 2003; Bhamidi et al., 2009). The modeling of the delay distributions is different from the learning of latent tree graphical models discussed in this paper. 3. A quartet is simply an unrooted binary tree on a set of four observed nodes. 4. The term consistent here is not the same as the estimation-theoretic one. Here, we say that a set of quartets is consistent if there exists a latent tree such that all quartets agree with the tree. 4 Learning Latent Tree Graphical Models 1.2 Paper Organization The rest of the paper is organized as follows. In Section 2, we introduce the notations and terminologies used in the paper. In Section 3, we introduce the notion of information distances which are used to reconstruct tree models. In the following two sections, we make two assumptions: Firstly, the true distribution is a latent tree and secondly, perfect knowledge of information distance of observed variables is available. We introduce recursive grouping in Section 4. This is followed by our second algorithm CLGrouping in Section 5. In Section 6, we relax the assumption that the information distances are known and develop sample based algorithms and at the same time provide sample complexity guarantees for recursive grouping and CLGrouping. We also discuss extensions of our algorithms for the case when the underlying model is not a tree and our goal is to learn an approximation to it using a latent tree model. We demonstrate the empirical performance of our algorithms in Section 7 and conclude the paper in Section 8. The Appendix includes proofs for the theorems presented in the paper. 2. Latent Tree Graphical Models 2.1 Mathematical Notation Let G = (W, E) be an undirected graph with vertex (or node) set W = {1, . . . , M } and edge set E ⊂ W 2 . Let nbd(i; G) and nbd[i; G] be the set of neighbors of node i and the closed neighborhood of i respectively, i.e., nbd[i; G] := nbd(i; G) ∪ {i}. For a tree T = (W, E), the set of leaf nodes (nodes with degree 1), the maximum degree, and the diameter are denoted by Leaf(T ), ∆(T ), and diam(T ) respectively. The path between two nodes i and j in a tree T = (W, E) is the set of edges connecting i and j and is denoted as Path((i, j); E). The depth of a node in a tree is the shortest distance (number of hops) to the leaf nodes in T . The parent of a node i that has depth s is the neighbor of i with depth s + 1. The set of child nodes of a node i with depth s is the set of neighbors of i with depth s − 1, which we denote as C(i). A set of nodes that share the same parent is called a sibling group. A family is the union of the siblings and the associated parent. A latent tree is a tree with node set W := V ∪ H, the union of a set of observed nodes V (with m = |V |), and a set of latent (or hidden) nodes H. The effective depth δ(T ; V ) (with respect to V ) is the maximum distance of a hidden node to its closest observed node, i.e., δ(T ; V ) := max min |Path((i, j); T )|. i∈H j∈V (1) 2.2 Graphical Models An undirected graphical model (Lauritzen, 1996) is a family of multivariate probability distributions that factorize according to a graph G = (W, E). More precisely, let X = (X1 , . . . , XM ) be a random vector, where each random variable Xi , which takes on values in an alphabet X , corresponds to variable at node i ∈ V . The set of edges E encodes the set of conditional independencies in the model. The random vector X is said to be Markov on G if for every i, the random variable Xi is conditionally independent of all other variables 5 Choi, Tan, Anandkumar, and Willsky given its neighbors, i.e, if p is the joint distribution5 of X, then p(xi |xnbd(i;G) ) = p(xi |x\i ), (2) where x\i denotes the set of all variables6 excluding xi . Eqn. (2) is known as the local Markov property. In this paper, we consider both discrete and Gaussian graphical models. For discrete models, the alphabet X = {1, . . . , K} is a finite set. For Gaussian graphical models, X = R and furthermore, without loss of generality, we assume that the mean is known to be the zero vector and hence, the joint distribution p(x) ∝ exp(− 12 xT Σ−1 x) depends only on the covariance matrix Σ. An important and tractable class of graphical models is the set of tree-structured graphical models, i.e., multivariate probability distributions that are Markov on an undirected tree T = (W, E). It is known from junction tree theory (Cowell et al., 1999) that the joint distribution p for such a model factorizes as p(x1 , . . . , xM ) = Y p(xi ) i∈W Y (i,j)∈E p(xi , xj ) . p(xi )p(xj ) (3) That is, the sets of marginal {p(xi ) : i ∈ W } and pairwise joints on the edges {p(xi , xj ) : (i, j) ∈ E} fully characterize the joint distribution of a tree-structured graphical model. A special class of a discrete tree-structured graphical models is the set of symmetric discrete distributions. This class of models is characterized by the fact that the pairs of variables (Xi , Xj ) on all the edges (i, j) ∈ E follow the conditional probability law: 1 − (K − 1)θij , if xi = xj , p(xi |xj ) = (4) θij , otherwise, and the marginal distribution of every variable in the tree is uniform, i.e., p(xi ) = 1/K for all xi ∈ X and for all i ∈ V ∪ H. The parameter θij ∈ (0, 1/K) in (4), which does not depend on the states xi , xj ∈ X , is known as the crossover probability. Let xn := {x(1) , . . . , x(n) } be a set of n i.i.d. samples drawn from a graphical model (distribution) p, Markov on a latent tree Tp = (W, Ep ), where W = V ∪ H. Each sample x(l) ∈ X M is a length-M vector. In our setup, the learner only has access to samples drawn from the observed node set V , and we denote this set of sub-vectors containing (l) (n) (1) only the elements in V , as xnV := {xV , . . . , xV }, where each observed sample xV ∈ X m is a length-m vector. Our algorithms learn latent tree structures using the information distances (defined in Section 3) between pairs of observed variables, which can be estimated from samples. 2.3 Minimal Tree Extensions Our ultimate goal is to recover the graphical model p, i.e., the latent tree structure and its parameters, given n i.i.d. samples of the observed variables xnV . However, in general, there 5. We abuse the term distribution to mean a probability mass function in the discrete case (density with respect to the counting measure) and a probability density function (density with respect to the Lebesgue measure) in the continuous case. 6. We will use the terms node, vertex and variable interchangeably in the sequel. 6 Learning Latent Tree Graphical Models can be multiple latent tree models which result in the same observed statistics, i.e., the same joint distribution pV of the observed variables. We consider the class of tree models where it is possible to recover the latent tree model uniquely and provide necessary conditions for structure identifiability, i.e., the identifiability of the edge set E. Firstly, we limit ourselves to the scenario where all the random variables (both observed and latent) take values on a common alphabet X . Thus, in the Gaussian case, each hidden and observed variable is a univariate Gaussian. In the discrete case, each variable takes on values in the same finite alphabet X . Note that the model may not be identifiable if some of the hidden variables are allowed to have arbitrary alphabets. As an example, consider a discrete latent tree model with binary observed variables (K = 2). A latent tree with the simplest structure (fewest number of nodes) is a tree in which all m observed binary variables are connected to one hidden variable. If we allow the hidden variable to take on 2m states, then the tree can describe all possible statistics among the m observed variables, i.e., the joint distribution pV can be arbitrary.7 A probability distribution pV (xV ) is said to be tree-decomposable if it is the marginal (of variables in V ) of a tree-structured graphical model p(xV , xH ). In this case, p (over variables in W ) is said to be a tree extension of pV (Pearl, 1988). A distribution p is said to have a redundant hidden node h ∈ H if we can remove h and the marginal on the set of visible nodes V remains as pV . The following conditions ensure that a latent tree does not include a redundant hidden node (Pearl, 1988): (C1) Each hidden variable has at least three neighbors (which can be either hidden or observed). Note that this ensures that all leaf nodes are observed (although not all observed nodes need to be leaves). (C2) Any two variables connected by an edge in the tree model are neither perfectly dependent nor independent. Figure 1(a) shows an example of a tree satisfying (C1). If (C2), which is a condition on parameters, is also satisfied, then the tree in Figure 1(a) is identifiable. The tree shown in Figure 1(b) does not satisfy (C1) because h4 and h5 have degrees less than 3. In fact, if we marginalize out the hidden variables h4 and h5 , then the resulting model has the same tree structure as in Figure 1(a). We assume throughout the paper that (C2) is satisfied for all probability distributions. Let T≥3 be the set of (latent) trees satisfying (C1). We refer to T≥3 as the set of minimal (or identifiable) latent trees. Minimal latent trees do not contain redundant hidden nodes. The distribution p (over W and Markov on some tree in T≥3 ) is said to be a minimal tree extension of pV . As illustrated in Figure 1, using marginalization operations, any non-minimal latent tree distribution can be reduced to a minimal latent tree model. Proposition 1 (Minimal Tree Extensions) (Pearl, 1988, Section 8.3) (i) For every tree-decomposable distribution pV , there exists a minimal tree extension p Markov on a tree T ∈ T≥3 , which is unique up to the renaming of the variables or their values. 7. This follows from a elementary parameter counting argument. 7 Choi, Tan, Anandkumar, and Willsky h1 h1 h2 1 2 3 4 h2 1 h3 5 2 3 6 6 h4 h3 4 5 h5 (b) (a) Figure 1: Examples of minimal latent trees. Shaded nodes are observed and unshaded nodes are hidden. (a) An identifiable tree. (b) A non-identifiable tree because h4 and h5 have degrees less than 3. (ii) For Gaussian and binary distributions, if pV is known exactly, then the minimal tree extension p can be recovered. (iii) The structure of T is uniquely determined by the pairwise distributions of observed variables p(xi , xj ) for all i, j ∈ V . 2.4 Consistency We now define the notion of consistency. In Section 6, we show that our latent tree learning algorithms are consistent. Definition 2 (Consistency) A latent tree reconstruction algorithm A is a map from the observed samples xnV to an estimated tree Tbn and an estimated tree-structured graphical model pbn . We say that a latent tree reconstruction algorithm A is structurally consistent if there exists a graph homomorphism8 h such that lim Pr(h(Tbn ) 6= Tp ) = 0. (5) lim Pr (D(p || pbn ) > ε) = 0, (6) n→∞ Furthermore, we say that A is risk consistent if to every ε > 0, n→∞ where D(p || pbn ) is the KL-divergence (Cover and Thomas, 2006) between the true distribution p and the estimated distribution pbn . In the following sections, we design structurally and risk consistent algorithms for (minimal) Gaussian and symmetric discrete latent tree models, defined in (4). Our algorithms use pairwise distributions between the observed nodes. However, for general discrete models, pairwise distributions between observed nodes are, in general, not sufficient to recover 8. A graph homomorphism is a mapping between graphs that respects their structure. More precisely, a graph homomorphism h from a graph G = (W, E) to a graph G′ = (V ′ , E ′ ), written h : G → G′ is a mapping h : V → V ′ such that (i, j) ∈ E implies that (h(i), h(j)) ∈ E ′ . 8 Learning Latent Tree Graphical Models the parameters (Chang and Hartigan, 1991). Therefore, we only prove structural consistency, as defined in (5), for general discrete latent tree models. For such distributions, we consider a two-step procedure for structure and parameter estimation: Firstly, we estimate the structure of the latent tree using the algorithms suggested in this paper. Subsequently, we use the Expectation Maximization (EM) algorithm (Dempster et al., 1977) to infer the parameters. Note that, as mentioned, risk consistency will not be guaranteed in this case. 3. Information Distances The proposed algorithms in this paper receive as inputs the set of so-called (exact or estimated) information distances, which are functions of the pairwise distributions. These quantities are defined in Section 3.1 for the two classes of tree-structured graphical models discussed in this paper, namely the Gaussian and discrete graphical models. We also show that the information distances have a particularly simple form for symmetric discrete distributions. In Section 3.2, we use the information distances to infer the relationships between the observed variables such as j is a child of i or i and j are siblings. 3.1 Definitions of Information Distances We define information distances for Gaussian and discrete distributions and show that these distances are additive for tree-structured graphical models. Recall that for two random variables Xi and Xj , the correlation coefficient is defined as ρij := p Cov(Xi , Xj ) . Var(Xi )Var(Xj ) (7) For Gaussian graphical models, the information distance associated with the pair of variables Xi and Xj is defined as: dij := − log |ρij |. (8) Intuitively, if the information distance dij is large, then Xi and Xj are weakly correlated and vice-versa. For discrete graphical models, let Jij denote the joint probability matrix between Xi and ij Xj (i.e., Jab = p(xi = a, xj = b), a, b ∈ X ). Also let Mi be the diagonal marginal probability i = p(x = a)). For discrete graphical models, the information distance matrix of Xi (i.e., Maa i associated with the pair of variables Xi and Xj is defined as (Lake, 1994): dij := − log √ | det Jij | det Mi det Mj . (9) Note that for binary variables, i.e., K = 2, the value of dij in (9) reduces to the expression in (8), i.e., the information distance is a function of the correlation coefficient, defined in (7), just as in the Gaussian case. For symmetric discrete distributions defined in (4), the information distance defined for discrete graphical models in (9) reduces to dij := −(K − 1) log(1 − Kθij ). 9 (10) Choi, Tan, Anandkumar, and Willsky Note that there is one-to-one correspondence between the information distances dij and model parameters for Gaussian distributions (parametrized by the correlation coefficient ρij ) in (8) and the symmetric discrete distributions (parametrized by the crossover probability θij ) in (10). This is, however, not true for general discrete distributions. Equipped with these definitions of information distances, assumption (C2) in Section 2.3 can be rewritten as the following: There exists constants 0 < l, u < ∞, such that (C2) l ≤ dij ≤ u, ∀ (i, j) ∈ Ep . (11) Proposition 3 (Additivity of Information Distances) The information distances dij defined in (8), (9), and (10) are additive tree metrics (Erdős et al., 1999). In other words, if the joint probability distribution p(x) is a tree-structured graphical model Markov on the tree Tp = (W, Ep ), then the information distances are additive on Tp : X dkl = dij , ∀k, l ∈ W. (12) (i,j)∈Path((k,l);Ep ) The property in (12) means that if each pair of vertices i, j ∈ W is assigned the weight dij , then Tp is a minimum spanning tree on W , denoted as MST(W ; D), where D is the information distance matrix with elements dij for all i, j ∈ V . It is straightforward to show that the information distances are additive for the Gaussian and symmetric discrete cases using the local Markov property of graphical models. For general discrete distributions with information distance as in (9), see Lake (1994) for the proof. In the rest of the paper, we map the parameters of Gaussian and discrete distributions to an information distance matrix D = [dij ] to unify the analyses for both cases. 3.2 Testing Inter-Node Relationships In this section, we use Proposition 3 to ascertain child-parent and sibling (cf. Section 2.1) relationships between the variables in a latent tree-structured graphical model. To do so, for any three variables i, j, k ∈ V , we define Φijk := dik − djk to be the difference between the information distances dik and djk . The following lemma suggests a simple procedure to identify the set of relationships between the nodes. Lemma 4 (Sibling Grouping) For distances dij for all i, j ∈ V on a tree T ∈ T≥3 , the following two properties on Φijk = dik − djk hold: (i) Φijk = dij for all k ∈ V \ {i, j} if and only if i is a leaf node and j is its parent. (ii) −dij < Φijk = Φijk′ < dij for all k, k′ ∈ V \ {i, j} if and only if both i and j are leaf nodes and they have the same parent, i.e., they belong to the same sibling group. The proof of the lemma uses Proposition 3 and is provided in Appendix A.1. Given Lemma 4, we can first determine all the values of Φijk for triples i, j, k ∈ V . Now we can determine the relationship between nodes i and j as follows: Fix the pair of nodes i, j ∈ V and consider all the other nodes k ∈ V \ {i, j}. Then, there are three possibilities for the set {Φijk : k ∈ V \ {i, j}}: 10 Learning Latent Tree Graphical Models e1 e1 e2 e3 e1 e2 e3 e1 e 2 e3 e4 e5 e6 e7 e4 e5 i e8 e8 e6 e7 i j e4 e 5 i e6 e7 k j e4 e5 i e8 j e1 e2 e3 k` k` k` e6 e7 (b) (c) e4 e5 j j e8 e8 e6 e7 k k (a) e2 e3 i (d) (e) Figure 2: Examples for each case in TestNodeRelationships. For each edge, ei represents the information distance associated with the edge. (a) Case 1: Φijk = −e8 = −dij for all k ∈ V \ {i, j}. (b) Case 2: Φijk = e6 − e7 6= dij = e6 + e7 for all k ∈ V \ {i, j} (c) Case 3a: Φijk = e4 + e2 + e3 − e7 6= Φijk′ = e4 − e2 − e3 − e7 . (d) Case 3b: Φijk = e4 + e5 6= Φijk′ = e4 − e5 . (e) Case 3c: Φijk = e5 6= Φijk′ = −e5 . 1. Φijk = dij for all k ∈ V \ {i, j}. Then, i is a leaf node and j is a parent of i. Similarly, if Φijk = −dij for all k ∈ V \ {i, j}, j is a leaf node and i is a parent of j. 2. Φijk is constant for all k ∈ V \ {i, j} but not equal to either dij or −dij . Then i and j are leaf nodes and they are siblings. 3. Φijk is not equal for all k ∈ V \ {i, j}. Then, there are three possibilities: Either (a) Nodes i and j are not siblings nor have a parent-child relationship or, (b) Nodes i and j are siblings but at least one of them is not a leaf or, (c) Nodes i and j have a parent-child relationship but the child is not a leaf. Thus, we have a simple test to determine the relationship between i and j and to ascertain whether i and j are leaf nodes. We call the above test TestNodeRelationships. See Figure 2 for examples. By running this test for all i and j, we can determine all the relationships among all pairs of observed variables. In the following section, we describe a recursive algorithm that is based on the above TestNodeRelationships procedure to reconstruct the entire latent tree model assuming that the true distance matrix D = [dij ] are known. In Section 5, we provide improved algorithms for the learning of latent trees again assuming that D is known. Subsequently, in Section 6, we develop algorithms for the consistent reconstruction of latent trees when information distances are unknown and we have to estimate them from the samples xnV . In addition, in Section 6.4 we discuss how to extend these algorithms for the case when pV is not necessarily tree-decomposable, i.e., the original graphical model is not assumed to be a latent tree. 4. Recursive Grouping Algorithm Given Information Distances This section is devoted to the development of the first algorithm for reconstructing latent tree models, recursive grouping (RG). At a high level, RG is a recursive procedure in which at each step, TestNodeRelationships is used to identify nodes that belong to the same family. Subsequently, RG introduces a parent node if a family of nodes (i.e., a sibling group) 11 Choi, Tan, Anandkumar, and Willsky does not contain an observed parent. This newly introduced parent node corresponds to a hidden node in the original unknown latent tree. Once such a parent (i.e., hidden) node h is introduced, the information distances from h to all other observed nodes can be computed. The inputs to RG are the vertex set V and the matrix of information distances D corresponding to a latent tree. The algorithm proceeds by recursively grouping nodes and adding hidden variables. In each iteration, the algorithm acts on a so-called active set of nodes Y , and in the process constructs a new active set Ynew for the next iteration.9 The steps are as follows: 1. Initialize by setting Y := V to be the set of observed variables. 2. Compute Φijk = dik − djk for all i, j, k ∈ Y . 3. Using the TestNodeRelationships procedure, define {Πl }L l=1 to be the coarsest partition10 of Y such that for every subset Πl (with |Πl | ≥ 2), any two nodes in Πl are either siblings which are leaf nodes or they have a parent-child relationship in which the child is a leaf. Note that for some l, Πl may S consist of a single node. Begin to construct the new active set by setting Ynew ← l:|Πl |=1 Πl . 4. For each l = 1, . . . , L with |Πl | ≥ 2, if Πl contains a parent node u, update Ynew ← Ynew ∪ {u}. Otherwise, introduce a new hidden node h, connect h (as a parent) to every node in Πl , and set Ynew ← Ynew ∪ {h}. 5. Update the active set: Yold ← Y and Y ← Ynew . 6. For each new hidden node h ∈ Y , compute the information distances dhl for all l ∈ Y using (13) and (14) described below. 7. If |Y | ≥ 3, return to step 2. Otherwise, if |Y | = 2, connect the two remaining nodes in Y with an edge then stop. If instead |Y | = 1, do nothing and stop. We now describe how to compute the information distances in Step 6 for each new hidden node h ∈ Y and all other active nodes l ∈ Y . Let i, j ∈ C(h) be two children of h, and let k ∈ Yold \ {i, j} be any other node in the previous active set. From Proposition 3, we have that dih − djh = dik − djk = Φijk and dih + djh = dij , from which we can recover the information distances between a previously active node i ∈ Yold and its new hidden parent h ∈ Y as follows: 1 (13) dih = (dij + Φijk ) . 2 For any other active node l ∈ Y , we can compute dhl using a child node i ∈ C(h) as follows: dil − dih , if l ∈ Yold , dhl = (14) dik − dih − dlk , otherwise, where k ∈ C(l). 9. Note that the current active set is also used (in Step 6) after the new active set has been defined. For clarity, we also introduce the quantity Yold in Steps 5 and 6. L 10. Recall that a partition P of a set Y is a collection of nonempty subsets {Πl ⊂ Y }L l=1 such that ∪l=1 Πl = Y and Πl ∩ Πl′ = ∅ for all l 6= l′ . A partition P is said to be coarser than another partition P ′ if every element of P ′ is a subset of some element of P . 12 Learning Latent Tree Graphical Models hh2 1 hh2 h3 h1 2 4 5 (a) 3 6 1 h1 2 4 5 1 3 h1 2 4 6 5 (c) (b) 3 6 h3 h2 h3 1 h1 2 4 5 3 6 (d) Figure 3: An illustrative example of RG. Solid nodes indicate the active set Y for each iteration. (a) Original latent tree. (b) Output after the first iteration of RG. Red dotted lines indicate the subsets Πl in the partition of Y . (c) Output after the second iteration of RG. Note that h3 , which was introduced in the first iteration, is an active node for the second iteration. Nodes 4,5, and 6 do not belong to the current active set and are represented in grey. (d) Output after the third iteration of RG, which is same as the original latent tree. Using equations (13) and (14), we can infer all the information distances dhl between a newly introduced hidden node h to all other active nodes l ∈ Y . Consequently, we have all the distances dij between all pairs of nodes in the active set Y . It can be shown that this algorithm recovers all minimal latent trees. The proof of the following theorem is provided in Appendix A.2. Theorem 5 (Correctness and Computational Complexity of RG) If Tp ∈ T≥3 and the matrix of information distances D (between nodes in V ) is available, then RG outputs the true latent tree Tp correctly in time O(diam(Tp )m3 ). We now use a concrete example to illustrate the steps involved in RG. In Figure 3(a), the original unknown latent tree is shown. In this tree, nodes 1, . . . , 6 are the observed nodes and h1 , h2 , h3 are the hidden nodes. We start by considering the set of observed nodes as active nodes Y := V = {1, . . . , 6}. Once Φijk are computed from the given distances dij , TestNodeRelationships is used to determine that Y is partitioned into four subsets: Π1 = {1}, Π2 = {2, 4}, Π3 = {5, 6}, Π4 = {3}. The subsets Π1 and Π4 contain only one node. The subset Π3 contains two siblings that are leaf nodes. The subset Π2 contains a parent node 2 and a child node 4, which is a leaf node. Since Π3 does not contain a parent, we introduce a new hidden node h1 and connect h1 to 5 and 6 as shown in Figure 3(b). The information distances d5h1 and d6h1 can be computed using (13), e.g., d5h1 = 21 (d56 + Φ561 ). The new active set is the union of all nodes in the single-node subsets, a parent node, and a new hidden node Ynew = {1, 2, 3, h1 }. Distances among the pairs of nodes in Ynew can be computed using (14) (e.g., d1h1 = d15 − d5h1 ). In the second iteration, we again use TestNodeRelationships to ascertain that Y can be partitioned into Π1 = {1, 2} and Π2 = {h1 , 3}. These two subsets do not have parents so h2 and h3 are added to Π1 and Π2 respectively. Parent nodes h2 and h3 are connected to their children in Π1 and Π2 as shown in Figure 3(c). Finally, we are left with the active set as Y = {h2 , h3 } and the algorithm terminates after h2 and h3 are connected by an edge. The hitherto unknown latent tree is fully reconstructed as shown in Figure 3(d). 13 Choi, Tan, Anandkumar, and Willsky A potential drawback of RG is that it involves multiple local operations, which may result in a high computational complexity. Indeed, from Theorem 5, the worst-case complexity is O(m4 ) which occurs when Tp , the true latent tree, is a hidden Markov model (HMM). This may be computationally prohibitive if m is large. In Section 5 we design an algorithm which uses a global pre-processing step to reduce the overall complexity substantially, especially for trees with large diameters (of which HMMs are extreme examples). 5. CLGrouping Algorithm Given Information Distances In this section, we present CLGrouping, an algorithm for reconstructing latent trees more efficiently than RG. As in Section 4, in this section, we assume that D is known; the extension to unknown D is discussed in Section 6.3. CLGrouping is a two-step procedure, the first of which is a global pre-processing step that involves the construction of a so-called Chow-Liu tree (Chow and Liu, 1968) over the set of observed nodes V . This step identifies nodes that do not belong to the same sibling group. In the second step, we complete the recovery of the latent tree by applying a distance-based latent tree reconstruction algorithm (such as RG or NJ) repeatedly on smaller subsets of nodes. We review the Chow-Liu algorithm in Section 5.1, relate the Chow-Liu tree to the true latent tree in Section 5.2, derive a simple transformation of the Chow-Liu tree to obtain the latent tree in Section 5.3 and propose CLGrouping in Section 5.4. For simplicity, we focus on the Gaussian distributions and the symmetric discrete distributions first, and discuss the extension to general discrete models in Section 5.5. 5.1 A Review of the Chow-Liu Algorithm In this section, we review the Chow-Liu tree reconstruction procedure. To do so, define T (V ) to be the set of trees with vertex set V and P(T (V )) to be the set of tree-structured graphical models whose graph has vertex set V , i.e., every q ∈ P(T (V )) factorizes as in (3). Given an arbitrary multivariate distribution pV (xV ), Chow and Liu (1968) considered the following KL-divergence minimization problem: pCL := argmin D(pV || q). (15) q∈P(T (V )) That is, among all the tree-structured graphical models with vertex set V , the distribution pCL is the closest one to pV in terms of the KL-divergence. By using the factorization property in (3), we can easily verify that pCL is Markov on the Chow-Liu tree TCL = (V, ECL ) which is given by the optimization problem:11 X I(Xi ; Xj ). (16) TCL = argmax T ∈T (V ) (i,j)∈T In (16), I(Xi ; Xj ) = D(p(xi , xj ) || p(xi ) p(xj )) is the mutual information (Cover and Thomas, 2006) between random variables Xi and Xj . The optimization in (16) is a max-weight spanning tree problem (Cormen et al., 2003) which can be solved efficiently in time O(m2 log m) 11. In (16) and the rest of the paper, we adopt the following simplifying notation; If T = (V, E) and if (i, j) ∈ E, we will also say that (i, j) ∈ T . 14 Learning Latent Tree Graphical Models using either Kruskal’s algorithm (Kruskal, 1956) or Prim’s algorithm (Prim, 1957). The edge weights for the max-weight spanning tree are precisely the mutual information quantities between random variables. Note that the parameters of pb in (15) are found by setting the pairwise distributions pCL (xi , xj ) on the edges to pV (xi , xj ), i.e., pCL (xi , xj ) = pV (xi , xj ) for all (i, j) ∈ ECL . We now relate the Chow-Liu tree on the observed nodes and the information distance matrix D. Lemma 6 (Correspondence between TCL and MST) If pV is a Gaussian distribution or a symmetric discrete distribution, then the Chow-Liu tree in (16) reduces to the minimum spanning tree (MST) where the edge weights are the information distances dij , i.e., X dij . (17) TCL = MST(V ; D) := argmin T ∈T (V ) (i,j)∈T Lemma 6, whose proof is omitted, follows because for Gaussian and symmetric discrete models, the mutual information12 I(Xi ; Xj ) is a monotonically decreasing function of the information distance dij .13 For other graphical models (e.g., non-symmetric discrete distributions), this relationship is not necessarily true. See Section 5.5 for a discussion. Note that when all nodes are observed (i.e., W = V ), Lemma 6 reduces to Proposition 3. 5.2 Relationship between the Latent Tree and the Chow-Liu Tree (MST) In this section, we relate MST(V ; D) in (17) to the original latent tree Tp . To relate the two trees, MST(V ; D) and Tp , we first introduce the notion of a surrogate node. Definition 7 (Surrogate Node) Given the latent tree Tp = (W, Ep ) and any node i ∈ W , the surrogate node of i with respect to V is defined as Sg(i; Tp , V ) := argmin dij . (18) j∈V Intuitively, the surrogate node of a hidden node h ∈ H is an observed node j ∈ V that is most strongly correlated to h. In other words, the information distance between h and j is the smallest. Note that if i ∈ V , then Sg(i; Tp , V ) = i since dii = 0. The map Sg(i; Tp , V ) is a many-to-one function, i.e., several nodes may have the same surrogate node, and its inverse is the inverse surrogate set of i denoted as Sg−1 (i; Tp , V ) := {h ∈ W : Sg(h; Tp , V ) = i}. (19) When the tree Tp and the observed vertex set V are understood from context, the surrogate node of h and the inverse surrogate set of i are abbreviated as Sg(h) and Sg−1 (i) respectively. We now relate the original latent tree Tp = (W, Ep ) to the Chow-Liu tree (also termed the MST) MST(V ; D) formed using the distance matrix D. 12. Note that, unlike information distances dij , the mutual information quantities I(Xi ; Xj ) do not form an additive metric on Tp . 13. For example, in the case of Gaussians, I(Xi ; Xj ) = − 12 log(1 − ρ2ij ) (Cover and Thomas, 2006). 15 Choi, Tan, Anandkumar, and Willsky Lemma 8 (Properties of the MST) The MST in (17) and surrogate nodes satisfy the following properties: (i) The surrogate nodes of any two neighboring nodes in Ep are neighbors in the MST, i.e., for all i, j ∈ W with Sg(i) 6= Sg(j), (i, j) ∈ Ep ⇒ (Sg(i), Sg(j)) ∈ MST(V ; D). (20) (ii) If j ∈ V and h ∈ Sg−1 (j), then every node along the path connecting j and h belongs to the inverse surrogate set Sg−1 (j). (iii) The maximum degree of the MST satisfies u ∆(MST(V ; D)) ≤ ∆(Tp )1+ l δ(Tp ;V ) , (21) where δ(Tp ; V ) is the effective depth defined in (1) and l, u are the bounds on the information distances on edges in Tp defined in (11). The proof of this result can be found in Appendix A.3. As a result of Lemma 8, the properties of MST(V ; D) can be expressed in terms of the original latent tree Tp . For example, in Figure 5(a), a latent tree is shown with its corresponding surrogacy relationships, and Figure 5(b) shows the corresponding MST over the observed nodes. The properties in Lemma 8(i-ii) can also be regarded as edge-contraction operations (Robinson and Foulds, 1981) in the original latent tree to obtain the MST. More precisely, an edge-contraction operation on an edge (j, h) ∈ V × H in the latent tree Tp is defined as the “shrinking” of (j, h) to a single node whose label is the observed node j. Thus, the edge (j, h) is “contracted” to a single node j. By using Lemma 8(i-ii), we observe that the Chow-Liu tree MST(V ; D) is formed by applying edge-contraction operations sequentially to each (j, h) pair for all h ∈ Sg−1 (j) ∩ H until all pairs have been contracted to a single node j. For example, the MST in Figure 5(b) is obtained by contracting edges (3, h3 ), (5, h2 ), and then (5, h1 ) in the latent tree in Figure 5(a). The properties in Lemma 8 can be used to design efficient algorithms based on transforming the MST to obtain the latent tree Tp . Note that the maximum degree of the MST, ∆(MST(V ; D)), is bounded by the maximum degree in the original latent tree. The quantity ∆(MST(V ; D)) determines the computational complexity of one of our proposed algorithms (CLGrouping) and it is small if the depth of the latent tree δ(Tp ; V ) is small and the information distances dij satisfy tight bounds (i.e., u/l is close to unity). The latter condition holds for (almost) homogeneous models in which all the information distances dij on the edges are almost equal. 5.3 Chow-Liu Blind Algorithm for a Subclass of Latent Trees In this section, we present a simple and intuitive transformation of the Chow-Liu tree that produces the original latent tree. However, this algorithm, called Chow-Liu Blind (or CLBlind), is applicable only to a subset of latent trees called blind latent tree-structured graphical models P(Tblind ). Equipped with the intuition from CLBlind, we generalize it in Section 5.4 to design the CLGrouping algorithm that produces the correct latent tree structure from the MST for all minimal latent tree models. If p ∈ P(Tblind ), then its structure Tp = (W, Ep ) and the distance matrix D satisfy the following properties: 16 Learning Latent Tree Graphical Models h1 1 h2 2 3 (a) 4 5 2 3 1 4 (b) 5 2 3 1 4 h2 2 5 1 (c) 3 4 h2 2 5 (d) 1 3 4 h1 5 1 h2 2 (e) 3 4 5 (f) Figure 4: An illustration of CLBlind. The shaded nodes are the observed nodes and the rest are hidden nodes. The dotted lines denote surrogate mappings for the hidden nodes. (a) Original latent tree, which belongs to the class of blind latent graphical models, (b) Chow-Liu tree over the observed nodes, (c) Node 3 is the input to the blind transformation, (d) Output after the blind transformation, (e) Node 2 is the input to the blind transformation, (f) Output after the blind transformation, which is same as the original latent tree. (i) The true latent tree Tp ∈ T≥3 and all the internal nodes14 are hidden, i.e., V = Leaf(Tp ). (ii) The surrogate node of (i.e., the observed node with the strongest correlation with) each hidden node is one of its children, i.e., Sg(h) ∈ C(h) for all h ∈ H. We now describe the CLBlind algorithm, which involves two main steps. Firstly, MST(V ; D) is constructed using the distance matrix D. Secondly, we apply the blind transformation of the Chow-Liu tree BlindTransform(MST(V ; D)), which proceeds as follows: 1. Identify the set of internal nodes in MST(V ; D). We perform an operation for each internal node as follows: 2. For internal node i, add a hidden node h to the tree. 3. Connect an edge between h and i (which now becomes a leaf node) and also connect edges between h and the neighbors of i in the current tree model. 4. Repeat steps 2 and 3 until all internal nodes have been operated on. See Figure 4 for an illustration of CLBlind. We use the adjective blind to describe the transformation BlindTransform(MST(V ; D)) since it does not depend on the distance matrix D but uses only the structure of the MST. The following theorem whose proof can be found in Appendix A.4 states the correctness result for CLBlind. Theorem 9 (Correctness and Computational Complexity of CLBlind) If p ∈ P(Tblind ) is a blind tree-structured graphical model Markov on Tp and the matrix of distances D is known, then CLBlind outputs the true latent tree Tp correctly in time O(m2 log m). The first condition on P(Tblind ) that all internal nodes are hidden is not uncommon in applications. For example, in phylogenetics, (DNA or amino acid) sequences of extant species 14. Recall that an internal node is one whose degree is greater than or equal to 2, i.e., a non-leaf. 17 Choi, Tan, Anandkumar, and Willsky at the leaves are observed, while the sequences of the extinct species are hidden (corresponding to the internal nodes), and the evolutionary (phylogenetic) tree is to be reconstructed. However, the second condition is more restrictive15 since it implies that each hidden node is connected to at least one observed node and that it is closer (i.e., more correlated) to one of its observed children compared to any other observed node. If the first constraint is satisfied but not the second, then the blind transformation BlindTransform(MST(V ; D)) does not overestimate the number of hidden variables in the latent tree (the proof follows from Lemma 8 and is omitted). Since the computational complexity of constructing the MST is O(m2 log m) where m = |V |, and the blind transformation is at most linear in m, the overall computational complexity is O(m2 log m). Thus, CLBlind is a computationally efficient procedure compared to RG, described in Section 4. 5.4 Chow-Liu Grouping Algorithm Even though CLBlind is computationally efficient, it only succeeds in recovering latent trees for a restricted subclass of minimal latent trees. In this section, we propose an efficient algorithm, called CLGrouping that reconstructs all minimal latent trees. We also illustrate CLGrouping using an example. CLGrouping uses the properties of the MST as described in Lemma 8. At a high-level, CLGrouping involves two distinct steps: Firstly, we construct the ChowLiu tree MST(V ; D) over the set of observed nodes V . Secondly, we apply RG or NJ to reconstruct a latent subtree over the closed neighborhoods of every internal node in MST(V ; D). If RG (respectively NJ) is used, we term the algorithm CLRG (respectively CLNJ). In the rest of the section, we only describe CLRG for concreteness since CLNJ proceeds along similar lines. Formally, CLRG proceeds as follows: 1. Construct the Chow-Liu tree MST(V ; D) as in (17). Set T = MST(V ; D). 2. Identify the set of internal nodes in MST(V ; D). 3. For each internal node i, let nbd[i; T ] be its closed neighborhood in T and let S = RG(nbd[i; T ], D) be the output of RG with nbd[i; T ] as the set of input nodes. 4. Replace the subtree over node set nbd[i; T ] in T with S. Denote the new tree as T . 5. Repeat steps 3 and 4 until all internal nodes have been operated on. Note that the only difference between the algorithm we just described and CLNJ is Step 3 in which the subroutine NJ replaces RG. Also, observe in Step 3 that RG is only applied to a small subset of nodes which have been identified in Step 1 as possible neighbors in the true latent tree. This reduces the computational complexity of CLRG compared to RG as seen in the following theorem whose proof is provided in Appendix A.5. Let |J| := |V \ Leaf(MST(V ; D))| < m be the number of internal nodes in the MST. Theorem 10 (Correctness and Computational Complexity of CLRG) If Tp ∈ T≥3 is a minimal latent tree and the matrix of information distances D is available, then CLRG outputs the true latent tree Tp correctly in time O(m2 log m + |J|∆3 (MST(V ; D))). 15. The second condition on P(Tblind ) holds when the tree is (almost) homogeneous. 18 Learning Latent Tree Graphical Models h1 h3 1 2 h2 1 3 4 5 2 4 6 1 1 3 5 3 2 h1 h1 3 1 h2 3 h1 h3 1 h2 h2 5 6 2 4 6 4 5 2 6 6 4 5 2 3 4 5 6 (a) (b) (c) (d) (e) (f) Figure 5: Illustration of CLRG. The shaded nodes are the observed nodes and the rest are hidden nodes. The dotted lines denote surrogate mappings for the hidden nodes so for example, node 3 is the surrogate of h3 . (a) The original latent tree, (b) The Chow-Liu tree (MST) over the observed nodes V , (c) The closed neighborhood of node 5 is the input to RG, (d) Output after the first RG procedure, (e) The closed neighborhood of node 3 is the input to the second iteration of RG, (f) Output after the second RG procedure, which is same as the original latent tree. Thus, the computational complexity of CLRG is low when the latent tree Tp has a small maximum degree and a small effective depth (such as the HMM) because (21) implies that ∆(MST(V ; D)) is also small. Indeed, we demonstrate in Section 7 that there is a significant speedup compared to applying RG over the entire observed node set V . We now illustrate CLRG using the example shown in Figure 5. The original minimal latent tree Tp = (W, E) is shown in Figure 5(a) with W = {1, 2, . . . , 6, h1 , h2 , h3 }. The set of observed nodes is V = {1, . . . , 6} and the set of hidden nodes is H = {h1 , h2 , h3 }. The Chow-Liu tree TCL = MST(V ; D) formed using the information distance matrix D is shown in Figure 5(b). Since nodes 3 and 5 are the only internal nodes in MST(V ; D), two RG operations will be executed on the closed neighborhoods of each of these two nodes. In the first iteration, the closed neighborhood of node 5 is the input to RG. This is shown in Figure 5(c) where nbd[4; MST(V ; D)] = {1, 3, 4, 5}, which is then replaced by the output of RG to obtain the tree shown in Figure 5(d). In the next iteration, RG is applied to the closed neighborhood of node 3 in the current tree nbd[3; T ] = {2, 3, 6, h1 } as shown in Figure 5(e). Note that nbd[3; T ] includes h1 ∈ H, which was introduced by RG in the previous iteration. This closed neighborhood is then replaced by the output of the second RG operation and the original latent tree Tp is obtained as shown in Figure 5(f). Observe that the trees obtained at each iteration of CLRG are related to the original latent tree in terms of edge-contraction operations (Robinson and Foulds, 1981), which were defined in Section 5.2. For example, the Chow-Liu tree in Figure 5(b) is obtained from the latent tree Tp in Figure 5(a) by sequentially contracting all edges connecting an observed node to its inverse surrogate set (cf. Lemma 8(ii)). Upon performing an iteration of RG, these contraction operations are inverted and hidden nodes are introduced. For example, in Figure 5(d), the hidden nodes h1 , h2 are introduced after performing RG on the closed neighborhood of node 4 on MST(V ; D). These newly introduced hidden nodes in fact, turn out to be the inverse surrogate set of node 4, i.e., Sg−1 (5) = {5, h1 , h2 }. This is not merely a coincidence and we prove in Appendix A.5 that at each iteration, the set of hidden nodes introduced corresponds to the inverse surrogate set of the internal node. 19 Choi, Tan, Anandkumar, and Willsky We conclude this section by emphasizing that CLGrouping (i.e., CLRG or CLNJ) has two primary advantages. Firstly, as demonstrated in Theorem 10, the structure of all minimal tree-structured graphical models can be recovered by CLGrouping in contrast to CLBlind. Secondly, it typically has much lower computational complexity compared to RG. 5.5 Extension to General Discrete Models For general (i.e., not symmetric) discrete models, the mutual information I(Xi ; Xj ) is in general not monotonic in the information distance dij , defined in (9).16 As a result, Lemma 6 does not hold, i.e., the Chow-Liu tree TCL is not necessarily the same as MST(V ; D). However, Lemma 8 does hold for all minimal latent tree models. Therefore, for general discrete models, we compute MST(V ; D) (instead of the Chow-Liu tree TCL with edge weights I(Xi ; Xj )), and apply RG or NJ to each internal node and its neighbors. This algorithm guarantees that the structure learned using CLGrouping is the same as Tp if the distance matrix D is available. 6. Sample-Based Algorithms for Learning Latent Tree Structures In Sections 4 and 5, we designed algorithms for the exact reconstruction of latent trees assuming that pV is a tree-decomposable distribution and the matrix of information distances D is available. In most (if not all) machine learning problems, the pairwise distributions p(xi , xj ) are unavailable. Consequently, D is also unavailable so RG, NJ and CLGrouping as stated in Sections 4 and 5 are not directly applicable. In this section, we consider extending RG, NJ and CLGrouping to the case when only samples xnV are available. We show how to modify the previously proposed algorithms to accommodate ML estimated distances and we also provide sample complexity results for relaxed versions of RG and CLGrouping. ML Estimation of Information Distances The canonical method for deterministic parameter estimation is via maximum-likelihood (ML) (Serfling, 1980). We focus on Gaussian and symmetric discrete distributions in this section. The generalization to general discrete models is straightforward. For Gaussians graphical models, we use ML to estimate the entries of the covariance matrix,17 i.e., n X (k) (k) b ij = 1 xi xj , Σ n k=1 ∀ i, j ∈ V. (22) b ij /(Σ b ii Σ b jj )1/2 . The The ML estimate of the correlation coefficient is defined as ρbij := Σ estimated information distance is then given by the analogue of (8), i.e., dbij = − log |b ρij |. For symmetric discrete distributions, we estimate the crossover probability θij via ML as18 n 1 X (k) (k) θbij = I xi 6= xj , n k=1 ∀ i, j ∈ V. (23) 16. The mutual information, however, is monotonic in dij for asymmetric binary discrete models. 17. Recall that we assume that the mean of the true random vector X is known and equals to the zero vector so we do not need to subtract the empirical mean in (22). 18. We use I{·} to denote the indicator function. 20 Learning Latent Tree Graphical Models The estimated information distance is given by the analogue of (10), i.e., dbij = −(K − 1) log(1 − K θbij ). For both classes of models, it can easily be verified from the Central Limit Theorem and continuity arguments (Serfling, 1980) that dbij − dij = Op (n−1/2 ), where n is the number of samples. This means that the estimates of the information distances are consistent with rate of convergence being n−1/2 . The m×m matrix of estimated information b = [dbij ]. distances is denoted as D 6.1 Relaxed Recursive Grouping (RG) Given Samples We now show how to relax the canonical RG algorithm described in Section 4 to handle the b is available. Recall that RG calls the TestNodeRelationships procedure case when only D recursively to ascertain child-parent and sibling relationships via equality tests Φijk = dik − djk (cf. Section 3.2). These equality constraints are, in general, not satisfied with the b ijk := dbik − dbjk , which are computed based on the estimated distance estimated differences Φ b in D. Besides, not all estimated distances are equally accurate. Longer distance estimates (i.e., lower correlation estimates) are less accurate for a given number of samples.19 As such, not all estimated distances can be used for testing inter-node relationships reliably. These observations motivate the following three modifications to the RG algorithm: b ijk is constant (across k). 1. Consider using a smaller subset of nodes to test whether Φ b ijk values. 2. Apply a threshold (inequality) test to the Φ 3. Improve on the robustness of the estimated distances dbih in (13) and (14) by averaging. We now describe each of these modifications in greater detail. Firstly, in the relaxed RG b ijk for those estimated distances dbij , dbik and dbjk that are algorithm, we only compute Φ below a prescribed threshold τ > 0 since longer distance estimates are unreliable. As such, for each pair of nodes (i, j) such that dbij < τ , associate the set n o Kij := k ∈ V \{i, j} : max{dbik , dbjk } < τ . (24) This is the subset of nodes in V whose estimated distances to i and j are less than τ . b ijk for all k ∈ Kij only. Compute Φ Secondly, instead of using equality tests in TestNodeRelationships to determine the relationship between nodes i and j, we relax this test and consider the statistic b ij := max Φ b ijk − min Φ b ijk Λ (25) b ij < ǫ, Λ (26) k∈Kij k∈Kij b ij in (25) is close to zero, then nodes i and j are likely to be in the same Intuitively, if Λ family. Thus, declare that nodes i, j ∈ V are in the same sibling group if 19. In fact, by using a large deviation result in Shen (2007, Theorem 1), we can formally show that a larger number of samples is required to get a good approximation of ρik if it is small compared to when ρik is large. 21 Choi, Tan, Anandkumar, and Willsky for another threshold ǫ > 0. Similarly, an observed node k is identified as a parent node if |dbik + dbkj − dbij | < ǫ for all i and j in the sibling group. Thirdly, in order to further improve on the quality of the distance estimate dbih of a newly introduced hidden node to observed nodes, we compute dbih using (13) with different pairs of j ∈ C(h) and k ∈ Kij , and take the average as follows:   X X 1 1 b ijk  .  Φ dbij + dbih = 2(|C(h)| − 1) |Kij | j∈C(h) (27) k∈Kij Similarly, for any other node k ∈ / C(h), we compute dbkh using all child nodes in C(h) and C(k) (if C(k) 6= ∅) as follows: ( 1 P b b if k ∈ V, i∈C(h) (dik − dih ), |C(h)| b P dkh = (28) 1 b b b (i,j)∈C(h)×C(k) (dij − dih − djk ), otherwise. |C(h)||C(k)| It is easy to verify that if dbih and dbkh are equal to dih and dkh respectively, then (27) and (28) reduce to (13) and (14) respectively. The following theorem shows that relaxed RG is consistent, and with appropriately chosen thresholds ǫ and τ , it has the sample complexity logarithmic in the number of observed variables. The proof follows from standard Chernoff bounds and is provided in Appendix A.6. Theorem 11 (Consistency and Sample Complexity of Relaxed RG) (i) Relaxed RG is structurally consistent for all Tp ∈ T≥3 . In addition, it is risk consistent for Gaussian and symmetric discrete distributions. (ii) Assume that the effective depth is δ(Tp ; V ) = O(1) b For every (i.e., constant in m) and relaxed RG is used to reconstruct the tree given D. η > 0, there exists thresholds ǫ, τ > 0 such that if √ n > C log(m/ 3 η) (29) for some constant C > 0, the error probability for structure reconstruction in (5) is bounded above by η. If, in addition, p is a Gaussian or symmetric discrete distribution and n > √ C ′ log(m/ 3 η), the error probability for distribution reconstruction in (6) is also bounded above by η. Thus, the sample complexity of relaxed RG, which is the number of samples required to achieve a desired level of accuracy, is logarithmic in m, the number of observed variables. As we observe from (29), the sample complexity for RG is logarithmic in m for shallow trees (i.e., trees where the effective depth is constant). This is in contrast to NJ where the sample complexity is super-polynomial in the number of observed nodes for the HMM (St. John et al., 2003; Lacey and Chang, 2006). RG with k-means Clustering In practice, if the number of samples is limited, the distance estimates dbij are noisy and it is difficult to select the threshold ǫ in Theorem 11 to identify sibling nodes reliably. In our 22 Learning Latent Tree Graphical Models experiments, we employ a modified version of the k-means clustering algorithm to cluster b ij , defined in (25), as a family. Recall that we test each Λ b ij a set of nodes with small Λ locally with a fixed threshold ǫ in (26). In contrast, the k-means algorithm provides a global scheme and circumvents the need to select the threshold ǫ. We adopt the silhouette method b ij to select optimal the number of clusters k. (Rousseeuw, 1987) with dissimilarity measure Λ 6.2 Relaxed Neighbor-Joining Given Samples In this section, we describe how NJ can be relaxed when the true distances are unavailable. We relax the NJ algorithm by using ML estimates of the distances dbij in place of unavailable distances dij . NJ typically assume that all observed nodes are at the leaves of the latent tree, so after learning the latent tree, we perform the following post-processing step: If there exists an edge (i, h) ∈ W × H with dbih < ǫ′ (for a given threshold ǫ′ > 0), then (i, h) is contracted to a single node whose label is i. The sample complexity of NJ is known to be O(exp(diam(Tp )) log m) (St. John et al., 2003) and thus does not scale well when the latent tree Tp has a large diameter. Comparisons between the sample complexities of other closely related latent tree learning algorithms are discussed in Atteson (1999); Erdős et al. (1999); Csűrös (2000) and St. John et al. (2003). 6.3 Relaxed CLGrouping Given Samples In this section, we discuss how to modify CLGrouping (CLRG and CLNG) when we only b The relaxed version of CLGrouping have access to the estimated information distance D. differs from CLGrouping in two main aspects. Firstly, we replace the edge weights in the construction of the MST in (17) with the estimated information distances dbij , i.e., b := argmin TbCL = MST(V ; D) X T ∈T (V ) (i,j)∈T dbij . (30) The procedure in (30) can be shown to be equivalent to the learning of the ML tree structure given samples xnV if pV is a Gaussian or symmetric discrete distribution.20 It has also been shown that the error probability of structure learning Pr(TbCL 6= TCL ) converges to zero exponentially fast in the number of samples n (Tan et al., 2009, 2010). Secondly, for CLRG (respectively CLNJ), we replace RG (respectively NJ) with the relaxed version of RG (respectively NJ). The sample complexity result of CLRG (and its proof) is similar to Theorem 11 and the proof is provided in Appendix A.7. Theorem 12 (Consistency and Sample Complexity of Relaxed CLRG) (i) Relaxed CLRG is structurally consistent for all Tp ∈ T≥3 . In addition, it is risk consistent for Gaussian and symmetric discrete distributions. (ii) Assume that the effective depth is δ(Tp ; V ) = O(1) (i.e., constant in m). Then the sample complexity of relaxed CLRG is logarithmic in m. After CLRG has been completed, as a final post-processing step, if we find that there exists an estimated distance dbih on an edge with i ∈ W and h ∈ H in the learned model 20. This follows from the observation that the ML search for the optimal structure is equivalent to the KL-divergence minimization problem in (15) with pV replaced by pbV , the empirical distribution of xn V. 23 Choi, Tan, Anandkumar, and Willsky which is smaller than some ǫ′ > 0 (which we specify in our experiments), then edge (i, h) is contracted to the single node i. This is similar to relaxed NJ and serves to contract all strong edges in the learned model. 6.4 Regularized CLGrouping for Learning Latent Tree Approximations For many practical applications, it is of interest to learn a latent tree that approximates the given empirical distribution. In general, introducing more hidden variables enables better fitting to the empirical distribution, but it increases the model complexity and may lead to overfitting. The Bayesian Information Criterion (Schwarz, 1978) provides a trade-off between model fitting and model complexity, and is defined as follows: κ(Tb) BIC(Tb) = log p(xnV ; Tb) − log n 2 (31) where Tb is a latent tree structure and κ(Tb) is the number of free parameters, which grows linearly with the number of hidden variables because Tb is a tree. Here, we describe regularized CLGrouping, in which we use the BIC in (31) to specify a stopping criterion on the number of hidden variables added. For each internal node and its neighbors in the Chow-Liu tree, we use relaxed NJ or RG to learn a latent subtree. Unlike in regular CLGrouping, before we integrate this subtree into our model, we compute its BIC score. Computing the BIC score requires estimating the maximum likelihood parameters for the models, so for general discrete distributions, we run the EM algorithm on the subtree to estimate the parameters.21 After we compute the BIC scores for all subtrees corresponding to all internal nodes in the Chow-Liu tree, we choose the subtree that results in the highest BIC score and incorporate that subtree into the current tree model. The BIC score can be computed efficiently on a tree model with a few hidden variables. Thus, for computational efficiency, each time a set of hidden nodes is added to the model, we generate samples of hidden nodes conditioned on the samples of observed nodes, and use these augmented samples to compute the BIC score approximately when we evaluate the next subtree to be integrated in the model. If none of the subtrees increases the BIC score (i.e., the current tree has the highest BIC score), the procedure stops and outputs the estimated latent tree. Alternatively, if we wish to learn a latent tree with a given number of hidden nodes, we can used the BIC-based procedure mentioned in the previous paragraph to learn subtrees until the desired number of hidden nodes is introduced. Depending on whether we use NJ or RG as the subroutine, we denote the specific regularized CLGrouping algorithm as regCLNJ or regCLRG. This approach of using an approximation of the BIC score has been commonly used to learn a graphical model with hidden variables (Elidan and Friedman, 2005; Zhang and Kočka, 2004). However, for these algorithms, the BIC score needs to be evaluated for a large subset of nodes, whereas in CLGrouping, the Chow-Liu tree among observed variables prunes out many subsets, so we need to evaluate BIC scores only for a small number of candidate subsets (the number of internal nodes in the Chow-Liu tree). 21. Note that for Gaussian and symmetric discrete distributions, the model parameters can be recovered from information distances directly using (8) or (10). 24 Learning Latent Tree Graphical Models (a) Double star (b) HMM (c) 5-complete Figure 6: Latent tree structures used in our simulations. 7. Experimental Results In this section, we compare the performances of various latent tree learning algorithms. We first show simulation results on synthetic datasets with known latent tree structures to demonstrate the consistency of our algorithms. We also analyze the performance of these algorithms when we change the underlying latent tree structures. Then, we show that our algorithms can approximate arbitrary multivariate probability distributions with latent trees by applying them to two real-world datasets, a monthly stock returns example and the 20 newsgroups dataset. 7.1 Simulations using Synthetic Datasets In order to analyze the performances of different tree reconstruction algorithms, we generate samples from known latent tree structures with varying sample sizes and apply reconstruction algorithms. We compare the neighbor-joining method (NJ) (Saitou and Nei, 1987) with recursive grouping (RG), Chow-Liu Neighbor Joining (CLNJ), and Chow-Liu Recursive Grouping (CLRG). Since the algorithms are given only samples of observed variables, we use the sample-based algorithms described in Section 6. For all our experiments, we use the same edge-contraction threshold ǫ′ = − log 0.9 (see Sections 6.2 and 6.3), and set τ in Section 6.1 to grow logarithmically with the number of samples. Figure 6 shows the three latent tree structures used in our simulations. The doublestar has 2 hidden and 80 observed nodes, the HMM has 78 hidden and 80 observed nodes, and the 5-complete tree has 25 hidden and 81 observed nodes including the root node. For simplicity, we present simulation results only on Gaussian models but note that the behavior on discrete models is similar. All correlation coefficients on the edges ρij were independently drawn from a uniform distribution supported on [0.2, 0.8]. The performance of each method is measured by averaging over 200 independent runs with different parameters. We use the following performance metrics to quantify the performance of each algorithm in Figure 7: 25 Choi, Tan, Anandkumar, and Willsky Struture Recovery Error Rate Robinson-Foulds Metric Error in Hidden Variables 70 50 CLRG 0.8 0 CLNJ 50 −0.5 40 −1 40 0.6 30 −1.5 30 −2 0.4 20 −2.5 20 0.2 log10 (KL-divergence) 0.5 NJ 60 Double Star 1 60 RG 1 −3 10 10 −3.5 0 100 1k 2k 10k 20k 100k 200k 0 100 1k 2k Number of Samples 10k 20k 100k 200k 1k 2k 10k 20k 100k 200k 1k 2k 10k 20k 100k 200k Number of Samples 0.5 60 180 0 50 160 −0.5 0.8 140 40 −1 120 HMM −4 100 Number of Samples 200 1 0 100 Number of Samples 0.6 100 −1.5 30 80 0.4 −2 20 60 −2.5 40 0.2 10 −3 20 0 1k 2k 10k 20k 0 1k 100k 200k 2k Number of Samples 10k 20k 100k 200k 0 1k 2k Number of Samples 10k 20k 100k 200k 2k Number of Samples 10k 20k 100k 200k Number of Samples 40 60 1 0.5 35 50 0 30 0.8 5-complete −3.5 1k 40 −0.5 25 0.6 30 −1 20 15 0.4 20 −1.5 10 0.2 10 −2 5 0 1k 2k 10k 20k Number of Samples 100k 200k 0 1k 2k 10k 20k 100k 200k Number of Samples 0 1k 2k 10k 20k Number of Samples 100k 200k −2.5 1k 2k 10k 20k 100k 200k Number of Samples Figure 7: Performance of RG, NJ, CLRG, and CLNJ for the latent trees shown in Figure 6. (i) Structure recovery error rate: This is the proportion of times that the proposed algorithm fails to recover the true latent tree structure. Note that this is a very strict measure since even a single wrong hidden node or misplaced edge results in an error for the entire structure. (ii) Robinson Foulds metric (Robinson and Foulds, 1981): This popular phylogenetic tree-distortion metric computes the number of graph transformations (edge contraction or expansion) needed to be applied to the estimated graph in order to get the correct structure. This metric quantifies the difference in the structures of the estimated and true models. (iii) Error in the number of hidden variables: We compute the average number of hidden variables introduced by each method and plot the absolute difference between the average estimated hidden variables and the number of hidden variables in the true structure. 26 Learning Latent Tree Graphical Models HMM 5-complete Double star RG 10.16 7.91 1.43 NJ 0.02 0.02 0.01 CLRG 0.10 0.26 0.76 CLNJ 0.05 0.06 0.20 Table 1: Average running time of each algorithm in seconds. (iv) KL-divergence D(pV || pbnV ): This is a measure of the distance between the estimated and the true models over the set of observed nodes V .22 We first note that from the structural error rate plots that the double star is the easiest structure to recover and the 5-complete tree is the hardest. In general, given the same number of observed variables, a latent tree with more hidden variables or larger effective depth (see Section 2) is more difficult to recover. For the double star, RG clearly outperforms all other methods. With only 1,000 samples, it recovers the true structure exactly in all 200 runs. On the other hand, CLGrouping performs significantly better than RG for the HMM. There are two reasons for such performance differences. Firstly, for Gaussian distributions, it was shown (Tan et al., 2010) that given the same number of variables and their samples, the Chow-Liu algorithm is most accurate for a chain and least accurate for a star. Since the Chow-Liu tree of a latent double star graph is close to a star, and the Chow-Liu tree of a latent HMM is close to a chain, the Chow-Liu tree tend to be more accurate for the HMM than for the double star. Secondly, the internal nodes in the Chow-Liu tree of the HMM tend to have small degrees, so we can apply RG or NJ to a very small neighborhood, which results in a significant improvement in both accuracy and computational complexity. Note that NJ is particularly poor at recovering the HMM structure. In fact, it has been shown that even if the number of samples grows polynomially with the number of observed variables (i.e., n = O(mB ) for any B > 0), it is insufficient for NJ to recover HMM structures (Lacey and Chang, 2006). The 5-complete tree has two layers of hidden nodes, making it very difficult to recover the exact structure using any method. CLNJ has the best structure recovery error rate and KL divergence, while CLRG has the smallest Robinson-Foulds metric. Table 1 shows the running time of each algorithm averaged over 200 runs and all sample sizes. All algorithms are implemented in MATLAB. As expected, we observe that CLRG is significantly faster than RG for HMM and 5-complete graphs. NJ is fastest, but CLNJ is also very efficient and leads to much more accurate reconstruction of latent trees. Based on the simulation results, we conclude that for a latent tree with a few hidden variables, RG is most accurate, and for a latent tree with a large diameter, CLNJ performs the best. A latent tree with multiple layers of hidden variables is more difficult to recover correctly using any method, and CLNJ and CLRG outperform NJ and RG. 22. Note that this is not the same quantity as in (6) because if the number of hidden variables is estimated incorrectly, D(p || pbn ) is infinite so we plot D(pV || pbn V ) instead. However, for Gaussian and symmetric discrete distributions, D(p || pbn ) converges to zero in probability since the number of hidden variables is estimated correctly asymptotically. 27 Choi, Tan, Anandkumar, and Willsky CL NJ RG CLNJ CLRG Log-Likelihood -13,321 -12,400 -14,042 -11,990 -12,879 BIC -13,547 -12,747 -14,300 -12,294 -13,174 # Hidden 0 45 12 29 26 # Parameters 84 129 96 113 110 Time (secs) 0.15 0.02 21.15 0.24 0.40 Table 2: Comparison of the log-likelihood, BIC, number of hidden variables introduced, number of parameters, and running time for the monthly stock returns example. BIC score -12,000 -12,500 -13,000 -13,500 -14,000 -14,500 CL NJ RG CLNJ CLRG Figure 8: Plot of BIC scores for the monthly stock returns example. 7.2 Monthly Stock Returns We apply our latent tree learning algorithms to model the dependency structure of monthly stock returns of 84 companies in the S&P 100 stock index.23 We use the samples of the monthly returns from 1990 to 2007. As shown in Table 2 and Figure 8, CLNJ achieves the highest log-likelihood and BIC scores. NJ introduces more hidden variables than CLNJ and has lower log-likelihoods, which implies that starting from a Chow-Liu tree helps to get a better latent tree approximation. Figure 11 shows the latent tree structure learned using the CLNJ method. Each observed node is labeled with the ticker of the company. Note that related companies are closely located on the tree. Many hidden nodes can be interpreted as industries or divisions. For example, h1 has Verizon, Sprint, and T-mobile as descendants, and can be interpreted as the telecom industry, and h3 correspond to the technology division with companies such as Microsoft, Apple, and IBM as descendants. Nodes h26 and h27 group commercial banks together, and h25 has all retail stores as child nodes. 7.3 20 Newsgroups with 100 Words For our last experiment, we apply our latent tree learning algorithms to the 20 Newsgroups dataset with 100 words.24 The dataset consists of 16,242 binary samples of 100 23. We disregard 16 companies that have been listed on S&P 100 only after 1990. 24. http://cs.nyu.edu/~ roweis/data/20news_w100.mat 28 Learning Latent Tree Graphical Models words, indicating whether each word appears in each posting or not. In addition to the Chow-Liu tree (CL), NJ, RG, CLNJ, and CLRG, we also compare the performances with the regCLNJ and regCLRG (described in Section 6.4), the latent cluster model (LCM) (Lazarsfeld and Henry, 1968), and BIN, which is a greedy algorithm for learning latent trees (Harmeling and Williams, 2010). Table 3 shows the performance of different algorithms, and Figure 9 plots the BIC score. We use the MATLAB code (a small part of it is implemented in C) provided by Harmeling and Williams (2010)25 to run LCM and BIN. Note that although LCM has only one hidden node, the hidden node has 16 states, resulting in many parameters. We also tried to run the algorithm by Chen et al. (2008), but their JAVA implementation on this dataset did not complete even after several days. For NJ, RG, CLNJ, and CLRG, we learned the structures using only information distances (defined in (9)) and then used the EM algorithm to fit the parameters. For regCLNJ and regCLRG, the model parameters are learned during the structure learning procedure by running the EM algorithm locally, and once the structure learning is over, we refine the parameters by running the EM algorithm for the entire latent tree. All methods are implemented in MATLAB except the E-step of the EM algorithm, which is implemented in C++. Despite having many parameters, the models learned via LCM have the best BIC score. However, it does not reveal any interesting structure and is computationally more expensive to learn. In addition, it may result in overfitting. In order to show this, we split the dataset randomly and use half as the training set and the other half as the test set. Table 4 shows the performance of applying the latent trees learned from the training set to the test set, and Figure 10 shows the log-likelihood on the training and the test sets. For LCM, the test log-likelihood drops significantly compared to the training log-likelihood, indicating that LCM is overfitting the training data. NJ, CLNJ, and CLRG achieve high log-likelihood scores on the test set. Although regCLNJ and regCLRG do not result in a better BIC score, they introduce fewer hidden variables, which is desirable if we wish to learn a latent tree with small computational complexity, or if we wish to discover a few hidden variables that are meaningful in explaining the dependencies of observed variables. Figure 12 shows the latent tree structure learned using regCLRG from the entire dataset. Many hidden variables in the tree can be interpreted as topics - h5 as sports, h9 as computer technology, h13 as medical, etc. Note that some words have multiple meanings and appear in different topics - e.g., program can be used in the phrase “space program” as well as “computer program”, and win may indicate the windows operating system or winning in sports games. 8. Conclusion In this paper, we proposed algorithms to learn a latent tree model from the information distances of observed variables. Our first algorithm, recursive grouping, identifies sibling and parent-child relationships and introduces hidden nodes recursively. Our second algorithm, CLGrouping, first learns the Chow-Liu tree among observed variables and then applies latent-tree-learning subroutines such as recursive grouping or neighbor joining locally to each internal node in the Chow-Liu tree and its neighbors. These algorithms are 25. http://people.kyb.tuebingen.mpg.de/harmeling/code/ltt-1.3.tar 29 Choi, Tan, Anandkumar, and Willsky CL LCM BIN NJ RG CLNJ CLRG regCLNJ regCLRG Log-Likelihood BIC Hidden Params -238,713 -223,096 -232,042 -230,575 -239,619 -230,858 -231,279 -235,326 -234,012 -239,677 -230,925 -233,952 -232,257 -240,875 -232,540 -232,738 -236,553 -235,229 0 1 98 74 30 74 51 27 26 199 1,615 394 347 259 347 301 253 251 Total 8.9 8,835.9 3022.6 1611.2 927.1 1479.6 1224.6 630.8 606.9 Time (s) Structure 3.3 30.8 2.7 3.1 449.7 493.0 EM 1608.2 896.4 1476.8 1224.6 181.1 113.9 Table 3: Comparison between various algorithms on the newsgroup set. BIC score -230,000 -232,000 -234,000 -236,000 -238,000 -240,000 -242,000 CL LCM BIN NJ RG CLNJ CLRG regCLNJ regCLRG Figure 9: The BIC scores of various algorithms on the newsgroup set. structurally consistent (and risk consistent as well in the case of Gaussian and discrete symmetric distributions), and have sample complexity logarithmic in the number of observed variables. Using simulations with synthetic datasets, we showed that RG performs well when the number of hidden variables is small, and that CLGrouping performs significantly better than other algorithms when there are many hidden variables in the latent tree. Using both Gaussian and discrete real-world data, we compared the performances of our algorithms to other EM-based approaches and the neighbor-joining method, and our algorithm show superior results in both accuracy (measured by KL-divergence and graph distance) and computational efficiency. In addition, we introduced regularized CLGrouping, which is useful in learning a latent tree approximation with a given number of hidden nodes. The MATLAB implementation of our algorithms can be downloaded from the project webpage http://people.csail.mit.edu/myungjin/latentTree.html. Acknowledgement The authors thank Prof. Sekhar Tatikonda (Yale University) for discussions and comments. This work was supported in part by AFOSR under Grant FA9550-08-1-1080 and in part by MURI under AFOSR Grant FA9550-06-1-0324. Vincent Tan is also supported by A*STAR, Singapore. 30 Learning Latent Tree Graphical Models CL LCM BIN NJ RG CLNJ CLRG regCLNJ regCLRG Train Log-Like BIC -119,013 -119,909 -112,746 -117,288 -117,172 -118,675 -115,319 -116,908 -118,280 -119,248 -115,372 -116,987 -115,565 -116,920 -117,723 -118,924 -116,980 -118,119 Test Log-Like BIC -120,107 -121,003 -116,884 -120,949 -117,957 -119,460 -116,011 -117,600 -119,181 -120,149 -116,036 -117,652 -116,199 -117,554 -118,606 -119,808 -117,652 -118,791 Hidden Params 0 1 78 77 8 80 51 34 27 199 1,009 334 353 215 359 301 267 253 Total 3.0 3,197.7 1,331.3 802.8 137.6 648.0 506.0 425.5 285.7 Time (s) Struct 1.3 7.6 1.5 1.7 251.3 236.5 EM 801.5 130.0 646.5 504.3 174.2 49.2 Table 4: Comparison between various algorithms on the newsgroup dataset with a train/test split. Log-likelihood -112,000 -113,000 Train Test -114,000 -115,000 -116,000 -117,000 -118,000 -119,000 -120,000 -121,000 CL LCM BIN NJ RG CLNJ CLRG regCLNJ regCLRG Figure 10: Train and test log-likelihood scores of various algorithms on the newsgroup dataset with a train/test split. Appendix A. Proofs A.1 Proof of Lemma 4: Sibling Grouping We prove statement (i) in Lemma 4 using (12) in Proposition 3. Statement (ii) follows along similar lines and its proof is omitted for brevity. If : From the additive property of information distances in (12), if i is a leaf node and j is its parent, dik = dij + djk and thus Φijk = dij for all k 6= i, j. Only If: Now assume that Φijk = dij for all k ∈ V \ {i, j}. In order to prove that i is a leaf node and j is its parent, assume to the contrary, that i and j are not connected with an edge, then there exists a node u 6= i, j on the path connecting i and j. If u ∈ V , then let k = u. Otherwise, let k be an observed node in the subtree away from i and j (see Figure 13(a)), which exists since Tp ∈ T≥3 . By the additive property of information distances in (12) and the assumption that all distances are positive, dij = diu + duj > diu − duj = dik − dkj = Φijk 31 (32) UTX h2 BA SLB h11 h3 h28 VZ h21 HPQ GE h1 S h13 MS C AXP h4 MER h27 AIG WB BK h24 h25 INTC MSFT h20 IBM XRX T h8 JPM h29 TXN CAT IP AA CBS HON DELL AAPL WY PEP MDT h9 UNH h26 SNS ORCL CL h6 ABT h7 h22 BAC DD RF BNI HNZ CPB SLE KO PG h16 AMGN CMCSA XOM RTN OXY CVS WMB AEP WFC USB h12 JNJ NSC GD h18 MMM MO CVX TYC TGT AVP HAL h23 WMT HD h19 h14 COP DIS F FDX DOW BHI MCD h15 EMC h10 h5 h17 WYE MRK BAX CI PFE BMY EXC ETR SO Choi, Tan, Anandkumar, and Willsky 32 Figure 11: Tree structure learned from monthly stock returns using CLNJ. SP100 h17 h3 h4 power h14 religion h20 war children orbit lunar earth satellite moon law state world rights human israel jews h8 solar technology mission god bible mars h1 gun h2 christian jesus space launch shuttle nasa health case food aids course evidence fact question program h9 h21 insurance version msg studies water h12 car files h25 dealer h15 cancer disease doctor patients engine honda games baseball h18 h11 image puck league season players driver card team fans h7 hockey win nhl won video graphics h16 h23 disk system h10 number memory data h22 scsi h26 h19 dos pc drive software problem help mac display computer science hit phone oil h5 h6 format vitamin windows bmw email ftp medicine h13 university server h24 research Learning Latent Tree Graphical Models 33 Figure 12: Tree structure learned from 20 newsgroup dataset using regCLRG. government president Choi, Tan, Anandkumar, and Willsky Vi\j j i u Sg(i) j Vj\i u k (a) h Sg(j) i k i u j l k (b) (c) j k (d) Figure 13: Shaded nodes indicate observed nodes and the rest indicate hidden nodes. (a),(b) Figures for Proof of Lemma 4. Dashed red line represent the subtrees away from i and j. (c) Figure for Proof of Lemma 8(i). (d) Figure for Proof of Lemma 8(iI) which is a contradiction. If i is not a leaf node in Tp , then there exist a node u 6= i, j such that (i, u) ∈ Ep . Let k = u if u ∈ V , otherwise, let k be an observed node in the subtree away from i and j (see Figure 13(b)). Then, Φijk = dik − djk = −dij < dij , (33) which is again a contradiction. Therefore, (i, j) ∈ Ep and i is a leaf node. A.2 Proof of Theorem 5: Correctness and Computational Complexity of RG The correctness of RG follows from the following observations: Firstly, at each iteration of RG, the sibling groups are identified correctly by Lemma 4. Since the new parent node added to a partition which does not contain an observed parent corresponds to a hidden node (in the original latent tree), a subforest of Tp is recovered at each iteration. Secondly, from Proposition 3, for all i, j in the active set Y , the information distances dij can be computed exactly with (13) and (14). Thirdly, since the information distances between nodes in the updated active set Y are known, they can now be regarded as observed nodes. Subforests of Tp are constructed at each iteration and when |Y | ≤ 2, the entire latent tree is recovered. The computational complexity follows from the fact there are a maximum of O(m3 ) differences Φijk = dik − djk that we have to compute at each iteration of RG. Furthermore, there are at most diam(Tp ) subsets in the coarsest partition (cf. step 3) of Y at the first iteration, and the number of subsets reduce at least by 2 from one iteration to the next due to the assumption that Tp ∈ T≥3 . This proves the claim that the computational complexity is upper bounded by O(diam(Tp )m3 ). A.3 Proof of Lemma 8: Properties of the MST (i) For an edge (i, j) ∈ Ep such that Sg(i) 6= Sg(j), let Vi\j ⊂ V and Vj\i ⊂ V denote observed nodes in the subtrees obtained by the removal of edge (i, j), where the former includes node i and excludes node j and vice versa (see Figure 13(c)). Using part (ii) of the lemma and the fat that Sg(i) 6= Sg(j), it can be shown that Sg(i) ∈ Vi\j and Sg(j) ∈ Vj\i . Since (i, j) lies on the unique path from k to l on Tp , for all observed nodes k ∈ Vi\j , l ∈ Vj\i , 34 Learning Latent Tree Graphical Models we have dkl = dki + dij + djl ≥ dSg(i),i + dij + dSg(j),j = dSg(i),Sg(j) , (34) where the inequality is from the definition of surrogacy and the final equality uses the fact that Sg(i) 6= Sg(j). By using the property of the MST that (Sg(i), Sg(j)) is the shortest edge from Vi\j to Vj\i , we have (20). (ii) First assume that we have a tie-breaking rule consistent across all hidden nodes so that if duh = dvh = mini∈V dih and duh′ = dvh′ = mini∈V dih′ then both h and h′ choose the same surrogate node. Let j ∈ V , h ∈ Sg−1 (j), and let u be a node on the path connecting h and j (see Figure 13(d)). Assume that Sg(u) = k 6= j. If duj > duk , then dhj = dhu + duj > dhu + duk = dhk , (35) which is a contradiction since j = Sg(h). If duj = duk , then dhj = dhk , which is again a contradiction to the consistent tie-breaking rule. Thus, the surrogate node of u is j. (iii) First we claim that u |Sg−1 (i)| ≤ ∆(Tp ) l δ(Tp ;V ) . (36) To prove this claim, let γ be the longest (worst-case) graph distance of any hidden node h ∈ H from its surrogate, i.e., γ := max |Path(h, Sg(h); Tp )|. (37) h∈H From the degree bound, for each i ∈ V , there are at most ∆(Tp )γ hidden nodes that are within the graph distance of γ,26 so |Sg−1 (i)| ≤ ∆(Tp )γ (38) for all i ∈ V . Let d∗ := maxh∈H dh,Sg(h) be the longest (worst-case) information distance between a hidden node and its surrogate. From the bounds on the information distances, lγ ≤ d∗ . In addition, for each h ∈ H, let z(h) := argminj∈V |Path((h, j); Tp )| be the observed node that is closest to h in graph distance. Then, by definition of the effective depth, dh,Sg(h) ≤ dh,z(h) ≤ uδ for all h ∈ H, and we have d∗ ≤ uδ. Since lγ ≤ d∗ ≤ uδ, we also have γ ≤ uδ/l. (39) Combining this result with (38) establishes the claim in (36). Now consider (a) (b) u ∆(MST(V ; D)) ≤ ∆(Tp ) max |Sg−1 (i)| ≤ ∆(Tp )1+ l δ(Tp ;V ) i∈V (40) where (a) is a result of the application of (20) and (b) results from (36). This completes the proof of the claim in (21) in Lemma 8. 26. The maximum size of the inverse surrogate set in (37) is attained by a ∆(Tp )-ary complete tree. 35 Choi, Tan, Anandkumar, and Willsky 2 v1 3 h1 h3 1 1 1 ð 5 6 4 h2 3 1 2 6 4 5 2 RG(A2,d) h1 h1 3 6 4 5 h3 1 h2 v2 v1 ð h2 H2 2 3 4 55 v2 T0 h2 A2 RG(A1,d) h1 H1 A1 T1 T1 T2 6 (b) 2 3 4 5 S1 6 (a) 2 v1 3 5 6 4 EC(Tp, V0) 1 1 ð h1 h1 3 H1 1 h2 2 6 3 4 5 2 6 4 5 h3 1 h2 v2 v1 h1 S2 ð h2 H2 2 3 4 5 v2 EC(Tp, V1) EC(Tp, V1) 6 EC(Tp, V2) (c) Figure 14: Figure for Proof of Theorem 10. (a) Original latent tree. (b) Illustration of CLGrouping. (c) Illustration of the trees constructed using edge contractions. A.4 Proof of Theorem 9: Correctness and Computational Complexity of CLBlind It suffices to show that the Chow-Liu tree MST(V ; d) is a transformation of the true latent tree Tp (with parameters such that p ∈ P(Tblind )) as follows: contract the edge connecting each hidden variable h with its surrogate node Sg(h) (one of its children and a leaf by assumption). Note that the blind transformation on the MST is merely the inverse mapping of the above. From (20), all the children of a hidden node h, except its surrogate Sg(h), are neighbors of its surrogate node Sg(h) in MST(V ; d). Moreover, these children of h which are not surrogates of any hidden nodes are leaf nodes in the MST. Similarly for two hidden nodes h1 , h2 ∈ H such that (h1 , h2 ) ∈ Ep , (Sg(h1 ), Sg(h2 )) ∈ MST(V ; d) from Lemma 8(i). Hence, CLBlind outputs the correct tree structure Tp . The computational complexity follows from the fact that the blind transformation is linear in the number of internal nodes, which is less than the number of observed nodes, and that learning the Chow-Liu tree takes O(m2 log m) operations. A.5 Proof of Theorem 10: Correctness and Computational Complexity of CLRG We first define some new notations. Notation: Let I := V \ Leaf(MST(V ; d)) be the set of internal nodes. Let v r ∈ I be the internal node visited at iteration r, and let H r be all hidden nodes in the inverse surrogate set Sg−1 (v r ), i.e., H r = Sg−1 (v r ) \ {v r }. Let Ar := nbd[v r ; T r−1 ], and hence Ar is the node set input to the recursive grouping routine at iteration r, and let RG(Ar , d) be the output latent tree learned by recursive grouping. Define T r as the tree output at the end 36 Learning Latent Tree Graphical Models of r iterations of CLGrouping. Let V r := {v r+1 , v r+2 , . . . , v |I| } be the set of internal nodes that have not yet been visited by CLGrouping at the end of r iterations. Let EC(Tp , V r ) be the tree constructed using edge contractions as follows: in the latent tree Tp , we contract edges corresponding to each node u ∈ V r and all hidden nodes in its inverse surrogate set Sg−1 (u). Let S r be a subtree of EC(Tp , V r ) spanning v r , H r and their neighbors. For example, in Figure 14, the original latent tree Tp is shown in Figure 14(a), and 0 T , T 1 , T 2 are shown in Figure 14(b). The set of internal nodes is I = {3, 5}. In the first iteration, v 1 = 5, A1 = {1, 3, 4, 5} and H 1 = {h1 , h2 }. In the second iteration, v 2 = 3, A2 = {2, 3, 6, h1 } and H 1 = {h3 }. V 0 = {3, 5}, V 1 = {3}, and V 2 = ∅, and in Figure 14(c), we show EC(Tp , V 0 ), EC(Tp , V 1 ), and EC(Tp , V 2 ). In EC(Tp , V 1 ), S 1 is the subtree spanning 5, h1 , h2 and their neighbors, i.e., {1, 3, 4, 5, h1 , h2 }. In EC(Tp , V 2 ), S 2 is the subtree spanning 3, h3 and their neighbors, i.e., {2, 3, 6, h1 , h3 }. Note that T 0 = EC(Tp , V 0 ), T 1 = EC(Tp , V 1 ), and T 2 = EC(Tp , V 2 ); we show below that this holds for all CLGrouping iterations in general. We prove the theorem by induction on the iterations r = 1, . . . , |I| of the CLGrouping algorithm. Induction Hypothesis: At the end of k iterations of CLGrouping, the tree obtained is T k = EC(Tp , V k ), ∀ k = 0, 1, . . . , |I|. (41) In words, the latent tree after k iterations of CLGrouping can be constructed by contracting each surrogate node in Tp that has not been visited by CLGrouping with its inverse surrogate set. Note that V |I| = ∅ and EC(Tp , V |I|) is equivalent to the original latent tree Tp . Thus, if the above induction in (41) holds, then the output of CLGrouping T |I| is the original latent tree. Base Step r = 0: The claim in (41) holds since V 0 = I and the input to the CLGrouping procedure is the Chow-Liu tree MST(V ; D), which is obtained by contracting all surrogate nodes and their inverse surrogate sets (see Section 5.2). Induction Step: Assume (41) is true for k = 1, . . . , r − 1. Now consider k = r. We first compare the two latent trees EC(Tp , V r ) and EC(Tp , V r−1 ). By the definition of EC, if we contract edges with v r and the hidden nodes in its inverse surrogate set H r on the tree EC(Tp , V r ), then we obtain EC(Tp , V r−1 ), which is equivalent to T r−1 by the induction assumption. Note that as shown in Figure 14, this transformation is local to the subtree S r : contracting v r with H r on EC(Tp , V r ) transforms S r into a star graph with v r at its center and the hidden nodes H r removed (contracted with v r ). Recall that the CLGrouping procedure replaces the induced subtree of Ar in T r−1 (which is precisely the star graph mentioned above by the induction hypothesis) with RG(Ar , d) to obtain T r . Thus, to prove that T r = EC(Tp , V r ), we only need to show that RG reverses the edge-contraction operations on v r and H r , that is, the subtree S r = RG(Ar , d). We first show that S r ∈ T≥3 , i.e., it is identifiable (minimal) when Ar is the set of visible nodes. This is because an edge contraction operation does not decrease the degree of any existing nodes. Since Tp ∈ T≥3 , all hidden nodes in EC(Tp , V r ) have degrees equal to or greater than 3, and since we are including all neighbors of H r in the subtree S r , we have S r ∈ T≥3 . By Theorem 5, RG reconstructs all latent trees in T≥3 and hence, S r = RG(Ar , d). 37 Choi, Tan, Anandkumar, and Willsky The computational complexity follows from the corresponding result in recursive grouping. The Chow-Liu tree can be constructed with O(m2 log m) complexity. The recursive b grouping procedure has complexity maxr |Ar |3 and maxr |Ar | ≤ ∆(MST(V ; d)). A.6 Proof of Theorem 11: Consistency and Sample Complexity of Relaxed RG (i) Structural consistency follows from Theorem 5 and the fact that the ML estimates of information distances dbij approach dij (in probability) for all i, j ∈ V as the number of samples tends to infinity. Risk consistency for Gaussian and symmetric discrete distributions follows from structural consistency. If the structure is correctly recovered, we can use the equations in (13) and (14) to infer the information distances. Since the distances are in one-to-one correspondence to the correlation coefficients and the crossover probability for Gaussian and symmetric discrete distribution respectively, the parameters are also consistent. This implies that the KL-divergence between p and pbn tends to zero (in probability) as the number of samples n tends to infinity. This completes the proof. (ii) The theorem follows by using the assumption that the effective depth δ = δ(Tp ; V ) is constant. Recall that τ > 0 is the threshold used in relaxed RG (see (24) in Section 6.1). Let the set of triples (i, j, k) whose pairwise information distances are less than τ apart be J , i.e., (i, j, k) ∈ J if and only if max{dij , djk , dki } < τ . Since we assume that the true information distances are uniformly bounded, there exist τ > 0 and some sufficiently small b ijk − Φijk | ≤ λ for all (i, j, k) ∈ J , then RG recovers the correct latent λ > 0 so that if |Φ structure. b ijk − Φijk | > λ}. We note that the probability of the Define the error event Eijk := {|Φ event Eijk decays exponentially fast, i.e., there exists Jijk > 0 such that for all n ∈ N, Pr(Eijk ) ≤ exp(−nJijk ). (42) The proof of (42) follows readily for Chernoff bounds (Hoeffding, 1958) and is omitted. The error probability associated to structure learning can be bounded as follows:   (a) [ X (b) Eijk  ≤ Pr(Eijk ) (43) Pr h(Tbn ) 6= Tp ≤ Pr  (i,j,k)∈J ≤ m 3 (i,j,k)∈J (c) max Pr(Eijk ) ≤ exp(3 log m) exp −n min Jijk , (i,j,k)∈J (i,j,k)∈J (44) where (a) follows from the fact that if the event {h(Tbn ) 6= Tp } occurs, then there is at least one sibling or parent-child relationship that is incorrect, which corresponds to the union of b ijk differs from Φijk by the events Eijk , i.e., there exists a triple (i, j, k) ∈ J is such that Φ more than λ. Inequality (b) follows from the union bound and (c) follows from (42). Because the information distances are uniformly bounded, there also exists a constant Jmin > 0 (independent of m) such that min(i,j,k)∈J Jijk ≥ Jmin for all m ∈ N. Hence √ for every η > 0, if the number of samples satisfies n > 3(log(m/ 3 η))/Jmin , the error probability is bounded above by η. Let C := 3/Jmin to complete the proof of the sample complexity result in (29). The proof for the logarithmic sample complexity of distribution 38 Learning Latent Tree Graphical Models reconstruction for Gaussian and symmetric discrete models follows from the logarithmic sample complexity result for structure learning and the fact that the information distances are in a one-to-one correspondence with the correlation coefficients (for Gaussian models) or crossover probabilities (for symmetric discrete models). A.7 Proof of Theorem 12: Consistency and Sample Complexity of Relaxed CLRG (i) Structural consistency of CLGrouping follows from structural consistency of RG (or NJ) and the consistency of the Chow-Liu algorithm. Risk consistency of CLGrouping for Gaussian or symmetric distributions follows from the structural consistency, and the proof is similar to the proof of Theorem 11(i). (ii) The input to the CLGrouping procedure TbCL is the Chow-Liu tree and has O(log m) sample complexity Tan et al. (2010), where m is the size of the tree. From Theorem 11, the recursive grouping procedure has O(log m) sample complexity (for appropriately chosen thresholds) when the input information distances are uniformly bounded. In any iteration of the CLGrouping, the information distances satisfy dij ≤ γu, where γ, defined in (37), is the worst-case graph distance of any hidden node from its surrogate. Since γ satisfies (39), dij ≤ u2 δ/l. If the effective depth δ = O(1) (as assumed), the distances dij = O(1) and the sample complexity is O(log m). References K. Atteson. The performance of neighbor-joining methods of phylogenetic reconstruction. Algorithmica, 25(2):251–278, 1999. H.-J. Bandelth and A. Dress. Reconstructing the shape of a tree from observed dissimilarity data. Adv. Appl. Math, 7:309–43, 1986. S. Bhamidi, R. Rajagopal, and S. Roch. Network delay inference from additive metrics. To appear in Random Structures and Algorithms, Arxiv preprint math/0604367, 2009. R. Castro, M. Coates, G. Liang, R. Nowak, and B. Yu. Network Tomography: Recent Developments. Stat. Science, 19:499–517, 2004. J. T. Chang and J. A. Hartigan. Reconstruction of evolutionary trees from pairwise distributions on current species. In Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface, pages 254–257, 1991. T. Chen, N. L. Zhang, and Y. Wang. Efficient model evaluation in the search based approach to latent structure discovery. In 4th European Workshop on Probabilistic Graphical Models, 2008. M. J. Choi, J. J. Lim, A. Torralba, and A. S. Willsky. Exploiting hierarchical context on a large database of object categories. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), San Francisco, CA, June 2010. 39 Choi, Tan, Anandkumar, and Willsky C. K. Chow and C. N. Liu. Approximating discrete probability distributions with dependence trees. IEEE Trans. on Information Theory, 3:462–467, 1968. T. Cormen, C. Leiserson, R. Rivest, and C. Stein. Introduction to Algorithms. McGraw-Hill Science/Engineering/Math, 2nd edition, 2003. T. M. Cover and J. A. Thomas. Elements of Information Theory. Wiley-Interscience, 2nd edition, 2006. R. G. Cowell, A. P. Dawid, S. L. Lauritzen, and D. J. Spiegelhalter. Probabilistic networks and expert systems. Statistics for Engineering and Information Science. Springer-Verlag, New York, 1999. M. Csűrös. Reconstructing Phylogenies in Markov Models of Sequence Evolution. PhD thesis, Yale University, 2000. C. Daskalakis, E. Mossel, and S. Roch. Optimal phylogenetic reconstruction. In STOC ’06: Proceedings of the thirty-eighth annual ACM symposium on Theory of computing, pages 159–168, 2006. A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum-likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, 39:1–38, 1977. R. Durbin, S. R. Eddy, A. Krogh, and G. Mitchison. Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge Univ. Press, 1999. G. Elidan and N. Friedman. Learning hidden variable networks: The information bottleneck approach. Journal of Machine Learning Research, 6:81–127, 2005. P. L. Erdős, L. A. Székely, M. A. Steel, and T. J. Warnow. A few logs suffice to build (almost) all trees: Part ii. Theoretical Computer Science, 221:153–184, 1999. J. Farris. Estimating phylogenetic trees from distance matrices. Amer. Natur., 106(951): 645–67, 1972. S. Harmeling and C. K. I. Williams. Greedy learning of binary latent trees. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2010. W. Hoeffding. Probability inequalities for sums of bounded random variables. Journal of the American Statistical Association, 58:13–30, 1958. D. Hsu, S.M. Kakade, and T. Zhang. A Spectral Algorithm for Learning Hidden Markov Models. In Intl. Conf. on Learning Theory (COLT), 2009. T. Jiang, P. E. Kearney, and M. Li. A polynomial-time approximation scheme for inferring evolutionary trees from quartet topologies and its application. SIAM J. Comput., 30(6): 194261, 2001. C. Kemp and J. B. Tenenbaum. The discovery of structural form. Proceedings of the National Academy of Science, 105(31):10687–10692, 2008. 40 Learning Latent Tree Graphical Models J. B. Kruskal. On the Shortest Spanning Subtree of a Graph and the Traveling Salesman Problem. Proceedings of the American Mathematical Society, 7(1), Feb 1956. M. R. Lacey and J. T. Chang. A signal-to-noise analysis of phylogeny estimation by neighbor-joining: Insufficiency of polynomial length sequences. Mathematical Biosciences, 199:188–215, 2006. J. A. Lake. Reconstructing evolutionary trees from dna and protein sequences: Parallnear distances. Proceedings of the National Academy of Science, 91:1455–1459, 1994. S. L. Lauritzen. Graphical models. Clarendon Press, 1996. P. F. Lazarsfeld and N.W. Henry. Latent structure analysis. Boston: Houghton Mifflin, 1968. D. G. Luenberger. Introduction to Dynamic Systems: Theory, Models, and Applications. Wiley, 1979. D. Parikh and T. H. Chen. Hierarchical Semantics of Objects (hSOs). In ICCV, pages 1–8, 2007. J. Pearl. Probabilistic Reasoning in Intelligent Systems: Network of Plausible inference. Morgan Kaufmann, 1988. R. C. Prim. Shortest connection networks and some generalizations. Bell System Technical Journal, 36, 1957. D. F. Robinson and L. R. Foulds. Comparison of Phylogenetic Trees. Mathematical Biosciences, 53:131–147, 1981. S. Roch. A short proof that phylogenetic tree reconstruction by maximum likelihood is hard. IEEE/ACM Trans. Comput. Biol. Bioinformatics, 3(1), 2006. P. J. Rousseeuw. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Computational and Applied Mathematics, 20:53–65, 1987. N. Saitou and M. Nei. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol, 4(4):406–25, Jul 1987. S. Sattath and A. Tversky. Additive similarity trees. Psychometrika, 42:319–45, 1977. G. Schwarz. Estimating the dimension of a model. Annals of Statistics, 6:461–464, 1978. R. J. Serfling. Approximation Theorems of Mathematical Statistics. Wiley-Interscience, Nov 1980. S. Shen. Large deviation for the empirical correlation coefficient of two Gaussian random variables. Acta Mathematica Scientia, 27(4):821–828, Oct 2007. K. St. John, T. Warnow, B. M. E. Moret, and L. Vawter. Performance study of phylogenetic methods: (unweighted) quartet methods and neighbor-joining. J. Algorithms, 48(1):173– 193, 2003. 41 Choi, Tan, Anandkumar, and Willsky M. Steel. The complexity of reconstructing trees from qualitative characters and subtrees. Journal of Classification, 9:91–116, 1992. V. Y. F. Tan, A. Anandkumar, L. Tong, and A. S. Willsky. A Large-Deviation Analysis for the Maximum Likelihood Learning of Tree Structures. In Proceedings of IEEE International Symposium on Information Theory, pages 1140 – 1144, Seoul, Jul 2009. V. Y. F. Tan, A. Anandkumar, and A. S. Willsky. Learning Gaussian tree models: Analysis of error exponents and extremal structures. IEEE Transactions on Signal Processing, 58 (5):2701–2714, May 2010. Y. Tsang, M. Coates, and R. D. Nowak. Network Delay Tomography. IEEE Trans. Signal Processing, 51:21252136, 2003. N. L. Zhang. Hierarchical Latent Class Models for Cluster Analysis. Journal of Machine Learning Research, 5:697–723, 2004. N. L. Zhang and T Kočka. Efficient learning of hierarchical latent class models. In ICTAI, 2004. 42

RELATED PAPERS

RELATED TOPICS

Log In

Learning Latent Tree Graphical Models

Learning Latent Tree Graphical Models

Related Papers

RELATED PAPERS

RELATED TOPICS