{'fraction_non_alphanumeric': 0.053523922970021835, 'fraction_numerical': 0.011316259678379988, 'mean_word_length': 3.9873267326732673, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Number of Samples Needed to Learn the Correct Structure of a Bayesian Network

Bayesian Networks (BNs) are useful tools giving a natural and compact representation of joint probability distributions. In many applications one needs to learn a Bayesian Network (BN) from data. In this context, it is important to understand the number of samples needed in order to guarantee a successful learning. Previous work have studied BNs sample complexity, yet it mainly focused on the requirement that the learned distribution will be close to the original distribution which generated the data. In this work, we study a different aspect of the learning, namely the number of samples needed in order to learn the correct structure of the network. We give both asymptotic results, valid in the large sample limit, and experimental results, demonstrating the learning behavior for feasible sample sizes. We show that structure learning is a more difficult task, compared to approximating the correct distribution, in the sense that it requires a much larger number of samples, regardless of the computational power available for the learner.


Introduction
Bayesian Networks (BNs) are generative models which are well suited for representing knowledge under uncertainty. Their compact representation, modularity and intuitive causal interpretation has made them very popular, and today they are used in various fields such as AI, Expert systems, Economics and Computational Biology. A BN is specified by two components: a qualitative part (structure), which represents conditional (in)dependencies between the random variables (r.v.s.), and a quantitative part (parametrization), representing the exact joint probability distribution of related r.v.s. in the network. Often there is no way of building BNs according to expert knowledge. In these cases one needs to learn a BN from data, i.e. from a sample of N realizations (values taken by the network's variables), and it is important to assess the sample size one needs in order to learn a BN that approximates the true one with given quality. The subject of sample complexity of BNs has drawn some attention in the past decade. The PAC framework (Valiant 1984) was used in (Dasgupta 1997) for BNs with known structure, with and without hidden variables, and by various authors (Abbeel et al. 2005;Friedman and Yakhini 1996;Hoffgen 1993), for discrete BNs with unknown structure. All of these works used the Kullblack-Leibler distance (Kullback and Leibler 1951) (also known as relative entropy) between the original and learned model's distributions to measure the approximation quality. A different criterion was presented in (Greiner et al. 1997), who measured the performance of the learned model (compared to the original) when answering queries (which depends, of course, on the query distribution). Our study is concentrated on a different aspect of the learning task, which is the number of samples needed in order to learn the correct network structure.
Intuition and practice indicate that learning the correct structure 1 is a much harder task than approximating the original distribution, and hence requires a larger number of samples. In (Dai et al. 1997) the relation between the number of causal errors (e.g. misidentified edges and edges directions) and sample size (as well as links strength) was studied using computer simulations for two specific learning algorithms. We derive rigorous upper and lower-bounds on the error probability, as well as experimental results, which do not depend on the learning algorithm used.
For simplicity, we concentrate on BNs with boolean r.v.s, yet our techniques can be applied to any discretenode BNs. We give bounds on the error probability, i.e. the probability of learning a wrong structure, for the Bayesian-Information-Criteria (BIC), which is a specific case of the Minimum-Description-Length (MDL) score. Our results can be easily modified to give the asymptotic behavior for MDL scores with other penalty terms. We obtain different bounds for structures which cannot represent the true distribution ('under-fitting'), and structures which do represent the true distribution yet do not have the minimal number of parameters ('over-fitting'). For these two tasks we use, respectively, large and moderate deviations results. While these results are correct in the large sample limit, we also give experimental results for a moderate number of samples. We address the problem of learning the exact correct BN structure (up to equivalence classes considerations), and do not consider the case of learning a structure that approximates the original one. A relaxation of this requirement (for example allowing a certain fraction of mislearned edges), leads to lower sample complexity. On the other hand, we assume that the learner is computationally unbounded and learns by scoring all possible models exhaustively and selecting the best scoring one. Obviously, the number of samples needed by real-life algorithms for structure learning may be higher.
In the next section we give basic definitions and results for BNs, MDL score, and relative entropy, which are used later. In section 3 we give asymptotic sample complexity bounds for learning the correct structure. In section 4 we present computer simulations demonstrating the behavior of the error probability for a moderate number of samples. We end with conclusions and future directions.

Preliminaries and Definitions
Let X 1 , .., X n be binary random variables (r.v.s.), with joint probability distribution P . In general, uppercase will denote random variables, and lowercase will denote their realizations. We will sometimes omit the latter, so for example P (x 1 , .., x n ) = P (X 1 = x 1 , .., X n = x n ). Formally, a Bayesian Network (BN) is a pair B = G, Θ . Here G = V, E is an n-vertex directed acyclic graph (DAG) , called the structure, which represents (in)dependencies between the X i 's, such that a variable is independent of its non-descenders, given its parents (Pearl 1988) where P a G (i) are the parents of X i in the graph G. The set of parents of X j excluding X i is denoted P a G (j \ i), and its values are denoted by pa G (j \ i). Keeping the other parent's values and setting X i = k will be denoted by P a G (j \ i : k), for k = 0, 1.
The set of probability distributions on m values is the m − 1-dimensional simplex denoted S m . The joint probability distribution associated with B is denoted P B , P B ∈ S 2 n , and satisfies: The graph dimension |G| = n i=1 2 |P a G (i)| is the number of parameters needed to specify P B when the structure is G. We assume that the correct BN belongs to a set of distributions which are "bounded-away from zero" with a bound γ > 0. This simply means (2) For any edge (i, j) ∈ E, its 'information content' in B is defined as: We also define IC B ≡ min (i,j)∈E IC B (i, j). The parameters γ and IC B (which depend on B) influence the sample complexity bounds we will obtain, and when they tend to zero, these bounds become large.
For a distribution P we denote by I(P ) the set of all conditional independence relations of the form X ⊥ Y |Z that hold in P , for X, Y, Z disjoint sets of r.v.s. For a graph G we denote by I(G) all the independence assertions implied by G. G is called an I-map for P if I(P ) ⊃ I(G). G is called a P-map for P if I(G) = I(P ). We denote by M(G) the set of all distributions for which G is an I-map, M(G) ⊂ S 2 n . We will assume that our samples are generated from a probability distribution for which G * is a P-map. The BN generating the samples is denoted B * = G * , Θ * , and the corresponding distribution is denoted P B * . G * is the unique P-map for P B * , up to graph-equivalent structures. Moreover, |G * | < |G| for any other I-map G of P . We refer to G * as the 'correct' structure for P B * , and our purpose is to recover it (more precisely its equivalence class) from the data.
A fully-connected DAG which is consistent with the ordering of G * is denoted C * . The set of graphs obtained from fully-connected DAGs by removing the edge (i, j) (or (j, i)) is denoted C ij , where the same notation is used for all such DAGs. We denote by N the number of realizations (also referred to as samples) in the learning set. The samples are denoted by j is the value of the r.v. X j in the i-th sample. We assume the samples are i.i.d., with x (i) ∼ P B * . For a distribution P we denote by P (N ) the joint product measure of N i.i.d. r.v.s., each distributed according to P . Thus, for example, One can use a sample of size N to estimate P B * by simply counting the number of occurrences of each value x ∈ {0, 1} n in the sample. The resulting distribution is called the sample distri-butionP N , given byP For any B = G, Θ , the log-likelihood of the data is: Given G, the maximal likelihood parametrizationΘ = Θ(G), is simply given by the sample probability, sô Θ i|pa G (i) =P N X i = 1|pa G (i) . Using simple algebra, the likelihood of the model G is given by: Where H P is the entropy or conditional entropy of a variable with respect to P . Similarly, I P denotes the mutual information (or conditional mutual information) of two variables. We denote P N,G the sample distribution assuming that the correct structure is G, which is the probability distribution associated with the BN G,Θ(G) . If G is the best scoring structure, then P N,G =P N .
A score is called asymptotically consistent if, when N → ∞, the 'correct' model will attain the highest score with probability approaching one. Here, the 'correct' model refers to the graph (or actually its equivalence class) with minimal number of parameters which is an I-map for P B * . It can be easily shown that the log-likelihood is not a consistent score, and it is thus not useful for comparing structures. Adding edges to the graph always improves the likelihood, which makes the complete graph the highest scoring one, regardless of the true BN generating the data. A common strategy to cope with this overfitting problem, is based on the Minimal Description Length (MDL) principle (Rissanen 1978). The M DL score 'penalizes' complex models, thus giving a trade-off between data-fitting and model complexity. It usually takes the form: where Ψ ≡ Ψ(N ) is a penalty function. We assume that Ψ satisfies: Under these assumptions, the MDL score is asymptotically consistent (Haughton 1988(Haughton , 1989. Of particular interest is the choice Ψ(N ) = 1 2 log N , since in this case the score in eq. 6 is known to be asymptotically equivalent to the Bayesian score (with any nowherevanishing prior), and is also termed the Bayesian Information Criterion (BIC) (Schwarz 1978).

Relative Entropy Properties
Working with the relative entropy distance measure often possesses technical difficulties, since it does not satisfy the requirements of a norm. Here we prove two results for the relative entropy, analogous to the symmetry and triangle inequality properties of a norm. Our results rely on the fact that the reference distribution is bounded away from zero, and involve constants which depend on this proximity to zero.
Let P, Q be two strictly positive discrete probability distributions over m values. Their relative entropy is: For a graph G and a probability P we denote D(G||P ) = inf Q∈M(G) D(Q||P ). We can always write Q uniquely as Q = P + V , where > 0 and V is a unit vector. Moreover, since P and Q are probability distributions, V also satisfies: Taking Taylor expansion of the relative entropy gives: We see that the relative entropy behaves locally as O( 2 ). Moreover, it is 'locally symmetric', in the sense that D(P ||P + V ) D(P + V ||P ) → 1 as → 0. We can also formulate the following uniform bounds, independent of V 's direction: Using Taylor expansion we get: Where convergence is guaranteed since γ < 1 2 . For D(Q||P ), just flip between P and Q, and since min i Q i ≥ γ 2 , the bound is multiplied by four.

Learning the Correct Structure
This section is devoted to the problem of identifying the correct network structure from data. This problem is also known in the statistics literature as 'model selection'. For exponential families, the problem was studied in (Haughton 1988), where asymptotic consistency was established, and in (Haughton 1989), where asymptotic results on the error-probability as a function of the number of samples were obtained. In (Geiger et al. 2001) it was shown that directed graphical models, with no hidden variables, are in fact curved exponential families, and thus Haughton's results are applicable for them. Nevertheless, the explicit constants appearing in Haughton's asymptotical analysis were not fully characterized. We will present these results in the context of learning BNs structure, and, when possible, give bounds on the constants governing the decay of the error probability.
The different graph structures, which are curved exponential families, can be thought of as Riemannian manifolds M(G) ⊂ S 2 n (Geiger et al. 2001) of different dimensions which are subsets of the simplex S 2 n , 'competing' for their fit to the data, which can be casted simply as their (Kullback-Leibler) distance from the sample probabilityP N ∈ S 2 n , where each such manifold 'pays' a penalty proportional to its dimension. While analyzing these models together seems complicated, our approach is studying the error probability of one model at a time, i.e. P (N ) Following (Haughton 1988(Haughton , 1989, we divide the models into two disjoint subsets: 1. Graphs G which are not I-maps for P B * (i.e. P B * / ∈ M(G)), and 2. graphs G which are I-maps for P B * , yet have higher dimension than G * , P B * ∈ M(G) and |G| > |G * |.

Graphs which are not I-maps for P B *
We begin by bounding the error in learning a graph which is not an I-map for P B * : Let G be another graph which is not an I-map for P B * . Then ∃c > 0 such that: If, in addition, |G| ≤ |G * |, then: Proof The first inequality follows directly from (Haughton 1989), proposition 2. For the second inequality, note that ifP N ∈ M(G), then LL N (G * ) ≤ LL N (G), and therefore S N (G * ) < S N (G). This event happens, according to the Sanov theorem 2 (Sanov 1957), with probability Ω(2 −N D(G||P B * ) ), which completes the proof.
Note that for a given P and G, it might be difficult to compute the relative entropy D(G||P ). However, the relative entropy in the other direction D(P ||G) is easily given by: (14) One can apply lemma 2, to get an explicit lower-bound on the error exponent in the above theorem. Unfortunately, theorem 1 does not give us an explicit upperbound on the error exponent. The rest of this section is devoted to finding such a bound. Following (Friedman and Yakhini 1996), we derive our upper-bound in two steps. First, we assume that the observed distribution is ideal, that isP N = P B * . Although this assumption is unrealistic, and many times even not feasible, it helps us to understand the trade-off between data-fitting and model complexity in the MDL score. Next, we study the effect of sampling noise on our bound. By using concentration of measure arguments (specifically Chernoff bounds), we show that with high probabilityP N is close to P B * , allowing us to get an upper-bound in the presence of sampling noise. Before deriving our results, we bring a useful relation for Pmaps. From definition, if G = V, E is a P-map for P , then I P X i , X j |P a G (j \ i) > 0 for any (i, j) ∈ E. For strictly positive BNs, the following stronger relation holds: Proposition 1 Let B = G, Θ be a BN with G = V, E a P-map for P B and Θ ∈ (0, 1). Then IC B > 0.
The positivity of IC B means that no edge in G is 'redundant' with respect to the distribution P B . The magnitude of IC B can be thought of as a measure of how much information on P B is captured in G, compared to lower-dimensional models. Intuitively, one expect that the higher this value is, the easier it will be to separate G from lower-dimensional models based on data sampled from P B . The bounds obtained in the next section show indeed such dependence on IC B , as well as a dependence on γ. Our purpose is not to achieve the tightest bounds possible, but merely demonstrate the dependence of the convergence rate on these two parameters.

Deriving Bounds in the Ideal Case
In the ideal case all our statistics are deterministic, and depend only on the sample size. We denote them with the superscript (I), thus, for example, the loglikelihood is denoted LL Proof Assume, negatively, that ∃G, S The likelihood in the ideal case is given by: The likelihood loss, when removing an edge (i, j), is: Which is positive by proposition 1. By using lemma 4 and eq. 18, we get: Proposition 2 In the ideal case:

Treatment of the Noisy Case
Note that a weaker form of lemma 4 is still valid in the noisy case, assuming that we require only S N (G * ) ≥ max G,|G|≤|G * | S N (G). Our method of proof is showing that the likelihood difference LL N (G * ) − LL N (C ij ) in the noisy case, is close, with high probability, to the ideal version shown in eq. 18. In order to show this proximity, we use a series of concentration lemmas: Lemma 5 Let B = G, Θ be a BN. Take a subset S ⊂ {X 1 , .., X n } of r.v.s. Then ∀α ∈ (0, γ n ): =s} be the r.v. counting the number of samples in which the value of the S variables was s. Then Y N ∼ Binomial N, P B (S) . Using Chernoff bounds we get for α ∈ (0, 1), Noting that Y N = NP N (S) completes the proof.
Lemma 6 Let B = G, Θ be a BN. Take a subset S ⊂ {X 1 , .., X n } of r.v.s. Then ∀α ∈ (0, 1): Proof The function f (x) = x log x is Lipschitz continuous in the interval [a, 1] with Lipschitz constant log a + 1. From lemma 5, the inequality P N (S) − P B (S) holds with the desired probability. In that case, obviously we have γ n − α ≤P N (S), P B (S) and from here eq. 22 follows immediately.
We now give concentration bounds on the empirical entropy and conditional entropy functions: Lemma 7 Let B = G, Θ be a BN. Take two disjoint subsets S, T ⊂ {X 1 , .., X n } of r.v.s. Then ∀α ∈ (0, γ n ): P And: P Proof For the first part, simply sum over all possible realizations of S in lemma 6, and apply the union bound, noting that P B (S) < 1. For the second part, use the chain rule for entropy H(Y |X) = H(X, Y ) − H(X). Apply the first part on the sets S and S ∪ T , and apply the union bound again to get the desired result.
The next lemma bounds the difference in the loglikelihood gap LL N (G * )−LL N (C ij ) between the ideal and the noisy case: Lemma 8 Let B = G, Θ be a BN. Take a subset S ⊂ {X 1 , .., X n } of r.v.s. Then ∀α ∈ (0, 1): Proof Recall, from eq. 5, that the log-likelihood can be written as LL N (G) = −N n j=1 HP N X j |P a G (j) . Apply lemma 7 with T = X j and S = P a G (j), for G = G * and ∀j = 1, .., n. This can be used, along with the union bound, to bound |LL N (G * ) − LL Combining the bounds by using the triangle inequality gives eq. 25.
We are now ready to prove the main result of this section, giving our asymptotic upper-bound on the error exponent: Proof Set α = min γ n 2 , IC B * 2 n+2 [n log(γ/2)+1] . From lemma 8 and eq. 18, we get: Assume that N is large enough such that Ψ(N ) . For such N 's, we get: By using the weaker form of lemma 4, and taking the union bound over all C ij 's, we get: Taking the logarithms of both sides and substituting the α value we have chosen completes the proof.

Graphs which are I-maps for P B *
In the previous section we have seen that if a graph G is not an I-map for P B * , then M(G) is 'bounded away' from P B * in the simplex S 2 n . This was used to relate the error probability to a large deviation event. When the graph G is an I-map for P B * , we have P B * ∈ M(G). The event of choosing an overparameterized I-map is a moderate deviation event, and the next theorem characterizes its asymptotic probability: Proof Recall that our models are in fact curved exponential families. Then the first statement is a direct consequence of proposition 3 in (Haughton 1989). The second is due to (Woodroofe 1978(Woodroofe , 1982. Note that for nested models, the asymptotic error depends only on G * and not on the specific parametrization determining P B * . Computer simulations we have performed did not find a dependence on the parametrization also for non-nested models. This is substantially different from the under-fitting error shown in the previous section, where we have seen that the error exponents depends on the parametrization. The constant d is proportional (also in Haughton's result) to the difference in the dimensions |G|−|G * |: The more over-parameterized G is, the faster the decay of the error probability.

Experimental Results
In the previous section we noted a qualitative difference between the over and under-fitting errors. While the under-fitting error decays exponentially fast with N , the over-fitting probability decays slower, as a power of N . As we show now, this difference is relevant mainly in the large N limit, while for small N 's, as one might expect intuitively, the situation is the opposite and under-fitting is more likely. We examined the errors in the learning process for a small BN with four r.v.s. There are 543 different DAGs on four nodes, divided into 185 equivalence classes. We chose G * to be the graph with the edge set E * = (1, 2), (1, 3), (2, 4), (3, 4) , and the parametrization to be Θ * : Θ * 1 = 0.1, Θ * 2|0 = 0.1, Θ * 2|1 = 0.3, Θ * 3|0 = 0.1, Θ * 3|1 = 0.3, Θ * 4|00 = 0.1, Θ * 4|01 = 0.3, Θ * 4|10 = 0.8, Θ * 4|11 = 0.2 . There are 4 equivalence classes (including G * ) which are I-maps for P B * , and 181 which under-fit it. We also took two specific structures, G 1 with E 1 = (1, 2), (1, 3), (1, 4), (2, 4), (3, 4) , which is an over-parameterized I-map for P B * , and G 2 with E 2 = ∅, which is not an I-map for P B * . Since the errors become rare events for large N , we had to use importance-sampling methods in order to estimate their probability. Rather than generating samples directly from P B * , we used other distributions Q, for which these events are more likely (for example, for G 2 we took Qs for which the r.v.s. are almost independent), and applied the appropriate correction needed to estimate the probability given P B * . The error probabilities are presented in figure 1. For N < 2000, the error probability is dominated by the under-fitting structure, while for larger N s the over-fitting probability is higher. Qualitatively similar results were also obtained for other networks we have started with and other choices of competing wrong graphs.  Figure 1: The log-error probabilities for learning two wrong graphs, shown as functions of the sample size for G 1 (in blue) which is an over-parameterized I-map for P B * , and for G 2 (in red) which is not an I-map for P B * . We used importance sampling and averaged over 30 different Q distributions, where 6000 samples were drawn from each. The G 2 log-error exhibits an excellent fit for a straight line. The G 1 error is fitted better to a straight line on a loglog plot (not shown), in accordance with its power-law asymptotic behavior.
We presented new lower and upper-bounds on the number of samples needed to ensure that the true and learned models share the same structure. When considering only graphs which are not I-maps for the correct distribution P B * , we have shown that the error probability decays exponentially with N . When allowing graphs with are I-maps for P B * but are overparameterized, we showed that the error probability decays as a power of N . While these results imply that in the large N limit the error is dominated by over-fitting, for small values of N , the opposite situation of under-fitting the true model is more likely. This observation is corroborated by the fact that in few small BNs we have examined, the error exponent (both in our bounds and from numerical simulations) for the under-fitting case was close to zero, indicating a slow (although exponential) decay for moderate sample sizes. Although the BIC score has the desired property of consistency, perhaps a different choice of the penalty function (e.g. smaller than 1 2 log N for small values of N and larger for large values of N ), may lead to a smaller error probability.
The bounds we have obtained relate to the error of choosing a specific wrong model G, when its score is higher than that of the 'correct' model G * . Typically, however, one needs to identify G * from a set of candidate models. Here one needs to look simultaneously on the (dependent) events S N (G * ) ≤ S N (G) for all possible wrong graphs G. A challenging future direction is trying to bound the error probability for learning any model from a set of wrong models together. We note that the number of possible candidate graphs grows super-exponentially with the network size n, even if we consider only bounded in-degree networks. Therefore, a simple union-bound argument may not give satisfactory results here, and one needs to turn into more sophisticated techniques, if one wants to study the dependence of sample complexity on the size of the network.