High-Order Markov Blanket Discovery via a k-Order Relaxation of the Faithfulness Assumption

# High-Order Markov Blanket Discovery via a k-Order Relaxation of the Faithfulness Assumption
Source: [https://arxiv.org/html/2607.26357](https://arxiv.org/html/2607.26357)
Ragavi KrishnamoorthyHybrid Intelligence Fraunhofer IAIS Sankt Augustin, Germany Nico PiatkowskiHybrid 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](https://github.com/lklee9/k-order-Markov-blanket)\.

## 1Introduction

The goal of the Markov blanket \(MB\)discovery problem is to infer the MBof some variableYYin the graph\\gr\\grgiven a dataset\\db\\dbsampled from the joint distribution\\pr\\pr\. 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\\gr\\gr\.

X1X\_\{1\}X2X\_\{2\}X3X\_\{3\}YYZ:=\[∑i=13Xi=1\]Z:=\\Bigl\[\\sum\_\{i=1\}^\{3\}X\_\{i\}=1\\Bigr\]\\pr​\(Y=Z∣\\rmX\)=0\.9\\pr\(Y=Z\\mid\\rmX\)=0\.9Figure 1:Recovering the MBof variableYYwhenYYis the output of a noisy Boolean function that returns11when only one ofX1X\_\{1\},X2X\_\{2\}, andX3X\_\{3\}is11\. Our algorithm that assumes a22\-order faithfulness relaxation is able to fully recover the MBof each variable—even with just100100samples\.Most existing approaches to MB discovery operate under two main assumptions\. The first being that any separations in\\gr\\grimply the same conditional independence in\\pr\\pr\. The second assumption—known as the faithfulness assumption—goes the opposite direction, that is any conditional independencies in\\pr\\primply the same separation in\\gr\\gr\.

Theoretically, BNswith faithfulness violations have Lebesgue measure zero\[meek1995,spirtes2001,boeken2025\]; however these “typicality” results depend on somewhat arbitrary choices, like the choice ofσ\\sigma\-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](https://arxiv.org/html/2607.26357#Thmexample1)\.

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](https://arxiv.org/html/2607.26357#S1.F1)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](https://arxiv.org/html/2607.26357#S5)\. A full walkthrough of empirical faithfulness violation can be found later in[Example˜2](https://arxiv.org/html/2607.26357#Thmexample2)under[Section˜3](https://arxiv.org/html/2607.26357#S3)\. 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](https://arxiv.org/html/2607.26357#S6.F4)under[Section˜6\.1](https://arxiv.org/html/2607.26357#S6.SS1)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 recently22\-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 the22\-adjacency faithfulness can find XOR\-type relations, it is unable to find higher\-order parity relationships between more than33variables\.

#### Contributions

Therefore, motivated by the considerations above, in[Section˜4](https://arxiv.org/html/2607.26357#S4)we propose a relaxation of faithfulness for parity\-like relationships betweenk\+2k\+2variables: thekk\-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](https://arxiv.org/html/2607.26357#S5)by implementing a proof\-of\-concept algorithm called thekk\-order Markov blanket \(kOMB\)\. We then empirically show in[Section˜6](https://arxiv.org/html/2607.26357#S6)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 between44variables at a time, and not just33\.

## 2Probabilistic Graphical Models and Independencies

Consider the set ofnnvariables\\verts=\{V1,…,Vn\}\\verts=\\\{V\_\{1\},\\ldots,V\_\{n\}\\\}with no other hidden or latent variables\. The conditional independencies between these variables—\\rmY​\\indep​\\rmX∣\\rmZ\\rmY\\indep\{\}\\rmX\\mid\\rmZfor disjoint variables sets\\rmY,\\rmX,\\rmZ⊆\\verts\\rmY,\\rmX,\\rmZ\\subseteq\\verts—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\\gr=\(\\verts,\\edges\)\\gr=\(\\verts,\\edges\)since we are agnostic on the type of graph used for representing these graphical conditional independencies\\indeps​\\gr\\indeps\{\\gr\}\.

Regardless of the graph type of\\gr\\gr, for disjoint subsets\\rmX,\\rmY,\\rmS⊆\\verts\\rmX,\\rmY,\\rmS\\subseteq\\verts, we will denote that\\rmX\\rmXand\\rmY\\rmYare separated in\\gr\\gras,\\rmX​\\indep​\\gr​\\rmY∣\\rmS\\rmX\\indep\{\\gr\}\\rmY\\mid\\rmS\. Furthermore, we denote the set of conditional independencies encoded by the graph\\gr\\gr—i\.e\. the independence model of\\gr\\gr—as,

⟨\\rmX,\\rmY∣\\rmS⟩∈\\indeps​\\gr⇔\\rmX​\\indep​\\gr​\\rmY∣\\rmS,\\braket\{\\rmX,\\rmY\\mid\\rmS\}\\in\\indeps\{\\gr\}\\iff\\rmX\\indep\{\\gr\}\\rmY\\mid\\rmS,\(1\)where⟨\\rmX,\\rmY∣\\rmS⟩\\braket\{\\rmX,\\rmY\\mid\\rmS\}is known as an independence statement\[sadeghi2017\]\.

In addition to the graphical structure\\gr\\gr, we also consider a joint distribution\\pr\\prover\\verts\\vertsthat might encode its own set of conditional independencies between the variables in\\verts\\verts,

⟨\\rmX,\\rmY∣\\rmS⟩∈\\indeps​\\pr⇔\\rmX​\\indep​\\pr​\\rmY∣\\rmS,\\braket\{\\rmX,\\rmY\\mid\\rmS\}\\in\\indeps\{\\pr\}\\iff\\rmX\\indep\{\\pr\}\\rmY\\mid\\rmS,\(2\)where\\rmX​\\indep​\\pr​\\rmY∣\\rmS\\rmX\\indep\{\\pr\}\\rmY\\mid\\rmSdenotes\\rmX\\rmXand\\rmY\\rmYare probabilistically conditionally independent given\\rmS\\rmS\.

So far,\\indeps​\\gr\\indeps\{\\gr\}and\\indeps​\\pr\\indeps\{\\pr\}are two separate collections of conditional independencies in the graph\\gr\\grand distribution\\pr\\prrespectively\. Our base assumption throughout this paper, the*Global Markov Property*, links the two by requiring that every conditional independence encoded by the graph\\gr\\gralso hold in the distribution\\pr\\pr\. This is a common assumption in the Markov blanket discovery literatureschluter2014,yu2020a,marx2021\.

###### Assumption 1\(Global Markov Property\)\.

Given a PGM\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\), we say that the distribution\\pr\\pris*Markov*to the graph\\gr\\grwhen for all disjoint subsets\\rmX,\\rmY,\\rmS⊆\\verts\\rmX,\\rmY,\\rmS\\subseteq\\verts\[hammersley1971,lauritzen2018\],

\\rmX​\\indep​\\gr​\\rmY∣\\rmS⟹\\rmX​\\indep​\\pr​\\rmY∣\\rmS\.\\rmX\\indep\{\\gr\}\\rmY\\mid\\rmS\\implies\\rmX\\indep\{\\pr\}\\rmY\\mid\\rmS\.\(3\)Furthermore, by[Equations˜1](https://arxiv.org/html/2607.26357#S2.E1)and[2](https://arxiv.org/html/2607.26357#S2.E2), the relation in[Equation˜3](https://arxiv.org/html/2607.26357#S2.E3)is equivalent to\\indeps​\\gr⊆\\indeps​\\pr\\indeps\{\\gr\}\\subseteq\\indeps\{\\pr\}\.

Despite their differences, the independence models of\\gr\\grand\\pr\\prare both known to obey a set of axioms called the \(semi\)\-graphoid axioms\[pearl1989\]\(see[Definition˜6](https://arxiv.org/html/2607.26357#Thmdefinition6)\)\. As mentioned in[Section˜1](https://arxiv.org/html/2607.26357#S1), we are interested in learning the graphical Markov blanketof someY∈\\vertsY\\in\\vertsdefined as follows\.

###### Definition 1\(Markov Blanket and Boundary\)\.

For a variableY∈\\vertsY\\in\\verts, the MBofYYin\\gr\\gris any subset\\rmS⊆\\rmV∖Y\\rmS\\subseteq\\rmV\\setminus 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:\\rmV∖Y≔\\rmV∖\{Y\}\\rmV\\setminus Y\\coloneq\\rmV\\setminus\\\{Y\\\}\.s\.t\.

Y​\\indep​\\gr​\\verts∖\(Y∪\\rmS\)∣\\rmS,Y\\indep\{\\gr\}\\verts\\setminus\(Y\\cup\\rmS\)\\mid\\rmS,\(4\)or in other words, every other variable in\\verts∖\(Y∪\\rmS\)\\verts\\setminus\(Y\\cup\\rmS\)is graphically separated fromYYby\\rmS\\rmS\. The Markov boundary of the variableYY—henceforth denoted as\\mbs​\\gr​\(Y\)\\mbs\{\\gr\}\(Y\)—is then defined as the minimal Markov blanket ofYY\[pearl1988\]\.

An equivalent definition for the Markov blanket and boundary ofYYin\\pr\\prcan be made by just substituting the graph separation term in[Equation˜4](https://arxiv.org/html/2607.26357#S2.E4)for conditional independence in\\pr\\pr\.

Therefore, it is more accurate to say that we wish to learn the Markov boundary ofYYand 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 variableY∈\\vertsY\\in\\vertsin\\gr\\grwith 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\\pr\\pr—to grow and shrink a candidate MBofYY\. However, the MBofYYin the joint distribution\\pr\\prmight not be the same as its MBin the graph\\gr\\gr\. Following the Global Markov Property in[Assumption˜1](https://arxiv.org/html/2607.26357#Thmassumption1), we at know that the MBofYYin\\pr\\pris a subset of the MBin\\gr\\gr—which this might prevent us from fully recovering the MBin\\gr\\gr\. Therefore, in order to ensure that the Markov blankets ofYYin both\\pr\\prand\\gr\\grare equal, a common assumption constraint\-based approaches make is that\\gr\\gris*faithful*to the joint distribution\\pr\\pras defined in[Definition˜2](https://arxiv.org/html/2607.26357#Thmdefinition2)\[spirtes2001\]\.

###### Definition 2\(Faithfulness\)\.

Given a PGM\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\), the graph\\gr\\gris*faithful*to the joint distribution\\pr\\prif for all disjoint subsets\\rmX,\\rmY,\\rmS⊆\\verts\\rmX,\\rmY,\\rmS\\subseteq\\verts\[spirtes2001\],

\\rmX​\\indep​\\pr​\\rmY∣\\rmS⟹\\rmX​\\indep​\\gr​\\rmY∣\\rmS\\rmX\\indep\{\\pr\}\\rmY\\mid\\rmS\\implies\\rmX\\indep\{\\gr\}\\rmY\\mid\\rmS\(5\)which is basically the reverse direction of the global Markov property in[Assumption˜1](https://arxiv.org/html/2607.26357#Thmassumption1)\.

Therefore, by assuming both the Global Markov Property faithfulness, the conditional independencies of both\\gr\\grand\\pr\\prthen coincide,\\indeps​\\gr=\\indeps​\\pr\\indeps\{\\gr\}=\\indeps\{\\pr\}\. However faithfulness can be violated in many ways and we will explore some of them in the following section\.

## 3Faithfulness Violations and Existing Relaxations

Probably the most well\-known example of a faithfulness violation is the \(noisy\) XOR functionY≔X1⊕X2Y\\coloneq X\_\{1\}\\oplus X\_\{2\}with the graphical structureX1→Y←X2X\_\{1\}\\rightarrow Y\\leftarrow X\_\{2\}\[inazumi2011,marx2021\]222The \(noisy\) XOR function is a special case of the \(noisy\) parity function whose faithfulness violation we detail in[Example1](https://arxiv.org/html/2607.26357#Thmexample1)\.\. A recent faithfulness relaxation for tackling XOR\-type relations that is of particular interest to us is22\-adjacency faithfulness\[marx2021\]\.

###### Definition 3\(2\-adjacency faithfulness\)\.

Given the PGM\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\),\\gr\\gris22\-adjacent faithful to\\pr\\prif for all adjacent variables\(X,Y\)∈\\edges\(X,Y\)\\in\\edges,∃\\rmX⊆\\mbs​\\gr​\(Y\)\\exists\\rmX\\subseteq\\mbs\{\\gr\}\(Y\)whereX∈\\rmXX\\in\\rmXs\.t\.∀X∈\\rmX:X​\\dep​\\pr​Y∣\\rmX∖X\\forall X\\in\\rmX:X\\dep\{\\pr\}Y\\mid\\rmX\\setminus Xand if\|\\rmX\|=2\|\\rmX\|=2,X1​\\dep​\\pr​X2∣YX\_\{1\}\\dep\{\\pr\}X\_\{2\}\\mid 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\_\{1\}X2X\_\{2\}X3X\_\{3\}YY∀i∈\[3\]\\forall i\\in\[3\]\\prXi​\(0\)\\pr\_\{X\_\{i\}\}\(0\)\\prXi​\(1\)\\pr\_\{X\_\{i\}\}\(1\)0\.50\.5Figure 2:Example of a noisy parity BNover44variables\.###### Example 1\(Noisy Parity Function\)\.

Let\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\)be a PGMwhere\\rmX=\{X1,X2,X3\}\\rmX=\\\{X\_\{1\},X\_\{2\},X\_\{3\}\\\},\\verts=Y∪\\rmX\\verts=Y\\cup\\rmX, and\\gr\\grhas the structure of the DAGin[Figure˜2](https://arxiv.org/html/2607.26357#S3.F2)\. EachX∈\\rmXX\\in\\rmXare independent of each other and has the distribution of a fair coin\. The distribution ofYYconditioned on\\rmX\\rmXis then\\pr​\(Y=f​\(\\rvx\)∣\\rmX=\\rvx\)=0\.9\\pr\(Y=f\(\\rvx\)\\mid\\rmX=\\rvx\)=0\.9, where,f​\(\\rvx\)≡∑i=1nxn\(mod2\)f\(\\rvx\)\\equiv\\sum\_\{i=1\}^\{n\}x\_\{n\}\\pmod\{2\}\. For any random variableX∈\\rmXX\\in\\rmX, we know that\\pr​\(Y,X\)=0\.25=0\.5×0\.5=\\pr​\(Y\)​\\pr​\(X\)\\pr\(Y,X\)=0\.25=0\.5\\times 0\.5=\\pr\(Y\)\\pr\(X\)for any valueYYandXXtakes\. Therefore it must the case thatY​\\indep​\\pr​XY\\indep\{\\pr\}X, but sinceYYandXXare not separable in\\gr\\gr,Y​\\dep​\\gr​XY\\dep\{\\gr\}X;\\gr\\gris not faithful to\\pr\\pr\. Additionally for any distinct random variablesXi,Xj∈\\rmXX\_\{i\},X\_\{j\}\\in\\rmX,\\pr​\(Y,Xi,Xj\)=1/8=0\.53=\\pr​\(Y\)​\\pr​\(Xi\)​\\pr​\(Xj\)\\pr\(Y,X\_\{i\},X\_\{j\}\)=1/8=0\.5^\{3\}=\\pr\(Y\)\\pr\(X\_\{i\}\)\\pr\(X\_\{j\}\)for any values ofY,Xi,XjY,X\_\{i\},X\_\{j\}\. Therefore we know thatY​\\indep​\\pr​\{Xi,Xj\}Y\\indep\{\\pr\}\\\{X\_\{i\},X\_\{j\}\\\},Xi​\\indep​\\pr​\{Y,Xj\}X\_\{i\}\\indep\{\\pr\}\\\{Y,X\_\{j\}\\\}, andXj​\\indep​\\pr​\{Xi,Y\}X\_\{j\}\\indep\{\\pr\}\\\{X\_\{i\},Y\\\}\. SinceYYis adjacent to allX∈\\rmXX\\in\\rmX, but there is no22\-cardinality subsets\\rmX′⊂\\rmX\\rmX^\{\\prime\}\\subset\\rmXthat makesY​\\dep​\\pr​\\rmX′Y\\dep\{\\pr\}\\rmX^\{\\prime\};\\gr\\gris not22\-adjacent faithful to\\pr\\prfollowing the definition inmarx2021\.

Although we now know of a BNthat violates both faithfulness and22\-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](https://arxiv.org/html/2607.26357#Thmexample2), 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](https://arxiv.org/html/2607.26357#Thmexample1), let\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\)be a PGMwith random variables\\rmX\\rmXandYYwith the DAGstructure in[Figure˜3](https://arxiv.org/html/2607.26357#S3.F3)\. The distribution ofYYconditioned on\\rmX\\rmXis then\\pr​\(Y=\[∑X∈\\rmXX=1\]\|\\rmX\)=0\.9\\pr\\bigl\(Y=\\bigl\[\\sum\_\{X\\in\\rmX\}X=1\\bigr\]\\bigm\|\\rmX\\bigr\)=0\.9, where\[⋅\]\[\\cdot\]is the Iverson bracket notation for the indicator function\. In other words,Y=1Y=1has a probability of0\.90\.9when only one of the variables in\\rmX\\rmXis11\. Otherwise,Y=0Y=0has a probability of0\.90\.9\.

The marginal distribution ofYYis\\pr​\(Y=1\)=3/8\\pr\(Y=1\)=3/8since only 3 value\-combinations of\\rmX\\rmXhas just a single11among them\. We then know that∀X∈\\rmX:\\pr​\(Y=1,X=1\)=1/8≠\(3/8\)×\(4/8\)=\\pr​\(Y=1\)​\\pr​\(X=1\)\\forall X\\in\\rmX:\\pr\(Y=1,X=1\)=1/8\\neq\(3/8\)\\times\(4/8\)=\\pr\(Y=1\)\\pr\(X=1\)\. Therefore,∀X∈\\rmX:Y​\\dep​\\pr​X\\forall X\\in\\rmX:Y\\dep\{\\pr\}Xand consequently,\\gr\\gr*is*faithful to\\pr\\pr\.

However, it is possible for a sample\\db\\dbof\\pr\\prto be pathological in the sense that, the empirical distribution from\\db\\dbinduces a independence betweenYYand someX∈\\rmXX\\in\\rmX\. For example, the sample in[Figure˜3](https://arxiv.org/html/2607.26357#S3.F3)causes the marginal distribution ofYYto bePhys\.Rev\.E​\(Y\)=1/2\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\(Y\)=1/2and the marginal distribution overYYandX1X\_\{1\}to bePhys\.Rev\.E​\(Y,X1\)=1/4=\(1/2\)×\(1/2\)=Phys\.Rev\.E​\(Y\)​Phys\.Rev\.E​\(X1\)\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\(Y,X\_\{1\}\)=1/4=\(1/2\)\\times\(1/2\)=\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\(Y\)\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\(X\_\{1\}\)for all possible values ofYYandX1X\_\{1\}\. Therefore under this empirical distributionPhys\.Rev\.E\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\},Y​\\indep​Phys\.Rev\.E​X1Y\\indep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\_\{1\}—which empirically violates the faithfulness assumption\.

Fortunately in[Example˜2](https://arxiv.org/html/2607.26357#Thmexample2), the DAGin[Figure˜3](https://arxiv.org/html/2607.26357#S3.F3)is still22\-adjacent faithful to the empirical distribution, therefore the full MBis still recoverable by the algorithm outlined inmarx2021\. However, this is not always the case and as we will see in[Section˜6\.1](https://arxiv.org/html/2607.26357#S6.SS1), sometimes considering higher\-order relationships between variables can help better overcome empirical faithfulness violations\.

X1X\_\{1\}X2X\_\{2\}X3X\_\{3\}YYg​\(\\rmX\):=\[∑X∈\\rmXX=1\]g\(\\rmX\):=\\Bigl\[\\sum\_\{X\\in\\rmX\}X=1\\Bigr\]\\pr​\(Y=g​\(\\rmX\)∣\\rmX\)=0\.9\\pr\\bigl\(Y=g\(\\rmX\)\\mid\\rmX\\bigr\)=0\.9Pathological𝒟\\mathcal\{D\}X1X\_\{1\}X2X\_\{2\}X3X\_\{3\}YY00000011010101101001101011001001Figure 3:Example of a noisy exactly\-1 BNover44variables with a pathological sample that causes the DAGto be unfaithful to the empirical distribution\.
## 4The k\-Order Faithfulness Relaxation

In this section we will introduce our generalisation of the faithfulness assumption in[Definition˜2](https://arxiv.org/html/2607.26357#Thmdefinition2)for finding conditional dependencies between two variables in\\verts\\vertsthat require consideringkkadditional variables to be found\. But before that, we first define in[Definition˜4](https://arxiv.org/html/2607.26357#Thmdefinition4)what it means for two distinct variablesY,X∈\\vertsY,X\\in\\vertsto only be conditionally dependent given somekk\-cardinality subset\\rmZ⊆\\subverts,\|\\rmZ\|=k\\rmZ\\subseteq\\subverts,\|\\rmZ\|=k\.

###### Definition 4\(kk\-Order Dependence\)\.

Given a PGM\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\), letY,X∈\\rmVY,X\\in\\rmVbe two distinct random variables and\\rmS,\\rmZ⊆\\subverts\\rmS,\\rmZ\\subseteq\\subvertsbe two disjoint subsets\. Then we denote the absence of any subset of\\rmS′⊆\\rmS\\rmS^\{\\prime\}\\subseteq\\rmSthat rendersYYandXXconditionally independent given\\rmZ∪\\rmS′\\rmZ\\cup\\rmS^\{\\prime\}as

\\kassoc​\\rmS​Y,X​\\rmZ:=∀\\rmS′⊆\\rmS:Y​\\dep​\\pr​X∣\\rmZ∪\\rmS′\.\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ\}:=\\forall\\rmS^\{\\prime\}\\subseteq\\rmS:Y\\dep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\.\(6\)Here we say thatYYandXXare\|\\rmZ\|\|\\rmZ\|\-order associated w\.r\.t\.\\rmZ\\rmZover the “separating” set\\rmS\\rmS\. Of course under this definition there might be some proper subset\\rmZ′⊂\\rmZ\\rmZ^\{\\prime\}\\subset\\rmZwhere\\kassoc​\\rmS​Y,X​\\rmZ′\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ^\{\\prime\}\}still holds\. Therefore we sayYYandXXare\|\\rmZ\|\|\\rmZ\|\-order*dependent*on\\rmZ\\rmZonly if this is not the case,

\\kdep​\\rmS​Y,X​\\rmZ:=\\kassoc​\\rmS​Y,X​\\rmZ∧∀\\rmZ′⊂\\rmZ:¬\\kassoc​\\rmS​Y,X​\\rmZ′\.\\displaystyle\\begin\{aligned\} &\\kdep\{\\rmS\}\{Y,X\}\{\\rmZ\}\\\\ &:=\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ\}\\wedge\\forall\\rmZ^\{\\prime\}\\subset\\rmZ:\\neg\\,\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ^\{\\prime\}\}\.\\end\{aligned\}\(7\)

With this definition ofkk\-order dependence, we are able to capture parity\-like relations overk\+2k\+2variables—such as the one in[Example˜1](https://arxiv.org/html/2607.26357#Thmexample1)\. Specifically we know that the noisy\-parity relationship in[Example˜1](https://arxiv.org/html/2607.26357#Thmexample1)is22\-order dependent since for distinct variablesA,B∈\{X1,X2,X3,Y\}A,B\\in\\\{X\_\{1\},X\_\{2\},X\_\{3\},Y\\\},\\rmZ\\rmZmust always contain the other two variables for\\kassoc​∅​A,B​\\rmZ\\kassoc\{\\emptyset\}\{A,B\}\{\\rmZ\}to be true\.

Further note that whenk=1k=1, our notion ofkk\-order dependence— similar to*22\-association*inmarx2021—is able to capture the relationships between*Unfaithful Triples*which includes XOR\-type relationships\.

There are two aspects of[Definition˜4](https://arxiv.org/html/2607.26357#Thmdefinition4)whose purpose might not immediately jump out\. The first being that in[Equation˜6](https://arxiv.org/html/2607.26357#S4.E6), we requireY​\\dep​\\pr​X∣\\rmZ∪\\rmS′Y\\dep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}to be true for all subsets\\rmS′⊆\\rmS\\rmS^\{\\prime\}\\subseteq\\rmS—instead of justY​\\dep​\\pr​X∣\\rmZ∪\\rmSY\\dep\{\\pr\}X\\mid\\rmZ\\cup\\rmS\. This requirement ensures thatYYandXXare dependent as long as they are conditioned on\\rmZ\\rmZ—forcing a distinction between the possible dependants\\rmZ\\rmZand the possible separators\\rmS\\rmS\. The other odd aspect of[Definition˜4](https://arxiv.org/html/2607.26357#Thmdefinition4)is that we define bothkk\-order association and dependence over some arbitrary subset\\rmS⊆\\subverts\\rmS\\subseteq\\subvertsinstead of just the entirety of\\subverts\\subverts\. 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](https://arxiv.org/html/2607.26357#S5)that\\methodcorrectly recovers the Markov boundary of some variableY∈\\vertsY\\in\\verts\.

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\\verts\\vertsbe a set of elements andk∈ℕ:k≤\|\\verts\|k\\in\\mathbb\{N\}:k\\leq\|\\verts\|be some natural number333We adopt the definition of the natural numbers in which0is included\.smaller than the cardinality of\\verts\\verts\. Then we use,

\\powereq​k​\\verts:=\{\\rmZ⊆\\verts:\|\\rmZ\|=k\},\\powereq\{k\}\{\\verts\}:=\\\{\\rmZ\\subseteq\\verts:\|\\rmZ\|=k\\\},\(8\)to denote the set of all subsets of\\verts\\vertswith cardinalitykkand,

\\powerleq​k​\\verts:=\{\\rmZ⊆\\verts:\|\\rmZ\|≤k\},\\powerleq\{k\}\{\\verts\}:=\\\{\\rmZ\\subseteq\\verts:\|\\rmZ\|\\leq k\\\},\(9\)to denote the set of all the subsets of\\verts\\vertswith cardinality less than or equal tokk\.

Before introducing our generalisation of the faithfulness assumption tokk\-order dependencies, we first show in[Definition˜5](https://arxiv.org/html/2607.26357#Thmdefinition5)thatkk\-order associations w\.r\.t\.\\rmZ\\rmZimplies the existence of akk\-order dependence w\.r\.t\. some subset of\\rmZ\\rmZ\.

\{restatable\}

\[kk\-Order Dependence Arises from Association\]lemmakdeptoassoc Given a PGM\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\), letY,X∈\\rmVY,X\\in\\rmVbe distinct random variables and\\rmS,\\rmZ⊂\\subverts\\rmS,\\rmZ\\subset\\subvertsbe disjoint subsets\. Then the following holds for allY,X,\\rmS,\\rmZY,X,\\rmS,\\rmZ:

\\kassoc​\\rmS​Y,X​\\rmZ⟹∃\\rmZ′⊆\\rmZ:\\kdep​\\rmS​Y,X​\\rmZ′\.\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ\}\\implies\\exists\\rmZ^\{\\prime\}\\subseteq\\rmZ:\\kdep\{\\rmS\}\{Y,X\}\{\\rmZ^\{\\prime\}\}\.\(10\)The proof of[Definition˜5](https://arxiv.org/html/2607.26357#Thmdefinition5)and all other Lemmas in the paper can be found in[Appendix˜B](https://arxiv.org/html/2607.26357#A2)\.

###### Assumption 2\(kk\-Order Faithfulness\)\.

Given the PGM\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\),\\gr\\griskk\-order faithful to\\pr\\prif∀\\rmS⊆\\verts∖\{Y,X\}\\forall\\rmS\\subseteq\\verts\\setminus\\\{Y,X\\\},

∀\\rmZ∈\\power≤k​\(\\verts∖\{Y,X\}\)​∃\\rmS′⊆\\rmS:Y​\\indep​\\pr​X∣\\rmZ∪\\rmS′⟹Y​\\indep​\\gr​X∣\\rmS\.\\displaystyle\\begin\{aligned\} &\\forall\\rmZ\\in\\power\{\\leq k\}\(\\verts\\setminus\\\{Y,X\\\}\)\\\>\\\>\\exists\\rmS^\{\\prime\}\\subseteq\\rmS:Y\\indep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\\\\ &\\implies Y\\indep\{\\gr\}X\\mid\\rmS\.\\end\{aligned\}\(11\)Its contrapositive then states that∀\\rmS⊆\\verts∖\{Y,X\}\\forall\\rmS\\subseteq\\verts\\setminus\\\{Y,X\\\},

Y​\\dep​\\gr​X∣\\rmS⟹∃\\rmZ∈\\power≤k​\(\\verts∖\{Y,X\}\):\\kassoc​\\rmS​Y,X​\\rmZ\.\\displaystyle\\begin\{aligned\} &Y\\dep\{\\gr\}X\\mid\\rmS\\\\ &\\implies\\exists\\rmZ\\in\\power\{\\leq k\}\(\\verts\\setminus\\\{Y,X\\\}\):\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ\}\.\\end\{aligned\}\(12\)However, the definition in[Equation˜12](https://arxiv.org/html/2607.26357#S4.E12)allowskkto be larger than it needs to be—i\.e\.kkis not strict\. Therefore, we say that\\gr\\gris*strictly*kk\-order faithful to\\pr\\prif∀\\rmS⊆\\verts∖\{Y,X\}\\forall\\rmS\\subseteq\\verts\\setminus\\\{Y,X\\\},

Y​\\dep​\\gr​X∣\\rmS⟹∃\\rmZ∈\\powereq​k′​\\verts∖\{Y,X\}:\\kdep​\\rmS​Y,X​\\rmZ,\\displaystyle\\begin\{aligned\} &Y\\dep\{\\gr\}X\\mid\\rmS\\\\ &\\implies\\exists\\rmZ\\in\\powereq\{k^\{\\prime\}\}\{\\verts\\setminus\\\{Y,X\\\}\}:\\kdep\{\\rmS\}\{Y,X\}\{\\rmZ\},\\end\{aligned\}\(13\)wherek′≤kk^\{\\prime\}\\leq k\. This definition in[Equation˜13](https://arxiv.org/html/2607.26357#S4.E13)logically follows from[Equation˜12](https://arxiv.org/html/2607.26357#S4.E12)as a result of[Definition˜5](https://arxiv.org/html/2607.26357#Thmdefinition5)\.

With this new generalisation of faithfulness, we will now present an initial algorithm,\\method, that makes use of ourkk\-order faithfulness assumption for MBsdiscovery\.

## 5kOMB: k\-Order Markov Blanket

Our algorithm for Markov blanket \(MB\)discovery,\\method, is a modified version of the Grow and Shrink \(GS\)algorithm bymargaritis1999\. In GS, a candidate MB\\rmS\\rmSis iteratively grown by testing the conditional independence betweenYYand variables that are not yet in\\rmS\\rmS—conditioned the current candidate MB\\rmS\\rmS\. Therefore, the separating set used in the conditional independence tests change as the GSalgorithm progresses\.

In\\method, something similar occurs where we will only consider subsets of the candidate MBwhen trying to find subsets that successfully causeYYandXXto be conditionally independent\. In other words, instead of considering the dependence betweenYYandXXover all the variables in\\verts\\verts,\\kdep​\\subverts​Y,X​\\rmZ\\kdep\{\\subverts\}\{Y,X\}\{\\rmZ\}; we will instead consider the association betweenYYandXXover just the current candidate MB\\rmS\\rmS,\\kassoc​\\rmS​Y,X​\\rmZ\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ\}\. As a first step, we show that akk\-order dependence betweenYYandXXwith no separating variables to consider implies thatYYiskk\-order associated with every one of the involved variables\.

\{restatable\}

\[Marginalkk\-Order Dependence Implies Inter\-Association\]lemmakdepinter Given a PGM\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\), letX,Y∈\\vertsX,Y\\in\\vertsbe distinct variables and\\rmZ⊆\\subverts\\rmZ\\subseteq\\subverts\. Then ifYYandXXare\|\\rmZ\|\|\\rmZ\|\-order dependent w\.r\.t\.\\rmZ\\rmZover the empty separating set,YYis\|\\rmZ\|\|\\rmZ\|\-order associated with every variable in\\rmZ∪X\\rmZ\\cup X,

\\kdep​∅​Y,X​\\rmZ⟹\\kinter​∅​Y,\\rmZ∪X,\\displaystyle\\kdep\{\\emptyset\}\{Y,X\}\{\\rmZ\}\\implies\\kinter\{\\emptyset\}\{Y,\\rmZ\\cup X\},\(14\)where\\kinter​\\rmS​Y,\\rmW:=∀W∈\\rmW:\\kassoc​\\rmS​Y,W​\\rmW∖W\.\\displaystyle\\begin\{aligned\} \\kinter\{\\rmS\}\{Y,\\rmW\}:=\\forall W\\in\\rmW:\\kassoc\{\\rmS\}\{Y,W\}\{\\rmW\\setminus W\}\.\\end\{aligned\}\(15\)When[Equation˜14](https://arxiv.org/html/2607.26357#S5.E14)is true, we say thatYYis\|\\rmZ\|\|\\rmZ\|\-order*inter\-associated*withX∪\\rmZX\\cup\\rmZ\.

[Section˜5](https://arxiv.org/html/2607.26357#S5)requires no faithfulness assumption—it holds for every distribution via the semi\-graphoid axioms—and generalises the22\-association property of unfaithful triples inmarx2021to arbitrary orders; the noisy\-parity[Example˜1](https://arxiv.org/html/2607.26357#Thmexample1)is exactly this case with\\rmZ=\{X2,X3\}\\rmZ=\\\{X\_\{2\},X\_\{3\}\\\}\. It motivates themake\_strictroutine of[Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2), which prunes a found dependence set towards an inter\-associated one\.

The last consideration we need to make for\\methodto be remotely feasible—especially in cases where the true MBofY∈\\vertsY\\in\\vertsin\\gr\\gris large—is that, we need to somehow limit the size of the conditioning set when conducting the conditional independence \(CI\)tests in\\method\. When testing ifYYandXXarekk\-order associated, we will generally use CItests of the form,∀\\rmS′⊆\\rmS:Y​\\indep​\\pr​X∣\\rmZ∪\\rmS′\\forall\\rmS^\{\\prime\}\\subseteq\\rmS:Y\\indep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}, where\\rmZ\\rmZis the “dependence” set that rendersYYandXXconditionally dependent\. Fortunately, assumingkk\-order faithfulness, we know that the size of\\rmZ\\rmZis bounded bykk\.

However the same cannot be said for the set\\rmS′\\rmS^\{\\prime\}as it can be as large as the candidate MB\\rmS\\rmS–which can be arbitrarily large depending on\\gr\\gr\. Therefore, the third and last assumption we will make is that ifYYandXXare conditionally independent given\\rmZ∪\\rmS\\rmZ\\cup\\rmS, then there must exist some subset with size less thanll,\\rmS′∈\\powerleq​l​\\rmS\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\}, such thatYYandXXare still conditionally independent given\\rmZ∪\\rmS′\\rmZ\\cup\\rmS^\{\\prime\}\.

###### Assumption 3\(ll\-Bounded Separator Assumption\)\.

For a PGM\\model=\(\\gr,\\pr\)\\model=\(\\gr,\\pr\), distinct variablesY,X∈\\vertsY,X\\in\\verts, and\\verts′≔\\subverts\\verts^\{\\prime\}\\coloneq\\subverts, let us first recall thekk\-order faithfulness assumption as defined in[Equation˜11](https://arxiv.org/html/2607.26357#S4.E11)∀\\rmS⊆\\subverts\\forall\\rmS\\subseteq\\subverts,

∀\\rmZ∈\\power≤k​\(\\verts′\)​∃\\rmS′⊆\\rmS:Y​\\indep​\\pr​X∣\\rmZ∪\\rmS′⟹Y​\\indep​\\gr​X∣\\rmS\.\\displaystyle\\begin\{aligned\} &\\forall\\rmZ\\in\\power\{\\leq k\}\(\\verts^\{\\prime\}\)\\\>\\\>\\exists\\rmS^\{\\prime\}\\subseteq\\rmS:Y\\indep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\\\\ &\\implies Y\\indep\{\\gr\}X\\mid\\rmS\.\\end\{aligned\}\([11](https://arxiv.org/html/2607.26357#S4.E11)\)Then, for some natural numberl∈ℕl\\in\\mathbb\{N\}, we define anll\-bounded separator version of thekk\-order faithfulness assumption∀\\rmS⊆\\subverts\\forall\\rmS\\subseteq\\subverts,

∀\\rmZ∈\\power≤k​\(\\verts′\)​∃\\rmS′∈\\powerleq​l​\\rmS:Y​\\indep​\\pr​X∣\\rmZ∪\\rmS′⇔∀\\rmZ∈\\power≤k\(\\verts′\)∃\\rmS′⊆\\rmS:Y\\indep\\prX∣\\rmZ∪\\rmS′⟹Y​\\indep​\\gr​X∣\\rmS,\\displaystyle\\begin\{aligned\} &\\forall\\rmZ\\in\\power\{\\leq k\}\(\\verts^\{\\prime\}\)\\\>\\\>\\exists\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\}:Y\\indep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\\\\ &\\iff\\forall\\rmZ\\in\\power\{\\leq k\}\(\\verts^\{\\prime\}\)\\\>\\\>\\exists\\rmS^\{\\prime\}\\subseteq\\rmS:Y\\indep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\\\\ &\\implies Y\\indep\{\\gr\}X\\mid\\rmS,\\end\{aligned\}\(16\)which essentially states that as long asYYandXXare conditionally independent given\\rmZ\\rmZand every≤l\\leq l\-cardinality subset of\\rmS\\rmS, thenYYandXXare conditionally independent given\\rmZ\\rmZand every possible subset of\\rmS\\rmS—regardless of cardinality\.

The contrapositive of this assumption is then

Y​\\dep​\\gr​X∣\\rmS⟹∃\\rmZ∈\\power≤k​\(\\verts′\)​∀\\rmS′⊆\\rmS:Y​\\dep​\\pr​X∣\\rmZ∪\\rmS′⇔∃\\rmZ∈\\power≤k\(\\verts′\):\\kassoc\\rmSY,X\\rmZ⇔∃\\rmZ∈\\powerleqk\\verts′:\\lassoc\\rmSY,X\\rmZ,\\displaystyle\\begin\{aligned\} &Y\\dep\{\\gr\}X\\mid\\rmS\\\\ &\\implies\\exists\\rmZ\\in\\power\{\\leq k\}\(\\verts^\{\\prime\}\)\\\>\\\>\\forall\\rmS^\{\\prime\}\\subseteq\\rmS:Y\\dep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\\\\ &\\iff\\exists\\rmZ\\in\\power\{\\leq k\}\(\\verts^\{\\prime\}\):\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ\}\\\\ &\\iff\\exists\\rmZ\\in\\powerleq\{k\}\{\\verts^\{\\prime\}\}:\\lassoc\{\\rmS\}\{Y,X\}\{\\rmZ\},\\end\{aligned\}\(17\)where\\lassoc​\\rmS​Y,X​\\rmZ≔∀\\rmS′∈\\powerleq​l​\\rmS:Y​\\dep​\\pr​X∣\\rmZ∪\\rmS′\\lassoc\{\\rmS\}\{Y,X\}\{\\rmZ\}\\coloneq\\forall\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\}:Y\\dep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\.

The concept of bounding the maximum size of the separators in a graph has been explored before\. For instance,soh2019defined an undirected graph\\gr\\grto be weaklyKK\-separable if for any two non\-adjacent vertices in\\gr\\gr, there exists a set of vertices with cardinality≤K\\leq Kthat separates the two vertices in\\gr\\gr\. However, instead of focusing on separation in the graph\\gr\\gr, ourll\-bounded separator assumption focuses on the cardinality of “separating” sets that causes conditional independencies between variables in\\verts\\verts\.

Thell\-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](https://arxiv.org/html/2607.26357#Thmassumption3)simply makes such a bound explicit\. Furthermore, as shown later in[Appendix˜F](https://arxiv.org/html/2607.26357#A6), the bound is mild on the networks used in[Section˜6\.2](https://arxiv.org/html/2607.26357#S6.SS2)\.

### 5\.1Theoretical Algorithm

Input:

\\db\\db,

\\model\\model,

YY,

α\\alpha,

kk,

ll
Output:

\\mbs​\\gr​\(Y\)\\mbs\{\\gr\}\(Y\)
\\rmS←∅\\rmS\\leftarrow\\emptyset;

//Candidate Markov Blanket

1repeat

2

\\verts′←\\verts∖\(Y∪\\rmS\)\\verts^\{\\prime\}\\leftarrow\\verts\\setminus\(Y\\cup\\rmS\);

3

\\rmX′←\\rmX^\{\\prime\}\\leftarrowfind\_inter\_dep\(

YY,

∅\\emptyset,

\\rmS\\rmS,

\\verts′\\verts^\{\\prime\},

kk,

ll\);

4

\\rmS←\\rmS∪\\rmX′\\rmS\\leftarrow\\rmS\\cup\\rmX^\{\\prime\}
5until*\\rmX′=∅\\rmX^\{\\prime\}=\\emptyset*;

6while*∃X∈\\rmS:find\_cond​\(Y,X,∅,\\rmS∖X,k,l\)=⊥\\exists X\\in\\rmS:\\textit\{find\\\_cond\}\(Y,X,\\emptyset,\\rmS\\setminus X,k,l\)=\\bot*do

7

\\rmS←\\rmS∖X\\rmS\\leftarrow\\rmS\\setminus X;

8

9return*\\rmS\\rmS*;

Algorithm 1kOMBAs a proof\-of\-concept, in this section we propose\\method—a modification of the Grow and Shrink \(GS\)bymargaritis1999to discover Markov blankets with higher\-order relationships under the following assumptions: 1\) the Global Markov Property from[Assumption˜1](https://arxiv.org/html/2607.26357#Thmassumption1), 2\)kk\-order faithfulness as in[Assumption˜2](https://arxiv.org/html/2607.26357#Thmassumption2), and 3\) thell\-bounded separator assumption in[Assumption˜3](https://arxiv.org/html/2607.26357#Thmassumption3)\.

The original GSalgorithm starts with an empty MB\\rmS=∅\\rmS=\\emptysetand iteratively adds a variableXXinto\\rmS\\rmSif the conditional independence \(CI\)testY​\\indep​\\pr​X∣\\rmSY\\indep\{\\pr\}X\\mid\\rmSindicates thatYYis conditionally independent toXXgiven\\rmS\\rmS\. The main difference between\\methodand the GSalgorithm then is: 1\) considering the addition of entire sets of variablesX∪\\rmX:\\rmX∈\\powerleq​k​\\verts∖Y,X∈\\verts∖\(Y∪\\rmS∪\\rmX\)X\\cup\\rmX:\\rmX\\in\\powerleq\{k\}\{\\verts\\setminus Y\},X\\in\\verts\\setminus\(Y\\cup\\rmS\\cup\\rmX\)at each iteration to find higher\-order parity\-type relationships between variables, and 2\) conducting CItests using just subsets of the current candidate MB\\rmS\\rmSin the condition to ensure that these CItests remain feasible even when the cardinality of\\rmS\\rmSis large\. As such, we break up\\methodinto 3 sub\-algorithms:

[Algorithm1](https://arxiv.org/html/2607.26357#algorithm1)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,

[Algorithm2](https://arxiv.org/html/2607.26357#algorithm2)is overall responsible for finding some subset of the variables not already in the candidate MB,\\rmX∈\\powerleq​k​\\verts∖\(X∪Y∪\\rmS\)\\rmX\\in\\powerleq\{k\}\{\\verts\\setminus\(X\\cup Y\\cup\\rmS\)\}, that might causeYYandXXto be conditionally dependent given\\rmX\\rmXand some subset\\rmX′∈\\powerleq​k−\|\\rmX\|​\\rmS\\rmX^\{\\prime\}\\in\\powerleq\{k\-\|\\rmX\|\}\{\\rmS\}of the current candidate MB, finally,

[Algorithm3](https://arxiv.org/html/2607.26357#algorithm3)is then tasked with finding such a\\rmX′∈\\powerleq​k−\|\\rmX\|​\\rmS\\rmX^\{\\prime\}\\in\\powerleq\{k\-\|\\rmX\|\}\{\\rmS\}that causesY​\\dep​\\pr​Y∣\\rmX∪\\rmX′Y\\dep\{\\pr\}Y\\mid\\rmX\\cup\\rmX^\{\\prime\}and checking if there are anyll\-bounded subsets of the current candidate MB,\\rmS′∈\\powerleq​l​\\rmS\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\}, such thatY​\\indep​\\pr​Y∣\\rmX∪\\rmX′∪\\rmS′Y\\indep\{\\pr\}Y\\mid\\rmX\\cup\\rmX^\{\\prime\}\\cup\\rmS^\{\\prime\}\.

Keeping this outline in mind, we now present the full algorithm for\\methodand show that it correctly discovers the MBof someY∈\\vertsY\\in\\vertsassuming that[Assumptions˜1](https://arxiv.org/html/2607.26357#Thmassumption1),[2](https://arxiv.org/html/2607.26357#Thmassumption2)and[3](https://arxiv.org/html/2607.26357#Thmassumption3)hold\. The proofs for any subsequent Theorems can be found in[Appendix˜C](https://arxiv.org/html/2607.26357#A3)\.

1Function*make\_strict\(YY,XX,\\rmZ\\rmZ,\\rmS\\rmS,kk,ll\)*:

2for*Z∈\\rmZZ\\in\\rmZ*do

3if*¬\\negno\_seps\(YY,ZZ,X∪\\rmZ∖ZX\\cup\\rmZ\\setminus Z,\\rmS\\rmS,ll\)*then

4remove

ZZfrom

\\rmZ\\rmZ;

5

6

7return*X∪\\rmZX\\cup\\rmZ*;

8

9Function*find\_inter\_dep\(YY,\\rmX\\rmX,\\rmS\\rmS,\\rmV′\\rmV^\{\\prime\},kk,ll\)*:

10if*\|\\rmX\|\>k\|\\rmX\|\>k*then

11return*∅\\emptyset*;

12

13for*X∈\\rmV∖\(Y∪\\rmS∪\\rmX\)X\\in\\rmV\\setminus\(Y\\cup\\rmS\\cup\\rmX\)*do

14

\\rmZ←\\rmZ\\leftarrowfind\_cond\(

YY,

XX,

\\rmX\\rmX,

\\rmS\\rmS,

kk,

ll\);

15if*\\rmZ≠⊥\\rmZ\\neq\\bot*then

16return*make\_strict\(YY,XX,\\rmZ\\rmZ,\\rmS\\rmS,kk,ll\)*;

17

18

19for*X∈\\rmV′∖\(Y∪\\rmS∪\\rmX\)X\\in\\rmV^\{\\prime\}\\setminus\(Y\\cup\\rmS\\cup\\rmX\)*do

20

\\rmV′←\\rmV′∖X\\rmV^\{\\prime\}\\leftarrow\\rmV^\{\\prime\}\\setminus X;

21

\\rmX′←\\rmX^\{\\prime\}\\leftarrowfind\_inter\_dep\(

YY,

\\rmX∪X\\rmX\\cup X,

\\rmS\\rmS,

\\rmV′\\rmV^\{\\prime\},

kk,

ll\);

22if*\\rmX′≠∅\\rmX^\{\\prime\}\\neq\\emptyset*then

23return*\\rmX′\\rmX^\{\\prime\}*;

24

25

26return*∅\\emptyset*;

27

Algorithm 2Association Mining\{restatable\}

\[Correctness of[Algorithm˜3](https://arxiv.org/html/2607.26357#algorithm3)\]theoremcitest Given an estimated distributionPhys\.Rev\.E\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}, distinct variablesY,X∈\\vertsY,X\\in\\verts, disjoint subsets\\rmX,\\rmS⊂\\subverts\\rmX,\\rmS\\subset\\subverts, and the natural numbersk,l∈ℕk,l\\in\\mathbb\{N\}; the functionfind\_condin[Algorithm˜3](https://arxiv.org/html/2607.26357#algorithm3)returns the set\\rmX′∪\\rmX\\rmX^\{\\prime\}\\cup\\rmXif and only if∃\\rmX′∈\\powerleq​k−\|\\rmX\|​\\rmS\\exists\\rmX^\{\\prime\}\\in\\powerleq\{k\-\|\\rmX\|\}\{\\rmS\}where,

∀\\rmS′∈\\powerleql\\rmS∖\\rmX′:Y\\depPhys\.Rev\.EX∣\\rmX∪\\rmX′∪\\rmS′\.\\forall\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmX^\{\\prime\}\}:Y\\dep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmX\\cup\\rmX^\{\\prime\}\\cup\\rmS^\{\\prime\}\.\(18\)Otherwisefind\_condreturns the failure value⊥\\bot, distinct from the empty set, which it may return on success when\\rmX′=\\rmX=∅\\rmX^\{\\prime\}=\\rmX=\\emptyset\.

\{restatable\}

\[Correctness of[Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2)\]theoremassocmine Given an estimated distributionPhys\.Rev\.E\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}, the target variableY∈\\vertsY\\in\\verts, the current candidate MB\\rmS⊂\\verts∖Y\\rmS\\subset\\verts\\setminus Y, the remaining variables\\verts′=\\verts∖\(Y∪\\rmS\)\\verts^\{\\prime\}=\\verts\\setminus\(Y\\cup\\rmS\), and the natural numbersk,l∈ℕk,l\\in\\mathbb\{N\}; the callfind\_inter\_dep\(Y,∅,\\rmS,\\verts′,k,l\)\(Y,\\emptyset,\\rmS,\\verts^\{\\prime\},k,l\)returns the empty set if and only if

∀X∈\\verts′​∀\\rmZ∈\\powerleq​k​\\verts∖\{Y,X\}​∃\\rmS′∈\\powerleq​l​\\rmS∖\\rmZ:\\displaystyle\\forall X\\in\\verts^\{\\prime\}\\\>\\\>\\forall\\rmZ\\in\\powerleq\{k\}\{\\verts\\setminus\\\{Y,X\\\}\}\\\>\\\>\\exists\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmZ\}:\(19\)Y​\\indep​Phys\.Rev\.E​X∣\\rmZ∪\\rmS′\.\\displaystyle\\qquad Y\\indep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\.Otherwise it returns a nonempty setX∪\\rmZRX\\cup\\rmZ\_\{R\}withX∈\\verts′X\\in\\verts^\{\\prime\}and\\rmZR⊆\\rmZ\\rmZ\_\{R\}\\subseteq\\rmZfor some\\rmZ∈\\powerleq​k​\\verts∖\{Y,X\}\\rmZ\\in\\powerleq\{k\}\{\\verts\\setminus\\\{Y,X\\\}\}satisfying the estimatedll\-bounded association∀\\rmS′∈\\powerleql\\rmS∖\\rmZ:Y\\depPhys\.Rev\.EX∣\\rmZ∪\\rmS′\\forall\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmZ\}:Y\\dep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\.

1Function*no\_seps\(YY,XX,\\rmZ\\rmZ,\\rmS\\rmS,ll\)*:

2for*\\rmS′∈\\powerleq​l​\\rmS∖\\rmZ\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmZ\}*do

3if*Y​\\indep​X∣\\rmZ∪\\rmS′Y\\indep\{\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}*then

4return*False*;

5

6

7return*True*;

8

9Function*find\_cond\(YY,XX,\\rmX\\rmX,\\rmS\\rmS,kk,ll\)*:

10for*\\rmX′∈\\powerleq​k−\|\\rmX\|​\\rmS\\rmX^\{\\prime\}\\in\\powerleq\{k\-\|\\rmX\|\}\{\\rmS\}*do

11

\\rmZ←\\rmX′∪\\rmX\\rmZ\\leftarrow\\rmX^\{\\prime\}\\cup\\rmX;

12if*Y​\\dep​X∣\\rmZY\\dep\{\}X\\mid\\rmZ*then

13if*no\_seps\(YY,XX,\\rmZ\\rmZ,\\rmS\\rmS,ll\)*then

14return*\\rmZ\\rmZ*;

15

16

17

18return*⊥\\bot*;

19

Algorithm 3Practical Conditional Independence Test\{restatable\}

\[Correctness of[Algorithm˜1](https://arxiv.org/html/2607.26357#algorithm1)\]theoremalgkomb Given an estimated distributionPhys\.Rev\.E\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}over variables\\verts\\verts, the target variable whose graphical Markov blanket we want to discoverY∈\\vertsY\\in\\verts, and the natural numbersk,l∈ℕk,l\\in\\mathbb\{N\}; the grow phase of\\methodthen grows the candidate MB\\rmS\\rmSsuch that,

\\mbs​\\gr​\(Y\)⊆\\rmS\.\\mbs\{\\gr\}\(Y\)\\subseteq\\rmS\.\(20\)If, in addition, every true blanket member admits an in\-blanket dependence witness—i\.e\. for everyX∈\\mbs​\\gr​\(Y\)X\\in\\mbs\{\\gr\}\(Y\)and every\\rmS\\rmSwith\\mbs​\\gr​\(Y\)⊆\\rmS⊆\\verts∖Y\\mbs\{\\gr\}\(Y\)\\subseteq\\rmS\\subseteq\\verts\\setminus Ythere is a witness\\rmZ∈\\powerleq​k​\\mbs​\\gr​\(Y\)∖X\\rmZ\\in\\powerleq\{k\}\{\\mbs\{\\gr\}\(Y\)\\setminus X\}with\\lassoc​\\rmS∖\(X∪\\rmZ\)​Y,X​\\rmZ\\lassoc\{\\rmS\\setminus\(X\\cup\\rmZ\)\}\{Y,X\}\{\\rmZ\}—then the shrink phase of\\methodremoves variables from\\rmS\\rmSuntil,

\\mbs​\\gr​\(Y\)=\\rmS,\\mbs\{\\gr\}\(Y\)=\\rmS,\(21\)where it subsequently terminates and returns\\rmS\\rmS\.

Beyond[Assumptions˜1](https://arxiv.org/html/2607.26357#Thmassumption1),[2](https://arxiv.org/html/2607.26357#Thmassumption2)and[3](https://arxiv.org/html/2607.26357#Thmassumption3), the recovery guarantee of[Algorithm˜3](https://arxiv.org/html/2607.26357#algorithm3)assumes the CIdecisions are*sufficiently accurate*—i\.e\. they reflect the true conditional independencies of\\pr\\pr\.

### 5\.2Computational Complexity of\\method

As\\methodis 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](https://arxiv.org/html/2607.26357#S5.SS2)with proof given in[Appendix˜D](https://arxiv.org/html/2607.26357#A4)\.

\{restatable\}

\[Computational Complexity of\\method\]corollarykombcomplexity Letn=\|\\verts\|n=\|\\verts\|denote the number of variables,NNthe sample size, andKKthe maximum size of any candidate MB\\rmS\\rmSencountered by[Algorithm˜1](https://arxiv.org/html/2607.26357#algorithm1)during execution\. Furthermore assume each CItest runs inO​\(N\)O\(N\)time\. Then, for a given orderkkand separator boundll,\\method\([Algorithm˜1](https://arxiv.org/html/2607.26357#algorithm1)\) terminates in worst\-case time

O​\(K​n​∑i=0k\(ni\)​\[∑i′=0k−i\(Ki′\)​∑j=0l\(K−i′j\)​N\]\),O\\\!\\left\(K\\,n\\sum\_\{i=0\}^\{k\}\\binom\{n\}\{i\}\\left\[\\sum\_\{i^\{\\prime\}=0\}^\{k\-i\}\\binom\{K\}\{i^\{\\prime\}\}\\sum\_\{j=0\}^\{l\}\\binom\{K\-i^\{\\prime\}\}\{j\}\\,N\\right\]\\right\),\(22\)where we adopt the convention that\(nm\)=0\\binom\{n\}\{m\}=0whenevern<mn<m\.

Treatingkkandllas constants,[Equation˜22](https://arxiv.org/html/2607.26357#S5.E22)is bounded byO​\(nk\+1​Kk\+l\+1​N\)O\(n^\{k\+1\}\\,K^\{k\+l\+1\}\\,N\)—polynomial innn,KK, andNN, and fork=0k=0it reduces to a scan that is linear innn, recovering the original GSalgorithm\. Crucially, the search for a dependence\-exposing subset and its separator is governed byKKrather thannn, a direct consequence of thell\-bounded separator assumption \([Assumption˜3](https://arxiv.org/html/2607.26357#Thmassumption3)\)\. The remaining exponential dependence onkkandllmotivates keeping them small\. The ablation study in[Appendix˜G](https://arxiv.org/html/2607.26357#A7)showcases how differing values ofkkandllaffects[Algorithm˜1](https://arxiv.org/html/2607.26357#algorithm1)runtime and its ability to find high\-order dependencies\. Generally, from the results in[Appendix˜G](https://arxiv.org/html/2607.26357#A7)for the benchmark datasets used, we found that settingk=lk=land havingk≤2k\\leq 2provides a decent trade\-off between runtime and finding high\-order dependencies\.

Table 1:Results on Synthetically Generated Data

## 6Experiments

In order to empirically test\\methodwe implemented the algorithm as aPythonpackage written inRust\. The implementation uses the G\-Test withα=0\.01\\alpha=0\.01as 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\\methodand code for the experiments can be found at[https://github\.com/lklee9/k\-order\-Markov\-blanket](https://github.com/lklee9/k-order-Markov-blanket)\. All experiments were run on an Apple M4 Mac mini\.

Throughout this section we will compare\\methodagainst 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 thepyCausalFSpython library\[yu2020a\]with a significance level of0\.010\.01for any statistical tests as well\.

When comparing between different MBdiscovery methods, we will use the F1 score defined as:F​1=2×𝑝𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛×𝑟𝑒𝑐𝑎𝑙𝑙/\(𝑝𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛\+𝑟𝑒𝑐𝑎𝑙𝑙\)F1=2\\times\\mathit\{precision\}\\times\\mathit\{recall\}/\(\\mathit\{precision\}\+\\mathit\{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\. ThereforeF​1=1F1=1represents the perfect precision and recall\.

![Refer to caption](https://arxiv.org/html/2607.26357v1/x1.png)Figure 4:Precision, recall, and F1 score of GSand \(\\method\-kk\-ll\) when recovering the MBsof the synthetic BNsin[Section˜6\.1](https://arxiv.org/html/2607.26357#S6.SS1), 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\.### 6\.1Experiments with Synthetic Bayesian Networks

To systematically test how\\methodbehaves in the face of true and empirical faithfulness violations, we will task\\methodat learning the MBof every variable for five different basic toy Bayesian networks\. These problems will have44variables—\\rmX=\{X1,X2,X3\}\\rmX=\\\{X\_\{1\},X\_\{2\},X\_\{3\}\\\}andYY—with a DAGthat has edges∀i∈\[3\]:Xi→Y\\forall i\\in\[3\]:X\_\{i\}\\rightarrow Y\. Each variableX∈\\rmXX\\in\\rmXis an independent Bernoulli random variable withp=0\.5p=0\.5\. The conditional distribution ofYYis then\\pr​\(Y=f​\(\\rmX\)∣\\rmX\)=0\.9\\pr\(Y=f\(\\rmX\)\\mid\\rmX\)=0\.9wheref​\(\\rmX\)f\(\\rmX\)is one of55different Boolean functions, one for each toy problem which we describe in[Appendix˜E](https://arxiv.org/html/2607.26357#A5)\.

We then sample100100samples from the55Bayesian networks1010different times, leading to10×510\\times 5samples of size100100—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 all44variables from the1010different samples separately\. We then averaged the F1 score over the10×410\\times 4runs and variables for each Bayesian network, resulting in[Table˜1](https://arxiv.org/html/2607.26357#S5.T1)\. Note that we only include the results for\\methodand Grow and Shrink \(GS\)because GSwas the best performing method out of all the previous methods mentioned in[Section˜6](https://arxiv.org/html/2607.26357#S6)for these synthetic datasets\.

From the results in[Table˜1](https://arxiv.org/html/2607.26357#S5.T1), we can observe that\\methodconsistently outperforms the other approaches at MBdiscovery\. More interestingly—even for the problems whose Bayesian networks are faithful—as we increasekkand look at higher order relationships,\\methodachieves a better F1 score and therefore is capable of better overcoming empirical unfaithfulness\. Furthermore the MBfor the parity problem is only recoverable whenk=2k=2, which confirms\\method’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](https://arxiv.org/html/2607.26357#S6.F4)traces the precision, recall, and F1 score as the sample size grows from5050to100100\. As expected only\\method\-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\\method, the higher the value ofkkused, the greater its performance in the low sample regime and the quicker it recovers the true MBsfor all variables\.

### 6\.2Experiments with Benchmark Bayesian Networks

To further explore how\\methodbehaves in more realistic scenarios, we then repeated the same procedure as[Section˜6\.1](https://arxiv.org/html/2607.26357#S6.SS1), but with some well\-known benchmark datasets444The datasets can be accessed at:[https://pages\.mtu\.edu/~lebrown/supplements/mmhc\_paper/mmhc\_index\.html](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 of50005000and we only tasked each method at finding the MBfor the1010variables with the largest neighborhoods in the Bayesian network\. Furthermore, we also omit some of the poorer performing methods in[Table˜2](https://arxiv.org/html/2607.26357#S6.T2)to conserve space\.

From the results in[Table˜2](https://arxiv.org/html/2607.26357#S6.T2)we can observe that\\methoddoes 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\\methodwith order22improves on order11, whereas on the high\-cardinality Barley and Mildew networks order11is best and order22is substantially worse—there the larger conditioning sets leave the higher\-order CItests underpowered at this sample size\. Runtime, in turn, grows with both the orderkkand the cardinality of the network\. Specifically, the baselines finish within roughly66s per dataset, while\\methodranges from comparable at order11on the lower\-cardinality networks \(Alarm1, Insurance\) to a few hundred times slower at order22on 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](https://arxiv.org/html/2607.26357#A7.T5)and[6](https://arxiv.org/html/2607.26357#A7.T6)under[Appendix˜G](https://arxiv.org/html/2607.26357#A7)\.

## 7Conclusion

We introduced akk\-order relaxation of faithfulness that explicitly accounts for parity\-type dependencies amongk\+2k\{\+\}2variables, and used it to derive a proof\-of\-concept Markov blanket \(MB\)discovery algorithm,\\method; with theoretical guarantees that\\methodrecovers the graphical MBunder the proposed assumptions\. Empirically,\\methodovercomes both true and empirical faithfulness violations on synthetic problems, and it performs competitively on commonly used benchmark datasets\.

As\\methodis 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\\method’s worst\-case cost—for example, through more intricate pruning rules that make\\methodscale better with larger values ofkk\. Secondly, although our ablation \([Appendix˜G](https://arxiv.org/html/2607.26357#A7)\) shows thatl=kl=kwithk≤2k\\leq 2is a robust default, principled*a priori*selection ofkkandllwithout expert knowledge remains open\. A practical alternative is to run\\methodfor several values ofkkandllwhile 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 AGraphoid Axioms

###### Definition 6\(Graphoid Axioms\)\.

For disjoint subsets\\rmX,\\rmY,\\rmZ,\\rmW⊆\\verts\\rmX,\\rmY,\\rmZ,\\rmW\\subseteq\\verts, a semi\-graphoid independence model\\indeps⋅\\indeps\{\\cdot\}obeys the following properties\[dawid1979,pearl1987\]:

1. 1\.Symmetry:\\rmX​\\indep​\\rmY∣\\rmZ⟹\\rmY​\\indep​\\rmX∣\\rmZ\\rmX\\indep\{\}\\rmY\\mid\\rmZ\\implies\\rmY\\indep\{\}\\rmX\\mid\\rmZ
2. 2\.Decomposition:\\rmX​\\indep​\\rmY∪\\rmW∣\\rmZ⟹\\rmX​\\indep​\\rmY∣\\rmZ\\rmX\\indep\{\}\\rmY\\cup\\rmW\\mid\\rmZ\\implies\\rmX\\indep\{\}\\rmY\\mid\\rmZ
3. 3\.Weak Union:\\rmX​\\indep​\\rmY∪\\rmW∣\\rmZ⟹\\rmX​\\indep​\\rmY∣\\rmW∪\\rmZ\\rmX\\indep\{\}\\rmY\\cup\\rmW\\mid\\rmZ\\implies\\rmX\\indep\{\}\\rmY\\mid\\rmW\\cup\\rmZ
4. 4\.Contraction:\(\\rmX​\\indep​\\rmY∣\\rmW∪\\rmZ\)∧\(\\rmX​\\indep​\\rmW∣\\rmZ\)⟹\\rmX​\\indep​\\rmY∪\\rmW∣\\rmZ\(\\rmX\\indep\{\}\\rmY\\mid\\rmW\\cup\\rmZ\)\\wedge\(\\rmX\\indep\{\}\\rmW\\mid\\rmZ\)\\implies\\rmX\\indep\{\}\\rmY\\cup\\rmW\\mid\\rmZ

The independence model is called a graphoid if it also obeys the following additional property\[pearl1985\]:

1. 5\.Intersection:\(\\rmX​\\indep​\\rmY∣\\rmZ∪\\rmW\)∧\(\\rmX​\\indep​\\rmW∣\\rmZ∪\\rmY\)⟹\\rmX​\\indep​\\rmY∪\\rmW∣\\rmZ\(\\rmX\\indep\{\}\\rmY\\mid\\rmZ\\cup\\rmW\)\\wedge\(\\rmX\\indep\{\}\\rmW\\mid\\rmZ\\cup\\rmY\)\\implies\\rmX\\indep\{\}\\rmY\\cup\\rmW\\mid\\rmZ

All graph independence models\\indeps​\\gr\\indeps\{\\gr\}are guaranteed to obey all the graphoid axioms\[lauritzen2018\]\. However, probabilistic independence models\\indeps​\\pr\\indeps\{\\pr\}are only guaranteed to obey the semi\-graphoid axioms\. However if the joint distribution is strictly positive, then\\indeps​\\pr\\indeps\{\\pr\}is guaranteed to be a graphoid as well\[pearl1988,pearl1989\]\.

## Appendix BProofs

###### Proof\.

Consider the family of association witnesses contained in\\rmZ\\rmZ,

𝒲≔\{\\rmW⊆\\rmZ:\\kassoc​\\rmS​Y,X​\\rmW\}\.\\mathcal\{W\}\\coloneq\\\{\\rmW\\subseteq\\rmZ:\\kassoc\{\\rmS\}\{Y,X\}\{\\rmW\}\\\}\.By the hypothesis\\kassoc​\\rmS​Y,X​\\rmZ\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ\}we have\\rmZ∈𝒲\\rmZ\\in\\mathcal\{W\}, so𝒲\\mathcal\{W\}is nonempty; and since\\rmZ\\rmZis finite, so is𝒲\\mathcal\{W\}\. Hence𝒲\\mathcal\{W\}contains an element\\rmZ′\\rmZ^\{\\prime\}of minimum cardinality\. Every proper subset\\rmZ∗⊊\\rmZ′\\rmZ^\{\*\}\\subsetneq\\rmZ^\{\\prime\}satisfies\|\\rmZ∗\|<\|\\rmZ′\|\|\\rmZ^\{\*\}\|<\|\\rmZ^\{\\prime\}\|and\\rmZ∗⊆\\rmZ\\rmZ^\{\*\}\\subseteq\\rmZ, so by the minimality of\\rmZ′\\rmZ^\{\\prime\}we have\\rmZ∗∉𝒲\\rmZ^\{\*\}\\notin\\mathcal\{W\}, i\.e\.¬\\kassoc​\\rmS​Y,X​\\rmZ∗\\neg\\,\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ^\{\*\}\}\. Together with\\kassoc​\\rmS​Y,X​\\rmZ′\\kassoc\{\\rmS\}\{Y,X\}\{\\rmZ^\{\\prime\}\}, this is by[Equation˜7](https://arxiv.org/html/2607.26357#S4.E7)exactly\\kdep​\\rmS​Y,X​\\rmZ′\\kdep\{\\rmS\}\{Y,X\}\{\\rmZ^\{\\prime\}\}\. As\\rmZ′⊆\\rmZ\\rmZ^\{\\prime\}\\subseteq\\rmZ, this establishes[Equation˜10](https://arxiv.org/html/2607.26357#S4.E10)\. ∎

###### Proof\.

Assume\\kdep​∅​Y,X​\\rmZ\\kdep\{\\emptyset\}\{Y,X\}\{\\rmZ\}\. For the variableXX, the inter\-association clause is\\kassoc​∅​Y,X​\\rmZ\\kassoc\{\\emptyset\}\{Y,X\}\{\\rmZ\}, i\.e\.Y​\\dep​\\pr​X∣\\rmZY\\dep\{\\pr\}X\\mid\\rmZ, which is exactly the first conjunct of[Equation˜7](https://arxiv.org/html/2607.26357#S4.E7)and hence immediate\.

Now suppose, towards a contradiction, that the clause fails for someZ∈\\rmZZ\\in\\rmZ, i\.e\.Y​\\indep​\\pr​Z∣X∪\(\\rmZ∖Z\)Y\\indep\{\\pr\}Z\\mid X\\cup\(\\rmZ\\setminus Z\)\. Since\\rmZ∖Z\\rmZ\\setminus Zis a proper subset of\\rmZ\\rmZand the separating set is empty, the second conjunct of[Equation˜7](https://arxiv.org/html/2607.26357#S4.E7)givesY​\\indep​\\pr​X∣\\rmZ∖ZY\\indep\{\\pr\}X\\mid\\rmZ\\setminus Z\. Applying the contraction and then weak union graphoid axioms of[Definition˜6](https://arxiv.org/html/2607.26357#Thmdefinition6),

\(Y​\\indep​\\pr​Z∣X∪\(\\rmZ∖Z\)\)∧\(Y​\\indep​\\pr​X∣\\rmZ∖Z\)\\displaystyle\(Y\\indep\{\\pr\}Z\\mid X\\cup\(\\rmZ\\setminus Z\)\)\\wedge\(Y\\indep\{\\pr\}X\\mid\\rmZ\\setminus Z\)\(23\)⟹Y​\\indep​\\pr​\{X,Z\}∣\\rmZ∖Z⟹Y​\\indep​\\pr​X∣\\rmZ,\\displaystyle\\implies Y\\indep\{\\pr\}\\\{X,Z\\\}\\mid\\rmZ\\setminus Z\\implies Y\\indep\{\\pr\}X\\mid\\rmZ,contradictingY​\\dep​\\pr​X∣\\rmZY\\dep\{\\pr\}X\\mid\\rmZ\. HenceY​\\dep​\\pr​Z∣X∪\(\\rmZ∖Z\)Y\\dep\{\\pr\}Z\\mid X\\cup\(\\rmZ\\setminus Z\)for everyZ∈\\rmZZ\\in\\rmZ, which together withY​\\dep​\\pr​X∣\\rmZY\\dep\{\\pr\}X\\mid\\rmZis precisely\\kinter​∅​Y,\\rmZ∪X\\kinter\{\\emptyset\}\{Y,\\rmZ\\cup X\}\. ∎

## Appendix CCorrectness of\\method

###### Proof\.

The functionfind\_condfirst iterates through all the subsets\\rmX′∈\\powerleq​k−\|\\rmX\|​\\rmS\\rmX^\{\\prime\}\\in\\powerleq\{k\-\|\\rmX\|\}\{\\rmS\}and for each\\rmX′\\rmX^\{\\prime\}whereY​\\dep​Phys\.Rev\.E​X∣\\rmX∪\\rmX′Y\\dep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmX\\cup\\rmX^\{\\prime\}\(line 9 infind\_cond\), it passes\\rmX′\\rmX^\{\\prime\}to the functionno\_seps\(line 10 infind\_cond\) to determine if there are any subsets\\rmS′∈\\powerleq​l​\\rmS∖\\rmX′\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmX^\{\\prime\}\}that results inY​\\indep​Phys\.Rev\.E​X∣\\rmX∪\\rmX′∪\\rmS′Y\\indep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmX\\cup\\rmX^\{\\prime\}\\cup\\rmS^\{\\prime\}\(line 4 inno\_seps\)\. As soon asno\_sepsfinds such a\\rmS′\\rmS^\{\\prime\}, it returnsFalse\(line 4 inno\_seps\), andfind\_condwill need to go to the next iteration\. Otherwise, ifno\_sepsis unable to find any such subset\\rmS′\\rmS^\{\\prime\}, it returnsTrue\(line 5 inno\_seps\), andfind\_condthen immediately returns\\rmZ=\\rmX′∪\\rmX\\rmZ=\\rmX^\{\\prime\}\\cup\\rmX\(line 11 infind\_cond\) as there are noll\-bounded subsets\\rmS′∈\\powerleq​l​\\rmS∖\\rmZ\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmZ\}that renderYYandXXconditionally independent given\\rmZ∪\\rmS′\\rmZ\\cup\\rmS^\{\\prime\}\. Here we use that the external set\\rmX\\rmXis disjoint from\\rmS\\rmS, so\\rmS∖\\rmZ=\\rmS∖\\rmX′\\rmS\\setminus\\rmZ=\\rmS\\setminus\\rmX^\{\\prime\}; and since∅∈\\powerleq​l​\\rmS∖\\rmX′\\emptyset\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmX^\{\\prime\}\}, the dependence checkY​\\dep​Phys\.Rev\.E​X∣\\rmZY\\dep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmZon line 9 is precisely the\\rmS′=∅\\rmS^\{\\prime\}=\\emptysetinstance tested byno\_seps, so a returned\\rmZ\\rmZsatisfies[Equation˜18](https://arxiv.org/html/2607.26357#S5.E18)in full\. However, iffind\_condcompletes all its iteration without returning early, then we know that all subset\\rmX′\\rmX^\{\\prime\}has somell\-bounded subset\\rmS′\\rmS^\{\\prime\}that will causeY​\\indep​Phys\.Rev\.E​X∣\\rmX′∪\\rmX∪\\rmS′Y\\indep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmX^\{\\prime\}\\cup\\rmX\\cup\\rmS^\{\\prime\}\. Thereforefind\_condreturns the failure value⊥\\botin this case \(line 12 infind\_cond\)\. ∎

###### Proof\.

First consider a call that returns a nonempty set\.find\_inter\_depdoes so only by returningmake\_strict\(Y,X,\\rmZ,\\rmS,k,l\)\(Y,X,\\rmZ,\\rmS,k,l\)after a callfind\_cond\(Y,X,\\rmX,\\rmS,k,l\)\(Y,X,\\rmX,\\rmS,k,l\)has returned some\\rmZ≠⊥\\rmZ\\neq\\bot, for someX∈\\verts∖\(Y∪\\rmS∪\\rmX\)X\\in\\verts\\setminus\(Y\\cup\\rmS\\cup\\rmX\)\(as\\rmS\\rmSis passed unchanged through every recursive call, any suchXXsatisfiesX∈\\verts∖\(Y∪\\rmS\)=\\verts′X\\in\\verts\\setminus\(Y\\cup\\rmS\)=\\verts^\{\\prime\}, matching the theorem’s candidate set\)\. By[Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2),\\rmZ=\\rmX∪\\rmX′\\rmZ=\\rmX\\cup\\rmX^\{\\prime\}for some\\rmX′∈\\powerleq​k−\|\\rmX\|​\\rmS\\rmX^\{\\prime\}\\in\\powerleq\{k\-\|\\rmX\|\}\{\\rmS\}such thatY​\\dep​Phys\.Rev\.E​X∣\\rmZ∪\\rmS′Y\\dep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}holds for all\\rmS′∈\\powerleq​l​\\rmS∖\\rmX′\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmX^\{\\prime\}\}\. As the dependence set\\rmX\\rmXis disjoint from\\rmS\\rmS, we have\\rmS∖\\rmX′=\\rmS∖\\rmZ\\rmS\\setminus\\rmX^\{\\prime\}=\\rmS\\setminus\\rmZ, soY​\\dep​Phys\.Rev\.E​X∣\\rmZ∪\\rmS′Y\\dep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}holds for all\\rmS′∈\\powerleq​l​\\rmS∖\\rmZ\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmZ\}with\\rmZ∈\\powerleq​k​\\verts∖\{Y,X\}\\rmZ\\in\\powerleq\{k\}\{\\verts\\setminus\\\{Y,X\\\}\}—the estimatedll\-bounded association claimed in the theorem\. Sincemake\_strictonly removes variables from\\rmZ\\rmZ, the returned set isX∪\\rmZRX\\cup\\rmZ\_\{R\}for some\\rmZR⊆\\rmZ\\rmZ\_\{R\}\\subseteq\\rmZ, as claimed\.

Next consider a call that returns the empty set\. Here every internalfind\_condcall must have failed, sofind\_inter\_deptraverses the full enumeration tree of the dependence sets\\rmX∈\\powerleq​k​\\verts∖\(Y∪\\rmS\)\\rmX\\in\\powerleq\{k\}\{\\verts\\setminus\(Y\\cup\\rmS\)\}\(lines 13–17 offind\_inter\_dep\), and for each visits every candidateX∈\\verts∖\(Y∪\\rmS∪\\rmX\)X\\in\\verts\\setminus\(Y\\cup\\rmS\\cup\\rmX\)throughfind\_cond\(Y,X,\\rmX,\\rmS,k,l\)\(Y,X,\\rmX,\\rmS,k,l\), which by[Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2)searches all\\rmX′∈\\powerleq​k−\|\\rmX\|​\\rmS\\rmX^\{\\prime\}\\in\\powerleq\{k\-\|\\rmX\|\}\{\\rmS\}\. Every conditioning set\\rmZ∈\\powerleq​k​\\verts∖\{Y,X\}\\rmZ\\in\\powerleq\{k\}\{\\verts\\setminus\\\{Y,X\\\}\}decomposes uniquely into its part outside and its part inside the candidate blanket,\\rmZ=\(\\rmZ∖\\rmS\)∪\(\\rmZ∩\\rmS\)=\\rmX∪\\rmX′\\rmZ=\(\\rmZ\\setminus\\rmS\)\\cup\(\\rmZ\\cap\\rmS\)=\\rmX\\cup\\rmX^\{\\prime\}, so this enumeration covers every pair\(X,\\rmZ\)\(X,\\rmZ\)withX∈\\verts′X\\in\\verts^\{\\prime\}and\\rmZ∈\\powerleq​k​\\verts∖\{Y,X\}\\rmZ\\in\\powerleq\{k\}\{\\verts\\setminus\\\{Y,X\\\}\}\. By[Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2), the call for\(X,\\rmX\)\(X,\\rmX\)returns⊥\\botexactly when every such\\rmX′\\rmX^\{\\prime\}admits a separator\\rmS′∈\\powerleq​l​\\rmS∖\\rmX′\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmX^\{\\prime\}\}withY​\\indep​Phys\.Rev\.E​X∣\\rmX∪\\rmX′∪\\rmS′Y\\indep\{\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}\}X\\mid\\rmX\\cup\\rmX^\{\\prime\}\\cup\\rmS^\{\\prime\}\. Thereforefind\_inter\_depreturns the empty set if and only if[Equation˜19](https://arxiv.org/html/2607.26357#S5.E19)holds\. The distinguished value⊥\\botlets the callers tell such a failure from a success that returns the empty conditioning set\. ∎

###### Proof\.

Throughout we assume that the CIdecisions are accurate, soPhys\.Rev\.E\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}and\\pr\\pragree on every tested statement\.

In the grow phase, every nonempty set returned byfind\_inter\_depcontains a variableX∉\\rmSX\\notin\\rmSby[Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2), so each iteration of the repeat loop strictly grows\\rmS\\rmSand the phase terminates\. Upon terminationfind\_inter\_depreturned the empty set, so by[Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2)the condition in[Equation˜19](https://arxiv.org/html/2607.26357#S5.E19)holds inPhys\.Rev\.E\{\\rm Phys\.\\penalty 10000\\ Rev\.\\penalty 10000\\ E\}, and by accuracy in\\pr\\pr\. Fix anyX∈\\verts∖\(\\rmS∪Y\)X\\in\\verts\\setminus\(\\rmS\\cup Y\)\. For every\\rmZ∈\\powerleq​k​\\verts∖\{Y,X\}\\rmZ\\in\\powerleq\{k\}\{\\verts\\setminus\\\{Y,X\\\}\},[Equation˜19](https://arxiv.org/html/2607.26357#S5.E19)supplies an\\rmS′∈\\powerleq​l​\\rmS∖\\rmZ⊆\\powerleq​l​\\rmS\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmZ\}\\subseteq\\powerleq\{l\}\{\\rmS\}withY​\\indep​\\pr​X∣\\rmZ∪\\rmS′Y\\indep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}; thell\-bounded separator assumption \([Assumption˜3](https://arxiv.org/html/2607.26357#Thmassumption3)\) then lifts this, for every such\\rmZ\\rmZ, to the unbounded statement that some\\rmS′⊆\\rmS\\rmS^\{\\prime\}\\subseteq\\rmSgivesY​\\indep​\\pr​X∣\\rmZ∪\\rmS′Y\\indep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}, which is the antecedent ofkk\-order faithfulness \([Equation˜11](https://arxiv.org/html/2607.26357#S4.E11)\)\. HenceY​\\indep​\\gr​X∣\\rmSY\\indep\{\\gr\}X\\mid\\rmSfor everyX∈\\verts∖\(\\rmS∪Y\)X\\in\\verts\\setminus\(\\rmS\\cup Y\)\. Since the parents and children ofYY—and, for an undirected\\gr\\gr, its neighbours—can never be separated fromYY, they must all lie in\\rmS\\rmS\. Any spouseXXshares a childC∈ch​\(Y\)⊆\\rmSC\\in\\text\{ch\}\(Y\)\\subseteq\\rmS, so the collider pathY→C←XY\\rightarrow C\\leftarrow Xis active given\\rmS\\rmS; henceY​\\dep​\\gr​X∣\\rmSY\\dep\{\\gr\}X\\mid\\rmS, which forcesX∈\\rmSX\\in\\rmSas well\. Therefore\\mbs​\\gr​\(Y\)⊆\\rmS\\mbs\{\\gr\}\(Y\)\\subseteq\\rmS, establishing[Equation˜20](https://arxiv.org/html/2607.26357#S5.E20)\.

During the shrink phase we maintain the invariant\\mbs​\\gr​\(Y\)⊆\\rmS\\mbs\{\\gr\}\(Y\)\\subseteq\\rmS, which holds on entry by[Equation˜20](https://arxiv.org/html/2607.26357#S5.E20)\. First, we show that no true member is ever removed\. ForX∈\\mbs​\\gr​\(Y\)X\\in\\mbs\{\\gr\}\(Y\), the in\-blanket witness hypothesis of[Algorithm˜3](https://arxiv.org/html/2607.26357#algorithm3)supplies a set\\rmZ∈\\powerleq​k​\\mbs​\\gr​\(Y\)∖X⊆\\powerleq​k​\\rmS∖X\\rmZ\\in\\powerleq\{k\}\{\\mbs\{\\gr\}\(Y\)\\setminus X\}\\subseteq\\powerleq\{k\}\{\\rmS\\setminus X\}such thatY​\\dep​\\pr​X∣\\rmZ∪\\rmS′Y\\dep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}holds for all\\rmS′∈\\powerleq​l​\(\\rmS∖X\)∖\\rmZ\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\(\\rmS\\setminus X\)\\setminus\\rmZ\}\. By accuracy and[Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2)\(with candidate blanket\\rmS∖X\\rmS\\setminus X, external dependence set\\rmX=∅\\rmX=\\emptyset, and\\rmX′=\\rmZ\\rmX^\{\\prime\}=\\rmZ\), the callfind\_cond\(Y,X,∅,\\rmS∖X,k,l\)\(Y,X,\\emptyset,\\rmS\\setminus X,k,l\)does not return⊥\\bot, soXXis never selected for removal \(a member whose only witness is\\rmZ=∅\\rmZ=\\emptysetmakesfind\_condreturn the empty set rather than⊥\\bot\)\.

Next, we show that every non\-member is removed once selected\. LetX∈\\rmS∖\\mbs​\\gr​\(Y\)X\\in\\rmS\\setminus\\mbs\{\\gr\}\(Y\); by the invariant\\mbs​\\gr​\(Y\)⊆\\rmS∖X\\mbs\{\\gr\}\(Y\)\\subseteq\\rmS\\setminus X\. We first observe that\\mbs​\\gr​\(Y\)\\mbs\{\\gr\}\(Y\)separatesYYfromXXeven after conditioning on any further set\\rmZ⊆\\verts∖\{Y,X\}\\rmZ\\subseteq\\verts\\setminus\\\{Y,X\\\}\. SinceX∉\\mbs​\\gr​\(Y\)X\\notin\\mbs\{\\gr\}\(Y\),XXis not adjacent toYY, so every path between them has at least one interior vertex; write such a path asY=N0∼N1∼⋯∼Nm=XY=N\_\{0\}\\sim N\_\{1\}\\sim\\cdots\\sim N\_\{m\}=Xwithm≥2m\\geq 2\. Its first interior vertexN1N\_\{1\}is adjacent toYYand hence lies in\\mbs​\\gr​\(Y\)\\mbs\{\\gr\}\(Y\)\. IfN1N\_\{1\}is a non\-collider on the path, it blocks the path\. If insteadN1N\_\{1\}is a colliderY→N1←N2Y\\rightarrow N\_\{1\}\\leftarrow N\_\{2\}, thenN1N\_\{1\}is a common child ofYYand its successorN2N\_\{2\}on the path, soN2N\_\{2\}is a spouse ofYYand henceN2∈\\mbs​\\gr​\(Y\)N\_\{2\}\\in\\mbs\{\\gr\}\(Y\)\. In particularN2≠XN\_\{2\}\\neq X\(asX∉\\mbs​\\gr​\(Y\)X\\notin\\mbs\{\\gr\}\(Y\)\), soN2N\_\{2\}is a genuine interior vertex; and since the edgeN2→N1N\_\{2\}\\rightarrow N\_\{1\}points out ofN2N\_\{2\}, it is a non\-collider on the path and blocks it\. For an undirected\\gr\\gr,N1N\_\{1\}is instead a neighbour ofYYand blocks the path directly\. As each path is thus blocked at a conditioned vertex of\\mbs​\\gr​\(Y\)\\mbs\{\\gr\}\(Y\)irrespective of\\rmZ\\rmZ\[pearl1988,lauritzen1996\],Y​\\indep​\\gr​X∣\\mbs​\\gr​\(Y\)∪\\rmZY\\indep\{\\gr\}X\\mid\\mbs\{\\gr\}\(Y\)\\cup\\rmZfor every\\rmZ∈\\powerleq​k​\\verts∖\{Y,X\}\\rmZ\\in\\powerleq\{k\}\{\\verts\\setminus\\\{Y,X\\\}\}\.

The Global Markov Property transfers this to\\pr\\prwith the separator\\rmS′≔\\mbs​\\gr​\(Y\)∖\\rmZ⊆\\rmS∖X\\rmS^\{\\prime\}\\coloneq\\mbs\{\\gr\}\(Y\)\\setminus\\rmZ\\subseteq\\rmS\\setminus X, so the unbounded side of[Assumption˜3](https://arxiv.org/html/2607.26357#Thmassumption3)—at the separating set\\rmS∖X\\rmS\\setminus X—holds for every such\\rmZ\\rmZ\. The biconditional then yields some\\rmS′∈\\powerleq​l​\\rmS∖X\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus X\}; discarding any overlap with\\rmZ\\rmZleaves\\rmZ∪\\rmS′\\rmZ\\cup\\rmS^\{\\prime\}unchanged and gives\\rmS′∈\\powerleq​l​\(\\rmS∖X\)∖\\rmZ\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\(\\rmS\\setminus X\)\\setminus\\rmZ\}withY​\\indep​\\pr​X∣\\rmZ∪\\rmS′Y\\indep\{\\pr\}X\\mid\\rmZ\\cup\\rmS^\{\\prime\}\. By accuracy and[Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2),find\_cond\(Y,X,∅,\\rmS∖X,k,l\)\(Y,X,\\emptyset,\\rmS\\setminus X,k,l\)returns⊥\\bot, soXXis removed once selected\.

Each iteration of the while loop removes one variable, so the phase terminates; it exits only when noX∈\\rmSX\\in\\rmSyields⊥\\bot, whence\\rmS∖\\mbs​\\gr​\(Y\)=∅\\rmS\\setminus\\mbs\{\\gr\}\(Y\)=\\emptysetand, with the invariant,\\rmS=\\mbs​\\gr​\(Y\)\\rmS=\\mbs\{\\gr\}\(Y\), establishing[Equation˜21](https://arxiv.org/html/2607.26357#S5.E21)\. ∎

## Appendix DComplexity of\\method

###### Proof\.

We chargeO​\(N\)O\(N\)to each CItest and bound the running time by analysing the three sub\-algorithms of[Section˜5](https://arxiv.org/html/2607.26357#S5)from the innermost outwards\.

First considerfind\_cond\([Algorithm˜3](https://arxiv.org/html/2607.26357#algorithm3)\) for a fixed candidate variableXXand external dependence set\\rmX\\rmXwith\|\\rmX\|=i\|\\rmX\|=i\. It iterates over every subset\\rmX′∈\\powerleq​k−i​\\rmS\\rmX^\{\\prime\}\\in\\powerleq\{k\-i\}\{\\rmS\}, of which there are at most∑i′=0k−i\(Ki′\)\\sum\_\{i^\{\\prime\}=0\}^\{k\-i\}\\binom\{K\}\{i^\{\\prime\}\}since\|\\rmS\|≤K\|\\rmS\|\\leq K\. For each\\rmX′\\rmX^\{\\prime\}it performs one CItest and, whenever that test indicates a dependence, callsno\_seps, which iterates over every separator\\rmS′∈\\powerleq​l​\\rmS∖\\rmX′\\rmS^\{\\prime\}\\in\\powerleq\{l\}\{\\rmS\\setminus\\rmX^\{\\prime\}\}and performs one further CItest per separator\. Writing\|\\rmX′\|=i′\|\\rmX^\{\\prime\}\|=i^\{\\prime\}, there are at most∑j=0l\(K−i′j\)\\sum\_\{j=0\}^\{l\}\\binom\{K\-i^\{\\prime\}\}\{j\}such separators, so a single call tofind\_condcosts

O​\(∑i′=0k−i\(Ki′\)​∑j=0l\(K−i′j\)​N\),O\\\!\\left\(\\sum\_\{i^\{\\prime\}=0\}^\{k\-i\}\\binom\{K\}\{i^\{\\prime\}\}\\sum\_\{j=0\}^\{l\}\\binom\{K\-i^\{\\prime\}\}\{j\}\\,N\\right\),\(24\)where the convention\(nm\)=0\\binom\{n\}\{m\}=0forn<mn<maccounts for the separator set\\rmS∖\\rmX′\\rmS\\setminus\\rmX^\{\\prime\}being exhausted\.

Next considerfind\_inter\_dep\([Algorithm˜2](https://arxiv.org/html/2607.26357#algorithm2)\)\. It explores the enumeration tree of\\verts∖Y\\verts\\setminus Yup to subsets of cardinalitykk\(lines 13–17 offind\_inter\_dep\), and hence visits at most∑i=0k\(ni\)\\sum\_\{i=0\}^\{k\}\\binom\{n\}\{i\}dependence sets\\rmX∈\\powerleq​k​\\verts∖Y\\rmX\\in\\powerleq\{k\}\{\\verts\\setminus Y\}\(a loose bound on the exact count∑i=0k\(n−1i\)\\sum\_\{i=0\}^\{k\}\\binom\{n\-1\}\{i\}, as\|\\verts∖Y\|=n−1\|\\verts\\setminus Y\|=n\-1\)\. For every such\\rmX\\rmXit loops over theO​\(n\)O\(n\)candidate variablesX∈\\verts∖\(Y∪\\rmS∪\\rmX\)X\\in\\verts\\setminus\(Y\\cup\\rmS\\cup\\rmX\)and invokesfind\_cond\. Multiplying these two factors by[Equation˜24](https://arxiv.org/html/2607.26357#A4.E24)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\\\!\\left\(n\\sum\_\{i=0\}^\{k\}\\binom\{n\}\{i\}\\left\[\\sum\_\{i^\{\\prime\}=0\}^\{k\-i\}\\binom\{K\}\{i^\{\\prime\}\}\\sum\_\{j=0\}^\{l\}\\binom\{K\-i^\{\\prime\}\}\{j\}\\,N\\right\]\\right\)\.\(25\)Each successful sweep additionally invokesmake\_strictexactly once; it iterates over the found dependence set\\rmZ\\rmZ—of size at mostmin⁡\(k,n−2\)\\min\(k,n\-2\), since\\rmZ⊆\\verts∖\{Y,X\}\\rmZ\\subseteq\\verts\\setminus\\\{Y,X\\\}and\|\\rmZ\|≤k\|\\rmZ\|\\leq k—and performs a singleno\_sepscall per element, at costO​\(min⁡\(k,n−2\)​∑j=0l\(Kj\)​N\)O\\\!\\left\(\\min\(k,n\-2\)\\sum\_\{j=0\}^\{l\}\\binom\{K\}\{j\}\\,N\\right\)\. This is dominated by the per\-sweepfind\_condtotal already counted in[Equation˜25](https://arxiv.org/html/2607.26357#A4.E25)—whose bracketed factor already contains the term∑j=0l\(Kj\)\\sum\_\{j=0\}^\{l\}\\binom\{K\}\{j\}\(ati′=0i^\{\\prime\}=0\) multiplied by the outer factorn​∑i=0k\(ni\)≥n≥min⁡\(k,n−2\)n\\sum\_\{i=0\}^\{k\}\\binom\{n\}\{i\}\\geq n\\geq\\min\(k,n\-2\)—and is therefore absorbed into the bound\.

Finally, consider\\methoditself \([Algorithm˜1](https://arxiv.org/html/2607.26357#algorithm1)\)\. Each successful call tofind\_inter\_depin the grow phase returns a nonempty setX∪\\rmZX\\cup\\rmZthat is added to\\rmS\\rmS, so the candidate Markov blanket grows by at least one variable per successful call; as\|\\rmS\|≤K\|\\rmS\|\\leq Kthroughout, the grow phase performs at mostO​\(K\)O\(K\)association\-mining sweeps\. The shrink phase re\-evaluates its guard after each removal, callingfind\_condat mostO​\(K2\)O\(K^\{2\}\)times over the at mostKKremovals; asK≤n−1K\\leq n\-1, this cost is dominated by that of the grow phase\. Multiplying[Equation˜25](https://arxiv.org/html/2607.26357#A4.E25)by theseO​\(K\)O\(K\)sweeps yields the bound in[Equation˜22](https://arxiv.org/html/2607.26357#S5.E22), completing the proof\. ∎

## Appendix EBoolean Functions for the Synthetic Bayesian Networks

parity:f​\(X1,…,X3\)=\(∑i=13Xi\)\(mod2\)f\(X\_\{1\},\\ldots,X\_\{3\}\)=\(\\sum\_\{i=1\}^\{3\}X\_\{i\}\)\\pmod\{2\}

ex\-ss:f​\(X1,…,X3;s\)=\[\(∑i=13Xi\)=s\]f\(X\_\{1\},\\ldots,X\_\{3\};s\)=\\bigl\[\(\\sum\_\{i=1\}^\{3\}X\_\{i\}\)=s\\bigr\]

and:f​\(X1,…,X3\)=\[\(∑i=13Xi\)=3\]f\(X\_\{1\},\\ldots,X\_\{3\}\)=\\bigl\[\(\\sum\_\{i=1\}^\{3\}X\_\{i\}\)=3\\bigr\]

or:f​\(X1,…,X3\)=\[\(∑i=13Xi\)≥1\]f\(X\_\{1\},\\ldots,X\_\{3\}\)=\\bigl\[\(\\sum\_\{i=1\}^\{3\}X\_\{i\}\)\\geq 1\\bigr\]

wheres∈\{1,2\}s\\in\\\{1,2\\\}\.

## Appendix FSeparator Sizes in the Benchmark Datasets

To gauge how restrictive thell\-bounded separator assumption \([Assumption˜3](https://arxiv.org/html/2607.26357#Thmassumption3)\) is in practice, we examine the benchmark networks of[Section˜6\.2](https://arxiv.org/html/2607.26357#S6.SS2)\. 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](https://arxiv.org/html/2607.26357#A6.T3)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 exceeds2\.192\.19and size\-≤3\\leq\\\!3separators cover at least83%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 size99; this is also the network on which\\methodis the most expensive and least accurate \([Section˜6\.2](https://arxiv.org/html/2607.26357#S6.SS2)\), illustrating that the difficulty of MBdiscovery tracks the separator sizes the assumption must accommodate\. Overall, a modest separator boundllalready 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\\leq 3” is the fraction whose minimal separator has size at most three\.
## Appendix GAblation

For completeness, we report the full results of the experiments in[Section˜6\.2](https://arxiv.org/html/2607.26357#S6.SS2)for every method and every\(k,l\)\(k,l\)configuration of\\method\.[Table˜5](https://arxiv.org/html/2607.26357#A7.T5)gives the performance of all methods, while[Table˜6](https://arxiv.org/html/2607.26357#A7.T6)gives their runtimes in seconds\. In both tables, “–” denotes a configuration that failed to complete within the allotted time budget\.

[Table˜4](https://arxiv.org/html/2607.26357#A7.T4)summarises the F1 score of\\methodacross all synthetic and benchmark datasets for every\(k,l\)\(k,l\)configuration\. F1 generally improves fromk=0k=0tok=2k=2—higher orders expose dependencies that lower orders miss—but the gains flatten or reverse beyondk=2k=2as the larger conditioning sets incur both higher runtime \([Table˜6](https://arxiv.org/html/2607.26357#A7.T6)\) and more error\-prone CItests\. Settingl=kl=kwithk≤2k\\leq 2therefore offers a robust default, the heuristic we adopt in[Section˜5\.2](https://arxiv.org/html/2607.26357#S5.SS2)\. For a fixedkk, raisingllbeyondkkdoes not consistently improve F1 \([Table˜4](https://arxiv.org/html/2607.26357#A7.T4)\), which motivates simply settingl=kl=krather than tuningllseparately\.

Table 4:F1 score across datasets for each\(k,l\)\(k,l\)configuration of\\method\.Table 5:Full results on benchmark datasets\.Table 6:Full runtimes \(in seconds\) on benchmark datasets\.