# A Decomposition Algorithm for Learning Bayesian Networks Based on Scoring Function

Authors: Mingmin Zhu; Sanyang Liu
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Learning Bayesian network (BN) structure from data is a typical NP-hard problem. But almost existing algorithms have the very high complexity when the number of variables is large. In order to solve this problem(s), we present an algorithm that integrates with a decomposition-based approach and a scoring-function-based approach for learning BN structures. Firstly, the proposed algorithm decomposes the moral graph of BN into its maximal prime subgraphs. Then it orientates the local edges in each subgraph by the K2-scoring greedy searching. The last step is combining directed subgraphs to obtain final BN structure. The theoretical and experimental results show that our algorithm can efficiently and accurately identify complex network structures from small data set.


---

## Body

## 1. Introduction

Bayesian networks (BNs), also known as belief networks, are becoming a popular tool for representing uncertainty in artificial intelligence. They have been applied to a wide range of tasks such as natural spoken dialog systems, vision recognition, expert systems, medical diagnosis, and genetic regulatory network inference [1–5]. A BN consists of two important components: a directed acyclic graph (DAG) representing the dependency structure among the variables and a conditional probability table (CPT) for each variable given its parent set. There has been a lot of work in the last ten years on the learning of BNs for both graph structure and probability parameters. However, learning the structure is harder and, arguably, more critical [6–8]. Most of these algorithms can be grouped into two different categories: constraint-based methods and search-and-score methods. The algorithms based on constraint generate a list of conditional independence (CI) relationships among the variables in the domain and attempt to find a network that represents these relationships as far as possible [9–11]. The number, complexity, and reliability of the required independence tests are the main concerns regarding this type of algorithms. The algorithms based on the score-and-search paradigm see the learning task as a combinatorial optimization problem, where a search method operates on a search space associated with BNs. They evaluate the degree of fitness between each element and the available data by a scoring function [12–14]. Among the search-and-score algorithms, the K2 algorithm [15] is commonly used which learns the network structure from data by requiring a prior ordering of nodes as input. However, searching for the best structure is difficult because the search space increases exponentially with the number of variables [16].In this paper, we propose a BN structure learning algorithm, which combined with the merits of the constraint-based method and the K2 algorithm. The proposed algorithm not only employs constraint knowledge to decompose the search space but also uses the K2 score as heuristic knowledge to induce the process of local greedy search. It uses the property that the local information of variables in each maximal prime sub-graph cannot be destroyed by decomposing the undirected independence graph. At the same time, the K2 algorithm can be used to learn the local structure of each undirected subgraph and obtain a group of directed subgraphs. Thus, by combining these directed subgraphs, we obtain the final structure. Theoretical results and large number of experiments show that the new algorithm is effective and efficient, especially in small data set. Moreover, the proposed algorithm uses maximal prime decomposition to identify the whole graph structure, which reduces the search space significantly and greatly enhances learning speed.The remainder of this paper is organized as follows. Section2 introduces notation and definitions. We describe our algorithm and its theoretical proofs in Section 3. Section 4 discusses how to construct the moral graph from observed data. Simulation studies are conducted to demonstrate the performance of our algorithm and existing algorithms in Section 5. Finally, in Section 6, we conclude and outline our future work. All proofs will be presented in the Appendix section.


## 2. Notation and Definitions

In this section, we provide some basic technical terminologies that are sufficient for understanding this paper.A BN is a tuple(G,P), where G=(V,E) is a directed acyclic graph (DAG) with nodes representing the random variables in V and P a joint probability distribution on V. A node Vi is called a parent of Vj if the directed edge Vi→Vj∈E. The set of all parents of Vj is denoted as pa(Vj). A path from Vi to Vj is a sequence [Vi=Vi0,Vi1,…,Vin=Vj] of distinct nodes such that Vi(k-1)→Vik or Vik→Vi(k-1) for k=1,2,…,n. We say that Vj is an ancestor of Vi and Vi is a descendant of Vj if there is a directed path from Vj to Vi in G. In addition, G and P must satisfy the Markov condition: every variable Vi∈V is independent of any subset of its nondescendant variables conditioned on the set of its parents, denoted by pa(Vi). We denote the conditional independence of the variable sets X and Y given Z in some distribution P with IndP(X;Y∣Z). A path ρ is said to be d-separated by a set Z in a DAG G, if and only if (1) ρ contains a “head-to-tail meeting”: Vi→Vj→Vk or a “tail-to-tail meeting”: Vi←Vj→Vk such that the middle node Vj is in Z, or (2) ρ contains a “head-to-head meeting”: Vi→Vj←Vk such that the middle node Vj is not in Z and no descendant of Vj is in Z. Two distinct sets X and Y of nodes are said to be d-separated by a set Z in G if Z d-separates every path from any node in X to any node in Y. We use DsepG(X;Y∣Z) to denote the assertion that X is d-separated from Y given Z in G.For a DAGG=(V,E), its moral graph Gm=(V,Em) is an undirected graph obtained by connecting every pair of nodes with a common child which are not already connected in G and then dropping the directionality of all directed edges. An undirected path ρ is said to be separated by a set Z in a moral graph Gm, if and only if ρ passes through Z. A BN (G,P) satisfies the faithfulness condition if the d-separations in G identify all and only the conditional independencies in P, that is, IndP(X;Y∣Z) if and only if DsepG(X;Y∣Z). We will drop the subscript G or P in the notation of conditional independence when the faithfulness condition holds. Let G(V) denote a graph consisting of a finite set V of nodes.Definition 2.1 (see [17]).
Suppose the triplet(V′,S,V′′) denote a partition of V where V=V′∪S∪V′′; if every path in Gm(V) between Vi∈V′ and Vj∈V′′ contains a node in S and Gm(S) is complete, then Gm is decomposable and (V′,S,V′′) is a decomposition of Gm into subgraphs Gm(V′∪S) and Gm(V′′∪S); otherwise Gm is prime. Moreover, S is called a complete separator of Gm with respect to (V′,S,V′′).Furthermore, an undirected graphGm is said to be decomposed maximally if Gm and all its sub-graphs can be decomposed recursively until all sub-graphs are not decomposable. These sub-graphs are the maximal prime sub-graphs of Gm. Transformation of a BN into its maximal prime sub-graphs is equivalent to recursively decomposing the moral graph of DAG. The most well-known algorithm is the MPD-JT algorithm [17] which first triangulates the moral graph Gm by adding a fill-in edge to every cycle whose length is greater than three and then identifies all the cliques of triangulation graph and arranges them as a junction tree (JT), and finally, recursively aggregates cliques connected by incomplete separators (in the moral graph). Figure 2 schematically illustrates the process of the MPD-JT algorithm.Example 2.2.
Consider the Asia network [18] in Figure 1(a). Its moral graph Gm is shown in Figure 1(b). By our definition, Gm is decomposable. Figure 1(c) provides one of its triangulation graphs of Gm and Figure 2(a) gives its corresponding junction tree. Since  Gm(BL) is an incomplete separator in the moral graph Gm, we aggregate cliques BLS and BEL according to the MPD-JT algorithm. The resulting graph is shown in Figure 2(b). From Figure 2(b), we can see that all the separators are complete in Gm, so, there are no cliques that need to be aggregated. We obtain five maximal prime sub-graphs (Figure 2(c)).Asia network, its moral graph, and one of its triangulation graphs.
(a)
(b)
(c)Figure 2
The process of decomposing the Asia moral graph into its maximal prime subgraphs.The K2 algorithm is a greedy search algorithm that learns the network structure from data. It attempts to select the network structure which maximizes the network's posterior probability, given the observed data. Its scoring function can be expressed as follows:(2.1)fK2(G,D)=logp(G)+∑i=1n‍[∑j=1qi‍[log((ri-1)!(Nij+ri-1)!)+∑k=1ri‍log(Nijk!)]].Nijk is the number of cases in D where Vi is in its kth state and its parents are in their jth state. Nij=∑k=1riNijk denotes the number of cases in D where pa(Vi) is in its jth state, ri denotes the number of states of variable Vi, and qi is the number of parent configurations of Vi. It is obvious that the search space is the main element influencing the performance of the K2 algorithm. In theory, the more constraint conditions, the smaller the search space of BNs. And then the search efficiency will be much higher. Hence, it is very necessary to reduce the search space. Considering the BN's own characteristic, we combine the constraint approach with maximal prime decomposition to learn a BN structure.


## 3. The Structure Learning Algorithm and Its Correctness

In this section, theoretical results are presented for learning the structure of BN. We show how the problem of learning the structure over the full set of all variables can be split into its subproblems. Below we first give two theorems based on which we propose the decomposition algorithm for structural learning of BNs.Theorem 3.1.
LetGm be a moral graph for a DAG G=(V,E). Then a subset of variables S⊂V d-separates Vi∈V from Vj∈V in G if and only if S separates Vi from Vj in Gm.By Theorem3.1, we can know that the condition of d-separation in a DAG is equivalent to the separation in its moral graph.Example 3.2.
Consider the Asia network in Figure1(a). By the definition of d-separation, we can see that A and X are d-separated by E, because the path A→T→E→X is blocked at E. Since there is a path A→T→E←L←S→B→D and E is a head-to-head node, A and D are not d-separated by E. Thus, we can conclude that in the moral graph of Asia network in Figure 1(b), A and X are separated by E and E does not separate A from D. In fact, it is obviously held in Figure 1(b).Theorem 3.3.
LetGm be a moral graph for a DAG G=(V,E) and let G1m,G2m,…,Gkm be maximal prime sub-graphs of Gm. For any variable Vi∈V, there exists at least one maximal prime sub-graph of Gm which contains Vi and Pa(Vi).A consequence of Theorem3.3 is that a problem of structural learning from a data set  can be decomposed into problems of local structural learning from a subdata set. The guarantee of such approach is that each subdata must contain sufficient data information of variable  and its parent set.Based on the previous analysis, an improved algorithm was proposed combined with the merits of previous two basic methods for BN structure learning. First, according to the observed data or domain knowledge, we construct the independent graph (the moral map) of the target DAG applying constraint-based approach and then decompose the independent graph. It is shown that each maximal prime sub-graph contains the sufficient information of local variables, so, the search space of the score function can be effectively reduced. Second, the local structure of each sub-graph is learned by using the score function and a directed acyclic graph is obtained for each sub-graph. Finally, we combine all these directed acyclic sub-graphs. Theoretical and experimental results show that the new algorithm is effective and reasonable. Now we formalize our algorithm in Algorithm1.Algorithm 1:Decomposition algorithm for learning BNs.
(1)Input: Data set D; Variable set V=(V1,V2,…,Vn); Node ordering ρ.
(2) Construct the moral graphGm=(V,Em) from the data set D.
(3) Decompose the moral graphGm into its maximal prime subgraphs G1m,G2m,…,Gkm.
(4) For each sub-graphGlm(l=1,2,…,k) of Gm, call the K2 algorithm with the local node
ordering ofGlm to construct a directed acyclic graph Gl(l=1,2,…,k).
(5) CombineG1,G2,…,Gk into a directed graph G=(V,E) where
V=V(G1)∪V(G2)∪⋯∪V(Gk) and E=E(G1)∪E(G2)∪⋯∪E(Gk).
(6)Output: A directed acyclic graph G.As shown in Algorithm1, the proposed algorithm first decomposed the entire variable set into its subsets. Then the local directed graph of each subset is recovered by K2 algorithm. Unlike the case of the X-G-Z algorithm [19], the final result of our algorithm is a DAG, not a partial directed graph, and the procedure of finding minimal d-separators is avoided. Furthermore, the computational complexity of Algorithm 1 is less than that of the K2 algorithm. In fact, the triangulation used to construct a junction tree from an undirected independence graph is the main cost of the MPD-JT algorithm. Although the problem of optimally triangulating an undirected graph is NP hard, suboptimal triangulation methods may be used provided that the obtained tree does not contain too large cliques. Two most well-known algorithms are the lexicographic search [20] and the maximum cardinality search [21], and their complexities are O(ne) and O(n+e), respectively, where n is the number of nodes and e is the number of edges in the graph. Thus, decomposition of  Gm is a computationally simple task compared to the exhaustive search using the K2 score. Let n denote the number of variables in V, and m the number of cases in data set D. In the worst case, the complexity of the K2 algorithm is O(mn4r)  [15] where r is the largest number of states of a variable in V. Suppose that Gm is decomposed into k subgraphs, where k≤n. Let n0 denote the number of variables in the largest subgraph, that is, n0=maxl|Glm|, where |Glm| denotes the number of variables in Glm. The complexity for constructing a directed acyclic subgraph Gl in Algorithm 1 is O(mn04r), and thus that of all directed subgraphs is O(kmn04r). Since n0 is usually much less than n, our algorithm is less computationally complex than the K2 algorithm. We now establish the correctness of our algorithm by showing that our final result is a DAG.Theorem 3.4.
Given a data setD and a variable (node) set V, then the graph G=(V,E) returned by Algorithm 1 is a DAG.


## 4. Construction of the Moral Graph from Observed Data

In general, the moral graph of an underlying BN can be obtained from the observed data based on the conditional independence tests. An edgeVi-Vj is included in the moral graph if and only if  Vi and Vj are not conditionally independent given all other variables. However, these require sufficient data for estimating the parameters and for improving the power of tests. To avoid testing high order conditional independence relations, we propose a Markov boundary discovery algorithm which is based on a subroutine MMPC [22]. The Markov boundary of a variable T, denoted as MB(T), is a minimal set of variables conditioned on which all other variables that are probabilistically independent of the target T. Furthermore, the set of parents, children, and spouses of T is its unique Markov boundary under the faithfulness condition [23]. The MMPC algorithm is sketched in the appendix section.According to Algorithm2, we can see that our algorithm starts with a two-phase approach. A shrinking phase attempts to remove the most irrelevant variables to T, followed by a growing phase that attempts to add as many dependent variables as possible. The growing phase is interleaved with the shrinking phase. Interleaving the two phases allows to eliminate some of the false positives in the current boundary as the algorithm progresses during the growing phase. Theoretically, the more the constraint knowledge obtained by CI tests, the smaller the search space, and the higher the searching efficiency. However, the results of higher-order CI tests may be unreliable. Steps 3 and 4 of Algorithm 2 only use order-0 and order-1 CI tests to reduce the search space whose number of CI tests is bound by O(n3). In Steps 6 and 7 of Algorithm 2, we only condition on subsets of sizes up to one instead of conditioning on all subsets of the PC¯∪SP¯. The order of the complexity is O(t·|PC|2) where PC is the largest set of parents and children over all variables in V and t is the number of variables in PC¯∪SP¯, t≪n. Thus, the total complexity of Algorithm 2 is O(n·t·|PC|2).Algorithm 2:Algorithm for discovering Markov boundary.
(1)Input: Data set D; Variable set V=(V1,V2,…,Vn); Target variable T.
(2) Initialization:PC¯=V∖T,SP¯=∅.
(3) Order-0 CI test: for each variableVi∈PC¯, if Ind(Vi;T) is hold, then PC¯=PC¯∖Vi,S(Vi)=∅.
(4) Order-1 CI test: for each variableVi∈PC¯, if there is a variable Vj∈PC¯∖Vi such that
Ind(Vi;T∣Vj), then PC¯=PC¯∖Vi,S(Vi)=S(Vi)∪Vj.
(5) Find superset of spouses: for each variableVi∈PC¯, if there is a variable
Vj∈V∖{T∪PC¯}, such that ¬Ind(Vj;T∣S(Vj)∪Vi), then SP¯=SP¯∪Vj.
(6) Find parents and children ofT: call the MMPC algorithm to get
PC(T)=MMPC(T,D(PC¯∪SP¯)). For each Vi∈PC¯∖PC(T), if T∈MMPC(Vi,D),
thenPC(T)=PC(T)∪Vi.
(7) Find spouses ofT: for each variable Vi∈PC(T), if there is a variable
Vj∈MMPC(Vi,D)∖{PC(T)∪T} and a subset Z⊂PC¯∪SP¯∖{T∪Vj}, such that
Ind(T;Vj∣Z) and ¬Ind(T;Vj∣Z∪Vi), then SP(T)=SP(T)∪Vj.
(8) ReturnMB(T)=PC(T)∪SP(T).
(9)Output: A Markov boundary MB(T) of T.Theorem 4.1.
Suppose  (G,P) satisfies the faithfulness condition, where  G=(V,E) is a DAG and P is a joint probability distribution of the variables in V. For each target variable T, Algorithm 2 returns the Markov boundary of  T.By Theorem4.1, a moral graph can be constructed from observed data.On the other hand,Gm can be constructed based on the prior or domain knowledge rather than conditional independence tests. The domain knowledge may be experts' prior knowledge of dependencies among variables, such as Markov chains, chain graphical models, and dynamic or temporal models. Based on the domain knowledge of dependencies, data patterns of databases can be represented as a collection of variable sets D={D1,D2,…,Dk}, in which variables contained in the same set may associate each other directly but variables contained in different sets associate each other through other variables. This means that two variables that are not contained in the same set are independent conditionally of all other variables. From the k data patterns, we get separately k undirected subgraphs. Combining them together, we obtain the undirected graph Gm.


## 5. Experimental Results

In this section, we present the experimental results carried out with our algorithm on two standard network data sets (Alarm [24] and Insurance [25]). The first one is a medical diagnostic system for patient monitoring. It consists of 37 nodes and 46 edges. The random variables in the Alarm network are discrete in nature and can take two, three, or four states (Figure 3). The second example is a network for estimating the expected claim costs for a car Insurance policyholder. It consists of 27 nodes and 52 edges connecting them. Our implementation is based on the Bayesian network toolbox written by Murphy [26] and the Causal Explorer System developed by Aliferis et al. [27]. The experimental platform was a personal computer with Pentium 4, 3.06 GHz CPU, 0.99 GB memory, and Windows XP.Figure 3
Alarm network.


### 5.1. Learning Networks with Node Ordering

In this subsection, we show simulations of the Alarm network when the node ordering is known. Although our algorithm is based on a constraint-based method and a search-and-score method, given that the result of network returned by Algorithm1 is a DAG, which has some similarities with the search-and-score methods, it is interesting to include one of these methods in comparison. We have selected the well-known K2 algorithm. We begin with a BN G0 which is completely specified in terms of structure and parameters, and we obtain a data set of a given size by sampling from G0. Then, using our algorithm and the K2 algorithm, respectively, we obtain two learned networks G1 and G2, which must be compared with the original network G0. More precisely, the experiments are carried out on different data set randomly selected from the database (100000 data points). The size of the data set is varied from 500 to 10000, and 10 replicates are done for each of the different network parameters and sample sizes. The average number of missing edges, extra edges, and reversed edges in the learned networks with respect to the original one is computed.We compare our methods with two different significance levels on Alarm network, that is,α=0.01 and. The results are shown in Table 1 for different sample sizes 500, 1000, 2000, 5000, and 10000. In the second row of table, four values in a bracket denote the number of missing edges, extra edges, reversed edges, and computation time, respectively, and other rows give values relative to the second row, which are obtained by dividing their real values by the values in the second row. A relative value larger than 1 denotes that its real value is larger than the corresponding value in the second row. From Table 1, the first thing that can be observed is that these results seem to confirm our intuition about the need to use Algorithm 1 with a smaller significance level α than those typically used for independence tests, since Algorithm 1 with the value α=0.01 offers better results than that with α=0.05. It also shows that the structure obtained by our algorithm (α=0.01) has the least number of missing edges, extra edges, and reversed edges. And further, our algorithm costs the least time.Table 1
Results relative to Algorithm1 with α=0.01 and α=0.05: missing edges, extra edges, reversed edges, and computation time.
Algorithm (levelα)
n
=
500
n
=
1000
n
=
2000
n
=
5000
n
=
10000
(2.4, 3.1, 1.2, 6.9951)
(1.6, 3.0, 1.1, 7.6500)
(1.2, 3.0, 1.0, 9.9092)
(0.8, 2.9, 0.9, 12.7819)
(0.6, 2.7, 0.7, 18.8523)
Algorithm1 (0.01)
(1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0000)
Algorithm1 (0.05)
(0.9, 2.0, 1.3, 0.9199)
(0.9, 1.8, 1.3, 0.9330)
(1.0, 1.4, 1.2, 1.1031)
(1.1, 1.1, 1.2, 1.1730)
(1.2, 1.0, 1.1, 1.1883)
K2 algorithm
(1.4, 2.2, 1.3, 1.5819)
(1.3, 1.9, 1.2, 1.7235)
(1.2, 1.5, 1.2, 2.1592)
(1.2, 1.4, 1.1, 2.7847)
(1.0, 1.2, 1.0, 2.9670)Table2 displays the K2 scores obtained for the original Alarm network, network returned by Algorithm 1 and network returned by K2 algorithm, respectively. It is easy to see that the larger the dataset size is, the closer the score to the original one. However, when the dataset size is smaller than 5000, the score returned by our algorithm is more approximate to the original one. At the same time, it can be found that our algorithm performs better than the K2 algorithm in terms of the running results on all data sets. Moreover, the advantage is very obvious when the data set is small; namely, the smaller the data set size, the more obvious the improvement.Table 2
Scores on Alarm for different dataset sizen: score of original network Sc (G0), score of network returned by Algorithm 1 Sc (G1), and score of network returned by K2 algorithm Sc (G2).
Score
n
=
500
n
=
1000
n
=
2000
n
=
5000
n
=
10000
Sc (G0)
-
5.2004
e
+
003
-
1.0000
e
+
004
-
1.9440
e
+
004
-
4.7737
e
+
004
-
9.4515
e
+
004
Sc (G1)
-
4.9831
e
+
003
-
1.0652
e
+
004
-
2.0502
e
+
004
-
4.7558
e
+
004
-
9.4514
e
+
004
Sc (G2)
-
6.8038
e
+
003
-
1.9824
e
+
004
-
3.0308
e
+
004
-
4.7826
e
+
004
-
9.4533
e
+
004


### 5.2. Learning Networks without Node Ordering

From the previous experimental results, we can see that our algorithm is capable of identifying structures that are close to optimal ones, given a prior ordering of the variables. However, this ordering information may not be available in real-life applications. Thus, in this subsection, Algorithm1 is extended to manage the problem in which the node ordering is not available. In fact, we only need to delete the directed cycles after combining all directed subgraphs at the expense of introducing a little complexity.Now we compare our extended method with the X-G-Z [19], RAI [28], and MMHC [22] algorithms. Similar with Section 5.1, the size of the data set is varied from 1000 to 10000, and 10 replicates are done for each of the different network parameters and sample sizes. We use two significance levels (α=0.01 and α=0.05) in the simulations. Unlike the situation of Section 5.1, algorithms that return a DAG are converted to the corresponding partial directed acyclic graph (PDAG) before comparing the quality of reconstruction. The PDAG is equivalent to its DAG. The average number of missing edges, extra edges, and reversed edges in the learned PDAG with respect to the underlying PDAG is computed.We summarize the simulation results in Tables3 and 4. In the second row of tables, five numbers in a bracket denote the number of missing edges, extra edges, reversed edges, the sum of the first three values, and computation time, respectively. Other rows give values relative to the second row, which are obtained by dividing their real values by the values in the second row. A relative value larger than 1 denotes that its real value is larger than the corresponding value in the second row. In each column, the best of the eight results are displayed in bold. From a general point of view, we can see that the X-G-Z algorithm obtains the least number of missing edges, the MMHC obtains the least number of extra edges, and Algorithm 1 has the least number of reversed edges. In terms of the sum of the three kinds of edges, the RAI and MMHC algorithms perform better than the X-G-Z algorithm and Algorithm 1 performs best. Although it can be seen that Algorithm 1 seems to have a similar performance with RAI in most cases, the running time on all data sets cost for Algorithm 1 is least. Moreover, the advantage is very obvious when the data set is large; namely, the bigger the sample size, the more obvious the improvement. The main reason is that Algorithm 1 uses lower-order CI tests and employs maximal prime decomposition technique to effectively decompose the search space which cut down many computings of statistic factors, scorings of the structures, and comparisons of the solutions and thus greatly enhances the time performance. Contrast Table 4 with Table 3, since the Insurance network has more variables than the Alarm network, it can be found that the running time of X-G-Z fast increases as the number of variable increases. However, Algorithm 1 is not sensitive to the increase of the variable capacity. In conclusion, our algorithm has a better overall performance compared to the other state-of-the-art algorithms.Table 3
Results relative to Algorithm1: missing edges, extra edges, reversed edges, and computation time for Alarm network.
Algorithm (levelα)
n
=
1000
n
=
2000
n
=
5000
n
=
10000
(10.0, 1.3, 10.2, 21.5, 8.9908)
(6.6, 1.2, 8.9, 16.74, 17.4378)
(4.0, 0.8, 5.4, 10.2, 28.9233)
(2.6, 0.7, 2.7, 6.0, 41.4688)
Algorithm1 (0.01)
(1.0, 1.0,1.0, 1.0, 1.0000)
(1.0, 1.0,1.0, 1.0, 1.0000)
(1.0, 1.0,1.0, 1.0, 1.0000)
(1.0, 1.0,1.0, 1.0, 1.0000)
Algorithm1 (0.05)
(0.8, 4.0, 1.8, 1.4, 1.0321)
(0.8, 3.6, 1.5, 1.3, 0.9389)
(0.8, 4.8, 1.6, 1.5, 1.1315)
(0.8, 4.9, 1.8, 1.7, 1.0356)
RAI (0.01)
(0.9, 1.6, 1.2, 1.1, 1.2879)
(1.0, 1.5, 1.1, 1.1, 1.6792)
(1.0, 1.8, 1.1, 1.1, 2.6955)
(1.0, 1.8, 1.2, 1.2, 4.4231)
RAI (0.05)
(0.8, 3.1, 1.5, 1.3, 1.4489)
(0.9, 2.8, 1.3, 1.3, 1.6977)
(0.9, 3.7, 1.4, 1.4, 3.6929)
(0.9, 3.7, 1.6, 1.5, 4.4982)
X-G-Z (0.01)
(0.8, 2.0, 1.1, 1.0, 5.0260)
(0.8, 2.4, 1.1, 1.1, 2.3356)
(0.9, 3.6, 1.1, 1.2, 1.9726)
(0.9, 4.3, 1.9, 1.7, 1.6813)
X-G-Z (0.05)
(0.8, 2.2, 1.0, 1.0, 3.5679)
(0.9, 2.7,1.0, 1.1, 3.8967)
(0.8, 3.8, 1.2, 1.2, 3.5876)
(0.7, 4.2, 1.3, 1.4, 1.8952)
MMHC (0.01)
(1.3,0.3, 1.8, 1.5, 1.4798)
(1.4,0.3, 1.2, 1.2, 2.1321)
(1.5,0.4, 1.2, 1.2, 3.8258)
(1.5,0.7, 1.4, 1.4, 6.7052)
MMHC (0.05)
(1.2, 0.5, 1.5, 1.3, 1.7269)
(1.3, 0.4, 1.1, 1.1, 2.2036)
(1.3, 0.5, 1.2, 1.3, 5.4004)
(1.3, 1.0, 1.1, 1.2, 6.8066)Table 4
Results relative to Algorithm1: missing edges, extra edges, reversed edges, and computation time for Insurance  network.
Algorithm (levelα)
n
=
1000
n
=
2000
n
=
5000
n
=
10000
(13.0, 2.4, 20.6, 36.0, 2.7031)
(10.0, 1.5, 21.5, 33, 4.1406)
(7.4, 1.3, 16.3, 25.0, 6.9063)
(6.7, 1.1, 12.2, 20.0, 11.6094)
Algorithm1 (0.01)
(1.0, 1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0, 1.0000)
Algorithm1 (0.05)
(1.2, 1.5, 1.0, 1.1, 1.1850)
(1.1, 1.3,1.0, 1.0, 1.0981)
(1.0, 1.0, 1.3, 1.1, 1.0158)
(1.0, 1.1, 1.1, 1.1, 1.0202)
RAI (0.01)
(1.0, 1.1, 1.0, 1.1, 1.2879)
(1.0, 1.1, 1.0, 1.0, 1.6792)
(1.0, 1.1, 1.1, 1.1, 2.6955)
(1.0, 1.0, 1.0, 1.1, 4.4231)
RAI (0.05)
(1.1, 1.3, 1.0, 1.1, 1.6636)
(1.1, 1.2,1.0, 1.0, 1.9857)
(1.0, 1.1, 1.2, 1.2, 3.3154)
(1.0, 1.1, 1.1, 1.1, 5.2848)
X-G-Z (0.01)
(1.0, 1.6, 1.0, 1.0, 7.2486)
(1.2, 2.3, 1.2, 1.3, 6.0151)
(1.3, 2.5, 1.3, 1.4, 7.8868)
(1.4, 2.9, 1.7, 1.7, 8.5505)
X-G-Z (0.05)
(1.1, 1.6,0.9, 1.0, 9.9249)
(1.2, 2.4, 1.1, 1.2, 8.1661)
(1.3, 2.7, 1.4, 1.4, 16.0270)
(1.3, 3.0, 1.9, 1.8, 11.6998)
MMHC (0.01)
(1.2,0.6, 1.1, 1.1, 1.4798)
(1.4,1.0, 1.1, 1.2, 2.1321)
(1.4, 1.1, 1.2, 1.2, 3.8258)
(1.3,1.0, 1.2, 1.2, 6.7052)
MMHC (0.05)
(1.3,0.6, 1.1, 1.2, 1.9827)
(1.4, 1.1, 1.2, 1.3, 2.5774)
(1.5, 1.2, 1.2, 1.3, 4.8484)
(1.4, 1.1, 1.2, 1.3, 8.1278)


## 5.1. Learning Networks with Node Ordering

In this subsection, we show simulations of the Alarm network when the node ordering is known. Although our algorithm is based on a constraint-based method and a search-and-score method, given that the result of network returned by Algorithm1 is a DAG, which has some similarities with the search-and-score methods, it is interesting to include one of these methods in comparison. We have selected the well-known K2 algorithm. We begin with a BN G0 which is completely specified in terms of structure and parameters, and we obtain a data set of a given size by sampling from G0. Then, using our algorithm and the K2 algorithm, respectively, we obtain two learned networks G1 and G2, which must be compared with the original network G0. More precisely, the experiments are carried out on different data set randomly selected from the database (100000 data points). The size of the data set is varied from 500 to 10000, and 10 replicates are done for each of the different network parameters and sample sizes. The average number of missing edges, extra edges, and reversed edges in the learned networks with respect to the original one is computed.We compare our methods with two different significance levels on Alarm network, that is,α=0.01 and. The results are shown in Table 1 for different sample sizes 500, 1000, 2000, 5000, and 10000. In the second row of table, four values in a bracket denote the number of missing edges, extra edges, reversed edges, and computation time, respectively, and other rows give values relative to the second row, which are obtained by dividing their real values by the values in the second row. A relative value larger than 1 denotes that its real value is larger than the corresponding value in the second row. From Table 1, the first thing that can be observed is that these results seem to confirm our intuition about the need to use Algorithm 1 with a smaller significance level α than those typically used for independence tests, since Algorithm 1 with the value α=0.01 offers better results than that with α=0.05. It also shows that the structure obtained by our algorithm (α=0.01) has the least number of missing edges, extra edges, and reversed edges. And further, our algorithm costs the least time.Table 1
Results relative to Algorithm1 with α=0.01 and α=0.05: missing edges, extra edges, reversed edges, and computation time.
Algorithm (levelα)
n
=
500
n
=
1000
n
=
2000
n
=
5000
n
=
10000
(2.4, 3.1, 1.2, 6.9951)
(1.6, 3.0, 1.1, 7.6500)
(1.2, 3.0, 1.0, 9.9092)
(0.8, 2.9, 0.9, 12.7819)
(0.6, 2.7, 0.7, 18.8523)
Algorithm1 (0.01)
(1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0000)
Algorithm1 (0.05)
(0.9, 2.0, 1.3, 0.9199)
(0.9, 1.8, 1.3, 0.9330)
(1.0, 1.4, 1.2, 1.1031)
(1.1, 1.1, 1.2, 1.1730)
(1.2, 1.0, 1.1, 1.1883)
K2 algorithm
(1.4, 2.2, 1.3, 1.5819)
(1.3, 1.9, 1.2, 1.7235)
(1.2, 1.5, 1.2, 2.1592)
(1.2, 1.4, 1.1, 2.7847)
(1.0, 1.2, 1.0, 2.9670)Table2 displays the K2 scores obtained for the original Alarm network, network returned by Algorithm 1 and network returned by K2 algorithm, respectively. It is easy to see that the larger the dataset size is, the closer the score to the original one. However, when the dataset size is smaller than 5000, the score returned by our algorithm is more approximate to the original one. At the same time, it can be found that our algorithm performs better than the K2 algorithm in terms of the running results on all data sets. Moreover, the advantage is very obvious when the data set is small; namely, the smaller the data set size, the more obvious the improvement.Table 2
Scores on Alarm for different dataset sizen: score of original network Sc (G0), score of network returned by Algorithm 1 Sc (G1), and score of network returned by K2 algorithm Sc (G2).
Score
n
=
500
n
=
1000
n
=
2000
n
=
5000
n
=
10000
Sc (G0)
-
5.2004
e
+
003
-
1.0000
e
+
004
-
1.9440
e
+
004
-
4.7737
e
+
004
-
9.4515
e
+
004
Sc (G1)
-
4.9831
e
+
003
-
1.0652
e
+
004
-
2.0502
e
+
004
-
4.7558
e
+
004
-
9.4514
e
+
004
Sc (G2)
-
6.8038
e
+
003
-
1.9824
e
+
004
-
3.0308
e
+
004
-
4.7826
e
+
004
-
9.4533
e
+
004


## 5.2. Learning Networks without Node Ordering

From the previous experimental results, we can see that our algorithm is capable of identifying structures that are close to optimal ones, given a prior ordering of the variables. However, this ordering information may not be available in real-life applications. Thus, in this subsection, Algorithm1 is extended to manage the problem in which the node ordering is not available. In fact, we only need to delete the directed cycles after combining all directed subgraphs at the expense of introducing a little complexity.Now we compare our extended method with the X-G-Z [19], RAI [28], and MMHC [22] algorithms. Similar with Section 5.1, the size of the data set is varied from 1000 to 10000, and 10 replicates are done for each of the different network parameters and sample sizes. We use two significance levels (α=0.01 and α=0.05) in the simulations. Unlike the situation of Section 5.1, algorithms that return a DAG are converted to the corresponding partial directed acyclic graph (PDAG) before comparing the quality of reconstruction. The PDAG is equivalent to its DAG. The average number of missing edges, extra edges, and reversed edges in the learned PDAG with respect to the underlying PDAG is computed.We summarize the simulation results in Tables3 and 4. In the second row of tables, five numbers in a bracket denote the number of missing edges, extra edges, reversed edges, the sum of the first three values, and computation time, respectively. Other rows give values relative to the second row, which are obtained by dividing their real values by the values in the second row. A relative value larger than 1 denotes that its real value is larger than the corresponding value in the second row. In each column, the best of the eight results are displayed in bold. From a general point of view, we can see that the X-G-Z algorithm obtains the least number of missing edges, the MMHC obtains the least number of extra edges, and Algorithm 1 has the least number of reversed edges. In terms of the sum of the three kinds of edges, the RAI and MMHC algorithms perform better than the X-G-Z algorithm and Algorithm 1 performs best. Although it can be seen that Algorithm 1 seems to have a similar performance with RAI in most cases, the running time on all data sets cost for Algorithm 1 is least. Moreover, the advantage is very obvious when the data set is large; namely, the bigger the sample size, the more obvious the improvement. The main reason is that Algorithm 1 uses lower-order CI tests and employs maximal prime decomposition technique to effectively decompose the search space which cut down many computings of statistic factors, scorings of the structures, and comparisons of the solutions and thus greatly enhances the time performance. Contrast Table 4 with Table 3, since the Insurance network has more variables than the Alarm network, it can be found that the running time of X-G-Z fast increases as the number of variable increases. However, Algorithm 1 is not sensitive to the increase of the variable capacity. In conclusion, our algorithm has a better overall performance compared to the other state-of-the-art algorithms.Table 3
Results relative to Algorithm1: missing edges, extra edges, reversed edges, and computation time for Alarm network.
Algorithm (levelα)
n
=
1000
n
=
2000
n
=
5000
n
=
10000
(10.0, 1.3, 10.2, 21.5, 8.9908)
(6.6, 1.2, 8.9, 16.74, 17.4378)
(4.0, 0.8, 5.4, 10.2, 28.9233)
(2.6, 0.7, 2.7, 6.0, 41.4688)
Algorithm1 (0.01)
(1.0, 1.0,1.0, 1.0, 1.0000)
(1.0, 1.0,1.0, 1.0, 1.0000)
(1.0, 1.0,1.0, 1.0, 1.0000)
(1.0, 1.0,1.0, 1.0, 1.0000)
Algorithm1 (0.05)
(0.8, 4.0, 1.8, 1.4, 1.0321)
(0.8, 3.6, 1.5, 1.3, 0.9389)
(0.8, 4.8, 1.6, 1.5, 1.1315)
(0.8, 4.9, 1.8, 1.7, 1.0356)
RAI (0.01)
(0.9, 1.6, 1.2, 1.1, 1.2879)
(1.0, 1.5, 1.1, 1.1, 1.6792)
(1.0, 1.8, 1.1, 1.1, 2.6955)
(1.0, 1.8, 1.2, 1.2, 4.4231)
RAI (0.05)
(0.8, 3.1, 1.5, 1.3, 1.4489)
(0.9, 2.8, 1.3, 1.3, 1.6977)
(0.9, 3.7, 1.4, 1.4, 3.6929)
(0.9, 3.7, 1.6, 1.5, 4.4982)
X-G-Z (0.01)
(0.8, 2.0, 1.1, 1.0, 5.0260)
(0.8, 2.4, 1.1, 1.1, 2.3356)
(0.9, 3.6, 1.1, 1.2, 1.9726)
(0.9, 4.3, 1.9, 1.7, 1.6813)
X-G-Z (0.05)
(0.8, 2.2, 1.0, 1.0, 3.5679)
(0.9, 2.7,1.0, 1.1, 3.8967)
(0.8, 3.8, 1.2, 1.2, 3.5876)
(0.7, 4.2, 1.3, 1.4, 1.8952)
MMHC (0.01)
(1.3,0.3, 1.8, 1.5, 1.4798)
(1.4,0.3, 1.2, 1.2, 2.1321)
(1.5,0.4, 1.2, 1.2, 3.8258)
(1.5,0.7, 1.4, 1.4, 6.7052)
MMHC (0.05)
(1.2, 0.5, 1.5, 1.3, 1.7269)
(1.3, 0.4, 1.1, 1.1, 2.2036)
(1.3, 0.5, 1.2, 1.3, 5.4004)
(1.3, 1.0, 1.1, 1.2, 6.8066)Table 4
Results relative to Algorithm1: missing edges, extra edges, reversed edges, and computation time for Insurance  network.
Algorithm (levelα)
n
=
1000
n
=
2000
n
=
5000
n
=
10000
(13.0, 2.4, 20.6, 36.0, 2.7031)
(10.0, 1.5, 21.5, 33, 4.1406)
(7.4, 1.3, 16.3, 25.0, 6.9063)
(6.7, 1.1, 12.2, 20.0, 11.6094)
Algorithm1 (0.01)
(1.0, 1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0, 1.0000)
(1.0, 1.0, 1.0, 1.0, 1.0000)
Algorithm1 (0.05)
(1.2, 1.5, 1.0, 1.1, 1.1850)
(1.1, 1.3,1.0, 1.0, 1.0981)
(1.0, 1.0, 1.3, 1.1, 1.0158)
(1.0, 1.1, 1.1, 1.1, 1.0202)
RAI (0.01)
(1.0, 1.1, 1.0, 1.1, 1.2879)
(1.0, 1.1, 1.0, 1.0, 1.6792)
(1.0, 1.1, 1.1, 1.1, 2.6955)
(1.0, 1.0, 1.0, 1.1, 4.4231)
RAI (0.05)
(1.1, 1.3, 1.0, 1.1, 1.6636)
(1.1, 1.2,1.0, 1.0, 1.9857)
(1.0, 1.1, 1.2, 1.2, 3.3154)
(1.0, 1.1, 1.1, 1.1, 5.2848)
X-G-Z (0.01)
(1.0, 1.6, 1.0, 1.0, 7.2486)
(1.2, 2.3, 1.2, 1.3, 6.0151)
(1.3, 2.5, 1.3, 1.4, 7.8868)
(1.4, 2.9, 1.7, 1.7, 8.5505)
X-G-Z (0.05)
(1.1, 1.6,0.9, 1.0, 9.9249)
(1.2, 2.4, 1.1, 1.2, 8.1661)
(1.3, 2.7, 1.4, 1.4, 16.0270)
(1.3, 3.0, 1.9, 1.8, 11.6998)
MMHC (0.01)
(1.2,0.6, 1.1, 1.1, 1.4798)
(1.4,1.0, 1.1, 1.2, 2.1321)
(1.4, 1.1, 1.2, 1.2, 3.8258)
(1.3,1.0, 1.2, 1.2, 6.7052)
MMHC (0.05)
(1.3,0.6, 1.1, 1.2, 1.9827)
(1.4, 1.1, 1.2, 1.3, 2.5774)
(1.5, 1.2, 1.2, 1.3, 4.8484)
(1.4, 1.1, 1.2, 1.3, 8.1278)


## 6. Conclusions

In this paper, we have given a more accurate characterization of moral graph and proposed a new algorithm for structural learning, which substantially improves on the K2 algorithm. We also extend our algorithm to manage networks without node ordering. Although the new algorithm depends on the quality of constructed moral graph, simulation studies illustrate that our method yields good results in a variety of situations, especially when the underlying data set is small.The results in this paper also raised a number of interesting questions for future research. We briefly comment on some of those questions here. First, maximal prime decomposition plays a key role in Algorithm1. Although decomposition of an undirected graph into its maximal prime sub-graphs has been discussed a lot, we believe that there is room for further improvements. Second, we have applied the K2 algorithm for the learning of local structures in Algorithm 1. It will be interesting to see whether there exists some alternative approach to serve the same purpose here. Finally, although we assume in this paper that the data are completely observed, missing data or data with latent variables may arise in practice. Generalization of the proposed algorithm to missing data or data with latent variables is of great research interest.



---

*Source: 974063-2012-10-31.xml*