Paper deep dive
High-Order Markov Blanket Discovery via a k-Order Relaxation of the Faithfulness Assumption
Loong Kuan Lee, Ragavi Krishnamoorthy, Nico Piatkowski
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 8/4/2026, 10:30:14 AM
Summary
This paper addresses the limitations of the faithfulness assumption in Markov Blanket (MB) discovery, which fails to capture higher-order dependencies like XOR and parity relations. The authors propose a 'k-order' relaxation of the faithfulness assumption and introduce the k-order Markov Blanket (kOMB) algorithm. Empirical results demonstrate that kOMB successfully recovers the Markov blanket under both true and empirical faithfulness violations, outperforming existing constraint-based methods.
Entities (8)
Relation Signals (6)
k-order Markov Blanket → proposesrelaxationof → Faithfulness Assumption
confidence 95% · We then propose a proof of concept algorithm called k-order Markov blanket (kOMB) that uses this relaxation for MB discovery.
k-order Markov Blanket → overcomes → Empirical Faithfulness Violations
confidence 92% · empirically show how kOMB can recover the MB of a variable under both true and empirical violations of faithfulness.
Faithfulness Assumption → violatedby → XOR Function
confidence 90% · this assumption can be violated by higher-order dependencies such as XOR and parity-type relations
Faithfulness Assumption → violatedby → Parity Function
confidence 90% · captures parity type relationships between k+2 variables
2-adjacency faithfulness → failstohandle → Parity Function
confidence 85% · However, although the 22-adjacency faithfulness can find XOR-type relations, it is unable to find higher-order parity relationships between more than 33 variables.
2-adjacency faithfulness → handles → XOR Function
confidence 85% · 2-adjacency faithfulness in particular is quite practical as by limiting faithfulness violations to “Unfaithful Triples”
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The problem of learning the graphical Markov blanket (MB) of a variable from data has applications in many areas such as structure learning for Bayesian networks and Markov random fields, causal discovery, and feature selection. However, a common assumption most methods make is that the conditional independencies in the distribution imply the same separation in the graphical structure -- also known as the faithfulness assumption. Unfortunately, this assumption can be violated by higher-order dependencies such as XOR and parity-type relations, and -- on finite samples -- by empirical violations that, in extreme cases, even induce spurious dependencies absent from the true distribution. Therefore, in this paper we propose a "k-order" relaxation of the faithfulness assumption that captures parity type relationships between k+2 variables. We then propose a proof of concept algorithm called k-order Markov blanket (kOMB) that uses this relaxation for MB discovery. Finally, we empirically show how kOMB can recover the MB of a variable under both true and empirical violations of faithfulness. Code available at: this https URL
Tags
Links
- Source: https://arxiv.org/abs/2607.26357v1
- Canonical: https://arxiv.org/abs/2607.26357v1
Trouble viewing inline? Open PDF directly →
Full Text
77,628 characters extracted from source content.
Expand or collapse full text
High-Order Markov Blanket Discovery via a k-Order Relaxation of the Faithfulness Assumption Loong Kuan Lee Hybrid Intelligence Fraunhofer IAIS Sankt Augustin, Germany Ragavi Krishnamoorthy Hybrid Intelligence Fraunhofer IAIS Sankt Augustin, Germany Nico Piatkowski Hybrid Intelligence Fraunhofer IAIS Sankt Augustin, Germany Abstract The problem of learning the graphical Markov blanket (MB) of a variable from data has applications in many areas such as structure learning for Bayesian networks and Markov random fields, causal discovery, and feature selection. However, a common assumption most methods make is that the conditional independencies in the distribution imply the same separation in the graphical structure—also known as the faithfulness assumption. Unfortunately, this assumption can be violated by higher-order dependencies such as XOR and parity-type relations, and—on finite samples—by empirical violations that, in extreme cases, even induce spurious dependencies absent from the true distribution. Therefore, in this paper we propose a “k-order” relaxation of the faithfulness assumption that captures parity type relationships between k+2 variables. We then propose a proof of concept algorithm called k-order Markov blanket (kOMB) that uses this relaxation for MB discovery. Finally, we empirically show how kOMB can recover the MB of a variable under both true and empirical violations of faithfulness. Code available at: https://github.com/lklee9/k-order-Markov-blanket. 1 Introduction The goal of the Markov blanket (MB)discovery problem is to infer the MBof some variable Y in the graph given a dataset sampled from the joint distribution . Originating from the Bayesian network and causal learning literature [margaritis1999, tsamardinos2003, ling2021, hu2024], the MB discovery problem has also been discussed and utilised in the feature selection literature as well [tsamardinos2003a, wang2020]. More recently, MBdiscovery has also been employed as a subroutine within distributed, divide-and-conquer approaches to large-scale BNstructure learning [dong2025a]. That said, we are interested in inferring MBsin both Bayesian network (BN)and Markov random field (MRF). Therefore, whenever possible we will try to keep our notation and exposition agnostic to the graph type of . X1X_1X2X_2X3X_3YYZ:=[∑i=13Xi=1]Z:= [ _i=1^3X_i=1 ] (Y=Z∣ )=0.9 (Y=Z )=0.9 Recall (%) Sample Size Method 100 1000 Grow and Shrink 28.33 100 0-OMB (ours) 28.33 100 1-OMB (ours) 77.50 100 2-OMB (ours) 100.0 100 Figure 1: Recovering the MBof variable Y when Y is the output of a noisy Boolean function that returns 11 when only one of X1X_1, X2X_2, and X3X_3 is 11. Our algorithm that assumes a 22-order faithfulness relaxation is able to fully recover the MBof each variable—even with just 100100 samples. Most existing approaches to MB discovery operate under two main assumptions. The first being that any separations in imply the same conditional independence in . The second assumption—known as the faithfulness assumption—goes the opposite direction, that is any conditional independencies in imply the same separation in . Theoretically, BNswith faithfulness violations have Lebesgue measure zero [meek1995, spirtes2001, boeken2025]; however these “typicality” results depend on somewhat arbitrary choices, like the choice of σ-ideal [boeken2025]. There are well-known examples of some—very simple—relations that violate faithfulness. The most famous of which is the (noisy) XOR function for violations involving three variables [inazumi2011, marx2021]; which is a special case of the (noisy) parity function whose higher-order faithfulness violations we detail in Example˜1. That said, even if unfaithful BNsand MRFsare truly almost non-existent in nature, a lack of sufficient samples can still cause an empirical faithfulness violation [uhler2013, lemeire2012, boeken2025]. For example, the small example in Figure˜1 is one such case where the true BNis faithful, but due to a finite sample, empirical faithfulness violations occur—which can be overcome by the algorithm we will propose later in Section˜5. A full walkthrough of empirical faithfulness violation can be found later in Example˜2 under Section˜3. Furthermore, in extreme cases, a limited sample size can cause MBdiscovery algorithms to falsely discover spurious dependencies absent from the true distribution. We later observe this phenomenon in Figure˜4 under Section˜6.1 with the low precision exhibited by the different MBdiscovery methods at low sample sizes. There are several weaker assumptions compared to faithfulness such as adjacency faithfulness [spirtes2001], Pearl-minimality [pearl1988, pearl2009], SGS-minimality [spirtes2001], frugality [forster2018], and more recently 22-adjacency faithfulness [marx2021]. 22-adjacency faithfulness in particular is quite practical as by limiting faithfulness violations to “Unfaithful Triples”, it yields an algorithm with a much smaller search space compared to frugality. However, although the 22-adjacency faithfulness can find XOR-type relations, it is unable to find higher-order parity relationships between more than 33 variables. Contributions Therefore, motivated by the considerations above, in Section˜4 we propose a relaxation of faithfulness for parity-like relationships between k+2k+2 variables: the k-order faithfulness assumption. This relaxation of faithfulness allows us to infer—from sample data—the graphical Markov blanketof a variable under mild assumptions. We show this in Section˜5 by implementing a proof-of-concept algorithm called the k-order Markov blanket (kOMB). We then empirically show in Section˜6 that kOMB can overcome both true and empirical faithfulness violations and that it outperforms some existing constraint-based methods on well-known benchmark datasets—where we find merit in exploring relationships between 44 variables at a time, and not just 33. 2 Probabilistic Graphical Models and Independencies Consider the set of n variables =V1,…,Vn =\V_1,…,V_n\ with no other hidden or latent variables. The conditional independencies between these variables— ∣ for disjoint variables sets , , ⊆ , , —can be represented by either a directed acyclic graph (DAG) [pearl1988] or an undirected graph (UG) lauritzen1996. That said, we will denote both graphs with the tuple =( , ) =( , ) since we are agnostic on the type of graph used for representing these graphical conditional independencies . Regardless of the graph type of , for disjoint subsets , , ⊆ , , , we will denote that and are separated in as, ∣ . Furthermore, we denote the set of conditional independencies encoded by the graph —i.e. the independence model of —as, ⟨ , ∣ ⟩∈ ⇔ ∣ , , ∈ , (1) where ⟨ , ∣ ⟩ , is known as an independence statement [sadeghi2017]. In addition to the graphical structure , we also consider a joint distribution over that might encode its own set of conditional independencies between the variables in , ⟨ , ∣ ⟩∈ ⇔ ∣ , , ∈ , (2) where ∣ denotes and are probabilistically conditionally independent given . So far, and are two separate collections of conditional independencies in the graph and distribution respectively. Our base assumption throughout this paper, the Global Markov Property, links the two by requiring that every conditional independence encoded by the graph also hold in the distribution . This is a common assumption in the Markov blanket discovery literature schluter2014, yu2020a, marx2021. Assumption 1 (Global Markov Property). Given a PGM =( , ) =( , ), we say that the distribution is Markov to the graph when for all disjoint subsets , , ⊆ , , [hammersley1971, lauritzen2018], ∣ ⟹ ∣ . . (3) Furthermore, by Equations˜1 and 2, the relation in Equation˜3 is equivalent to ⊆ . Despite their differences, the independence models of and are both known to obey a set of axioms called the (semi)-graphoid axioms [pearl1989] (see Definition˜6). As mentioned in Section˜1, we are interested in learning the graphical Markov blanketof some Y∈ ∈ defined as follows. Definition 1 (Markov Blanket and Boundary). For a variable Y∈ ∈ , the MBof Y in is any subset ⊆ ∖Y Y111We allow set operations between variable sets and single variables to implicitly mean the operation between the set and the atomic set of the single variable—for example: ∖Y≔ ∖Y Y \Y\. s.t. Y ∖(Y∪ )∣ ,Y (Y∪ ) , (4) or in other words, every other variable in ∖(Y∪ ) (Y∪ ) is graphically separated from Y by . The Markov boundary of the variable Y—henceforth denoted as (Y) (Y)—is then defined as the minimal Markov blanket of Y [pearl1988]. An equivalent definition for the Markov blanket and boundary of Y in can be made by just substituting the graph separation term in Equation˜4 for conditional independence in . Therefore, it is more accurate to say that we wish to learn the Markov boundary of Y and not just its Markov blanket. However, since colloquially the term Markov blanket is commonly used to implicitly refer to the Markov boundary, we shall use the term Markov blanket in such a fashion as well. Several methods exist for discovering the Markov blanket (MB)of a target variable Y∈ ∈ in with one of the larger class of methods being constraint-based approaches [margaritis1999, tsamardinos2003, schluter2014, ling2019]. These approaches iteratively test conditional independence statements—using data sampled from —to grow and shrink a candidate MBof Y. However, the MBof Y in the joint distribution might not be the same as its MBin the graph . Following the Global Markov Property in Assumption˜1, we at know that the MBof Y in is a subset of the MBin —which this might prevent us from fully recovering the MBin . Therefore, in order to ensure that the Markov blankets of Y in both and are equal, a common assumption constraint-based approaches make is that is faithful to the joint distribution as defined in Definition˜2 [spirtes2001]. Definition 2 (Faithfulness). Given a PGM =( , ) =( , ), the graph is faithful to the joint distribution if for all disjoint subsets , , ⊆ , , [spirtes2001], ∣ ⟹ ∣ (5) which is basically the reverse direction of the global Markov property in Assumption˜1. Therefore, by assuming both the Global Markov Property faithfulness, the conditional independencies of both and then coincide, = = . However faithfulness can be violated in many ways and we will explore some of them in the following section. 3 Faithfulness Violations and Existing Relaxations Probably the most well-known example of a faithfulness violation is the (noisy) XOR function Y≔X1⊕X2Y X_1 X_2 with the graphical structure X1→Y←X2X_1→ Y← X_2 [inazumi2011, marx2021]222The (noisy) XOR function is a special case of the (noisy) parity function whose faithfulness violation we detail in Example 1.. A recent faithfulness relaxation for tackling XOR-type relations that is of particular interest to us is 22-adjacency faithfulness [marx2021]. Definition 3 (2-adjacency faithfulness). Given the PGM =( , ) =( , ), is 22-adjacent faithful to if for all adjacent variables (X,Y)∈ (X,Y)∈ , ∃ ⊆ (Y)∃ (Y) where X∈ ∈ s.t. ∀X∈ :X Y∣ ∖X∀ X∈ :X Y X and if | |=2| |=2, X1 X2∣YX_1 X_2 Y. That said XOR-type relationships only involve three variables. Since we are interested in higher-order unfaithfulness, we will use the noisy parity function as a primary example of the type of faithfulness violation we are interested in. X1X_1X2X_2X3X_3Y ∀i∈[3]∀ i∈[3] (0) _X_i(0) (1) _X_i(1) 0.5 0.5 X1X_1 X2X_2 X3X_3 (1) _Y(1) 0 0 0 0.1 0 0 1 0.9 0 1 0 0.9 0 1 1 0.1 1 0 0 0.9 1 0 1 0.1 1 1 0 0.1 1 1 1 0.9 Figure 2: Example of a noisy parity BNover 44 variables. Example 1 (Noisy Parity Function). Let =( , ) =( , ) be a PGMwhere =X1,X2,X3 =\X_1,X_2,X_3\, =Y∪ =Y∪ , and has the structure of the DAGin Figure˜2. Each X∈ ∈ are independent of each other and has the distribution of a fair coin. The distribution of Y conditioned on is then (Y=f( )∣ = )=0.9 (Y=f( ) = )=0.9, where, f( )≡∑i=1nxn(mod2)f( )≡ _i=1^nx_n 2. For any random variable X∈ ∈ , we know that (Y,X)=0.25=0.5×0.5= (Y) (X) (Y,X)=0.25=0.5× 0.5= (Y) (X) for any value Y and X takes. Therefore it must the case that Y XY X, but since Y and X are not separable in , Y XY X; is not faithful to . Additionally for any distinct random variables Xi,Xj∈ _i,X_j∈ , (Y,Xi,Xj)=1/8=0.53= (Y) (Xi) (Xj) (Y,X_i,X_j)=1/8=0.5^3= (Y) (X_i) (X_j) for any values of Y,Xi,XjY,X_i,X_j. Therefore we know that Y Xi,XjY \X_i,X_j\, Xi Y,XjX_i \Y,X_j\, and Xj Xi,YX_j \X_i,Y\. Since Y is adjacent to all X∈ ∈ , but there is no 22-cardinality subsets ′⊂ ⊂ that makes Y ′Y ; is not 22-adjacent faithful to following the definition in marx2021. Although we now know of a BNthat violates both faithfulness and 22-adjacency faithfulness, these types of relationships between variables might very well not be common in nature. In fact, it has been shown that unfaithful BNshas Lebesgue measure zero [meek1995, spirtes2001, boeken2025]. However, empirical violations of faithfulness due to limited samples can be quite common [uhler2013, lemeire2012, boeken2025]. Therefore in Example˜2, we will present a scenario where, even though the true BNis faithful; due to sampling error, the DAGof the BNis not faithful to the empirical distribution obtained from the sample. Example 2. Similar to Example˜1, let =( , ) =( , ) be a PGMwith random variables and Y with the DAGstructure in Figure˜3. The distribution of Y conditioned on is then (Y=[∑X∈ =1]| )=0.9 (Y= [ _X∈ X=1 ] | )=0.9, where [⋅][·] is the Iverson bracket notation for the indicator function. In other words, Y=1Y=1 has a probability of 0.90.9 when only one of the variables in is 11. Otherwise, Y=0Y=0 has a probability of 0.90.9. The marginal distribution of Y is (Y=1)=3/8 (Y=1)=3/8 since only 3 value-combinations of has just a single 11 among them. We then know that ∀X∈ : (Y=1,X=1)=1/8≠(3/8)×(4/8)= (Y=1) (X=1)∀ X∈ : (Y=1,X=1)=1/8≠(3/8)×(4/8)= (Y=1) (X=1). Therefore, ∀X∈ :Y X∀ X∈ :Y X and consequently, is faithful to . However, it is possible for a sample of to be pathological in the sense that, the empirical distribution from induces a independence between Y and some X∈ ∈ . For example, the sample in Figure˜3 causes the marginal distribution of Y to be Phys.Rev.E(Y)=1/2 Phys. 10000\ Rev. 10000\ E(Y)=1/2 and the marginal distribution over Y and X1X_1 to be Phys.Rev.E(Y,X1)=1/4=(1/2)×(1/2)=Phys.Rev.E(Y)Phys.Rev.E(X1) Phys. 10000\ Rev. 10000\ E(Y,X_1)=1/4=(1/2)×(1/2)= Phys. 10000\ Rev. 10000\ E(Y) Phys. 10000\ Rev. 10000\ E(X_1) for all possible values of Y and X1X_1. Therefore under this empirical distribution Phys.Rev.E Phys. 10000\ Rev. 10000\ E, Y Phys.Rev.EX1Y Phys. 10000\ Rev. 10000\ EX_1—which empirically violates the faithfulness assumption. Fortunately in Example˜2, the DAGin Figure˜3 is still 22-adjacent faithful to the empirical distribution, therefore the full MBis still recoverable by the algorithm outlined in marx2021. However, this is not always the case and as we will see in Section˜6.1, sometimes considering higher-order relationships between variables can help better overcome empirical faithfulness violations. X1X_1X2X_2X3X_3YYg( ):=[∑X∈ =1]g( ):= [ _X∈ X=1 ] (Y=g( )∣ )=0.9 (Y=g( ) )=0.9 Pathological D X1X_1 X2X_2 X3X_3 Y 0 0 0 0 0 0 1 1 0 1 0 1 0 1 1 0 1 0 0 1 1 0 1 0 1 1 0 0 1 0 0 1 Figure 3: Example of a noisy exactly-1 BNover 44 variables with a pathological sample that causes the DAGto be unfaithful to the empirical distribution. 4 The k-Order Faithfulness Relaxation In this section we will introduce our generalisation of the faithfulness assumption in Definition˜2 for finding conditional dependencies between two variables in that require considering k additional variables to be found. But before that, we first define in Definition˜4 what it means for two distinct variables Y,X∈ ,X∈ to only be conditionally dependent given some k-cardinality subset ⊆ ,| |=k ,| |=k. Definition 4 (k-Order Dependence). Given a PGM =( , ) =( , ), let Y,X∈ ,X∈ be two distinct random variables and , ⊆ , be two disjoint subsets. Then we denote the absence of any subset of ′⊆ that renders Y and X conditionally independent given ∪ ′ ∪ as Y,X :=∀ ′⊆ :Y X∣ ∪ ′. Y,X :=∀ :Y X ∪ . (6) Here we say that Y and X are | || |-order associated w.r.t. over the “separating” set . Of course under this definition there might be some proper subset ′⊂ ⊂ where Y,X ′ Y,X still holds. Therefore we say Y and X are | || |-order dependent on only if this is not the case, Y,X := Y,X ∧∀ ′⊂ :¬ Y,X ′. aligned & Y,X \\ &:= Y,X ∀ ⊂ : \, Y,X . aligned (7) With this definition of k-order dependence, we are able to capture parity-like relations over k+2k+2 variables—such as the one in Example˜1. Specifically we know that the noisy-parity relationship in Example˜1 is 22-order dependent since for distinct variables A,B∈X1,X2,X3,YA,B∈\X_1,X_2,X_3,Y\, must always contain the other two variables for ∅A,B A,B to be true. Further note that when k=1k=1, our notion of k-order dependence— similar to 22-association in marx2021—is able to capture the relationships between Unfaithful Triples which includes XOR-type relationships. There are two aspects of Definition˜4 whose purpose might not immediately jump out. The first being that in Equation˜6, we require Y X∣ ∪ ′Y X ∪ to be true for all subsets ′⊆ —instead of just Y X∣ ∪ X ∪ . This requirement ensures that Y and X are dependent as long as they are conditioned on —forcing a distinction between the possible dependants and the possible separators . The other odd aspect of Definition˜4 is that we define both k-order association and dependence over some arbitrary subset ⊆ instead of just the entirety of . This allows us to reason about associations relative to the candidate Markov blanket maintained by our algorithm, which is central to proving in Section˜5 that correctly recovers the Markov boundary of some variable Y∈ ∈ . Since henceforth we need to reason about subsets with a restricted cardinality; before continuing we will first define some notational shorthand for dealing with such sets. Definition 5 (Power Sets of Limited Cardinality). Let be a set of elements and k∈ℕ:k≤| |k :k≤| | be some natural number333We adopt the definition of the natural numbers in which 0 is included. smaller than the cardinality of . Then we use, k := ⊆ :| |=k, k :=\ :| |=k\, (8) to denote the set of all subsets of with cardinality k and, k := ⊆ :| |≤k, k :=\ :| |≤ k\, (9) to denote the set of all the subsets of with cardinality less than or equal to k. Before introducing our generalisation of the faithfulness assumption to k-order dependencies, we first show in Definition˜5 that k-order associations w.r.t. implies the existence of a k-order dependence w.r.t. some subset of . restatable [k-Order Dependence Arises from Association]lemmakdeptoassoc Given a PGM =( , ) =( , ), let Y,X∈ ,X∈ be distinct random variables and , ⊂ , ⊂ be disjoint subsets. Then the following holds for all Y,X, , ,X, , : Y,X ⟹∃ ′⊆ : Y,X ′. Y,X ∃ : Y,X . (10) The proof of Definition˜5 and all other Lemmas in the paper can be found in Appendix˜B. Assumption 2 (k-Order Faithfulness). Given the PGM =( , ) =( , ), is k-order faithful to if ∀ ⊆ ∖Y,X∀ \Y,X\, ∀ ∈ ≤k( ∖Y,X)∃ ′⊆ :Y X∣ ∪ ′⟹Y X∣ . aligned &∀ ∈ ≤ k( \Y,X\)\>\>∃ :Y X ∪ \\ & Y X . aligned (11) Its contrapositive then states that ∀ ⊆ ∖Y,X∀ \Y,X\, Y X∣ ⟹∃ ∈ ≤k( ∖Y,X): Y,X . aligned &Y X \\ & ∃ ∈ ≤ k( \Y,X\): Y,X . aligned (12) However, the definition in Equation˜12 allows k to be larger than it needs to be—i.e. k is not strict. Therefore, we say that is strictly k-order faithful to if ∀ ⊆ ∖Y,X∀ \Y,X\, Y X∣ ⟹∃ ∈ k′ ∖Y,X: Y,X , aligned &Y X \\ & ∃ ∈ k \Y,X\: Y,X , aligned (13) where k′≤k ≤ k. This definition in Equation˜13 logically follows from Equation˜12 as a result of Definition˜5. With this new generalisation of faithfulness, we will now present an initial algorithm, , that makes use of our k-order faithfulness assumption for MBsdiscovery. 5 kOMB: k-Order Markov Blanket Our algorithm for Markov blanket (MB)discovery, , is a modified version of the Grow and Shrink (GS)algorithm by margaritis1999. In GS, a candidate MB is iteratively grown by testing the conditional independence between Y and variables that are not yet in —conditioned the current candidate MB . Therefore, the separating set used in the conditional independence tests change as the GSalgorithm progresses. In , something similar occurs where we will only consider subsets of the candidate MBwhen trying to find subsets that successfully cause Y and X to be conditionally independent. In other words, instead of considering the dependence between Y and X over all the variables in , Y,X Y,X ; we will instead consider the association between Y and X over just the current candidate MB , Y,X Y,X . As a first step, we show that a k-order dependence between Y and X with no separating variables to consider implies that Y is k-order associated with every one of the involved variables. restatable [Marginal k-Order Dependence Implies Inter-Association]lemmakdepinter Given a PGM =( , ) =( , ), let X,Y∈ ,Y∈ be distinct variables and ⊆ . Then if Y and X are | || |-order dependent w.r.t. over the empty separating set, Y is | || |-order associated with every variable in ∪X ∪ X, ∅Y,X ⟹ ∅Y, ∪X, Y,X Y, ∪ X, (14) where Y, :=∀W∈ : Y,W ∖W. aligned Y, :=∀ W∈ : Y,W W. aligned (15) When Equation˜14 is true, we say that Y is | || |-order inter-associated with X∪ ∪ . Section˜5 requires no faithfulness assumption—it holds for every distribution via the semi-graphoid axioms—and generalises the 22-association property of unfaithful triples in marx2021 to arbitrary orders; the noisy-parity Example˜1 is exactly this case with =X2,X3 =\X_2,X_3\. It motivates the make_strict routine of Algorithm˜2, which prunes a found dependence set towards an inter-associated one. The last consideration we need to make for to be remotely feasible—especially in cases where the true MBof Y∈ ∈ in is large—is that, we need to somehow limit the size of the conditioning set when conducting the conditional independence (CI)tests in . When testing if Y and X are k-order associated, we will generally use CItests of the form, ∀ ′⊆ :Y X∣ ∪ ′∀ :Y X ∪ , where is the “dependence” set that renders Y and X conditionally dependent. Fortunately, assuming k-order faithfulness, we know that the size of is bounded by k. However the same cannot be said for the set ′ as it can be as large as the candidate MB –which can be arbitrarily large depending on . Therefore, the third and last assumption we will make is that if Y and X are conditionally independent given ∪ ∪ , then there must exist some subset with size less than l, ′∈ l ∈ l , such that Y and X are still conditionally independent given ∪ ′ ∪ . Assumption 3 (l-Bounded Separator Assumption). For a PGM =( , ) =( , ), distinct variables Y,X∈ ,X∈ , and ′≔ , let us first recall the k-order faithfulness assumption as defined in Equation˜11 ∀ ⊆ ∀ , ∀ ∈ ≤k( ′)∃ ′⊆ :Y X∣ ∪ ′⟹Y X∣ . aligned &∀ ∈ ≤ k( )\>\>∃ :Y X ∪ \\ & Y X . aligned (11) Then, for some natural number l∈ℕl , we define an l-bounded separator version of the k-order faithfulness assumption ∀ ⊆ ∀ , ∀ ∈ ≤k( ′)∃ ′∈ l :Y X∣ ∪ ′⇔∀ ∈ ≤k( ′)∃ ′⊆ :Y ∣ ∪ ′⟹Y X∣ , aligned &∀ ∈ ≤ k( )\>\>∃ ∈ l :Y X ∪ \\ & ∀ ∈ ≤ k( )\>\>∃ :Y X ∪ \\ & Y X , aligned (16) which essentially states that as long as Y and X are conditionally independent given and every ≤l≤ l-cardinality subset of , then Y and X are conditionally independent given and every possible subset of —regardless of cardinality. The contrapositive of this assumption is then Y X∣ ⟹∃ ∈ ≤k( ′)∀ ′⊆ :Y X∣ ∪ ′⇔∃ ∈ ≤k( ′): ,X ⇔∃ ∈ ′: ,X , aligned &Y X \\ & ∃ ∈ ≤ k( )\>\>∀ :Y X ∪ \\ & ∃ ∈ ≤ k( ): Y,X \\ & ∃ ∈ k : Y,X , aligned (17) where Y,X ≔∀ ′∈ l :Y X∣ ∪ ′ Y,X ∀ ∈ l :Y X ∪ . The concept of bounding the maximum size of the separators in a graph has been explored before. For instance, soh2019 defined an undirected graph to be weakly K-separable if for any two non-adjacent vertices in , there exists a set of vertices with cardinality ≤K≤ K that separates the two vertices in . However, instead of focusing on separation in the graph , our l-bounded separator assumption focuses on the cardinality of “separating” sets that causes conditional independencies between variables in . The l-bounded separator assumption is also analogous to assumptions that existing constraint-based MBdiscovery methods make about the size of the conditioning sets needed to decide conditional independence reliably from finite samples [margaritis1999, tsamardinos2003, aliferis2003, tsamardinos2003b]. Assumption˜3 simply makes such a bound explicit. Furthermore, as shown later in Appendix˜F, the bound is mild on the networks used in Section˜6.2. 5.1 Theoretical Algorithm Input: , , Y, α, k, l Output: (Y) (Y) ←∅ ← ; // Candidate Markov Blanket 1 repeat 2 ′← ∖(Y∪ ) ← (Y∪ ); 3 ′← ← find_inter_dep(Y, ∅ , , ′ , k, l); 4 ← ∪ ′ ← ∪ 5until ′=∅ = ; 6while ∃X∈ :find_cond(Y,X,∅, ∖X,k,l)=⊥∃ X∈ :find\_cond(Y,X, , X,k,l)= do 7 ← ∖X ← X; 8 9return ; Algorithm 1 kOMB As a proof-of-concept, in this section we propose —a modification of the Grow and Shrink (GS)by margaritis1999 to discover Markov blankets with higher-order relationships under the following assumptions: 1) the Global Markov Property from Assumption˜1, 2) k-order faithfulness as in Assumption˜2, and 3) the l-bounded separator assumption in Assumption˜3. The original GSalgorithm starts with an empty MB =∅ = and iteratively adds a variable X into if the conditional independence (CI)test Y X∣ X indicates that Y is conditionally independent to X given . The main difference between and the GSalgorithm then is: 1) considering the addition of entire sets of variables X∪ : ∈ k ∖Y,X∈ ∖(Y∪ ∪ )X∪ : ∈ k Y,X∈ (Y∪ ∪ ) at each iteration to find higher-order parity-type relationships between variables, and 2) conducting CItests using just subsets of the current candidate MB in the condition to ensure that these CItests remain feasible even when the cardinality of is large. As such, we break up into 3 sub-algorithms: Algorithm 1 is just the entry-point function that implements the overall logic of the GSalgorithm, with the main modifications and additions located in the other two algorithms, Algorithm 2 is overall responsible for finding some subset of the variables not already in the candidate MB, ∈ k ∖(X∪Y∪ ) ∈ k (X∪ Y∪ ), that might cause Y and X to be conditionally dependent given and some subset ′∈ k−| | ∈ k-| | of the current candidate MB, finally, Algorithm 3 is then tasked with finding such a ′∈ k−| | ∈ k-| | that causes Y Y∣ ∪ ′Y Y ∪ and checking if there are any l-bounded subsets of the current candidate MB, ′∈ l ∈ l , such that Y Y∣ ∪ ′∪ ′Y Y ∪ ∪ . Keeping this outline in mind, we now present the full algorithm for and show that it correctly discovers the MBof some Y∈ ∈ assuming that Assumptions˜1, 2 and 3 hold. The proofs for any subsequent Theorems can be found in Appendix˜C. 1Function make_strict(Y, X, , , k, l): 2 for Z∈ ∈ do 3 if ¬ no_seps(Y, Z, X∪ ∖ZX∪ Z, , l) then 4 remove Z from ; 5 6 7 return X∪ ∪ ; 8 9Function find_inter_dep(Y, , , ′ , k, l): 10 if | |>k| |>k then 11 return ∅ ; 12 13 for X∈ ∖(Y∪ ∪ )X∈ (Y∪ ∪ ) do 14 ← _cond (Y, X, , , k, l); 15 if ≠⊥ ≠ then 16 return make_strict(Y, X, , , k, l); 17 18 19 for X∈ ′∖(Y∪ ∪ )X∈ (Y∪ ∪ ) do 20 ′← ′∖X ← X; 21 ′← ← find_inter_dep( Y, ∪X ∪ X, , ′ , k, l); 22 if ′≠∅ ≠ then 23 return ′ ; 24 25 26 return ∅ ; 27 Algorithm 2 Association Mining restatable [Correctness of Algorithm˜3]theoremcitest Given an estimated distribution Phys.Rev.E Phys. 10000\ Rev. 10000\ E, distinct variables Y,X∈ ,X∈ , disjoint subsets , ⊂ , ⊂ , and the natural numbers k,l∈ℕk,l ; the function find_cond in Algorithm˜3 returns the set ′∪ ∪ if and only if ∃ ′∈ k−| | ∃ ∈ k-| | where, ∀ ′∈ ∖ ′:Y .Rev.EX∣ ∪ ′∪ ′.∀ ∈ l :Y Phys. 10000\ Rev. 10000\ EX ∪ ∪ . (18) Otherwise find_cond returns the failure value ⊥ , distinct from the empty set, which it may return on success when ′= =∅ = = . restatable [Correctness of Algorithm˜2]theoremassocmine Given an estimated distribution Phys.Rev.E Phys. 10000\ Rev. 10000\ E, the target variable Y∈ ∈ , the current candidate MB ⊂ ∖Y ⊂ Y, the remaining variables ′= ∖(Y∪ ) = (Y∪ ), and the natural numbers k,l∈ℕk,l ; the call find_inter_dep(Y,∅, , ′,k,l)(Y, , , ,k,l) returns the empty set if and only if ∀X∈ ′∀ ∈ k ∖Y,X∃ ′∈ l ∖ : ∀ X∈ \>\>∀ ∈ k \Y,X\\>\>∃ ∈ l : (19) Y Phys.Rev.EX∣ ∪ ′. Y Phys. 10000\ Rev. 10000\ EX ∪ . Otherwise it returns a nonempty set X∪ ∪ _R with X∈ ′X∈ and ⊆ _R for some ∈ k ∖Y,X ∈ k \Y,X\ satisfying the estimated l-bounded association ∀ ′∈ ∖ :Y .Rev.EX∣ ∪ ′∀ ∈ l :Y Phys. 10000\ Rev. 10000\ EX ∪ . 1Function no_seps(Y, X, , , l): 2 for ′∈ l ∖ ∈ l do 3 if Y X∣ ∪ ′Y X ∪ then 4 return False; 5 6 7 return True; 8 9Function find_cond(Y, X, , , k, l): 10 for ′∈ k−| | ∈ k-| | do 11 ← ′∪ ← ∪ ; 12 if Y X∣ X then 13 if no_seps(Y, X, , , l) then 14 return ; 15 16 17 18 return ⊥ ; 19 Algorithm 3 Practical Conditional Independence Test restatable [Correctness of Algorithm˜1]theoremalgkomb Given an estimated distribution Phys.Rev.E Phys. 10000\ Rev. 10000\ E over variables , the target variable whose graphical Markov blanket we want to discover Y∈ ∈ , and the natural numbers k,l∈ℕk,l ; the grow phase of then grows the candidate MB such that, (Y)⊆ . (Y) . (20) If, in addition, every true blanket member admits an in-blanket dependence witness—i.e. for every X∈ (Y)X∈ (Y) and every with (Y)⊆ ⊆ ∖Y (Y) Y there is a witness ∈ k (Y)∖X ∈ k (Y) X with ∖(X∪ )Y,X (X∪ )Y,X —then the shrink phase of removes variables from until, (Y)= , (Y)= , (21) where it subsequently terminates and returns . Beyond Assumptions˜1, 2 and 3, the recovery guarantee of Algorithm˜3 assumes the CIdecisions are sufficiently accurate—i.e. they reflect the true conditional independencies of . 5.2 Computational Complexity of As is a proof-of-concept, we prioritised correctness and clarity over efficiency. It is nonetheless instructive to characterise its worst-case time complexity, which we state in Section˜5.2 with proof given in Appendix˜D. restatable [Computational Complexity of ]corollarykombcomplexity Let n=| |n=| | denote the number of variables, N the sample size, and K the maximum size of any candidate MB encountered by Algorithm˜1 during execution. Furthermore assume each CItest runs in O(N)O(N) time. Then, for a given order k and separator bound l, (Algorithm˜1) terminates in worst-case time O(Kn∑i=0k(ni)[∑i′=0k−i(Ki′)∑j=0l(K−i′j)N]),O\! (K\,n _i=0^k ni [ _i =0^k-i Ki _j=0^l K-i j\,N ] ), (22) where we adopt the convention that (nm)=0 nm=0 whenever n<mn<m. Treating k and l as constants, Equation˜22 is bounded by O(nk+1Kk+l+1N)O(n^k+1\,K^k+l+1\,N)—polynomial in n, K, and N, and for k=0k=0 it reduces to a scan that is linear in n, recovering the original GSalgorithm. Crucially, the search for a dependence-exposing subset and its separator is governed by K rather than n, a direct consequence of the l-bounded separator assumption (Assumption˜3). The remaining exponential dependence on k and l motivates keeping them small. The ablation study in Appendix˜G showcases how differing values of k and l affects Algorithm˜1 runtime and its ability to find high-order dependencies. Generally, from the results in Appendix˜G for the benchmark datasets used, we found that setting k=lk=l and having k≤2k≤ 2 provides a decent trade-off between runtime and finding high-order dependencies. Table 1: Results on Synthetically Generated Data Method and ex-1 ex-2 or parity GS 0.410 0.313 0.288 0.365 0.025 ± 0.405 ± 0.419 ± 0.407 ± 0.374 ± 0.110 kOMB 0.410 0.283 0.268 0.350 0.025 (k=0,l=3)(k=0,l=3) ± 0.405 ± 0.372 ± 0.374 ± 0.351 ± 0.110 kOMB 0.618 0.675 0.675 0.560 0.050 (k=1,l=3)(k=1,l=3) ± 0.374 ± 0.474 ± 0.474 ± 0.435 ± 0.221 kOMB 0.618 1.000 1.000 0.710 1.000 (k=2,l=3)(k=2,l=3) ± 0.374 ± 0.000 ± 0.000 ± 0.384 ± 0.000 6 Experiments In order to empirically test we implemented the algorithm as a Python package written in Rust. The implementation uses the G-Test with α=0.01α=0.01 as its conditional independence test of choice when possible and falls back to the SCI test [marx2019] when G-test is too weak for the test being conducted—i.e. when the average frequency for each cell is <5<5. A public repo for our implementation of and code for the experiments can be found at https://github.com/lklee9/k-order-Markov-blanket. All experiments were run on an Apple M4 Mac mini. Throughout this section we will compare against existing methods for Markov blanket (MB)discovery. Specifically we will compare against the following eight MBdiscovery methods: GS [margaritis1999], IAMB [tsamardinos2003], HITON [aliferis2003], MMMB [tsamardinos2003b], PCMB [pena2007], LRH [liu2016], STMB [gao2017a], BAMB [ling2019]. We will use the implementation of these methods found in the pyCausalFS python library [yu2020a] with a significance level of 0.010.01 for any statistical tests as well. When comparing between different MBdiscovery methods, we will use the F1 score defined as: F1=2×/(+)F1=2×precision×recall/(precision+recall), which is the harmonic mean of precision and recall. Here precision denotes the ratio between the number of true positives in the inferred MBand the size of the inferred MB. On the other hand, Recall denotes the ratio between the number of true positives in the inferred MBand the size of the true MB. Therefore F1=1F1=1 represents the perfect precision and recall. Figure 4: Precision, recall, and F1 score of GSand ( -k-l) when recovering the MBsof the synthetic BNsin Section˜6.1, over sample size. Each row corresponds to a metric and each column to a BN. Table 2: F1 score and runtime (in seconds) on benchmark datasets. Best F1 and fastest runtime per dataset in bold. F1 Score Runtime (s) Method Alarm1 Barley Insurance Mildew Alarm1 Barley Insurance Mildew BAMB 0.6742 0.3385 0.6588 0.4945 3.517 1.149 5.615 0.797 ± 0.1771 ± 0.1226 ± 0.0797 ± 0.1084 ± 4.305 ± 0.681 ± 5.358 ± 0.392 GS 0.3434 0.2113 0.4615 0.2867 4.875 1.575 4.433 1.035 ± 0.1503 ± 0.0807 ± 0.1477 ± 0.1581 ± 5.183 ± 0.606 ± 3.792 ± 0.341 HITON_MB 0.7636 0.3400 0.6746 0.4918 2.560 1.027 3.514 0.604 ± 0.1861 ± 0.1284 ± 0.1594 ± 0.1487 ± 1.277 ± 0.298 ± 2.525 ± 0.218 MMMB 0.7692 0.3298 0.6978 0.4861 2.724 1.090 3.508 0.553 ± 0.1874 ± 0.1274 ± 0.1405 ± 0.1810 ± 1.278 ± 0.374 ± 2.525 ± 0.128 kOMB 0.7804 0.5158 0.7136 0.6220 5.258 58.190 4.856 17.698 (k=1,l=1)(k=1,l=1) ± 0.1145 ± 0.1964 ± 0.1328 ± 0.1595 ± 2.921 ± 31.619 ± 2.041 ± 10.205 kOMB 0.8209 0.3128 0.7336 0.4360 18.144 335.821 14.493 121.909 (k=2,l=2)(k=2,l=2) ± 0.1914 ± 0.1582 ± 0.1624 ± 0.1394 ± 7.566 ± 135.966 ± 7.280 ± 33.153 6.1 Experiments with Synthetic Bayesian Networks To systematically test how behaves in the face of true and empirical faithfulness violations, we will task at learning the MBof every variable for five different basic toy Bayesian networks. These problems will have 44 variables— =X1,X2,X3 =\X_1,X_2,X_3\ and Y—with a DAGthat has edges ∀i∈[3]:Xi→Y∀ i∈[3]:X_i→ Y. Each variable X∈ ∈ is an independent Bernoulli random variable with p=0.5p=0.5. The conditional distribution of Y is then (Y=f( )∣ )=0.9 (Y=f( ) )=0.9 where f( )f( ) is one of 55 different Boolean functions, one for each toy problem which we describe in Appendix˜E. We then sample 100100 samples from the 55 Bayesian networks 1010 different times, leading to 10×510× 5 samples of size 100100—hence the size of each sample is much larger than the cardinality of the domain over all variables, 100>>24100>>2^4. We then tasked each method to infer—for each Bayesian network—the MBof all 44 variables from the 1010 different samples separately. We then averaged the F1 score over the 10×410× 4 runs and variables for each Bayesian network, resulting in Table˜1. Note that we only include the results for and Grow and Shrink (GS)because GSwas the best performing method out of all the previous methods mentioned in Section˜6 for these synthetic datasets. From the results in Table˜1, we can observe that consistently outperforms the other approaches at MBdiscovery. More interestingly—even for the problems whose Bayesian networks are faithful—as we increase k and look at higher order relationships, achieves a better F1 score and therefore is capable of better overcoming empirical unfaithfulness. Furthermore the MBfor the parity problem is only recoverable when k=2k=2, which confirms ’s ability to find these higher-order parity-like relations. In order to better understand the source of these empirical unfaithfulness in the low sample size regime, Figure˜4 traces the precision, recall, and F1 score as the sample size grows from 5050 to 100100. As expected only -2-2 is capable of recovering the higher-order dependence of the parity network regardless of sample size. From the other faithful networks, we can observe that at small sample size, all methods exhibit lower recall and precision which indicates dependencies missed by the MBdiscovery methods and spurious dependencies erroneously detected by said methods respectively. That said, we can also observe that for , the higher the value of k used, the greater its performance in the low sample regime and the quicker it recovers the true MBsfor all variables. 6.2 Experiments with Benchmark Bayesian Networks To further explore how behaves in more realistic scenarios, we then repeated the same procedure as Section˜6.1, but with some well-known benchmark datasets444The datasets can be accessed at: https://pages.mtu.edu/~lebrown/supplements/mmhc_paper/mmhc_index.html in the MBdiscovery literature. The main difference in the experiments in this section is that we used samples of 50005000 and we only tasked each method at finding the MBfor the 1010 variables with the largest neighborhoods in the Bayesian network. Furthermore, we also omit some of the poorer performing methods in Table˜2 to conserve space. From the results in Table˜2 we can observe that does generally outperform the other existing approaches to MBdiscovery. Whether a higher order helps, however, depends on the network’s cardinality. On the lower-cardinality Alarm1 and Insurance networks with order 22 improves on order 11, whereas on the high-cardinality Barley and Mildew networks order 11 is best and order 22 is substantially worse—there the larger conditioning sets leave the higher-order CItests underpowered at this sample size. Runtime, in turn, grows with both the order k and the cardinality of the network. Specifically, the baselines finish within roughly 66 s per dataset, while ranges from comparable at order 11 on the lower-cardinality networks (Alarm1, Insurance) to a few hundred times slower at order 22 on the high-cardinality networks (Barley, Mildew). The full results and runtimes for every method and every (k,l)(k,l) configuration are given in Tables˜5 and 6 under Appendix˜G. 7 Conclusion We introduced a k-order relaxation of faithfulness that explicitly accounts for parity-type dependencies among k+2k+2 variables, and used it to derive a proof-of-concept Markov blanket (MB)discovery algorithm, ; with theoretical guarantees that recovers the graphical MBunder the proposed assumptions. Empirically, overcomes both true and empirical faithfulness violations on synthetic problems, and it performs competitively on commonly used benchmark datasets. As is a proof-of-concept demonstrating that exploiting higher-order dependencies can improve MBdiscovery, an important direction for future work is the development of approximate algorithms that uncover these higher-order dependencies without ’s worst-case cost—for example, through more intricate pruning rules that make scale better with larger values of k. Secondly, although our ablation (Appendix˜G) shows that l=kl=k with k≤2k≤ 2 is a robust default, principled a priori selection of k and l without expert knowledge remains open. A practical alternative is to run for several values of k and l while memoising results across runs to speed up the total runtime. Finally, our evaluation compares against established constraint-based MBdiscovery methods using the G-test (with an SCI fallback) for conditional independence. Therefore, utilising more recent alternative CItests is a further avenue for future work. Acknowledgements.This research has been funded by the Federal Ministry of Research, Technology and Space of Germany and the state of North Rhine-Westphalia as part of the Lamarr Institute for Machine Learning and Artificial Intelligence. References High-Order Markov Blanket Discovery via a k-Order Relaxation of the Faithfulness Assumption (Supplementary Material) Appendix A Graphoid Axioms Definition 6 (Graphoid Axioms). For disjoint subsets , , , ⊆ , , , , a semi-graphoid independence model ⋅ · obeys the following properties [dawid1979, pearl1987]: 1. Symmetry: ∣ ⟹ ∣ 2. Decomposition: ∪ ∣ ⟹ ∣ ∪ 3. Weak Union: ∪ ∣ ⟹ ∣ ∪ ∪ ∪ 4. Contraction: ( ∣ ∪ )∧( ∣ )⟹ ∪ ∣ ( ∪ ) ( ) ∪ The independence model is called a graphoid if it also obeys the following additional property [pearl1985]: 5. Intersection: ( ∣ ∪ )∧( ∣ ∪ )⟹ ∪ ∣ ( ∪ ) ( ∪ ) ∪ All graph independence models are guaranteed to obey all the graphoid axioms [lauritzen2018]. However, probabilistic independence models are only guaranteed to obey the semi-graphoid axioms. However if the joint distribution is strictly positive, then is guaranteed to be a graphoid as well [pearl1988, pearl1989]. Appendix B Proofs * Proof. Consider the family of association witnesses contained in , ≔ ⊆ : Y,X .W \ : Y,X \. By the hypothesis Y,X Y,X we have ∈ , so W is nonempty; and since is finite, so is W. Hence W contains an element ′ of minimum cardinality. Every proper subset ∗⊊ ′ ^* satisfies | ∗|<| ′|| ^*|<| | and ∗⊆ ^* , so by the minimality of ′ we have ∗∉ ^* , i.e. ¬ Y,X ∗ \, Y,X ^*. Together with Y,X ′ Y,X , this is by Equation˜7 exactly Y,X ′ Y,X . As ′⊆ , this establishes Equation˜10. ∎ * Proof. Assume ∅Y,X Y,X . For the variable X, the inter-association clause is ∅Y,X Y,X , i.e. Y X∣ X , which is exactly the first conjunct of Equation˜7 and hence immediate. Now suppose, towards a contradiction, that the clause fails for some Z∈ ∈ , i.e. Y Z∣X∪( ∖Z)Y Z X∪( Z). Since ∖Z Z is a proper subset of and the separating set is empty, the second conjunct of Equation˜7 gives Y X∣ ∖ZY X Z. Applying the contraction and then weak union graphoid axioms of Definition˜6, (Y Z∣X∪( ∖Z))∧(Y X∣ ∖Z) (Y Z X∪( Z)) (Y X Z) (23) ⟹Y X,Z∣ ∖Z⟹Y X∣ , Y \X,Z\ Z Y X , contradicting Y X∣ X . Hence Y Z∣X∪( ∖Z)Y Z X∪( Z) for every Z∈ ∈ , which together with Y X∣ X is precisely ∅Y, ∪X Y, ∪ X. ∎ Appendix C Correctness of * Proof. The function find_cond first iterates through all the subsets ′∈ k−| | ∈ k-| | and for each ′ where Y Phys.Rev.EX∣ ∪ ′Y Phys. 10000\ Rev. 10000\ EX ∪ (line 9 in find_cond), it passes ′ to the function no_seps (line 10 in find_cond) to determine if there are any subsets ′∈ l ∖ ′ ∈ l that results in Y Phys.Rev.EX∣ ∪ ′∪ ′Y Phys. 10000\ Rev. 10000\ EX ∪ ∪ (line 4 in no_seps). As soon as no_seps finds such a ′ , it returns False (line 4 in no_seps), and find_cond will need to go to the next iteration. Otherwise, if no_seps is unable to find any such subset ′ , it returns True (line 5 in no_seps), and find_cond then immediately returns = ′∪ = ∪ (line 11 in find_cond) as there are no l-bounded subsets ′∈ l ∖ ∈ l that render Y and X conditionally independent given ∪ ′ ∪ . Here we use that the external set is disjoint from , so ∖ = ∖ ′ = ; and since ∅∈ l ∖ ′ ∈ l , the dependence check Y Phys.Rev.EX∣ Phys. 10000\ Rev. 10000\ EX on line 9 is precisely the ′=∅ = instance tested by no_seps, so a returned satisfies Equation˜18 in full. However, if find_cond completes all its iteration without returning early, then we know that all subset ′ has some l-bounded subset ′ that will cause Y Phys.Rev.EX∣ ′∪ ∪ ′Y Phys. 10000\ Rev. 10000\ EX ∪ ∪ . Therefore find_cond returns the failure value ⊥ in this case (line 12 in find_cond). ∎ * Proof. First consider a call that returns a nonempty set. find_inter_dep does so only by returning make_strict(Y,X, , ,k,l)(Y,X, , ,k,l) after a call find_cond(Y,X, , ,k,l)(Y,X, , ,k,l) has returned some ≠⊥ ≠ , for some X∈ ∖(Y∪ ∪ )X∈ (Y∪ ∪ ) (as is passed unchanged through every recursive call, any such X satisfies X∈ ∖(Y∪ )= ′X∈ (Y∪ )= , matching the theorem’s candidate set). By Algorithm˜2, = ∪ ′ = ∪ for some ′∈ k−| | ∈ k-| | such that Y Phys.Rev.EX∣ ∪ ′Y Phys. 10000\ Rev. 10000\ EX ∪ holds for all ′∈ l ∖ ′ ∈ l . As the dependence set is disjoint from , we have ∖ ′= ∖ = , so Y Phys.Rev.EX∣ ∪ ′Y Phys. 10000\ Rev. 10000\ EX ∪ holds for all ′∈ l ∖ ∈ l with ∈ k ∖Y,X ∈ k \Y,X\—the estimated l-bounded association claimed in the theorem. Since make_strict only removes variables from , the returned set is X∪ ∪ _R for some ⊆ _R , as claimed. Next consider a call that returns the empty set. Here every internal find_cond call must have failed, so find_inter_dep traverses the full enumeration tree of the dependence sets ∈ k ∖(Y∪ ) ∈ k (Y∪ ) (lines 13–17 of find_inter_dep), and for each visits every candidate X∈ ∖(Y∪ ∪ )X∈ (Y∪ ∪ ) through find_cond(Y,X, , ,k,l)(Y,X, , ,k,l), which by Algorithm˜2 searches all ′∈ k−| | ∈ k-| | . Every conditioning set ∈ k ∖Y,X ∈ k \Y,X\ decomposes uniquely into its part outside and its part inside the candidate blanket, =( ∖ )∪( ∩ )= ∪ ′ =( )∪( ∩ )= ∪ , so this enumeration covers every pair (X, )(X, ) with X∈ ′X∈ and ∈ k ∖Y,X ∈ k \Y,X\. By Algorithm˜2, the call for (X, )(X, ) returns ⊥ exactly when every such ′ admits a separator ′∈ l ∖ ′ ∈ l with Y Phys.Rev.EX∣ ∪ ′∪ ′Y Phys. 10000\ Rev. 10000\ EX ∪ ∪ . Therefore find_inter_dep returns the empty set if and only if Equation˜19 holds. The distinguished value ⊥ lets the callers tell such a failure from a success that returns the empty conditioning set. ∎ * Proof. Throughout we assume that the CIdecisions are accurate, so Phys.Rev.E Phys. 10000\ Rev. 10000\ E and agree on every tested statement. In the grow phase, every nonempty set returned by find_inter_dep contains a variable X∉ ∉ by Algorithm˜2, so each iteration of the repeat loop strictly grows and the phase terminates. Upon termination find_inter_dep returned the empty set, so by Algorithm˜2 the condition in Equation˜19 holds in Phys.Rev.E Phys. 10000\ Rev. 10000\ E, and by accuracy in . Fix any X∈ ∖( ∪Y)X∈ ( ∪ Y). For every ∈ k ∖Y,X ∈ k \Y,X\, Equation˜19 supplies an ′∈ l ∖ ⊆ l ∈ l l with Y X∣ ∪ ′Y X ∪ ; the l-bounded separator assumption (Assumption˜3) then lifts this, for every such , to the unbounded statement that some ′⊆ gives Y X∣ ∪ ′Y X ∪ , which is the antecedent of k-order faithfulness (Equation˜11). Hence Y X∣ X for every X∈ ∖( ∪Y)X∈ ( ∪ Y). Since the parents and children of Y—and, for an undirected , its neighbours—can never be separated from Y, they must all lie in . Any spouse X shares a child C∈ch(Y)⊆ (Y) , so the collider path Y→C←XY→ C← X is active given ; hence Y X∣ X , which forces X∈ ∈ as well. Therefore (Y)⊆ (Y) , establishing Equation˜20. During the shrink phase we maintain the invariant (Y)⊆ (Y) , which holds on entry by Equation˜20. First, we show that no true member is ever removed. For X∈ (Y)X∈ (Y), the in-blanket witness hypothesis of Algorithm˜3 supplies a set ∈ k (Y)∖X⊆ k ∖X ∈ k (Y) X k X such that Y X∣ ∪ ′Y X ∪ holds for all ′∈ l( ∖X)∖ ∈ l( X) . By accuracy and Algorithm˜2 (with candidate blanket ∖X X, external dependence set =∅ = , and ′= = ), the call find_cond(Y,X,∅, ∖X,k,l)(Y,X, , X,k,l) does not return ⊥ , so X is never selected for removal (a member whose only witness is =∅ = makes find_cond return the empty set rather than ⊥ ). Next, we show that every non-member is removed once selected. Let X∈ ∖ (Y)X∈ (Y); by the invariant (Y)⊆ ∖X (Y) X. We first observe that (Y) (Y) separates Y from X even after conditioning on any further set ⊆ ∖Y,X \Y,X\. Since X∉ (Y)X∉ (Y), X is not adjacent to Y, so every path between them has at least one interior vertex; write such a path as Y=N0∼N1∼⋯∼Nm=XY=N_0 N_1 ·s N_m=X with m≥2m≥ 2. Its first interior vertex N1N_1 is adjacent to Y and hence lies in (Y) (Y). If N1N_1 is a non-collider on the path, it blocks the path. If instead N1N_1 is a collider Y→N1←N2Y→ N_1← N_2, then N1N_1 is a common child of Y and its successor N2N_2 on the path, so N2N_2 is a spouse of Y and hence N2∈ (Y)N_2∈ (Y). In particular N2≠XN_2≠ X (as X∉ (Y)X∉ (Y)), so N2N_2 is a genuine interior vertex; and since the edge N2→N1N_2→ N_1 points out of N2N_2, it is a non-collider on the path and blocks it. For an undirected , N1N_1 is instead a neighbour of Y and blocks the path directly. As each path is thus blocked at a conditioned vertex of (Y) (Y) irrespective of [pearl1988, lauritzen1996], Y X∣ (Y)∪ X (Y)∪ for every ∈ k ∖Y,X ∈ k \Y,X\. The Global Markov Property transfers this to with the separator ′≔ (Y)∖ ⊆ ∖X (Y) X, so the unbounded side of Assumption˜3—at the separating set ∖X X—holds for every such . The biconditional then yields some ′∈ l ∖X ∈ l X; discarding any overlap with leaves ∪ ′ ∪ unchanged and gives ′∈ l( ∖X)∖ ∈ l( X) with Y X∣ ∪ ′Y X ∪ . By accuracy and Algorithm˜2, find_cond(Y,X,∅, ∖X,k,l)(Y,X, , X,k,l) returns ⊥ , so X is removed once selected. Each iteration of the while loop removes one variable, so the phase terminates; it exits only when no X∈ ∈ yields ⊥ , whence ∖ (Y)=∅ (Y)= and, with the invariant, = (Y) = (Y), establishing Equation˜21. ∎ Appendix D Complexity of * Proof. We charge O(N)O(N) to each CItest and bound the running time by analysing the three sub-algorithms of Section˜5 from the innermost outwards. First consider find_cond (Algorithm˜3) for a fixed candidate variable X and external dependence set with | |=i| |=i. It iterates over every subset ′∈ k−i ∈ k-i , of which there are at most ∑i′=0k−i(Ki′) _i =0^k-i Ki since | |≤K| |≤ K. For each ′ it performs one CItest and, whenever that test indicates a dependence, calls no_seps, which iterates over every separator ′∈ l ∖ ′ ∈ l and performs one further CItest per separator. Writing | ′|=i′| |=i , there are at most ∑j=0l(K−i′j) _j=0^l K-i j such separators, so a single call to find_cond costs O(∑i′=0k−i(Ki′)∑j=0l(K−i′j)N),O\! ( _i =0^k-i Ki _j=0^l K-i j\,N ), (24) where the convention (nm)=0 nm=0 for n<mn<m accounts for the separator set ∖ ′ being exhausted. Next consider find_inter_dep (Algorithm˜2). It explores the enumeration tree of ∖Y Y up to subsets of cardinality k (lines 13–17 of find_inter_dep), and hence visits at most ∑i=0k(ni) _i=0^k ni dependence sets ∈ k ∖Y ∈ k Y (a loose bound on the exact count ∑i=0k(n−1i) _i=0^k n-1i, as | ∖Y|=n−1| Y|=n-1). For every such it loops over the O(n)O(n) candidate variables X∈ ∖(Y∪ ∪ )X∈ (Y∪ ∪ ) and invokes find_cond. Multiplying these two factors by Equation˜24 bounds the cost of a single association-mining sweep by O(n∑i=0k(ni)[∑i′=0k−i(Ki′)∑j=0l(K−i′j)N]).O\! (n _i=0^k ni [ _i =0^k-i Ki _j=0^l K-i j\,N ] ). (25) Each successful sweep additionally invokes make_strict exactly once; it iterates over the found dependence set —of size at most min(k,n−2) (k,n-2), since ⊆ ∖Y,X \Y,X\ and | |≤k| |≤ k—and performs a single no_seps call per element, at cost O(min(k,n−2)∑j=0l(Kj)N)O\! ( (k,n-2) _j=0^l Kj\,N ). This is dominated by the per-sweep find_cond total already counted in Equation˜25—whose bracketed factor already contains the term ∑j=0l(Kj) _j=0^l Kj (at i′=0i =0) multiplied by the outer factor n∑i=0k(ni)≥n≥min(k,n−2)n _i=0^k ni≥ n≥ (k,n-2)—and is therefore absorbed into the bound. Finally, consider itself (Algorithm˜1). Each successful call to find_inter_dep in the grow phase returns a nonempty set X∪ ∪ that is added to , so the candidate Markov blanket grows by at least one variable per successful call; as | |≤K| |≤ K throughout, the grow phase performs at most O(K)O(K) association-mining sweeps. The shrink phase re-evaluates its guard after each removal, calling find_cond at most O(K2)O(K^2) times over the at most K removals; as K≤n−1K≤ n-1, this cost is dominated by that of the grow phase. Multiplying Equation˜25 by these O(K)O(K) sweeps yields the bound in Equation˜22, completing the proof. ∎ Appendix E Boolean Functions for the Synthetic Bayesian Networks parity: f(X1,…,X3)=(∑i=13Xi)(mod2)f(X_1,…,X_3)=( _i=1^3X_i) 2 ex-s: f(X1,…,X3;s)=[(∑i=13Xi)=s]f(X_1,…,X_3;s)= [( _i=1^3X_i)=s ] and: f(X1,…,X3)=[(∑i=13Xi)=3]f(X_1,…,X_3)= [( _i=1^3X_i)=3 ] or: f(X1,…,X3)=[(∑i=13Xi)≥1]f(X_1,…,X_3)= [( _i=1^3X_i)≥ 1 ] where s∈1,2s∈\1,2\. Appendix F Separator Sizes in the Benchmark Datasets To gauge how restrictive the l-bounded separator assumption (Assumption˜3) is in practice, we examine the benchmark networks of Section˜6.2. We use every variable in turn as a target and, for every pair formed by a target and a variable outside its Markov blanket, compute a minimal d-separator drawn from the true Markov blanket. Table˜3 reports the mean and worst-case (maximum) separator size and the fraction of pairs whose separator has size at most three. Separators are small across all four networks—the mean never exceeds 2.192.19 and size-≤3≤\!3 separators cover at least 83%83\% of pairs—and no pair had to be skipped, i.e. a separator within the true Markov blanket always existed. The bound is loosest on Barley, where the largest separator reaches size 99; this is also the network on which is the most expensive and least accurate (Section˜6.2), illustrating that the difficulty of MBdiscovery tracks the separator sizes the assumption must accommodate. Overall, a modest separator bound l already suffices to capture most of the conditional independencies needed to recover the Markov blanket on these networks. Table 3: Minimal d-separator sizes restricted to the true Markov blanket, computed using every variable of each benchmark network as a target. “Targets” is the number of targets (all variables) and “Pairs” the number of target–non-blanket pairs; “Mean sep. size” and “Max” are the mean and maximum separator size; and “Cover ≤3≤ 3” is the fraction whose minimal separator has size at most three. Network Targets Pairs Mean sep. size Max Cover ≤3≤ 3 Alarm1 37 1202 0.621 3 1.000 Barley 48 2004 1.739 9 0.835 Insurance 27 562 2.190 7 0.925 Mildew 35 1030 1.065 6 0.976 Appendix G Ablation For completeness, we report the full results of the experiments in Section˜6.2 for every method and every (k,l)(k,l) configuration of . Table˜5 gives the performance of all methods, while Table˜6 gives their runtimes in seconds. In both tables, “–” denotes a configuration that failed to complete within the allotted time budget. Table˜4 summarises the F1 score of across all synthetic and benchmark datasets for every (k,l)(k,l) configuration. F1 generally improves from k=0k=0 to k=2k=2—higher orders expose dependencies that lower orders miss—but the gains flatten or reverse beyond k=2k=2 as the larger conditioning sets incur both higher runtime (Table˜6) and more error-prone CItests. Setting l=kl=k with k≤2k≤ 2 therefore offers a robust default, the heuristic we adopt in Section˜5.2. For a fixed k, raising l beyond k does not consistently improve F1 (Table˜4), which motivates simply setting l=kl=k rather than tuning l separately. Table 4: F1 score across datasets for each (k,l)(k,l) configuration of . Method parity ex-1 ex-2 and or Alarm1 Barley Insurance Mildew kOMB 0.0250 0.2825 0.2675 0.4100 0.3500 – – – – (k=0,l=1)(k=0,l=1) ± 0.1104 ± 0.3720 ± 0.3737 ± 0.4050 ± 0.3508 – – – – kOMB 0.0250 0.2825 0.2675 0.4100 0.3500 – – – – (k=0,l=2)(k=0,l=2) ± 0.1104 ± 0.3720 ± 0.3737 ± 0.4050 ± 0.3508 – – – – kOMB 0.0250 0.2825 0.2675 0.4100 0.3500 – – – – (k=0,l=3)(k=0,l=3) ± 0.1104 ± 0.3720 ± 0.3737 ± 0.4050 ± 0.3508 – – – – kOMB 0.0500 0.6750 0.6750 0.6175 0.5600 0.7804 0.5158 0.7136 0.6220 (k=1,l=1)(k=1,l=1) ± 0.2207 ± 0.4743 ± 0.4743 ± 0.3741 ± 0.4349 ± 0.1145 ± 0.1964 ± 0.1328 ± 0.1595 kOMB 0.0500 0.6750 0.6750 0.6175 0.5600 0.8187 0.3128 0.7336 0.4360 (k=1,l=2)(k=1,l=2) ± 0.2207 ± 0.4743 ± 0.4743 ± 0.3741 ± 0.4349 ± 0.1944 ± 0.1582 ± 0.1624 ± 0.1394 kOMB 0.0500 0.6750 0.6750 0.6175 0.5600 0.7923 0.2498 0.7076 0.3726 (k=1,l=3)(k=1,l=3) ± 0.2207 ± 0.4743 ± 0.4743 ± 0.3741 ± 0.4349 ± 0.1885 ± 0.1140 ± 0.1304 ± 0.1558 kOMB 1.0000 1.0000 1.0000 0.6175 0.7100 0.7775 0.4634 0.7249 0.5059 (k=2,l=1)(k=2,l=1) ± 0.0000 ± 0.0000 ± 0.0000 ± 0.3741 ± 0.3842 ± 0.1511 ± 0.1774 ± 0.1464 ± 0.1458 kOMB 1.0000 1.0000 1.0000 0.6175 0.7100 0.8209 0.3128 0.7336 0.4360 (k=2,l=2)(k=2,l=2) ± 0.0000 ± 0.0000 ± 0.0000 ± 0.3741 ± 0.3842 ± 0.1914 ± 0.1582 ± 0.1624 ± 0.1394 kOMB 1.0000 1.0000 1.0000 0.6175 0.7100 0.7967 0.2498 0.7058 0.3726 (k=2,l=3)(k=2,l=3) ± 0.0000 ± 0.0000 ± 0.0000 ± 0.3741 ± 0.3842 ± 0.1828 ± 0.1140 ± 0.1311 ± 0.1558 kOMB – – – – – 0.7721 – 0.7227 0.5243 (k=3,l=1)(k=3,l=1) – – – – – ± 0.1599 – ± 0.1552 ± 0.1075 kOMB – – – – – 0.8001 – 0.7288 0.4194 (k=3,l=2)(k=3,l=2) – – – – – ± 0.1948 – ± 0.1471 ± 0.1414 kOMB – – – – – 0.7967 – 0.7058 0.3726 (k=3,l=3)(k=3,l=3) – – – – – ± 0.1828 – ± 0.1311 ± 0.1558 Table 5: Full results on benchmark datasets. Method Alarm1 Barley Insurance Mildew BAMB 0.6742 0.3385 0.6588 0.4945 ± 0.1771 ± 0.1226 ± 0.0797 ± 0.1084 GS 0.3434 0.2113 0.4615 0.2867 ± 0.1503 ± 0.0807 ± 0.1477 ± 0.1581 HITON_MB 0.7636 0.3400 0.6746 0.4918 ± 0.1861 ± 0.1284 ± 0.1594 ± 0.1487 IAMB 0.6931 0.3373 0.6293 0.4677 ± 0.1827 ± 0.1223 ± 0.0778 ± 0.1224 LRH 0.6813 0.3373 0.6281 0.4674 ± 0.1777 ± 0.1223 ± 0.0782 ± 0.1230 MMMB 0.7692 0.3298 0.6978 0.4861 ± 0.1874 ± 0.1274 ± 0.1405 ± 0.1810 PCMB 0.7371 0.2284 0.6299 0.4156 ± 0.2127 ± 0.1199 ± 0.1291 ± 0.1445 STMB 0.6598 0.2145 0.5648 0.2734 ± 0.1627 ± 0.1281 ± 0.1324 ± 0.1765 kOMB 0.7804 0.5158 0.7136 0.6220 (k=1,l=1)(k=1,l=1) ± 0.1145 ± 0.1964 ± 0.1328 ± 0.1595 kOMB 0.8187 0.3128 0.7336 0.4360 (k=1,l=2)(k=1,l=2) ± 0.1944 ± 0.1582 ± 0.1624 ± 0.1394 kOMB 0.7923 0.2498 0.7076 0.3726 (k=1,l=3)(k=1,l=3) ± 0.1885 ± 0.1140 ± 0.1304 ± 0.1558 kOMB 0.7775 0.4588 0.7249 0.5059 (k=2,l=1)(k=2,l=1) ± 0.1511 ± 0.1825 ± 0.1464 ± 0.1458 kOMB 0.8209 0.3128 0.7336 0.4360 (k=2,l=2)(k=2,l=2) ± 0.1914 ± 0.1582 ± 0.1624 ± 0.1394 kOMB 0.7967 0.2498 0.7058 0.3726 (k=2,l=3)(k=2,l=3) ± 0.1828 ± 0.1140 ± 0.1311 ± 0.1558 kOMB 0.7721 0.0000 0.7227 0.0419 (k=3,l=1)(k=3,l=1) ± 0.1599 ± 0.0000 ± 0.1552 ± 0.1458 kOMB 0.8001 0.0000 0.7288 0.3607 (k=3,l=2)(k=3,l=2) ± 0.1948 ± 0.0000 ± 0.1471 ± 0.1964 kOMB 0.7967 0.0000 0.7058 0.3726 (k=3,l=3)(k=3,l=3) ± 0.1828 ± 0.0000 ± 0.1311 ± 0.1558 Table 6: Full runtimes (in seconds) on benchmark datasets. Method Alarm1 Barley Insurance Mildew BAMB 3.517 1.149 5.615 0.797 ± 4.305 ± 0.681 ± 5.358 ± 0.392 GS 4.875 1.575 4.433 1.035 ± 5.183 ± 0.606 ± 3.792 ± 0.341 HITON_MB 2.560 1.027 3.514 0.604 ± 1.277 ± 0.298 ± 2.525 ± 0.218 IAMB 2.527 1.104 1.650 0.943 ± 1.710 ± 0.496 ± 0.967 ± 0.420 LRH 89.877 5.071 45.401 1.833 ± 93.285 ± 2.496 ± 24.653 ± 0.651 MMMB 2.724 1.090 3.508 0.553 ± 1.278 ± 0.374 ± 2.525 ± 0.128 PCMB 15.515 4.982 27.213 3.066 ± 10.203 ± 1.872 ± 21.695 ± 1.486 STMB 5.305 1.895 14.658 0.561 ± 5.855 ± 2.367 ± 33.396 ± 0.524 kOMB 5.258 58.190 4.856 17.698 (k=1,l=1)(k=1,l=1) ± 2.921 ± 31.619 ± 2.041 ± 10.205 kOMB 3.974 12.369 3.828 7.562 (k=1,l=2)(k=1,l=2) ± 2.468 ± 7.709 ± 2.235 ± 2.301 kOMB 4.111 9.639 3.171 7.384 (k=1,l=3)(k=1,l=3) ± 2.791 ± 3.794 ± 1.707 ± 2.127 kOMB 26.092 788.531 19.901 311.324 (k=2,l=1)(k=2,l=1) ± 14.668 ± 453.059 ± 7.226 ± 167.079 kOMB 18.144 335.821 14.493 121.909 (k=2,l=2)(k=2,l=2) ± 7.566 ± 135.966 ± 7.280 ± 33.153 kOMB 18.016 290.512 10.424 121.146 (k=2,l=3)(k=2,l=3) ± 9.135 ± 85.391 ± 4.290 ± 32.239 kOMB 174.900 1800.000 96.149 1761.246 (k=3,l=1)(k=3,l=1) ± 83.504 ± 0.000 ± 46.913 ± 188.357 kOMB 132.619 1800.000 58.350 1534.465 (k=3,l=2)(k=3,l=2) ± 46.566 ± 0.000 ± 21.745 ± 300.944 kOMB 127.861 1800.000 53.452 1356.514 (k=3,l=3)(k=3,l=3) ± 49.652 ± 0.000 ± 18.285 ± 304.778