# Learning the Finer Things: Bayesian Structure Learning at the Instantiation
Level

Chase Yakaboski, Eugene Santos, Jr

###### Abstract

Successful machine learning methods require a trade-off between memorization
and generalization. Too much memorization and the model cannot generalize to
unobserved examples. Too much over-generalization and we risk under-fitting
the data. While we commonly measure their performance through cross validation
and accuracy metrics, how should these algorithms cope in domains that are
extremely under-determined where accuracy is always unsatisfactory? We present
a novel probabilistic graphical model structure learning approach that can
learn, generalize and explain in these elusive domains by operating at the
random variable instantiation level. Using Minimum Description Length (MDL)
analysis, we propose a new decomposition of the learning problem over all
training exemplars, fusing together minimal entropy inferences to construct a
final knowledge base. By leveraging Bayesian Knowledge Bases (BKBs), a
framework that operates at the instantiation level and inherently subsumes
Bayesian Networks (BNs), we develop both a theoretical MDL score and
associated structure learning algorithm that demonstrates significant
improvements over learned BNs on 40 benchmark datasets. Further, our algorithm
incorporates recent off-the-shelf DAG learning techniques enabling tractable
results even on large problems. We then demonstrate the utility of our
approach in a significantly under-determined domain by learning gene
regulatory networks on breast cancer gene mutational data available from The
Cancer Genome Atlas (TCGA).

##  1 Introduction

††Github URL: https://github.com/di2ag/pybkb

Since popularization by Pearl (1986), learning Bayesian Networks (BNs) has
solidified into a steadfast research area for 40 years. It has become an
important paradigm for modeling and reasoning under uncertainty and has seen
applications from stock market prediction (Malagrino, Roman, and Monteiro
2018) and medical diagnosis (Shih, Choi, and Darwiche 2018) to Gene Regulatory
Networks (GRNs) (Sauta et al. 2020). Despite Bayesian Network Structure
Learning (BNSL) being NP-hard (Chickering, Heckerman, and Meek 2004) and even
simpler structures like polytrees being NP-hard(er) (Dasgupta 1999), new
constraints (Grüttemeier and Komusiewicz 2020), improvements (Trösser, de
Givry, and Katsirelos 2021), and scalings (Scanagatta et al. 2015) are
presented at major AI conferences every year. This is because BNs and
affliated structures like Markov (Koller and Friedman 2009) and Dependency
Networks (DNs) (Heckerman et al. 2000) offer a quality that other methods such
as deep learning do not; explainability (Došilović, Brčić, and Hlupić 2018;
Burkart and Huber 2021).

The optimization that occurs in Probabilistic Graphical Model (PGM) structure
learning is

$\displaystyle G^{*}=\operatorname*{argmax}_{G}F(G,D)$
$\displaystyle\text{subject to }G\in\mathbb{\Omega}$

where $D$ is a database, $G$ is a graph structure such as a BN, $F$ is a
scoring function that yields the goodness of fit of the structure $G$, and
$\mathbb{\Omega}$ is the set of allowed structures for $G$; for BNs this would
be the space of all possible Directed Acyclic Graphs (DAGs).

Scoring functions are essential to the structure learning problem and should
have a theoretical justification in information theory or otherwise. For
instance, the most common scoring functions such as Bayesian Information
Criteria (BIC) (Schwarz 1978), Minimum Description Length (MDL) (Rissanen
1978), and Akaike Information Criterion (Akaike 1974) are all based on
information theoretic criteria or can be viewed from this perspective. While
we will spend part of this paper in theoretically justifying our model scoring
approach, our goal is _not_ to present a better scoring function. Instead, our
goal is to illustrate that no matter the scoring function or learning
algorithm, an over-generalization is encountered when modeling at the Random
Variable (RV) level.

By operating at the RV level, models force a complete distribution, as is the
case with BNs. While a complete distribution is often desired, this has an
unintended over-generalization consequence, particularly in under-determined
domains. This phenomenon even occurs in deep learning systems, and is
generally referred to as fooling (Szegedy et al. 2014; Nguyen, Yosinski, and
Clune 2015; Kardan and Stanley 2018). However, we will limit our scope to PGMs
as our end goal is to analyze and/or hypothesize structural dependency
relationships. Given this goal, such over-generalization could yield non-
optimal structures, biasing analysis and derived hypotheses leading to
misguided conclusions. To illustrate this over-generalization and provide
intuition for learning at the RV instantiation level, we provide a motivating
example taken from real-world data.

##### Motivating Example

It is well known in cancer research that the genes TP53 and TTN have somatic
mutations that affect chemotherapy responses (Xue et al. 2021). To demonstrate
a real-world effect of BN over-generalization, we learned a simple BN for this
interaction over the TCGA (Tomczak, Czerwińska, and Wiznerowicz 2015)
mutational dataset as seen in Figure 1(a). This BN encodes four possible
worlds represented by distinctly styled arrows in Figure 1(b). For this
example we have reduced the state space of each gene to just mutated or not
mutated. Assume our goal is to minimize the entropy or uncertainty of each
world or explanation. Then the conditional entropy of the model is the sum
over each world’s conditional entropy which is inherently direction dependent.
Since there exists many possible world edge configurations (RV instantiation
dependencies), there may exist a better set of edges than those induced by the
BN. Figure 1(c) shows this is true and illustrates the best collection of
minimal entropy inferences for this example.

(a)

(b)

(c)

Figure 1: (a) A simple learned BN over the TCGA gene mutation dataset using
GOBNILP (Cussens, Haws, and Studeny 2015) where the variable states are either
mutated or not mutated. (b) A graph of all corresponding worlds represented in
(a) delineated by different line styles. (c) A better orientation of intra-
world dependency relationships that lead to a lower total conditional entropy.
All values are conditional entropies calculated from the TCGA gene mutational
dataset.

##### Contributions

To address the over-generalization described we develop a structure learning
algorithm leveraging the Bayesian Knowledge Base (BKB) framework as it
inherently operates on the RV instantiation level. We accomplish this by
detailing a necessary scoring metric to rank BKB models based on an MDL
analysis and show theoretically that our MDL score takes over-generalization
into account. Leveraging this theoretical result, we then develop our BKB
Structure Learning (BKBSL) algorithm to minimize MDL and demonstrate
empirically both competitive accuracy and better data fit compared to learned
BNs. Further, we show that our algorithm can utilize existing optimization
frameworks for DAG learning bringing BKBSL into the realm of these well
studied off-the-shelf methods. Lastly, we conclude by utilizing a learned BKB
to explain possible gene associations over TCGA breast cancer data.

##  2 Related Work and Preliminaries

As the MDL principle will be our guiding force in both theoretical and
empirical analysis, we provide a brief review of its applications to directed
PGMs, e.g. Bayesian Networks, as these models are most applicable to our study
of BKBs. Lam and Bacchus (1994) first presented an MDL learning approach for
BNs based on a heuristic search method seeking to spend equal time between
simple and more complex BN structures. This was accomplished by extending Chow
and Liu’s (1968) result on recovering polytrees to general BNs via Kullback-
Leibler cross entropy minimization allowing them to develop a weighting
function over nodes. Their approach demonstrated that minimizing the MDL of
BNs performs an intuitive trade-off between model accuracy and complexity. In
their work, they also eluded to a potential subjectivity in choosing a model
encoding strategy leading to research for improved MDL scores for BNs (Yun and
Keong 2004; Drugan and Wiering 2010). Hansen and Yu (2001) detail a complete
review of various MDL formulations.

Empirical evaluation of MDL as a scoring function for BN learning has also
been well studied. Yang and Chang (2002) analyzed the performance of five
scoring functions: uniform prior score metric (UPSM), conditional uniform
prior score metric (CUPSM), Dirichlet prior score metric (DPSM), likelihood-
equivalence Bayesian Dirichlet score metric (BDeu), and MDL. They showed that
MDL was able to correctly identify ground-truth network structures from a
variety of possible candidates, yielding the highest discrimination ability.
Liu et al (2012) also performed empirical BN learning analysis over different
scoring function, namely: MDL, Akaike’s information criterion (AIC), BDeu, and
factorized normalized maximum likelihood (fNML). Their approach tested the
recovery accuracy of each scoring method over various gold standard networks
as compared to the random networks used by Yang and Chang. Their results
confirm the utility of MDL as it performed best in recovering the optimal
networks when sufficient data was given.

To our knowledge there has been no work in structure learning on the RV
instantiation level, likely due to the desire to learn complete distributions.
Further, we have limited our comparisons to BNs as they are a predominant
model in literature and provide a comparison to judge empirical results.

###  2.1 Bayesian Networks

A BN is a DAG $G=(V,E)$ that represents the factorized joint probability
distribution over random variables $X=(X_{1},\dots,X_{n})$ of the form:

$\text{Pr}(X)=\prod_{i}^{n}P(X_{i}\mid\Pi(X_{i}))$ (1)

such that $\Pi(X_{i})$, or more concisely $\pi_{i}$, are the structural
parents of the random variable $X_{i}$ according to $G$ and each node
$V_{i}\in V$ correspond directly to a random variable $X_{i}$ and $n$ is the
number of random variables in $X$. As the BN MDL formulation is well known in
the literature, we point the reader to Appendix D.1 for a review.

###  2.2 Bayesian Knowledge Bases

Santos, Jr. and Santos (1999) developed the BKB framework to generalize BNs to
the random variable instantiation level and to offer knowledge engineers a
semantically sound and intuitive knowledge representation. BKBs unify “if-
then” rules with probability theory to provide several modeling benefits
compared to BNs, fuzzy logics, etc. First, BKBs do not require complete
accessibility and can even leverage incomplete information to develop
potentially more representative models. Second, since BKBs operate at the
instantiation level, they can handle various types of cyclic knowledge.
Lastly, BKBs have both robust tuning algorithms to assist in model
correction/validation (Santos, Gu, and Santos 2013; Yakaboski and Santos 2018)
and information fusion algorithms to incorporate knowledge efficiently and
soundly from disparate knowledge sources (Santos, Wilkinson, and Santos 2011;
Yakaboski and Santos 2021).

BKBs consist of two components: instantiation nodes (I-nodes) which represent
instantiations of random variables of the form $X_{i}=x_{ik}$ where $k$ is the
$k$-th state of $X_{i}$, and support nodes (S-nodes) that represent the
conditional probabilities between I-node relationships of the form
$X_{i}=x_{ik}\rightarrow q=0.87\rightarrow X_{j}=x_{jl}$. The collection of
these (in)dependencies describe the BKB correlation graph. For a precise
definition and a graphical depiction see Appendix D.

While BKBs can handle incompleteness and various forms of cyclicity, it is
necessary that all S-nodes in a BKB obey mutual exclusivity (mutex) and
associated probability semantics. Mutex guarantees that mutually exclusive
events cannot be true at the same time. Concretely, we say that two sets of
I-nodes, $I_{1}$ and $I_{2}$, are mutex if there exists an I-node
$X_{i}=x_{ik_{1}}\in I_{1}$ and $X_{i}=x_{ik_{2}}\in I_{2}$ such that
$k_{1}\neq k_{2}$. We say that two S-nodes are mutex if the two I-node parent
sets of each S-node are mutex.

BKBs use a weight function to map S-node weights to associated conditional
probabilities associated with the conditional dependence relationship
described by the S-node. This weight function is analogous to the conditional
probability tables (CPTs) of BNs except that BKBs do not require complete
information. This weight function along with the correlation graph defines the
BKB.

##### Bayesian Knowledge Base

A Bayesian Knowledge Base (BKB) is a tuple $K=(G,w)$ where $G$ is a
correlation graph $G=(I\cup S,E)$ where $I$ and $S$ are the sets of I- and
S-nodes, $E\subset\\{I\times S\\}\cup\\{S\times I\\}$, and
$w:S\rightarrow[0,1]$ is a weight function, such that the following properties
are met:

  1. 1.

$\forall q\in S$, the set of all incident I-nodes to $q$, denoted $\Pi(q)$,
contains at most one instantiation of each random variable.

  2. 2.

For distinct S-nodes $q_{1},q_{2}\in S$ that support the same I-node, denoted
$Head_{G}(q_{i})$, the sets $\Pi(q_{1})$ and $\Pi(q_{2})$ must be mutex.

  3. 3.

For any $Q\subseteq S$ such that (i) $Head_{G}(q_{1})$ and $Head_{G}(q_{2})$
are mutex, and (ii) $\Pi(q_{1})$ and $\Pi(q_{2})$ are not mutex, then $\forall
q_{1},q_{2}\in Q,\sum_{q\in Q}w(q)\leq 1$.

The probability of an inference (world) in a BKB, denoted $\tau$, is the
product of all S-node weights $w(q)$ consistent with that world (Appendix
D.2). The joint probability distribution of a BKB is the sum of all inference
probabilities consistent with an evidence set $\Theta$ represented as
$I_{\Theta}$ and given by

$P(\Theta)=\sum_{\tau\in I_{\Theta}}\prod_{q\in\tau}w(q)$ (2)

##  3 The BKB Minimum Description Length

To construct a BKB structure learning algorithm, we first need to define a
scoring function that can rank BKB structures as well as provide a theoretical
justification for its utility. For these reasons we focus our attention on a
Minimum Description Length (MDL) score as it is a well studied tenet of
learning theory (Rissanen 1978). The idea is that the best model given data
should minimize both (1) the encoding length of the model and (2) the length
of encoding the data given the model. This is akin to applying Occam’s Razor
to model selection, i.e., choose the simplest model that still describes the
data reasonably well.

###  3.1 Encoding the BKB

The minimum length needed to encode a BKB is related directly to the number of
S-nodes modeled by the BKB. The encoding of each S-node will contain a
probability value and a parent set of I-nodes. For a problem with $n$ RVs such
that each RV can have $r_{i}$ number of instantiations, to encode all I-nodes
we would need $\log_{2}(m)$ number of bits where $m=\prod_{i}r_{i}$. The
general BKB MDL is

$\sum_{q\in
S}\bigg{(}(|\Pi(q)|+1)\log_{2}(m)+\delta\bigg{)}-N\sum_{\tau}p_{\tau}\log_{2}(q_{\tau})$
(3)

where $\delta$ is the number of bits needed to store the probability value and
the first term is the BKB model encoding length. From Equation 3 it is clear
that we incur a modeling cost for the BKB based on the finer granularity of
the model. This cost is derived in Appendix C.1. Therefore, if we know that
the distribution factorizes to a BN, there is no reason to use a BKB. However,
rarely in practice do we know the true distribution of the data and in most
cases the data we learn from is incomplete. This eludes to a natural rationale
for using a BKB that will be theoretically justified.

###  3.2 Encoding the data with a BKB

The task of encoding dataset $D$ is centered on learning a joint distribution
over random variables $X=(X_{1},\dots,X_{n})$. As we are focused on the
discrete case, each variable $X_{i}$ will have $r_{i}$ discrete states and
every unique choice of random variable instantiations defines a complete world
$\tau$ of the ground-truth distribution of the data and is assigned an
associated probability value, $p_{\tau}$.

We will make several standard statistical assumptions. We assume that each
data instance in $D$ is a complete world such that each data instance
specifies a value for every random variable in $X$111Only for simplicity do we
make the assumption that a data instance must be complete, i.e., not contain
any missing values.. We assume that each data instance is a result of
independent random trails and expect that each world would appear in the
dataset with a frequency $\approx Np_{\tau}$.

The main theme of the MDL metric is to efficiently compress $N$ data instances
as a binary string such that we minimize the string’s encoding length. There
are many different coding and compression algorithms we can use, but since we
simply care about comparing different MDLs as a metric we can limit our focus
to just symbol codes (MacKay 2005). Symbols codes assign a unique binary
string to each world in the dataset. For instance, the world
$\\{X_{1}=1,X_{2}=1,X_{3}=0\\}$ might be assigned the symbol 0001. We can then
encode the dataset as just the concatenation of all these world symbols and
judge our compression effectiveness by calculating the length of this final
binary string.

Research in information theory has proved that for symbol codes it is possible
to minimize the length of this encoded data string by leveraging the
probability/frequency of each world symbol occurring. Specifically, Huffman’s
algorithm generates optimal Huffman codes (MacKay 2005), i.e., optimal symbol
code mappings, that yield a minimum length. The key intuition for Huffman
codes is that we should give shorter code words to more frequently occurring
worlds and longer ones to less probable worlds. Lam and Bacchus (1994) proved
that the encoding length of the data is a monotonically increasing function of
the cross-entropy between the distribution defined by the model and the true
distribution.

Therefore, if we have a true distribution $P$ and a model distribution $Q$
over the same set of worlds, $\tau_{1},\dots,\tau_{t}$, where each world
$\tau_{i}$ is assigned the probability $p_{i}$ by $P$ and $q_{i}$ by $Q$, then
the cross entropy between these distributions is

$C(P,Q)=\sum_{i=1}^{t}p_{i}\log_{2}\frac{p_{i}}{q_{i}}=\sum_{i=1}^{t}p_{i}(\log_{2}p_{i}-\log_{2}q_{i})$
(4)

Calculating Equation 4 is not appealing as the number of worlds is
combinatorial in the number of variables. Chow and Lui (1968) developed a
famous simplification of Equation 4 as just a local computation over low-order
marginals when $Q$ has a tree factorization. Lam and Bucchus (1994) extended
their result to models that have a general DAG structure. Their main result
concludes that $C(P,Q)$ is a monotonically decreasing function of the sum of
each random variable’s mutual information with their parents,
$I(X_{i};\Pi(X_{i}))$. Their exact result is restated in Appendix Theorem 3.

Further generalizing these results to the instantiation level, we can now show
that an optimally learned BKB can encode the distribution as well or better
than a BN. Consider again the fundamental MDL problem of learning the dataset
encoding. Equation 4 says that we need to minimize the total cross-entropy
between $P$ and $Q$. However, in terms of data encoding, we only need to
minimize the difference between each $p_{i}$ and $q_{i}$ for unique worlds
that are actually in $D$. In this sense, our database encoding doesn’t care
about the worlds that aren’t in the dataset, but for which a model like a BN
naturally defines and generalizes. Therefore, BKBs handling of incompleteness
gives us an opportunity to more tightly fit the data we actually know.
Consider the following cross-entropy

$\displaystyle C(P,Q)=\sum_{i=1}^{d}$ $\displaystyle
p_{i}(\log_{2}p_{i}-\log_{2}q_{i})$
$\displaystyle+\sum_{i=d+1}^{t}p_{i}(\log_{2}p_{i}-\log_{2}q_{i})$ (5)

where worlds $\tau_{1},\dots,\tau_{d}$ are represented by the unique exemplars
in $D$ that we hope to encode, i.e.,
$\\{\tau_{i},\dots,\tau_{d}\\}=\\{d_{1},d_{2},d_{3},\dots\\}_{\neq}\subseteq
D$, and worlds $\tau_{d+1},\dots,\tau_{t}$ are worlds that our model induces.
In terms of encoding length we can narrow our focus to only considering worlds
present in $D$. As Lam and Bacchus (1994) proved the encoding length of the
data is a monotonically increasing function cross-entropy, it is trivial to
prove the following corollary.

######  Corollary 1.

Let’s Define the cross-entropy
$C_{D}(P,Q)=\sum_{i=1}^{d}p_{i}(\log_{2}p_{i}-\log_{2}q_{i})$ between two
distributions $P$ and $Q$ over the same set of worlds
$\tau_{1},\dots,\tau_{d}$ s.t. these worlds must be included in a dataset $D$.
Then the encoding length of the data $D$ is monotonically increasing function
of $C_{D}$.

Combining Corollary 1 with Lam and Bacchus mutual information theorem
(Appendix D Theorem 3) we arrive at our main theoretical result.

######  Theorem 1.

$C_{D}(P,Q)$ is a monotonically decreasing function of

$\sum_{\tau\in
D_{\neq}}\sum_{i=1}^{n}p(x_{i\tau},\pi_{i\tau})\log_{2}\frac{p(x_{i\tau},\pi_{i\tau})}{p(x_{i\tau})p(\pi_{i\tau})}$
(6)

where $x_{i\tau}$ is the instantiation of random variable $X_{i}$ determined
by data instance $\tau$, $\pi_{i\tau}$ is the parent set instantiation of
random variable $X_{i}$ governed by $\tau$, and $D_{\neq}$ is the set of
unique data instances (worlds) represented in $D$. Therefore, $C_{D}(P,Q)$
will be minimized when Equation 6 is maximized.

We leave the proof of this theorem to Appendix B.1 as the intuition is fairly
straightforward from Lam and Bacchus’ theorem (Appendix D Theorem 3). We have
established that the encoding length of the data is a increasing function
solely of $C_{D}$ and that maximizing Equation 6 minimizes $C_{D}$ and thereby
the encoding length. With these results, we can deduce the existence of a
theoretical BKB that will have an equal to or better data encoding length than
the induced worlds of a BN given $D$.

######  Theorem 2.

Given a BN $G$ learned over a dataset $D$ as to maximize the weight
$W_{G}=\sum_{i}I(X_{i};\Pi(X_{i}))$ given the structure $G$, there exists a
BKB $K$ with a weight $W_{K}$ according to Equation 6 such that $W_{K}\geq
W_{G}$.

We defer a detailed proof of Theorem 2 to Appendix B.2 and provide a more
concise proof sketch. The key insight is that any BN $G$ can be transformed
into a corresponding BKB $K_{G}$ as BKBs subsume BNs. We are only interested
in the data encoding length which can now be calculated over $K_{G}$ by
summing over the complete worlds represented in $G$. Consider a single random
variable $X_{i}$ and it’s associated parent set $\Pi(X_{i})$. For each
instantiation of $X_{i}$ and $\Pi(X_{i})$ there will be an associated S-node
created in $K_{G}$ with an instantiated weight according to Equation 6. Since
the choice of parent set for each random variable instantiation in
$\Pi(X_{i})\cup\\{X_{i}\\}$ is governed by $G$, we don’t consider other S-node
configurations of the same instantiated random variable set that may have
greater weight. The BN structure constrains our S-node structures. Therefore,
if we start with a optimal BN, transform it into a BKB $K_{G}$, and analyze
every permutation of each S-node’s possible configurations taking the
permutation that maximizes the instantiated weight, we will end up with a BKB
$K$ with the same number of S-nodes that has a total weight equal to or
greater than the BN representation $K_{G}$. This result allows us to also
state the following corollary based on the fact that
$I(X_{i};\Pi(X_{i}))=H(X_{i})-H(X_{i}|\Pi(X_{i}))$.

######  Corollary 2.

Since $I(X_{i};Pi(X_{i}))=H(X_{i})-H(X_{i}|\Pi(X_{i}))\geq 0$. We can maximize
Equation 6 by minimizing the instantiated conditional entropy
$H(x_{i\tau}|\pi_{i\tau})=p(x_{i\tau},\pi_{i\tau})\log_{2}\frac{p(x_{i\tau},\pi_{i\tau})}{p(\pi_{i\tau})}$.

##  4 BKB Structure Learning

Theorem 1 dictates that for every random variable instantiation $x_{i\tau}$ in
a data instance (world) $\tau\in D_{\neq}$ where $D_{\neq}$ is the set of
unique data instances in $D$, we should assign an instantiated parent set
$\pi_{i\tau}$ such that the instantiated conditional entropy is minimized
according to Corollary 2. The key insight of our structure learning approach
is that we can decompose our learning over the worlds represented in the data.
In each world, we will have at most a single instantiation of each RV and our
goal is to select a set of S-nodes with a structure that minimizes
instantiated conditional entropy for that world. We can view each world in the
data as a separate complete inference which form an acyclic subgraph of their
respective encompassing BKB. A precise definition of a BKB inference can be
found in Appendix D.

Our structure learning algorithm reduces to finding a directed acyclic
inference graph for each world that minimizes
$\sum_{\tau}\sum_{i}H(x_{i\tau}|\pi_{i\tau})=p(x_{i\tau},\pi_{i\tau})\log_{2}\frac{p(x_{i\tau},\pi_{i\tau})}{p(\pi_{i\tau})}$.
Further, we can use any off-the-shelf DAG learning algorithm to accomplish
this step so far as our scoring function inherently minimizes instantiated
conditional entropy and BKB encoding length. There has been significant
advancements in field of BN and DAG learning and we make no attempt in
covering all such procedures. Instead we will focus on the state-of-the-art
exact BN (DAG) solver GOBNILP (Cussens 2012; Cussens, Haws, and Studeny 2015).

Upon learning each minimal entropy inference, we then need a method for
merging this knowledge together that is semantically sound. A standard union
type operation will not generally be supported as the unioned BKB would likely
incur many mutex violations as seen in Appendix Figure 4a. Instead, we can
employ a well-studied BKB fusion (Santos, Wilkinson, and Santos 2011;
Yakaboski and Santos 2021) algorithm that supports the fusion of an arbitrary
number of BKB Fragments (BKFs) by attaching source I-nodes to every S-node
corresponding to the data instance from which the inference graph originated.
A graphical example of this approach is depicted in Appendix Figure 4b along
with additional information regarding BKB fusion in Appendix D.3. This
procedure ensures that that no mutual exclusion violations are present in the
fused BKB maintaining a consistent probability distribution over the data and
leading to model generalization. Appendix D.3 provides a detailed explanation
of generalization in fused BKBs.

Aside from forming a mutual exclusive BKB, fusion also presents us with
another degree of freedom during learning. If each data instance was generated
by an i.i.d. process, it is natural to assume a normalized probability over
all source nodes. However, many processes do not generate truly i.i.d or
representative samples. Therefore, if we view these source S-nodes as
reliabilities that can be tuned, we may be able to correct errors in higher
order inference calculations that arise due to under-fitting or over-
generalization. We leave such analysis to future work. Combining each of the
steps presented so far, we outline our general BKB structure learning
procedure in Algorithm 1.

Algorithm 1 BKB Structure Learning

Input: Dataset $D$, Source Reliabilities $R$, DAG learning algorithm $f$ and
hyperparameters $\Theta$

1: $K\leftarrow\emptyset$

2: for $\tau\in D_{\neq}$ do

3: $G_{\tau}\leftarrow f(\tau,R,\Theta)$

4: $K\leftarrow K\cup\\{G_{\tau}\\}$

5: end for

6: return $\text{BKB-Fusion}(K,R)$

##  5 Empirical Results

To demonstrate the utility of both our proposed algorithm as well as our
learned BKB models we conducted 40 experiments on benchmark datasets comparing
BKBSL and BNSL in terms of MDL and complexity performance. We then conducted
22 cross validation classification experiments to compare accuracy performance
with learned BNs as well as a use-case studying the under-determined
bioinformatics domain of structural dependency analysis among single-
nucleotide polymorphism (SNPs) in breast cancer.

###  5.1 Benchmark Analysis

When comparing MDL between our learned BKBs and BNs, we are only concerned
with comparing the data encoding length, as the model encoding length is only
used to penalize more complex models. Our data MDL results in Appendix Table 1
demonstrate that a BKB learned using our BKBSL algorithm finds a tighter data
fit than the best BN in all 40 dataset. Intuitively, this is because the BN
must generalize away from potentially good instantiated scores in favor of the
entire random variable score.

Our MDL experiments also demonstrates a practical strength of BKBSL over BNSL
related to the number of calls to a joint probability calculator or estimator.
In order to calculate the necessary scores for an exact DAG learning algorithm
like GOBNILP, we needed to calculate empirical joint probabilities from each
dataset. For all experiments we tracked the number of unique calls to this
function by our BKBSL algorithm and traditional BNSL algorithm. Since BNSL
operates at the RV level, it had to calculate all joint probabilities governed
by a given parent set configuration. However, BKBSL did not need to calculate
the full CPTs as it operates at the RV instantiation level and decomposes over
each data instance, reducing the number of calls to this calculator. We feel
that this is a more representative complexity performance metric as it is
agnostic to equipment configurations. This effect is detailed in Appendix
Table 1, and we can see strong correlation between performance savings over
BNs and the number of features (Pearson $r=-0.5994$, $p$-value
$=$4.363\text{\times}{10}^{-5}$$) as well as the number of I-nodes
($r=-0.4916$, $p$-value $=0.0013$) in the dataset. All learned BKBs and BNs
are hosted on Github and can be viewed within in a Jupyter Notebook for easier
visualization.

To finalize our BKBSL benchmark analysis we performed accuracy comparisons
between our learned BKBs and traditionally learned BNs using GOBNILP and
standard MDL scoring. We performed a 10-fold classification cross validation
on a subset of only 22 datasets due to the increased learning and reasoning
time incurred by running cross validation analysis. We can see from Appendix
Table 2 that our BKBSL algorithm is very competitive with BNs in terms of
precision, recall and F1-score. Further, our BKBSL models even beat BNs in
$63\%$ of cases in terms of precision with greater degradation of performance
in terms of recall and F1-score. The alternative hypothesis that either
traditionally learned BNs or our learned BKBs will outperform each other in
all accuracy categories (Chi2 Statistic $\chi=1.0$, $p$-value $=0.3173$) is
_not_ statistically significant. Therefore, we fail to reject the null that
learned BNs or BKBs perform better in these cases owing to approximately equal
total performance.

This raises the question: Why does our learned BKB perform better in some
datasets and not in others? While no feature of the datasets provided any
statistically significant predictor of superior performance and leaving more
in-depth analysis to future work, we do hypothesize an explanation. It is a
well-known problem that real world datasets are often unfaithful to DAGs, e.g.
BNs, due to the existence of multiple information equivalent Markov Boundaries
(MBs) (Statnikov, Lemeir, and Aliferis 2013; Wang and Wang 2020). Since our
BKBSL focuses on learning an optimal structure for every unique exemplar
$\tau$, we can view each learned BKF as an equivalent inference from a
hypothetical BN whose dependency structure matches that of the associated BKF.
As we are only concerned with the specific instantiations of $\tau$, the
hypothetical BN and BKF will yield the same probability for this world as
their parameters are governed by the same dataset. As our prediction task is
to determine the most likely state of a response variable $Y$ given a complete
set of evidence $E$, e.g., $y^{*}=\operatorname*{argmax}_{y}Q(Y=y\mid E)$,
then the closer our joint probability $Q(Y=y,E)$ is to the true data
distribution $P(Y=y,E)$ the more accurate our classifier. This is due to the
fact when comparing all conditional probabilities $Q(Y=y_{i}\mid
E)=\frac{Q(Y=y_{i},E)}{Q(E)}$ the denominator cancels out and we are only
concerned with the accuracy of $Q(Y=y_{i},E)$. If we imagine our learned BKFs
deriving from various hypothetical BNs each with uniquely induced MBs for
every RV, our fused BKB essentially incorporates a multiplicity of MBs choices
for each of these hypothetical BNs and selects the best performing choice for
every world of interest, i.e., prediction class given evidence. We hypothesize
that our BKBSL will then perform better on datasets that induce more
information equivalent MBs since a BN must select only one and our BKBSL can
incorporate multiple in its predictions. Whereas in datasets with fewer MBs,
our BKBSL performance may degrade due to overfitting. We intend to study this
area further as it may yield clear indications about when to use BKBSL over
BNSL in practice.

###  5.2 Gene Regulatory Network (GRN) Application in Breast Cancer

We applied our approach to somatic mutation profiles of breast cancer cases in
TCGA (Tomczak, Czerwińska, and Wiznerowicz 2015) to study whether our learned
model could still offer utility in this extremely under-determined domain.
Since prediction accuracy would not be a reliable metric of success in this
dataset, we focused our analysis on hypothesizing potentially significant
mutational interactions in cancer. However, if we are to trust any structural
hypotheses generated by our approach, we need to ensure the model captures two
fundamental biological concepts: (1) We can extract two- or three-hit
interactions that are supported in the literature (Knudson 1971, 2001;
Segditsas et al. 2009), and (2) we can identify and (possibly) handle genomic
instability (Bai et al. 2014; Croessmann et al. 2017).

(a)

(b)

Figure 2: (a) RV level graph of learned BKB over TCGA breast cancer subgraphed
on TP53 SNP relationships with HER genes. (b) Genomic instability evidence
from the RV level cycle between PIK3CA SNP and general TP53 gene variable. To
capture multiple levels of granularity, we included features related to exact
positional SNPs as well as including gene level features and general variant
classifications. For details about our naming conventions and feature
selection process see Appendix A.

Given the well-regarded two- and three-hit hypotheses for understanding the
role of genetic mutations in cancer development, a model attempting to
describe a mutational GRN should be able to capture this concept. The premise
behind the two- or three-hit hypotheses is that because many cancers are
driven by mutations in various tumor suppressor genes and these loss-of-
function mutations are recessive in nature, in general, at least two mutations
in these genes are needed to develop cancer. There are certain tumor
suppressors genes such as TP53 (Kandoth et al. 2013; Muller and Vousden 2014)
that are common in all cancer sub-types and are likely the first hit for non-
hereditary cancers. Looking at the dependence relationship subgraph in our
learned BKB between a first hit tumor suppressor gene such as TP53 and a
second hit tumor suppressor gene related to breast cancer such as a HER1 or
HER2 (Osborne, Wilson, and Tripathy 2004), we should observe a directed
dependency relationship from TP53 to HER1 or HER2. Figure 2(a) shows we
observe this relationship adding to the biological validity of our model.

It is well known that cancer is also driven by genomic instability (Negrini,
Gorgoulis, and Halazonetis 2010), the increased tendency for gene mutations to
accumulate, disrupting various biological pathways and in turn causing more
mutations. An artifact that we discovered in our BKB network analysis and
illustrative of our first example is that in some cases there exist cyclic
relationships on the random variable level. This is due to separate inferences
learning opposing dependency relationships due to the additional context of
their respective inference. This effect can be seen in Figure 2(b). Here we
have a cyclic random variable relationship between TP53 and the SNP located at
chromosome 17 position 179234297 of the PIK3CA gene. TP53 and PIK3CA are known
drivers of genomic instability (Bai et al. 2014), affecting the mutational
states of many other genes in their inference contexts. Since our inference
level parent set limit was set to one, our algorithm cannot reliably extract
mutual dependency relationship between TP53 and PIK3CA, thereby causing
different inferences to have different directionalities. This cannot be
captured by a BN. Further, this result is supported by in-vitro research
regarding the joint effect of PIK3CA and TP53 on genomic instability in breast
cancer (Croessmann et al. 2017) and is expected from this model given that
both genes drive many other downstream mutations. Overall we found 12 of these
cyclic relationships in our learned BKB. Of these 12 associations we found
literature support for four of them, namely, PIK3CA and TP53 (Croessmann et
al. 2017), OBSCN and TP53 (Sjöblom et al. 2006; Perry et al. 2013), MAP3K1 and
PIK3CA (Avivar-Valderas et al. 2018), CAD and PIK3CA (Wu et al. 2014).

##  6 Limitations and Future Work

The primary limitation of our approach is that we need to learn a BKB fragment
for every instance in the training dataset. While we have some complexity
benefits from making fewer calls to a given joint probability calculator, this
benefit could be neutralized due to the requirement of having to run the DAG
learning algorithm multiple times. For instance, once all scores are
calculated, BNSL only requires a single run of the respective DAG learning
algorithm, whereas our BKBSL approach requires running this algorithm $N$
times. We leave to future work the exploration of this trade off while
hypothesizing that there may be a DAG learning formulation in which all
fragments can be learned in a single pass by reusing/sharing local
scores/structures between fragments.

While our BKBSL algorithm seems to generalize well to unobserved examples in
our benchmark experiments, we still see instances with significant accuracy
degradation. It is a largely unanswered question as to why machine learning
algorithms perform better or worse on particular datasets (Pham and
Triantaphyllou 2008) and we have detailed a possible explanation in Section
5.1 to be explored in future work. As mentioned in Section 4, we could also
address accuracy degradation by tuning the source node reliabilities in our
model. Such an approach yields another degree of freedom to adjust the model
and also may highlight the importance/significance of individual data
instances in relation to over model accuracy. We also leave this direction for
future research.

##  7 Conclusions

We have presented a new approach for performing Bayesian structure learning at
the random variable instantiation level by leveraging Bayesian Knowledge
Bases. We have detailed a theoretical justification for our algorithm and
learned model as being the fusion of minimal entropy inferences or
explanations over each training exemplar. We demonstrated empirically that our
BKBSL algorithm finds a superior BKB than an equivalent BN (scored as BKB) on
40 benchmark datasets using our MDL formulation. Further, we demonstrated the
practical utility of our approach by presenting statistically competitive
accuracy with learned BNs over 22 benchmark datasets using a 10-fold cross
validation. This provides evidence that our algorithm adequately generalizes
to unseen data based on known knowledge rather than over-generalizing to force
a complete distribution. Lastly, we conducted a structural analysis over a
gene regulatory network learned from breast cancer mutation data taken from
TCGA. This analysis resulted in finding dependency relationships that matched
biological intuition and also revealed associations that are well known in the
bioinformatics community.

## Acknowledgements

This research was funded in part by the National Center for Advancing
Translational Sciences (NCATS), National Institutes of Health, through the
Biomedical Data Translator program (NIH award OT2TR003436), Air Force Office
of Scientific Research Grant No. FA9550-20-1-0032 and DURIP Grant No.
N00014-15-1-2514. Special thanks to Joseph Gormley and Eugene Hinderer for
their comments and encouragement.

## References

  * Akaike (1974) Akaike, H. 1974.  A new look at the statistical model identification.  _IEEE Transactions on Automatic Control_ , 19(6): 716–723.
  * Alcalá-Fdez et al. (2011) Alcalá-Fdez, J.; Fernández, A.; Luengo, J.; Derrac, J.; García, S.; Sánchez, L.; and Herrera, F. 2011.  Keel data-mining software tool: data set repository, integration of algorithms and experimental analysis framework.  _Journal of Multiple-Valued Logic & Soft Computing_, 17.
  * Avivar-Valderas et al. (2018) Avivar-Valderas, A.; McEwen, R.; Taheri-Ghahfarokhi, A.; Carnevalli, L. S.; Hardaker, E. L.; Maresca, M.; Hudson, K.; Harrington, E. A.; and Cruzalegui, F. 2018.  Functional significance of co-occurring mutations in PIK3CA and MAP3K1 in breast cancer.  _Oncotarget_ , 9(30): 21444–21458.
  * Bai et al. (2014) Bai, X.; Zhang, E.; Ye, H.; Nandakumar, V.; Wang, Z.; Chen, L.; Tang, C.; Li, J.; Li, H.; Zhang, W.; Han, W.; Lou, F.; Zhang, D.; Sun, H.; Dong, H.; Zhang, G.; Liu, Z.; Dong, Z.; Guo, B.; Yan, H.; Yan, C.; Wang, L.; Su, Z.; Li, Y.; Jones, L.; Huang, X. F.; Chen, S.-Y.; and Gao, J. 2014.  PIK3CA and TP53 Gene Mutations in Human Breast Cancer Tumors Frequently Detected by Ion Torrent DNA Sequencing.  _PLOS ONE_ , 9(6): e99306.
  * Bouckaert (1995) Bouckaert, R. R. 1995.  _Bayesian belief networks: from construction to inference_.  PhD Thesis, University of Utrecht.
  * Burkart and Huber (2021) Burkart, N.; and Huber, M. F. 2021.  A Survey on the Explainability of Supervised Machine Learning.  _Journal of Artificial Intelligence Research_ , 70: 245–317.
  * Chickering, Heckerman, and Meek (2004) Chickering, D. M.; Heckerman, D.; and Meek, C. 2004.  Large-Sample Learning of Bayesian Networks is NP-Hard.  _Journal of Machine Learning Research_ , 5: 44.
  * Chow and Liu (1968) Chow, C.; and Liu, C. 1968.  Approximating discrete probability distributions with dependence trees.  _IEEE Transactions on Information Theory_ , 14(3): 462–467.
  * Croessmann et al. (2017) Croessmann, S.; Wong, H. Y.; Zabransky, D. J.; Chu, D.; Rosen, D. M.; Cidado, J.; Cochran, R. L.; Dalton, W. B.; Erlanger, B.; Cravero, K.; Button, B.; Kyker-Snowman, K.; Hurley, P. J.; Lauring, J.; and Park, B. H. 2017.  PIK3CA mutations and TP53 alterations cooperate to increase cancerous phenotypes and tumor heterogeneity.  _Breast Cancer Research and Treatment_ , 162(3): 451–464.
  * Cussens (2012) Cussens, J. 2012.  Bayesian network learning with cutting planes.  _arXiv:1202.3713_.
  * Cussens, Haws, and Studeny (2015) Cussens, J.; Haws, D.; and Studeny, M. 2015.  Polyhedral aspects of score equivalence in Bayesian network structure learning.  _arXiv:1503.00829_.
  * Dasgupta (1999) Dasgupta, S. 1999.  Learning Polytrees.  In _Fifteenth Conference on Uncertainty in Artificial Intelligence_ , 134–141.  ArXiv: 1301.6688.
  * Došilović, Brčić, and Hlupić (2018) Došilović, F. K.; Brčić, M.; and Hlupić, N. 2018.  Explainable artificial intelligence: A survey.  In _2018 41st International Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO)_ , 0210–0215.
  * Drugan and Wiering (2010) Drugan, M. M.; and Wiering, M. A. 2010.  Feature selection for Bayesian network classifiers using the MDL-FS score.  _International journal of approximate reasoning_ , 51(6): 695–717.
  * Fayyad and Irani (1993) Fayyad, U.; and Irani, K. 1993.  Multi-interval discretization of continuous-valued attributes for classification learning.
  * Grüttemeier and Komusiewicz (2020) Grüttemeier, N.; and Komusiewicz, C. 2020.  Learning Bayesian Networks Under Sparsity Constraints: A Parameterized Complexity Analysis.  In _Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence_ , 4245–4251. Yokohama, Japan.  ISBN 978-0-9992411-6-5.
  * Hansen and Yu (2001) Hansen, M. H.; and Yu, B. 2001.  Model selection and the principle of minimum description length.  _Journal of the American Statistical Association_ , 96(454): 746–774.
  * Heckerman et al. (2000) Heckerman, D.; Chickering, D. M.; Meek, C.; Rounthwaite, R.; and Kadie, C. 2000\.  Dependency networks for inference, collaborative filtering, and data visualization.  _Journal of Machine Learning Research_ , 1(Oct): 49–75.
  * Jr and Santos (1999) Jr, E. S.; and Santos, E. S. 1999.  A framework for building knowledge-bases under uncertainty.  _J. Exp. Theor. Artif. Intell._ , 11(2): 265–286.
  * Kandoth et al. (2013) Kandoth, C.; McLellan, M. D.; Vandin, F.; Ye, K.; Niu, B.; Lu, C.; Xie, M.; Zhang, Q.; McMichael, J. F.; Wyczalkowski, M. A.; Leiserson, M. D. M.; Miller, C. A.; Welch, J. S.; Walter, M. J.; Wendl, M. C.; Ley, T. J.; Wilson, R. K.; Raphael, B. J.; and Ding, L. 2013.  Mutational landscape and significance across 12 major cancer types.  _Nature_ , 502(7471): 333–339.
  * Kardan and Stanley (2018) Kardan, N.; and Stanley, K. O. 2018.  Fitted Learning: Models with Awareness of their Limits.  _arXiv:1609.02226_.
  * Knudson (1971) Knudson, A. G. 1971.  Mutation and cancer: statistical study of retinoblastoma.  _Proceedings of the National Academy of Sciences_ , 68(4): 820–823.
  * Knudson (2001) Knudson, A. G. 2001.  Two genetic hits (more or less) to cancer.  _Nature Reviews Cancer_ , 1(2): 157–162.
  * Koller and Friedman (2009) Koller, D.; and Friedman, N. 2009.  _Probabilistic graphical models: principles and techniques_.  MIT press.
  * Lam and Bacchus (1994) Lam, W.; and Bacchus, F. 1994.  Learning Bayesian belief networks: An approach based on the MDL principle.  _Computational Intelligence_ , 10(3): 269–293.
  * Liu, Malone, and Yuan (2012) Liu, Z.; Malone, B.; and Yuan, C. 2012.  Empirical evaluation of scoring functions for Bayesian network model selection.  _BMC Bioinformatics_ , 13(S15): S14.
  * MacKay (2005) MacKay, D. J. C. 2005.  _Information Theory, Inference, and Learning Algorithms_.  Cambridge University Press.
  * Malagrino, Roman, and Monteiro (2018) Malagrino, L. S.; Roman, N. T.; and Monteiro, A. M. 2018.  Forecasting stock market index daily direction: A Bayesian Network approach.  _Expert Systems with Applications_ , 105: 11–22.
  * Moritz et al. (2018) Moritz, P.; Nishihara, R.; Wang, S.; Tumanov, A.; Liaw, R.; Liang, E.; Elibol, M.; Yang, Z.; Paul, W.; Jordan, M. I.; and others. 2018.  Ray: A distributed framework for emerging ${$AI$}$ applications.  In _13th USENIX Symposium on Operating Systems Design and Implementation (OSDI 18)_ , 561–577.
  * Muller and Vousden (2014) Muller, P.; and Vousden, K. 2014.  Mutant p53 in Cancer: New Functions and Therapeutic Opportunities.  _Cancer Cell_ , 25(3): 304–317.
  * Negrini, Gorgoulis, and Halazonetis (2010) Negrini, S.; Gorgoulis, V. G.; and Halazonetis, T. D. 2010.  Genomic instability — an evolving hallmark of cancer.  _Nature Reviews Molecular Cell Biology_ , 11(3): 220–228.
  * Nguyen, Yosinski, and Clune (2015) Nguyen, A.; Yosinski, J.; and Clune, J. 2015.  Deep Neural Networks are Easily Fooled: High Confidence Predictions for Unrecognizable Images.
  * Osborne, Wilson, and Tripathy (2004) Osborne, C.; Wilson, P.; and Tripathy, D. 2004.  Oncogenes and Tumor Suppressor Genes in Breast Cancer: Potential Diagnostic and Therapeutic Applications.  _The Oncologist_ , 9(4): 361–377.
  * Pearl (1986) Pearl, J. 1986.  Fusion, propagation, and structuring in belief networks.  _Artificial Intelligence_ , 29(3): 241–288.
  * Perry et al. (2013) Perry, N. A.; Ackermann, M. A.; Shriver, M.; Hu, L.-Y. R.; and Kontrogianni-Konstantopoulos, A. 2013.  Obscurins: Unassuming Giants Enter the Spotlight.  _IUBMB life_ , 65(6): 479–486.
  * Pham and Triantaphyllou (2008) Pham, H. N. A.; and Triantaphyllou, E. 2008.  The Impact of Overfitting and Overgeneralization on the Classification Accuracy in Data Mining.  In Maimon, O.; and Rokach, L., eds., _Soft Computing for Knowledge Discovery and Data Mining_ , 391–431. Boston, MA: Springer US.  ISBN 978-0-387-69935-6.
  * Rissanen (1978) Rissanen, J. 1978.  Modeling by shortest data description.  _Automatica_ , 14(5): 465–471.
  * Santos, Gu, and Santos (2013) Santos, E.; Gu, Q.; and Santos, E. E. 2013.  Bayesian knowledge base tuning.  _International Journal of Approximate Reasoning_ , 54(8): 1000–1012.
  * Santos, Wilkinson, and Santos (2011) Santos, E.; Wilkinson, J. T.; and Santos, E. E. 2011.  Fusing multiple Bayesian knowledge sources.  _International Journal of Approximate Reasoning_ , 52(7): 935–947.
  * Sauta et al. (2020) Sauta, E.; Demartini, A.; Vitali, F.; Riva, A.; and Bellazzi, R. 2020.  A Bayesian data fusion based approach for learning genome-wide transcriptional regulatory networks.  _BMC bioinformatics_ , 21: 1–28.  Publisher: Springer.
  * Scanagatta et al. (2015) Scanagatta, M.; de Campos, C. P.; Zaffalon, M.; and Corani, G. 2015.  Learning Bayesian Networks with Thousands of Variables.  _NIPS_.
  * Schwarz (1978) Schwarz, G. 1978.  Estimating the Dimension of a Model.  _The Annals of Statistics_ , 6(2): 461–464.
  * Segditsas et al. (2009) Segditsas, S.; Rowan, A. J.; Howarth, K.; Jones, A.; Leedham, S.; Wright, N. A.; Gorman, P.; Chambers, W.; Domingo, E.; Roylance, R. R.; Sawyer, E. J.; Sieber, O. M.; and Tomlinson, I. P. M. 2009.  APC and the three-hit hypothesis.  _Oncogene_ , 28(1): 146–155.
  * Shih, Choi, and Darwiche (2018) Shih, A.; Choi, A.; and Darwiche, A. 2018.  Formal Veriﬁcation of Bayesian Network Classiﬁers.  _Proceedings of Machine Learning Research_ , 72: 427–438.
  * Sjöblom et al. (2006) Sjöblom, T.; Jones, S.; Wood, L. D.; Parsons, D. W.; Lin, J.; Barber, T. D.; Mandelker, D.; Leary, R. J.; Ptak, J.; Silliman, N.; Szabo, S.; Buckhaults, P.; Farrell, C.; Meeh, P.; Markowitz, S. D.; Willis, J.; Dawson, D.; Willson, J. K. V.; Gazdar, A. F.; Hartigan, J.; Wu, L.; Liu, C.; Parmigiani, G.; Park, B. H.; Bachman, K. E.; Papadopoulos, N.; Vogelstein, B.; Kinzler, K. W.; and Velculescu, V. E. 2006.  The consensus coding sequences of human breast and colorectal cancers.  _Science (New York, N.Y.)_ , 314(5797): 268–274.
  * Statnikov, Lemeir, and Aliferis (2013) Statnikov, A.; Lemeir, J.; and Aliferis, C. F. 2013.  Algorithms for discovery of multiple Markov boundaries.  _The Journal of Machine Learning Research_ , 14(1): 499–566.
  * Suzuki (1993) Suzuki, J. 1993.  A Construction of Bayesian Networks from Databases Based on an MDL Principle.  In _Uncertainty in Artificial Intelligence_ , 266–273. Elsevier.  ISBN 978-1-4832-1451-1.
  * Suzuki (1999) Suzuki, J. 1999.  Learning Bayesian Belief Networks Based on the Minimum Description Length Principle: Basic Properties.  _IEICE Transactions on Fundamentals_ , E82.
  * Szegedy et al. (2014) Szegedy, C.; Zaremba, W.; Sutskever, I.; Bruna, J.; Erhan, D.; Goodfellow, I.; and Fergus, R. 2014.  Intriguing properties of neural networks.  _arXiv:1312.6199_.
  * Tomczak, Czerwińska, and Wiznerowicz (2015) Tomczak, K.; Czerwińska, P.; and Wiznerowicz, M. 2015.  The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge.  _Contemporary oncology_ , 19(1A): A68.
  * Trösser, de Givry, and Katsirelos (2021) Trösser, F.; de Givry, S.; and Katsirelos, G. 2021.  Improved Acyclicity Reasoning for Bayesian Network Structure Learning with Constraint Programming.  In _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence_ , 4250–4257. Montreal, Canada: International Joint Conferences on Artificial Intelligence Organization.  ISBN 978-0-9992411-9-6.
  * Tweedie et al. (2021) Tweedie, S.; Braschi, B.; Gray, K.; Jones, T. E.; Seal, R. L.; Yates, B.; and Bruford, E. A. 2021.  Genenames. org: the HGNC and VGNC resources in 2021.  _Nucleic acids research_ , 49(D1): D939–D946.
  * Wang and Wang (2020) Wang, Y.; and Wang, L. 2020.  Causal inference in degenerate systems: An impossibility result.  In _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , 3383–3392. PMLR.
  * Wu et al. (2014) Wu, X.; Renuse, S.; Sahasrabuddhe, N. A.; Zahari, M. S.; Chaerkady, R.; Kim, M.-S.; Nirujogi, R. S.; Mohseni, M.; Kumar, P.; Raju, R.; Zhong, J.; Yang, J.; Neiswinger, J.; Jeong, J.-S.; Newman, R.; Powers, M. A.; Somani, B. L.; Gabrielson, E.; Sukumar, S.; Stearns, V.; Qian, J.; Zhu, H.; Vogelstein, B.; Park, B. H.; and Pandey, A. 2014.  Activation of diverse signaling pathways by oncogenic PIK3CA mutations.  _Nature communications_ , 5: 4961.
  * Xue et al. (2021) Xue, D.; Lin, H.; Lin, L.; Wei, Q.; Yang, S.; and Chen, X. 2021.  TTN/TP53 mutation might act as the predictor for chemotherapy response in lung adenocarcinoma and lung squamous carcinoma patients.  _Translational Cancer Research_ , 10(3).
  * Yakaboski and Santos (2018) Yakaboski, C.; and Santos, E. 2018.  Bayesian Knowledge Base Distance-Based Tuning.  In _2018 IEEE/WIC/ACM International Conference on Web Intelligence, WI 2018, Santiago, Chile, December 3-6, 2018_ , 64–72. IEEE Computer Society.
  * Yakaboski and Santos (2021) Yakaboski, C.; and Santos, E. 2021.  Efficient Reasoning upon Fusion of Many Data Sources.  In Bell, E.; and Keshtkar, F., eds., _Proceedings of the Thirty-Fourth International Florida Artificial Intelligence Research Society Conference, North Miami Beach, Florida, USA, May 17-19, 2021_.
  * Yang and Chang (2002) Yang, S.; and Chang, K.-C. 2002.  Comparison of score metrics for Bayesian network learning.  In _IEEE Transactions on Systems, Man, and Cybernetics \- Part A: Systems and Humans_ , volume 32, 419–428.
  * Yun and Keong (2004) Yun, Z.; and Keong, K. 2004.  Improved MDL score for learning of Bayesian networks.  In _Proceedings of the International Conference on Artificial Intelligence in Science and Technology, AISAT_ , 98–103.

## Appendices

##  Appendix A Experimental Results and Design

Table 1: Data MDL results for BKBs learned with our BKB learning algorithm with GOBNILP DAG learning backend as well as BNs learned over same dataset with standard MDL scoring and GOBNILP. Also reported are the number of unique calls to the joint probability calculator for both BKB and BN learning. Dataset | BKB MDL | BN MDL | BKB Calls | BN Calls | Calls Difference
---|---|---|---|---|---
newthyroid | -73 | -462 | 1257 | 7200 | -5943
car | -994 | -5976 | 20369 | 40000 | -19631
banana | -2047 | -5929 | 103 | 135 | -32
poker | -976563 | -5263517 | 945 | 946 | -1
heart | -249 | -1572 | 1666695 | 26262720 | -24596025
saheart | -353 | -1866 | 48005 | 104976 | -56971
wine | -134 | -1103 | 1018988 | 39249360 | -38230372
mammographic | -300 | -2126 | 2735 | 6750 | -4015
pima | -628 | -3059 | 28309 | 43740 | -15431
housevotes | -494 | -1866 | 3209612 | 9746883 | -6537271
tic-tac-toe | -1007 | -4888 | 250986 | 786432 | -535446
winequality-red | -1087 | -7423 | 663679 | 4815936 | -4152257
ecoli | -199 | -875 | 9044 | 108864 | -99820
contraceptive | -1060 | -6083 | 194374 | 1080000 | -885626
yeast | -391 | -4624 | 47816 | 380160 | -332344
nursery | -7223 | -59544 | 307592 | 691200 | -383608
titanic | -957 | -2586 | 102 | 108 | -6
hayes-roth | -87 | -403 | 928 | 2000 | -1072
iris | -69 | -394 | 322 | 1024 | -702
bupa | -334 | -760 | 1248 | 2187 | -939
haberman | -104 | -384 | 80 | 81 | -1
winequality-white | -2294 | -24162 | 268095 | 624312 | -356217
glass | -155 | -841 | 42539 | 680400 | -637861
flare | -394 | -3407 | 275011 | 9933408 | -9658397
phoneme | -1265 | -13925 | 12389 | 55440 | -43051
bands | -290 | -1575 | 1038847 | 3064209 | -2025362
australian | -507 | -3793 | 1304495 | 6118368 | -4813873
tae | -65 | -263 | 540 | 972 | -432
post-operative | -61 | -319 | 16202 | 184320 | -168118
vowel | -387 | -4745 | 958224 | 13288429 | -12330205
balance | -505 | -1503 | 312 | 324 | -12
zoo | -59 | -650 | 1646540 | 37950312 | -36303772
hepatitis | -82 | -529 | 3634610 | 12986769 | -9352159
crx | -242 | -3788 | 1252575 | 7047840 | -5795265
cleveland | -333 | -1761 | 1536035 | 8398578 | -6862543
lymphography | -109 | -1083 | 2586093 | 33062432 | -30476339
breast | -122 | -1154 | 116161 | 5225472 | -5109311
led7digit | -122 | -1599 | 10204 | 24057 | -13853
kr-vs-k | -7629 | -78926 | 80933 | 146875 | -65942
monk-2 | -229 | -1548 | 6696 | 8640 | -1944
Table 2: Accuracy results on subset of benchmark datasets for both our learned BKBs and BNs using a 10-fold cross validation reported in terms of precision, recall and F1-score. | Precision |  | Recall |  | F1-Score |  | # Failed
---|---|---|---|---|---|---|---
Dataset | BKB | BN | BKB | BN | BKB | BN | BKB
newthyroid | 0.900 | 0.687 | 0.892 | 0.809 | 0.886 | 0.738 | 2
hayes-roth | 0.572 | 0.839 | 0.481 | 0.838 | 0.448 | 0.837 | 0
iris | 0.948 | 0.940 | 0.947 | 0.940 | 0.947 | 0.940 | 0
tae | 0.262 | 0.143 | 0.397 | 0.212 | 0.315 | 0.171 | 0
haberman | 0.744 | 0.739 | 0.742 | 0.758 | 0.644 | 0.743 | 0
banana | 0.704 | 0.734 | 0.698 | 0.730 | 0.688 | 0.725 | 0
titanic | 0.458 | 0.770 | 0.677 | 0.776 | 0.547 | 0.762 | 0
led7digit | 0.608 | 0.016 | 0.644 | 0.070 | 0.616 | 0.026 | 0
monk-2 | 0.974 | 0.870 | 0.972 | 0.861 | 0.972 | 0.859 | 0
balance | 0.669 | 0.522 | 0.726 | 0.566 | 0.697 | 0.543 | 0
mammographic | 0.838 | 0.841 | 0.831 | 0.841 | 0.830 | 0.841 | 0
bupa | 0.550 | 0.630 | 0.516 | 0.629 | 0.513 | 0.629 | 0
phoneme | 0.767 | 0.809 | 0.779 | 0.772 | 0.766 | 0.781 | 1
car | 0.752 | 0.490 | 0.277 | 0.700 | 0.187 | 0.577 | 0
post-operative | 0.586 | 0.508 | 0.476 | 0.713 | 0.501 | 0.593 | 3
ecoli | 0.595 | 0.181 | 0.725 | 0.426 | 0.651 | 0.254 | 1
yeast | 0.213 | 0.097 | 0.354 | 0.312 | 0.265 | 0.148 | 0
pima | 0.757 | 0.689 | 0.672 | 0.701 | 0.677 | 0.691 | 0
glass | 0.215 | 0.126 | 0.430 | 0.355 | 0.278 | 0.186 | 0
saheart | 0.730 | 0.669 | 0.416 | 0.684 | 0.321 | 0.670 | 0
tic-tac-toe | 0.589 | 0.683 | 0.652 | 0.697 | 0.554 | 0.683 | 9
contraceptive | 0.051 | 0.331 | 0.226 | 0.417 | 0.083 | 0.281 | 0
Table 3: Benchmark dataset information including the maximum parent set limit used in our experiments for both BKB and BN learning. Dataset | Parent Set Limit | # Features | # I-nodes | # Examples
---|---|---|---|---
newthyroid | 5 | 6 | 21 | 215
car | 6 | 7 | 25 | 1728
banana | 2 | 3 | 14 | 5300
poker | 1 | 11 | 45 | 1025009
heart | 10 | 14 | 35 | 270
saheart | 9 | 10 | 22 | 462
wine | 8 | 14 | 40 | 178
mammographic | 5 | 6 | 21 | 830
pima | 8 | 9 | 21 | 768
housevotes | 7 | 17 | 34 | 232
tic-tac-toe | 9 | 10 | 29 | 958
winequality-red | 10 | 12 | 33 | 1599
ecoli | 7 | 8 | 29 | 336
contraceptive | 9 | 10 | 31 | 1473
yeast | 8 | 9 | 32 | 1484
nursery | 8 | 9 | 32 | 12960
titanic | 3 | 4 | 9 | 2201
hayes-roth | 4 | 5 | 18 | 160
iris | 4 | 5 | 15 | 150
bupa | 6 | 7 | 14 | 345
haberman | 3 | 4 | 8 | 306
winequality-white | 4 | 12 | 45 | 4898
glass | 9 | 10 | 30 | 214
flare | 6 | 12 | 47 | 1066
phoneme | 5 | 6 | 34 | 5404
bands | 5 | 20 | 40 | 365
australian | 6 | 15 | 40 | 690
tae | 5 | 6 | 13 | 151
post-operative | 8 | 9 | 26 | 87
vowel | 5 | 14 | 60 | 990
balance | 4 | 5 | 11 | 625
zoo | 7 | 17 | 43 | 101
hepatitis | 6 | 20 | 40 | 80
crx | 5 | 16 | 54 | 653
cleveland | 10 | 14 | 31 | 297
lymphography | 5 | 19 | 63 | 148
breast | 9 | 10 | 43 | 277
led7digit | 7 | 8 | 24 | 500
kr-vs-k | 3 | 7 | 58 | 28056
monk-2 | 6 | 7 | 19 | 432

###  A.1 Benchmark Experiments

All benchmark datasets were taken directly from the KEEL data repository
(Alcalá-Fdez et al. 2011) and contained no missing values in order to
facilitate direct BNSL comparisons. All preprocessed KEEL datasets used in our
analyses are hosted on
Zenodo222https://zenodo.org/record/6584753#.Yo˙PK1RBxPZ. All 10-fold cross
validation experiments analysis used splits already implemented in the KEEL
repository for each dataset. All experiments were performed on cluster of 15
Dell PowerEdge M630 blade servers each with 32 processors and 500GB of RAM and
leveraging the Ray (Moritz et al. 2018) library to distribute all
calculations. While distributed, we still utilized non-optimized python code
for all GOBNILP score calculations limiting the number of parent sets that
could be explored, particularly on larger problems. In order to arrive at each
parent set limit in Table 3, we iterated parent set limits starting from one
to $n_{i}$ where $n_{i}$ is the parent set limit for each dataset $i$ that was
able to complete all necessary score calculations within 12 hours. For all MDL
and accuracy analysis we learned BKBs and BNs on their respective highest
parent set limit $n_{i}$ scores calculated for each dataset in order to
explore the greatest possible DAG space available. We calculated all joint
probabilities based on raw dataset counts without employing any estimation
techniques. All continuous valued features where discretized into at least two
splits using the Entropy MDL method (Fayyad and Irani 1993).

###  A.2 TCGA Breast Cancer Experiments

As there are over 500 features in our TCGA dataset we will not list them here
and instead point the reader to the pre-processed dataset available on
Zenodo333https://zenodo.org/record/6584753#.Yo˙PK1RBxPZ. Features annotated as
{Gene Name}-{Chromosome #}-{Position #} are SNP mutations with states
represented as a tuple (Reference Allele, Tumor Allele 1, Tumor Allele 2) and
features represented as just gene names (HGNC Symbols(Tweedie et al. 2021))
have states corresponding to the variant mutation classification that occurred
anywhere in that gene, i.e., missense mutation, nonsense mutation, etc.

As many mutations only occurred in a single case, we used a feature selection
process to filter out many mutations that would not have been statistically
significant. Our feature selection process involved a threshold for how many
cases a particular mutation and state appeared. For the SNP features we set an
extremely low inclusion threshold of 2, resulting in 72 unique SNP features,
and for whole gene features we set a threshold of 10, resulting in 428
features. Our filtered dataset contained 939 breast cancer cases and our
parent set limit was set to 1 in order to control complexity for analysis.

##  Appendix B Proofs

###  B.1 Theorem 1

$C_{D}(P,Q)$ is a monotonically decreasing function of

$\sum_{\tau\in
D_{\neq}}\sum_{i=1}^{n}p(x_{i\tau},\pi_{i\tau})\log_{2}\frac{p(x_{i\tau},\pi_{i\tau})}{p(x_{i\tau})p(\pi_{i\tau})}$
(7)

where $x_{i\tau}$ is the instantiation of random variable $X_{i}$ determined
by data instance $\tau$ and $\pi_{i\tau}$ is the parent set instantiation of
random variable $X_{i}$ governed by $\tau$.

###### Proof.

We can write $C_{D}(P,Q)=\sum_{\tau\in D_{\neq}}C_{\tau}(P,Q)$ where

$\displaystyle C_{\tau}(P,Q)$ $\displaystyle=p_{\tau}\log
p_{\tau}-p_{\tau}\log q_{\tau}$ (8) $\displaystyle=p_{\tau}\log
p_{\tau}-p_{\tau}\sum_{i=1}^{n}\log p(x_{i\tau}\mid\pi_{i\tau})$ (9)
$\displaystyle=p_{\tau}\log
p_{\tau}-p_{\tau}\sum_{i=1}^{n}\log\bigg{(}\frac{p(x_{i\tau},\pi_{i\tau})}{p(x_{i\tau})p(\pi_{i\tau})}\bigg{)}-p_{\tau}\sum_{i=1}^{n}\log
p(x_{i\tau})$ (10)

Now, consider that $p(x_{i\tau})=p_{\tau}+\tilde{p}(x_{i\tau})$, such that
$\tilde{p}(x_{i\tau})=P(X_{1}\neq
x_{1\tau},\dots,X_{i}=x_{i\tau},\dots,X_{n}\neq x_{n\tau})$ and similarly,
$p(x_{i\tau},\pi_{i\tau})=p_{\tau}+\tilde{p}(x_{i\tau},\pi_{i\tau})$. Then,

$-p_{\tau}\sum_{i=1}^{n}\log p(x_{ij})=-\sum_{i=1}^{n}p(x_{i\tau})\log
p(x_{i\tau})+\sum_{i=1}^{n}\tilde{p}(x_{i\tau})\log\tilde{p}(x_{i\tau})$ (11)

and

$\displaystyle-
p_{\tau}\sum_{i=1}^{n}\log\bigg{(}\frac{p(x_{i\tau},\pi_{i\tau})}{p(x_{i\tau})p(\pi_{i\tau})}\bigg{)}=-\sum_{i=1}^{n}$
$\displaystyle
p(x_{i\tau},\pi_{i_{\tau}})\log\bigg{(}\frac{p(x_{i\tau},\pi_{i\tau})}{p(x_{i\tau})p(\pi_{i\tau})}\bigg{)}$
(12)
$\displaystyle+\sum_{i=1}^{n}\tilde{p}(x_{i\tau},\pi_{i_{\tau}})\log\bigg{(}\frac{p(x_{i\tau},\pi_{i\tau})}{p(x_{i\tau})p(\pi_{i\tau})}\bigg{)}$
(13)

Further we can define the instantiated entropy as
$H_{1}(x_{ij})=-p(x_{ij})\log p(x_{ij})$ to yield

$\displaystyle C_{\tau}(P,Q)=-\sum_{i=1}^{n}$ $\displaystyle
p(x_{i\tau},\pi_{i_{\tau}})\log\bigg{(}\frac{p(x_{i\tau},\pi_{i\tau})}{p(x_{i\tau})p(\pi_{i\tau})}\bigg{)}+\sum_{i=1}^{n}\tilde{p}(x_{i\tau},\pi_{i_{\tau}})\log\bigg{(}\frac{p(x_{i\tau},\pi_{i\tau})}{p(x_{i\tau})p(\pi_{i\tau})}\bigg{)}$
(14)
$\displaystyle+\sum_{i=1}^{n}H_{1}(x_{i\tau})+\sum_{i=1}^{n}\tilde{p}(x_{i\tau})\log\tilde{p}(x_{i\tau})-H_{1}(\tau)$
(15)

The last three terms are functions of the data and not the underlying
structure of the model and thereby a constant when comparing networks. In all
worlds, it must be the case that
$p(x_{i\tau},\pi_{i_{\tau}})\geq\tilde{p}(x_{i\tau},\pi_{i_{\tau}})$.
Therefore, $C_{D}(P,Q)$ is a monotonically decreasing function of Equation 7.
∎

###  B.2 Theorem 2

Given a BN $G$ learned over a dataset $D$ as to maximize the weight $W_{G}$
given by Equation 30 and corresponding weight $W_{G}^{\prime}$ given by
Equation 6, there exists a BKB $K$ with a weight $W_{K}$ according to Equation
6 such that $W_{K}\geq W_{G}^{\prime}$.

###### Proof.

Assume we have a BN $G$ with weight $W_{G}$ according to Equation 30 and
$W_{G}^{\prime}$ according to Equation 6. We can transform $G$ into a BKB
$K_{G}$ such that an S-node is created corresponding to every entry in $G$’s
conditional probability table (CPT). To complete these proofs we need to
leverage a few definitions.

##### Definition B.2.1

A set of S-nodes $B_{i}$ is said to be a causal rule set (CRS) for the random
variable $X_{i}$ if $B_{i}$ contains all S-nodes pointing to the
instantiations of $X_{i}$.

##### Definition B.2.2

For every set of I-nodes $\alpha$ over unique random variables there is a set
of possible S-nodes $\psi$ such that every S-node $q\in\psi$ corresponds to a
permutation over $\alpha$. In other worlds, this set is all the unique S-nodes
that can be formed from $\alpha$ using all instantiations in $\alpha$.
Therefore, the cardinality of $|\psi|=|\alpha|$.

For the BKB $K_{G}$ the random variables associated with the parents of every
S-node in each CRS $B_{i}$ will be $\Pi_{G}(X_{i})$. Since a weight according
to Equation 6 is assigned to every S-node, we can also view each CRS as having
a weight $W(B_{i})=\sum_{q\in B_{i}}W(q)$.

##### Proposition

Starting from the BKB $K_{G}$ yielding the CRSs $B_{1},B_{2},\dots,B_{N}$ we
can find a new set of CRSs
$B_{1}^{\prime},B_{2}^{\prime},\dots,B_{N}^{\prime}$ that have a total weight
$\geq W(K_{G})$.

Now our proof reduces to proving the above proposition. Consider an S-node
$q\in B_{i}$ and the associated set of I-nodes $\alpha$ corresponding to $q$.
We know from the above definitions that there is a set $\psi$ of all possible
S-node permutations of the I-node set $\alpha$. Therefore, there must exist an
S-node $q^{\prime}\in\psi$ such that $W(q^{\prime})\geq W(q)$. In the case
that $q^{\prime}\neq q$, we can remove $q$ from $B_{i}$ and add $q^{\prime}$
to the CRS associated with $Head(q^{\prime})$ which will be different from
$Head(q)$. Performing this action will increase the weight of the BKB by a
$\Delta(W(q),W(q^{\prime}))\geq 0$. Therefore, we will always find a new BKB
from a BN that has a weight equal to or greater than that of the original BN.
∎

##  Appendix C Derivations

###  C.1 Derivation of BKB Modeling Cost

First, let’s define model MDL, denoted MMDL, for BKBs and BNs respectively

$MMDL(G_{BN})=\sum_{i=1}^{n}\bigg{(}|\Pi(X_{i})|\log(n)+\delta(r_{i}-1)\prod_{X_{j}\in\Pi(X_{i})}r_{j}\bigg{)}$
(16) $MMDL(G_{BKB})=\sum_{i=1}^{n}\sum_{q\in
CRS(X_{i})}\bigg{(}|\Pi(q)|\log(m)+\delta\bigg{)}$ (17)

where $CRS(X_{i})$ denotes the Causal Rule Set for random variable $X_{i}$ and
all logarithms in this derivation are base 2. Further, we know that $G_{BKB}$
will be complete and is constructed from the $G_{BN}$, and thereby $\forall
q\in CVS(X_{i})\rightarrow K=|\Pi_{BN}(X_{i})|=|\Pi_{BKB}(q)|$. Since we know
that $|\Pi_{BKB}(q)|$ will be a constant for all $q$’s in the same CRS, then
we can simply the summation and just multiply by the number of different
parent set instantiations of $X_{i}$. Therefore, let’s define
$R=\prod_{X_{j}\in\Pi(X_{i})}r_{j}$ yielding the simplified expressions

$MMDL(G_{BN})=\sum_{i=1}^{n}\bigg{(}K\log(n)+\delta(r_{i}-1)R\bigg{)}$ (18)
$MMDL(G_{BKB})=\sum_{i=1}^{n}\bigg{(}K\log(m)+\delta\bigg{)}R$ (19)

Then the cost for converting the BN to a complete BKB is

$\displaystyle C(G_{BKB},G_{BN})$ $\displaystyle=MMDL(G_{BKB})-MMDL(G_{BN})$
(20)
$\displaystyle=\sum_{i=1}^{n}\bigg{(}K\log(m)+\delta\bigg{)}R-\bigg{(}K\log(n)+(2-r_{i})\delta
R\bigg{)}$ (21)
$\displaystyle=\sum_{i=1}^{n}K\big{(}\log(m)R-\log(n)\big{)}+(2-r_{i})\delta
R$ (22)
$\displaystyle=\sum_{i=1}^{n}K\log\bigg{(}\frac{m^{R}}{n}\bigg{)}+(2-r_{i})\delta
R$ (23)

Now let’s define $m^{\prime}=r_{max}^{n}$ where $r_{max}$ the the maximum
number of instantiations of any variable in $X$. This yields

$\displaystyle C(G_{BKB},G_{BN})$
$\displaystyle\leq\sum_{i=1}^{n}K\log\bigg{(}\frac{m^{\prime
R}}{n}\bigg{)}+(2-r_{i})\delta R$ (24)
$\displaystyle=\sum_{i=1}^{n}K\big{(}Rn\log(r_{max})-\log(n)\big{)}+(2-r_{i})\delta
R$ (25)
$\displaystyle=\sum_{i=1}^{n}|\Pi(X_{i})|\big{(}Rn\log(r_{max})-\log(n)\big{)}+(2-r_{i})\delta
R$ (26)

##  Appendix D Other Background

###  D.1 Bayesian Networks

A BN is a Directed Acyclic Graph (DAG) $G=(V,E)$ that represents the
factorized joint probability distribution over random variables
$X=(X_{1},\dots,X_{n})$ of the form:

$\text{Pr}(X)=\prod_{i}^{n}P(X_{i}\mid\Pi(X_{i}))$ (27)

such that $\Pi(X_{i})$, or more concisely $\pi_{i}$, are the parents of the
random variable $X_{i}$ according to $G$ and each node $V_{i}\in V$ correspond
directly to each random variable $X_{i}$.

If a dataset $D$ and the BN structure $G$ are given, then parameters of a BN,
or conditional probabilities, can be readily calculated via Maximum Likelihood
Estimation (MLE) and resolve to simply taking counts over the dataset
(Bouckaert 1995), yielding the log-likelihood:

$LL_{D}(G)=\sum_{i=1}^{n}\sum_{j=1}^{c_{i}}\sum_{k=1}^{r_{i}}N_{ijk}\log\bigg{(}\frac{N_{ijk}}{N_{ij}}\bigg{)}$
(28)

where $r_{i}$ is the number of instantiations available to $X_{i}$, and
$c_{i}=\prod_{u_{j\in\pi_{i}}}r_{j}$ which is the number of possible parent
set $\pi_{i}$ configurations. This leads to the two predominant formulations
of the BN MDL from Lam and Bacchus (Lam and Bacchus 1994) and Suzuki (Suzuki
1999, 1993), namely

$\sum_{i=1}^{n}[k_{i}\log_{2}(n)+\delta(r_{i}-1)c_{i}]-N\sum_{\tau}p_{\tau}\log_{2}q_{\tau}\leq
H(D|G)+\frac{K(G)}{2}\log_{2}(N)$ (29)

where the left side of the inequality describes the MDL characterized by Lam
and Bacchus and the right corresponds to the derivation by Suzuki. Further,
$k_{i}$ is the number of parents of $X_{i}$, $\delta$ is the number of bits
required to encode a probability value, $\tau$ is an atomic world encoded by
the distribution $\text{Pr}(\cdot)$, $H(D|G)$ is the conditional empirical
entropy of the data given the model $G$, and $K(G)$ is the number of
parameters in the BN which is equal to the total number of conditional
probability values.

We can view both of these formulations as equivalent but taking different
perspectives. When calculating the model encoding term, Lam and Bucchus’s
formulation derives from the arithmetic coding literature and Suzuki leverages
a ”universal” or Lempel-Ziv coding paradigm (MacKay 2005).

###  D.2 Bayesian Knowledge Bases

Figure 3: Example of a BKB. The greyed I-nodes and S-nodes represent an
example inference (world) $\tau$ in the BKB.

In this section we will provide some useful definitions for better
understanding BKBs as well as a description of BKB fusion. A graphical
depiction of an example BKB can be seen in Figure 3. In this example we can
see the way in which the conditional (in)dependency relationships are defined
on the random variable instantiation level as opposed to the random variable
level. Figure 3 exemplifies how BKBs can encode cyclic dependencies on the
random variable level such as between random variables $C$ and $D$ as well as
forming an incomplete distribution in which all conditional probabilities
between random variables need not be specified.

#### Definitions

##### Correlation Graph

A correlation-graph is a directed graph $G=(I\cup S,E)$ s.t. $I\cap
S=\emptyset$, $E=\\{I\times S\\}\cup\\{I\times S\\}$, and $\forall q\in
S\exists\alpha\in I$, such that $(q,\alpha)\in E$, where $\alpha\in I$
represent the I-nodes and $q\in S$ represent the S-nodes.

##### BKB Inference

An inference is a tuple $\tau=(I^{\prime}\cup S^{\prime},E^{\prime})$ which
forms a subgraph over $K$ where $I^{\prime}\subseteq I$, $S^{\prime}\subseteq
S$, and $E^{\prime}\subseteq E$ and has a weight $w(\tau)=\prod_{a\in
S^{\prime}}w(q)$ that has the following properties:

  1. 1.

Every I-node $\alpha\in I^{\prime}$ is well-supported s.t. $\forall\alpha$
$\exists$ $(q,\alpha)inE^{\prime}$, i.e., all I-nodes have an incoming S-node,

  2. 2.

Every S-node $q\in S^{\prime}$is well-founded s.t. $\forall(\alpha,q)\in
E,(\alpha,q)\in E^{\prime}$, i.e., all incoming edges to an S-node included in
$\tau$ have all incoming edges associated with governing correlation graph,

  3. 3.

All S-nodes are well-defined s.t. $\forall q$ $\exists$ $(q,\alpha)\in
E^{\prime}$, i.e., there are no S-nodes without a head I-node,

  4. 4.

The subgraph representing $\tau$ is acyclic,

  5. 5.

Each random variable has at most one instantiation.

###  D.3 BKB Fusion

(a)

(b)

Figure 4: (a) Example of two BKB fragments that violate the mutual exclusion
property of a BKB. This can be seen by viewing the parent set of each S-node
incident on the $X_{1}=x_{1}$ and noticing that these two sets are not mutex.
Therefore, a simple union of these S-nodes would not be possible as it would
violate our BKB semantics. (b) The BKB fusion solution that creates two
distinct source I-nodes for $X_{1}=x_{1}$ thereby making the two parent sets
mutually exclusive.

The core idea behind BKB fusion is to take the union of all input BKB
Fragments (BKFs), which are also valid BKBs, and attach a source I-node to
every incident S-node on random variable $X_{i}$, denoted $S_{X_{i}}$, with an
instantiation equaling the designation of the source BKF, denoted $\sigma$,
from which the S-node originated. Further, BKB fusion takes into account the
reliability or probability of correctness of each source such that more
correct or reliable sources will have greater weight. Intuitively, adding
source I-nodes in this way transforms a rule described by an S-node such as
“If $A$ is true, then $B$ is true with probability $p$” into “If source
$\sigma$ is correct and $A$ is true, then $B$ is true with probability $p$.” A
graphical depiction of this procedure can be seen in Appendix Figure 4(b) and
an exact algorithm follows.

Given a set of $n$ BKFs $\\{K_{1},K_{2},\dots,K_{n}\\}$ where
$K_{i}=(G_{i},w_{i},\sigma_{i},r(\sigma_{i}))$, $G_{i}=(I_{i},S_{i},E_{i})$,
$\sigma_{i}$ is the name of the BKF source, and $r(\sigma_{i})$ is the
reliability index of that source. Then the output of Algorithm 2 is a new BKB,
$K^{\prime}=(G^{\prime},w^{\prime})$ with
$G^{\prime}=(I^{\prime},S^{\prime},E^{\prime})$ that is a fusion of these
fragments. For an I-node $\alpha$ in some fragment, let us denote $R_{\alpha}$
as the random variable corresponding to the instantiation $\alpha$. For a
source node $s=(S_{R_{\alpha}}=\sigma_{i})$ let the reliability of that $s$ be
$r(s)=r(\sigma_{i})$. A more detailed discussion can be found in the work of
Santos et al (Santos, Wilkinson, and Santos 2011).

Figure 5: An image depicting BKB generalization through the fusion of four BKB
fragments learned in the context of four different data exemplars. The circles
above represent the random variables $A$, $B$, and $C$ each having 3 states.
The four learned BKFs are denoted $\tau_{i}$ and the dependency structure for
each BKF in $\tau_{i}$ is designated by their respective colored arrows. The
bottom portion of the figure shows four other world instantations that could
be observed in unseen data and it is clear that our BKB can generalize with
fusion to three of the four unseen worlds.

##### BKB Generalization through Fusion

To better clarify the mechanism by which BKB fusion allows our BKBSL algorithm
to generalize to unseen data we leverage the visualization in Figure 5. Figure
5 depicts an example collection of three random variables, $A$, $B$, and $C$
each having three discrete states 1, 2, and 3. This random variable
specification yields a total of $3^{3}=27$ unique worlds that can be observed
from any given dataset. Therefore, let’s assume we observe four unique data
exemplars (worlds) denoted $\tau_{1}=\\{A=1,B=1,C=1\\}$,
$\tau_{2}=\\{A=2,B=2,C=2\\}$, $\tau_{3}=\\{A=3,B=2,C=1\\}$, and
$\tau_{4}=\\{A=1,B=3,C=3\\}$ and are illustrated by the colored boxes in
Figure 5. Next, from these four worlds we learn four corresponding BKFs
depicted by the colored arrows, where each colored arrow represents a S-node.
By fusing these four fragments, our fused BKB allows us to generalize to
unseen examples (the bottom illustration of Figure 5) by combining S-nodes
from different fragments to form a complete inference, weighted by each
sources’ reliability which, in the case our BKBSL implementation, are all the
same. Therefore, using BKB fusion we can generalize to at least 3 more worlds
that we have never been observed and cannot generalize to the world
$\\{A=1,B=2,C=1\\}$ as there is no way to combine S-nodes from these four BKFs
that yield a valid complete inference.

Algorithm 2 BKB-Fusion

Inputs: BKFs $K_{1},K_{2},\dots,K_{n}$ and Source Reliabilities $R$

1: Let $G^{\prime}=(I^{\prime},S^{\prime},E^{\prime})$ be an empty correlation
graph

2: $I^{\prime}\leftarrow\cup_{i}=1^{n}I_{i}$

3: $S^{\prime}\leftarrow\cup_{i}=1^{n}S_{i}$

4: $E^{\prime}\leftarrow\cup_{i}=1^{n}E_{i}$

5: for $K_{i}\in\\{K_{1},K_{2},\dots,K_{n}\\}$ do

6: for all S-nodes $q\in S_{i}$ do

7: Let $\alpha\leftarrow Head(q)$

8: Let the source I-node for $q$ be $s=(S_{R_{\alpha}}=\sigma_{i})$

9: Add $s$ to $I^{\prime}$ and add new S-node $q_{s}$ to $S^{\prime}$

10: Add edges $q_{s}\rightarrow s$ and $s\rightarrow q$ to $E^{\prime}$

11: Let $w^{\prime}(q)\leftarrow w(q)$

12: end for

13: end for

14: for all source variables $S_{R_{\alpha}}$ do

15: Let $\Lambda=\\{s\mid s\text{ is a source node which is a state of
}S_{R_{\alpha}}\\}$

16: $\rho\leftarrow\sum_{s\in\Lambda}r(s)$

17: for $s\in\Lambda$, let $q_{s}$ be the S-node such that $q_{s}\rightarrow
s\in E^{\prime}$ do

18: Let $w^{\prime}(q)\leftarrow\frac{r(s)}{\rho}$

19: end for

20: end for

21: return $K^{\prime}=(G^{\prime},w^{\prime})$

###  D.4 Additional Theorems

######  Theorem 3.

$C(P,Q)$ is a monotonically decreasing function of

$\sum_{i=1,\Pi(X_{i})\neq\emptyset}^{n}W(X_{i},\Pi(X_{i}))$ (30)

where

$W(X_{i},\Pi(X_{i}))=\sum_{j=1}^{r_{i}}\sum_{k=1}^{c_{i}}p(x_{ij},\pi_{ijk})\log_{2}\frac{p(x_{ij},\pi_{ijk})}{p(x_{ij})p(\pi_{ijk})}$
(31)

Therefore, $C(P,Q)$ will be minimized if and only if the sum is maximized.