Entropy and the Kullback–Leibler Divergence for Bayesian Networks: Computational Complexity and Efficient Implementation

Abstract

Bayesian networks (BNs) are a foundational model in machine learning and causal inference. Their graphical structure can handle high-dimensional problems, divide them into a sparse collection of smaller ones, underlies Judea Pearl’s causality, and determines their explainability and interpretability. Despite their popularity, there are almost no resources in the literature on how to compute Shannon’s entropy and the Kullback–Leibler (KL) divergence for BNs under their most common distributional assumptions. In this paper, we provide computationally efficient algorithms for both by leveraging BNs’ graphical structure, and we illustrate them with a complete set of numerical examples. In the process, we show it is possible to reduce the computational complexity of KL from cubic to quadratic for Gaussian BNs.

1. Introduction

Bayesian networks [1] (BNs) have played a central role in machine learning research since the early days of the field as expert systems [2,3], graphical models [4,5], dynamic and latent variables models [6], and as the foundation of causal discovery [7] and causal inference [8]. They have also found applications as diverse as comorbidities in clinical psychology [9], the genetics of COVID-19 [10], the Sustainable Development Goals of the United Nations [11], railway disruptions [12] and industry 4.0 [13].

Machine learning, however, has evolved to include a variety of other models and reformulated them into a very general information-theoretic framework. The central quantities of this framework are Shannon’s entropy and the Kullback–Leibler divergence. Learning models from data relies crucially on the former to measure the amount of information captured by the model (or its complement, the amount of information lost in the residuals) and on the latter as the loss function we want to minimise. For instance, we can construct variational inference [14], the Expectation-Maximisation algorithm [15], Expectation Propagation [16] and various dimensionality reduction approaches such as t-SNE [17] and UMAP [18] using only these two quantities. We can also reformulate classical maximum-likelihood and Bayesian approaches to the same effect, from logistic regression to kernel methods to boosting [19,20].

Therefore, the lack of literature on how to compute the entropy of a BN and the Kullback–Leibler divergence between two BNs is surprising. While both are mentioned in Koller and Friedman [5] and discussed at a theoretical level in Moral et al. [21] for discrete BNs, no resources are available on any other type of BN. Furthermore, no numerical examples of how to compute them are available even for discrete BNs. We fill this gap in the literature by:

• Deriving efficient formulations of Shannon’s entropy and the Kullback–Leibler divergence for Gaussian BNs and conditional linear Gaussian BNs.
• Exploring the computational complexity of both for all common types of BNs.
• Providing step-by-step numeric examples for all computations and all common types of BNs.

Our aim is to make apparent how both quantities are computed in their closed-form exact expressions and what is the associated computational cost.

The common alternative is to estimate both Shannon’s entropy and the Kullback–Leibler divergence empirically using Monte Carlo sampling. Admittedly, this approach is simple to implement for all types of BNs. However, it has two crucial drawbacks:

• Using asymptotic estimates voids the theoretical properties of many machine learning algorithms: Expectation-Maximisation is not guaranteed to converge [5], for instance.
• The number of samples required to estimate the Kullback–Leibler divergence accurately on the tails of the global distribution of both BNs is also an issue [22], especially when we need to evaluate it repeatedly as part of some machine learning algorithm. The same is true, although to a lesser extent, for Shannon’s entropy as well. In general, the rate of convergence to the true posterior in Monte Carlo particle filters is proportional to the number of variables squared [23].

Therefore, efficiently computing the exact value of Shannon’s entropy and the Kullback–Leibler divergence is a valuable research endeavour with a practical impact on BN use in machine learning. To help its development, we implemented the methods proposed in the paper in our bnlearn R package [24].

The remainder of the paper is structured as follows. In Section 2, we provide the basic definitions, properties and notation of BNs. In Section 3, we revisit the most common distributional assumptions in the BN literature: discrete BNs (Section 3.1), Gaussian BNs (Section 3.2) and conditional linear Gaussian BNs (Section 3.3). We also briefly discuss exact and approximate inferences for these types of BNs in Section 3.4 to introduce some key concepts for later use. In Section 4, we discuss how we can compute Shannon’s entropy and the Kullback–Leibler divergence for each type of BN. We conclude the paper by summarising and discussing the relevance of these foundational results in Section 5. Appendix A summarises all the computational complexity results from earlier sections, and Appendix B contains additional examples we omitted from the main text for brevity.

2. Bayesian Networks

Bayesian networks (BNs) are a class of probabilistic graphical models defined over a set of random variables $\mathbf{X}$$=\{X_1,\dots ,X_N\}$, each describing some quantity of interest, that are associated with the nodes of a directed acyclic graph (DAG) $\mathcal{G}$. Arcs in $\mathcal{G}$ express direct dependence relationships between the variables in $\mathbf{X}$, with graphical separation in $\mathcal{G}$ implying conditional independence in probability. As a result, $\mathcal{G}$ induces the factorisation

$$P\left(\mathbf{X}\mid \mathcal{G},\Theta \right)=\prod _{i=1}^{N}P\left(X_i\mid {\Pi }_{X_i},{\Theta }_{X_i}\right),$$

(1)

in which the global distribution (of $\mathbf{X}$, with parameters $\Theta$) decomposes into one local distribution for each $X_i$ (with parameters ${\Theta }_{X_i}$, ${\bigcup }_{\mathbf{X}}{\Theta }_{X_i}=\Theta$) conditional on its parents ${\Pi }_{X_i}$.

This factorisation is as effective at reducing the computational burden of working with BNs as the DAG underlying the BN is sparse, meaning that each node $X_i$ has a small number of parents ($\left|{\Pi }_{X_i}\right|<c$, usually with $c\in [2,5]$). For instance, learning BNs from data is only feasible in practice if this holds. The task of learning a BN $\mathcal{B}=(\mathcal{G},\Theta )$ from a data set $\mathcal{D}$ containing n observations comprises two steps:

$$\underset{\mathrm{learning}}{\underbrace{P\left(\mathcal{G},\Theta \mid \mathcal{D}\right)}}=\underset{\mathrm{structure}\mathrm{learning}}{\underbrace{P\left(\mathcal{G}\mid \mathcal{D}\right)}}·\underset{\mathrm{parameter}\mathrm{learning}}{\underbrace{P\left(\Theta \mid \mathcal{G},\mathcal{D}\right)}}.$$

If we assume that parameters in different local distributions are independent [25], we can perform parameter learning independently for each node. Each $X_i\mid {\Pi }_{X_i}$ will have a low-dimensional parameter space ${\Theta }_{X_i}$, making parameter learning computationally efficient. On the other hand, structure learning is well known to be both NP-hard [26] and NP-complete [27], even under unrealistically favourable conditions such as the availability of an independence and inference oracle [28]. However, if $\mathcal{G}$ is sparse, heuristic learning algorithms have been shown to run in quadratic time [29]. Exact learning algorithms, which have optimality guarantees that heuristic algorithms lack, retain their exponential complexity but become feasible for small problems because sparsity allows for tight bounds on goodness-of-fit scores and the efficient pruning of the space of the DAGs [30,31,32].

3. Common Distributional Assumptions for Bayesian Networks

While there are many possible choices for the distribution of $\mathbf{X}$ in principle, the literature has focused on three cases.

3.1. Discrete BNs

Discrete BNs [25] assume that both $\mathbf{X}$ and the $X_i$ are multinomial random variables (The literature sometimes denotes discrete BNs as “dBNs” or “DBNs”; we do not do that in this paper to avoid confusion with dynamic BNs, which are also commonly denoted as “dBNs”). Local distributions take the form

$$\begin{array}{cccc}\hfill & X_i\mid {\Pi }_{X_i}\sim \mathit{Mul}\left({\pi }_{ik\mid j}\right),\hfill & \hfill & {\pi }_{ik\mid j}=P\left(X_i=k\mid {\Pi }_{X_i}=j\right);\hfill \end{array}$$

their parameters are the conditional probabilities of $X_i$ given each configuration of the values of its parents, usually represented as a conditional probability table (CPT) for each $X_i$. The ${\pi }_{ik\mid j}$ can be estimated from data via the sufficient statistic $\{n_{ijk},$ $i=1,\dots N;$ $j=1,\dots ,q_i;$ $k=1,\dots ,r_i\}$, the corresponding counts tallied from $\{X_i,{\Pi }_{X_i}\}$ using maximum likelihood, Bayesian or shrinkage estimators as described in Koller and Friedman [5] and Hausser and Strimmer [33].

The global distribution takes the form of an N-dimensional probability table with one dimension for each variable. Assuming that each $X_i$ takes at most l values, the table will contain $\left|\mathit{Val}\left(\mathbf{X}\right)\right|=O\left(l^N\right)$ cells, where $\mathit{Val}(·)$ denotes the possible (configurations of the) values of its argument. As a result, it is impractical to use for medium and large BNs. Following standard practices from categorical data analysis [34], we can produce the CPT for each $X_i$ from the global distribution by marginalising (that is, summing over) all the variables other than $\{X_i,{\Pi }_{X_i}\}$ and then normalising over each configuration of ${\Pi }_{X_i}$. Conversely, we can compose the global distribution from the local distributions of the $X_i$ by multiplying the appropriate set of conditional probabilities. The computational complexity of the composition is $O\left(N{l}^N\right)$ because applying (1) for each of the $l^N$ cells yields

$$P\left(\mathbf{X}=\mathbf{x}\right)=\prod _{i=1}^{N}P\left(X_i=x_i\mid {\Pi }_{X_i}={\mathbf{x}}_{{\Pi }_{X_i}}\right),$$

which involves N multiplications. As for the decomposition, for each node, we:

• Sum over $N-|{\Pi }_{X_i}|-1$ variables to produce the joint probability table for $\{X_i,{\Pi }_{X_i}\}$, which contains $O\left(l^{|{\Pi }_{X_i}|+1}\right)$ cells. The value of each cell is the sum of $O\left(l^{N-|{\Pi }_{X_i}|-1}\right)$ probabilities.
• Normalise the columns of the joint probability table for $\{X_i,{\Pi }_{X_i}\}$ over each of the $O\left(l^{|{\Pi }_{X_i}|}\right)$ configurations of values of ${\Pi }_{X_i}$, which involves summing O(l) probabilities and dividing them by their total.

The resulting computational complexity is

$$\underset{\mathrm{marginalisation}}{\underbrace{O(l^{|{\Pi }_{X_i}|+1}·l^{N-|{\Pi }_{X_i}|-1})}}+\underset{\mathrm{normalisation}}{\underbrace{O(l·l^{|{\Pi }_{X_i}|})}}=O(l^N+l^{|{\Pi }_{X_i}|+1})$$

(2)

for each node and $O(N{l}^N+l{\sum }_{i=1}^N{l}^{|{\Pi }_{X_i}|})$ for the whole BN.

Example 1 

(Composing and decomposing a discrete BN). For reasons of space, this example is presented as Example A1 in Appendix B.

3.2. Gaussian BNs

Gaussian BNs [35] (GBNs) model $\mathbf{X}$ with a multivariate normal random variable N(${\mathit{\mu }}_{\mathcal{B}}$, ${\Sigma }_{\mathcal{B}})$ and assume that the $X_i$ are univariate normals linked by linear dependencies,

$$X_i\mid {\Pi }_{X_i}\sim N({\mu }_{X_i}+{\Pi }_{X_i}{\mathit{\beta }}_{X_i},{\sigma }_{X_i}^2),$$

(3)

which can be equivalently written as linear regression models of the form

$$\begin{array}{cccc}\hfill & X_i={\mu }_{X_i}+{\Pi }_{X_i}{\mathit{\beta }}_{X_i}+{\epsilon }_{X_i},\hfill & \hfill & {\epsilon }_{X_i}\sim N(0,{\sigma }_{X_i}^2).\hfill \end{array}$$

(4)

The parameters in (3) and (4) are the regression coefficients ${\mathit{\beta }}_{X_i}$ associated with the parents ${\Pi }_{X_i}$, an intercept term ${\mu }_{X_i}$ and the variance ${\sigma }_{X_i}^2$. They are usually estimated by maximum likelihood, but Bayesian and regularised estimators are available as well [1].

The link between the parameterisation of the global distribution of a GBN and that of its local distributions is detailed in Pourahmadi [36]. We summarise it here for later use.

• Composing the global distribution. We can create an $N\times N$ lower triangular matrix $C_{\mathcal{B}}$ from the regression coefficients in the local distributions such that $C_{\mathcal{B}}{C}_{\mathcal{B}}^{\mathrm{T}}$ gives ${\Sigma }_{\mathcal{B}}$ after rearranging rows and columns. In particular, we:
  • Arrange the nodes of $\mathcal{B}$ in the (partial) topological ordering induced by $\mathcal{G}$, denoted $X_{\left(i\right)},i=1,\dots ,N$.
  • The ith row of $C_{\mathcal{B}}$ (denoted $C_{\mathcal{B}}[i;•$ ], i = 1, …, N) is associated with X_(i). We compute its elements from the parameters of X_(i) | ${\Pi }_{X_i}$ as
    $\begin{array}{cccccc}\hfill & C_{\mathcal{B}}[i;i]=\sqrt{{\sigma }_{X_{\left(i\right)}}^2}\hfill & \hfill & \mathrm{and}\hfill & \hfill & C_{\mathcal{B}}[i;•\hfill \end{array}]={\mathit{\beta }}_{X_i}{C}_{\mathcal{B}}[{\Pi }_{X_{\left(i\right)}};•],$ where $C_{\mathcal{B}}[{\Pi }_{X_{\left(i\right)}};•$ ] are the rows of $C_{\mathcal{B}}$ that correspond to the parents of X_(i). The rows of $C_{\mathcal{B}}$ are filled following the topological ordering of the BN.
  • Compute ${{\displaystyle \sum ^{\sim }}}_{\mathcal{B}}=C_{\mathcal{B}}{C}_{\mathcal{B}}^{\mathrm{T}}$.
  • Rearrange the rows and columns of ${{\displaystyle \sum ^{\sim }}}_{\mathcal{B}}$ to obtain ${\sum }_{\mathcal{B}}$.
  Intuitively, we construct $C_{\mathcal{B}}$ by propagating the node variances along the paths in $\mathcal{G}$ while combining them with the regression coefficients, which are functions of the correlations between adjacent nodes. As a result, $C_{\mathcal{B}}{C}_{\mathcal{B}}^{\mathrm{T}}$ gives ${\Sigma }_{\mathcal{B}}$ after rearranging the rows and columns to follow the original ordering of the nodes.
  The elements of the mean vector ${\mathit{\mu }}_{\mathcal{B}}$ are similarly computed as $E\left(X_{\left(i\right)}\right)={\Pi }_{X_{\left(i\right)}}{\mathit{\beta }}_{X_i}$ iterating over the variables in topological order.
• Decomposing the global distribution. Conversely, we can derive the matrix $C_{\mathcal{B}}$ from ${\Sigma }_{\mathcal{B}}$ by reordering its rows and columns to follow the topological ordering of the variables in $\mathcal{G}$ and computing its Cholesky decomposition. Then
  $$R={\mathrm{I}}_N-diag\left(C_{\mathcal{B}}\right)C_{\mathcal{B}}^{-1},$$

  contains the regression coefficients ${\mathit{\beta }}_{X_i}$ in the elements corresponding to $X_{\left(i\right)},{\Pi }_{X_{\left(i\right)}}$ (Here $diag\left(C_{\mathcal{B}}\right)$ is a diagonal matrix with the same diagonal elements as $C_{\mathcal{B}}$ and ${\mathrm{I}}_N$ is the identity matrix.) Finally, we compute the intercepts ${\mu }_{X_i}$ as ${\mathit{\mu }}_{\mathcal{B}}-R{\mathit{\mu }}_{\mathcal{B}}$ by reversing the equations we used to construct ${\mathit{\mu }}_{\mathcal{B}}$ above.

The computational complexity of composing the global distribution is bound by the matrix multiplication $C_{\mathcal{B}}{C}_{\mathcal{B}}^{\mathrm{T}}$, which is $O\left(N^3\right)$; if we assume that $\mathcal{G}$ is sparse as in Scutari et al. [29], the number of arcs is bound by some $cN$, computing the ${\mathit{\mu }}_{\mathcal{B}}$ takes O(N) operations. The complexity of decomposing the global distribution is also $O\left(N^3\right)$ because both inverting $C_{\mathcal{B}}$ and multiplying the result by $diag\left(C_{\mathcal{B}}\right)$ are $O\left(N^3\right)$.

Example 2 

(Composing and decomposing a GBN). Consider the GBN $\mathcal{B}$ from Figure 1 top. The topological ordering of the variables defined by $\mathcal{B}$ is $\{\{X_1,X_2\},X_4,X_3\}$, so

$$C_{\mathcal{B}}=\begin{array}{c}\\ X_1\\ X_2\\ X_4\\ X_3\end{array}\begin{array}{cccccc}& X_1& X_2& X_4& X_3& \\ (& \begin{array}{c}0.894\\ 0\\ 1.341\\ 1.610\end{array}& \begin{array}{c}0\\ 0.774\\ 2.014\\ 2.416\end{array}& \begin{array}{c}0\\ 0\\ 1.049\\ 1.258\end{array}& \begin{array}{c}0\\ 0\\ 0\\ 0.948\end{array}& )\end{array}$$

where the diagonal elements are

$$\begin{array}{cccccccc}\hfill & C_{\mathcal{B}}[X_1;X_1]=\sqrt{0.8},\hfill & \hfill & C_{\mathcal{B}}[X_2;X_2]=\sqrt{0.6},\hfill & \hfill & C_{\mathcal{B}}[X_4;X_4]=\sqrt{1.1},\hfill & \hfill & C_{\mathcal{B}}[X_3;X_3]=\sqrt{0.9};\hfill \end{array}$$

and the elements below the diagonal are taken from the corresponding cells of

$$C_{\mathcal{B}}[X_4;•]=\left(\begin{array}{cc}1.5& 2.6\end{array}\right)\left(\begin{array}{cccc}0.894& 0& 0& 0\\ 0& 0.774& 0& 0\end{array}\right),$$

$$C_{\mathcal{B}}[X_3;•]=(1.2\left)\right(1.3412.0141.0490).$$

Computing $C_{\mathcal{B}}{C}_{\mathcal{B}}^{\mathrm{T}}$ gives

$${\tilde{\Sigma }}_{\mathcal{B}}=\begin{array}{c}\\ X_1\\ X_2\\ X_4\\ X_3\end{array}\begin{array}{cccccc}& X_1& X_2& X_4& X_3& \\ (& \begin{array}{c}0.800\\ 0\\ 1.200\\ 1.440\end{array}& \begin{array}{c}0\\ 0.600\\ 1.560\\ 1.872\end{array}& \begin{array}{c}1.200\\ 1.560\\ 6.956\\ 8.347\end{array}& \begin{array}{c}1.440\\ 1.872\\ 8.347\\ 10.916\end{array}& )\end{array}$$

and reordering the rows and columns of ${\tilde{\Sigma }}_{\mathcal{B}}$ gives

$${\Sigma }_{\mathcal{B}}=\begin{array}{c}\\ X_1\\ X_2\\ X_3\\ X_4\end{array}\begin{array}{cccccc}& X_1& X_2& X_3& X_4& \\ (& \begin{array}{c}0.800\\ 0\\ 1.440\\ 1.200\end{array}& \begin{array}{c}0\\ 0.600\\ 1.872\\ 1.560\end{array}& \begin{array}{c}1.440\\ 1.872\\ 10.916\\ 8.347\end{array}& \begin{array}{c}1.200\\ 1.560\\ 8.347\\ 6.956\end{array}& )\end{array}$$

The elements of the corresponding expectation vector ${\mathit{\mu }}_{\mathcal{B}}$ are then

$$\begin{array}{cc}\hfill E\left(X_1\right)& =2.400,\hfill \\ \hfill E\left(X_2\right)& =1.800,\hfill \\ \hfill E\left(X_4\right)& =0.2+1.5E\left(X_1\right)+2.6E\left(X_2\right)=8.480,\hfill \\ \hfill E\left(X_3\right)& =2.1+1.2E\left(X_4\right)=12.276.\hfill \end{array}$$

Starting from ${\Sigma }_{\mathcal{B}}$, we can reorder its rows and columns to obtain ${\tilde{\Sigma }}_{\mathcal{B}}$. The Cholesky decomposition of ${\tilde{\Sigma }}_{\mathcal{B}}$ is $C_{\mathcal{B}}$. Then

$$\begin{array}{cccc}\hfill & {\sigma }_{X_1}^2=C_{\mathcal{B}}{[X_1;X_1]}^2=0.8,\hfill & \hfill & {\sigma }_{X_2}^2=C_{\mathcal{B}}{[X_2;X_2]}^2=0.6,\hfill \\ \hfill & {\sigma }_{X_3}^2=C_{\mathcal{B}}{[X_3;X_3]}^2=0.9,\hfill & \hfill & {\sigma }_{X_4}^2=C_{\mathcal{B}}{[X_4;X_4]}^2=0.11.\hfill \end{array}$$

The coefficients ${\mathit{\beta }}_{X_i}$ of the local distributions are available from

$$\begin{array}{cc}\hfill R& ={\mathrm{I}}_N-\underset{diag\left(C_{\mathcal{B}}\right)}{\underbrace{\left[\begin{array}{cccc}0.894& 0& 0& 0\\ 0& 0.774& 0& 0\\ 0& 0& 1.049& 0\\ 0& 0& 0& 0.948\end{array}\right]}}\underset{C_{\mathcal{B}}^{-1}}{\underbrace{\left[\begin{array}{cccc}1.118& 0& 0& 0\\ 0& 1.291& 0& 0\\ -1.430& -2.479& 0.953& 0\\ 0& 0& -1.265& 1.054\end{array}\right]}}\hfill \\ \hfill & =\begin{array}{c}\\ X_1\\ X_2\\ X_4\\ X_3\end{array}\begin{array}{cccccc}& X_1& X_2& X_4& X_3& \\ (& \begin{array}{c}0\\ 0\\ 1.500\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 2.600\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\\ 1.200\end{array}& \begin{array}{c}0\\ 0\\ 0\\ 0\end{array}& )\end{array}\hfill \end{array}$$

where we can read $R_{X_4,X_1}=1.5={\beta }_{X_4,X_1}$, $R_{X_4,X_2}=2.6={\beta }_{X_4,X_2}$, $R_{X_3,X_4}=1.2={\beta }_{X_3,X_4}$.

We can read the standard errors of $X_1$, $X_2$, $X_3$ and $X_4$ directly from the diagonal elements of $C_{\mathcal{B}}$, and we can compute the intercepts from ${\mathit{\mu }}_{\mathcal{B}}-R{\mathit{\mu }}_{\mathcal{B}}$ which amounts to

$$\begin{array}{cc}\hfill {\mu }_{X_1}& =E\left(X_1\right)=2.400,\hfill \\ \hfill {\mu }_{X_2}& =E\left(X_2\right)=1.800,\hfill \\ \hfill {\mu }_{X_4}& =E\left(X_4\right)-E\left(X_1\right){\beta }_{X_4,X_1}-E\left(X_2\right){\beta }_{X_4,X_2}=0.200,\hfill \\ \hfill {\mu }_{X_3}& =E\left(X_3\right)-E\left(X_4\right){\beta }_{X_3,X_4}=2.100.\hfill \end{array}$$

3.3. Conditional Linear Gaussian BNs

Finally, conditional linear Gaussian BNs [37] (CLGBNs) subsume discrete BNs and GBNs as particular cases by combining discrete and continuous random variables in a mixture model. If we denote the former with ${\mathbf{X}}_D$ and the latter with ${\mathbf{X}}_G$, so that $\mathbf{X}={\mathbf{X}}_D\cup {\mathbf{X}}_G$, then:

• Discrete $X_i\in {\mathbf{X}}_D$ are only allowed to have discrete parents (denoted ${\Delta }_{X_i}$), and are assumed to follow a multinomial distribution parameterised with CPTs. We can estimate their parameters in the same way as those in a discrete BN.
• Continuous $X_i\in {\mathbf{X}}_G$ are allowed to have both discrete and continuous parents (denoted ${\Gamma }_{X_i}$, ${\Delta }_{X_i}\cup {\Gamma }_{X_i}={\Pi }_{X_i}$). Their local distributions are

  $$X_i\mid {\Pi }_{X_i}\sim N\left({\mu }_{X_i,{\delta }_{X_i}}+{\Gamma }_{X_i}{\mathit{\beta }}_{X_i,{\mathit{\delta }}_{X_i}},{\sigma }_{X_i,{\delta }_{X_i}}^2\right),$$

  which is equivalent to a mixture of linear regressions against the continuous parents with one component for each configuration ${\delta }_{X_i}\in \mathit{Val}\left({\Delta }_{X_i}\right)$ of the discrete parents:

  $$\begin{array}{cccc}\hfill & X_i={\mu }_{X_i,{\delta }_{X_i}}+{\Gamma }_{X_i}{\mathit{\beta }}_{X_i,{\delta }_{X_i}}+{\epsilon }_{X_i,{\delta }_{X_i}},\hfill & \hfill & {\epsilon }_{X_i,{\delta }_{X_i}}\sim N\left(0,{\sigma }_{X_i,{\delta }_{X_i}}^2\right).\hfill \end{array}$$
  If $X_i$ has no discrete parents, the mixture reverts to a single linear regression like that in (4). The parameters of these local distributions are usually estimated by maximum likelihood like those in a GBN; we have used hierarchical regressions with random effects in our recent work [38] for this purpose as well. Bayesian and regularised estimators are also an option [5].

If the CLGBN comprises $\left|{\mathbf{X}}_D\right|=M$ discrete nodes and $\left|{\mathbf{X}}_G\right|=N-M$ continuous nodes, these distributional assumptions imply the partial topological ordering

$$\underset{\mathrm{discrete}\mathrm{nodes}}{\underbrace{\left\{X_{\left(1\right)},\dots ,X_{\left(M\right)}\right\}}},\underset{\mathrm{continuous}\mathrm{nodes}}{\underbrace{\left\{X_{(M+1)},\dots ,X_{\left(N\right)}\right\}}}.$$

(5)

The discrete nodes jointly follow a multinomial distribution, effectively forming a discrete BN. The continuous nodes jointly follow a multivariate normal distribution, parameterised as a GBN, for each configuration of the discrete nodes. Therefore, the global distribution is a Gaussian mixture in which the discrete nodes identify the components, and the continuous nodes determine their distribution. The practical link between the global and local distributions follows directly from Section 3.1 and Section 3.2.

Example 3 

(Composing and decomposing a CLGBN). For reasons of space, this example is presented as Example A2 in Appendix B.

The complexity of composing and decomposing the global distribution is then

$$\underset{\mathrm{convert}\mathrm{between}\mathrm{CPTs}\mathrm{and}\mathrm{component}\mathrm{probabilities}}{\underbrace{O\left(M{l}^M\right).}}+\underset{\left(\mathrm{de}\right)\mathrm{compose}\mathrm{the}\mathrm{distinct}\mathrm{component}\mathrm{distributions}}{\underbrace{O\left({(N-M)}^3{l}^{\mathit{Val}\left({\mathbf{\Delta }}^{}\right)}\right)}}$$

where ${\mathbf{\Delta }}^{}={\bigcup }_{X_i\in {\mathbf{X}}_G}{\Delta }_{X_i}$ are the discrete parents of the continuous nodes.

3.4. Inference

For BNs, inference broadly denotes obtaining the conditional distribution of a subset of variables conditional on a second subset of variables. Following older terminology from expert systems [2], this is called formulating a query in which we ask the BN about the probability of an event of interest after observing some evidence. In conditional probability queries, the event of interest is the probability of one or more events in (or the whole distribution of) some variables of interest conditional on the values assumed by the evidence variables. In maximum a posteriori (“most probable explanation”) queries, we condition the values of the evidence variables to predict those of the event variables.

All inference computations on BNs are completely automated by exact and approximate algorithms, which we will briefly describe here. We refer the interested reader to the more detailed treatment in Castillo et al. [2] and Koller and Friedman [5].

Exact inference algorithms use local computations to compute the value of the query. The seminal works of Lauritzen and Spiegelhalter [39], Lauritzen and Wermuth [37] and Lauritzen and Jensen [40] describe how to transform a discrete BN or a (CL)GBN into a junction tree as a preliminary step before using belief propagation. Cowell [41] uses elimination trees for the same purpose in CLGBNs. (A junction tree is an undirected tree whose nodes are the cliques in the moral graph constructed from the BN and their intersections. A clique is the maximal subset of nodes such that every two nodes in the subset are adjacent).

Namasivayam et al. [42] give the computational complexity of constructing the junction tree from a discrete BN as $O(Nw+w{l}^{w}N)$ where w is the maximum number of nodes in a clique and, as before, l is the maximum number of values that a variable can take. We take the complexity of belief propagation to be $O(Nw{l}^w+|\Theta |)$, as stated in Lauritzen and Spiegelhalter [39] (“The global propagation is no worse than the initialisation [of the junction tree]”). This is confirmed by Pennock [43] and Namasivayam and Prasanna [44].

As for GBNs, we can also perform exact inference through their global distribution because the latter has only $O(N^2+N)$ parameters. The computational complexity of this approach is $O\left(N^3\right)$ because of the cost of composing the global distribution, which we derived in Section 3.2. However, all the operations involved are linear, making it possible to leverage specialised hardware such as GPUs and TPUs to the best effect. Koller and Friedman [5] (Section 14.2.1) note that “inference in linear Gaussian networks is linear in the number of cliques, and at most cubic in the size of the largest clique” when using junction trees and belief propagation. Therefore, junction trees may be significantly faster for GBNs when $w\ll N$. However, the correctness and convergence of belief propagation in GBNs require a set of sufficient conditions that have been studied comprehensively by Malioutov et al. [45]. Using the global distribution directly always produces correct results.

Approximate inference algorithms use Monte Carlo simulations to sample from the global distribution of $\mathbf{X}$ through the local distributions and estimate the answer queries by computing the appropriate summary statistics on the particles they generate. Therefore, they mirror the Monte Carlo and Markov chain Monte Carlo approaches in the literature: rejection sampling, importance sampling, and sequential Monte Carlo among others. Two state-of-the-art examples are the adaptive importance sampling (AIS-BN) scheme [46] and the evidence pre-propagation importance sampling (EPIS-BN) [47].

4. Shannon Entropy and Kullback–Leibler Divergence

The general definition of Shannon entropy for the probability distribution P of $\mathbf{X}$ is

$$H\left(P\right)=E_P(-logP\left(\mathbf{X}\right))=-{\int }_{\mathit{Val}\left(\mathbf{X}\right)}P\left(\mathbf{x}\right)logP\left(\mathbf{x}\right)d\mathbf{x}.$$

(6)

The Kullback–Leibler divergence between two distributions P and Q for the same random variables $\mathbf{X}$ is defined as

$$KL\left(P\parallel Q\right)=E_{P\left(\mathbf{X}\right)}\left(-log\frac{P\left(\mathbf{X}\right)}{Q\left(\mathbf{X}\right)}\right)=-{\int }_{\mathit{Val}\left(\mathbf{X}\right)}P\left(\mathbf{x}\right)log\frac{P\left(\mathbf{x}\right)}{Q\left(\mathbf{x}\right)}d\mathbf{x}.$$

(7)

They are linked as follows:

$$\underset{KL\left(P\left(\mathbf{X}\right)\parallel Q\left(\mathbf{X}\right)\right)}{\underbrace{E_{P\left(\mathbf{X}\right)}\left(-log\frac{P\left(\mathbf{X}\right)}{Q\left(\mathbf{X}\right)}\right)}}=\underset{H\left(P\left(\mathbf{X}\right)\right)}{\underbrace{E_{P\left(\mathbf{X}\right)}(-logP\left(\mathbf{X}\right))}}+\underset{H\left(P\left(\mathbf{X}\right),Q\left(\mathbf{X}\right)\right)}{\underbrace{E_{P\left(\mathbf{X}\right)}\left(logQ\left(\mathbf{X}\right)\right)}}$$

(8)

where $H\left(P\left(\mathbf{X}\right),Q\left(\mathbf{X}\right)\right)$ is the cross-entropy between $P\left(\mathbf{X}\right)$ and $Q\left(\mathbf{X}\right)$. For the many properties of these quantities, we refer the reader to Cover and Thomas [48] and Csiszár and Shields [49]. Their use and interpretation are covered in depth (and breadth!) in Murphy [19,20] for general machine learning and in Koller and Friedman [5] for BNs.

For a BN $\mathcal{B}$ encoding the probability distribution of $\mathbf{X}$, (6) decomposes into

$$H\left(\mathcal{B}\right)=\sum _{i=1}^{N}H\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}\right)$$

where ${\Pi }_{X_i}^{\mathcal{B}}$ are the parents of $X_i$ in $\mathcal{B}$. While this decomposition looks similar to (1), we see that its terms are not necessarily orthogonal, unlike the local distributions.

As for (7), we cannot simply write

$$KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)=\sum _{i=1}^{N}KL\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}\parallel X_i\mid {\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}\right)$$

because, in the general case, the nodes $X_i$ have different parents in $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$. This issue impacts the complexity of computing Kullback–Leibler divergences in different ways depending on the type of BN.

4.1. Discrete BNs

For discrete BNs, $H\left(\mathcal{B}\right)$ does not decompose into orthogonal components. As pointed out in Koller and Friedman [5] (Section 8.4.12),

$$\begin{array}{c}H\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}\right)=\sum _{j=1}^{q_i}P\left({\Pi }_{X_i}^{\mathcal{B}}=j\right)H\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}=j\right)\mathrm{where}\hfill \\ \hspace{1em}H\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}=j\right)=-\sum _{k=1}^{r_i}{\pi }_{ik\mid j}\left(\mathcal{B}\right)log{\pi }_{ik\mid j}\left(\mathcal{B}\right).\hfill \end{array}$$

(9)

If we estimated the conditional probabilities ${\pi }_{ik\mid j}\left(\mathcal{B}\right)$ from data, the $P\left({\Pi }_{X_i}^{\mathcal{B}}=j\right)$ are already available as the normalising constants of the individual conditional distributions $\{{\pi }_{ik\mid j}\left(\mathcal{B}\right),j=1,\dots ,q_i\}$ in the local distribution of $X_i$. In this case, the complexity of computing $H\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}\right)$ is linear in the number of parameters: $O\left(\right|\Theta \left|\right)={\sum }_{i=1}^{N}O\left(\right|{\Theta }_{X_i}\left|\right)$.

In the general case, we need exact inference to compute the probabilities $P\left({\Pi }_{X_i}^{\mathcal{B}}=j\right)$. Fortunately, they can be readily extracted from the junction tree derived from $\mathcal{B}$ as follows:

• Identify a clique containing both $X_i$ and ${\Pi }_{X_i}^{\mathcal{B}}$. Such a clique is guaranteed to exist by the family preservation property [5] (Definition 10.1).
• Compute the marginal distribution of ${\Pi }_{X_i}^{\mathcal{B}}$ by summing over the remaining variables in the clique.

Combining the computational complexity of constructing the junction tree from Section 3.4 and that of marginalisation, which is at most $O\left(l^{w-1}\right)$ for each node as in (2), we have

$$\begin{array}{c}\underset{\mathrm{create}\mathrm{the}\mathrm{junction}\mathrm{tree}}{\underbrace{O(Nw+w{l}^{w}N)}}+\underset{\mathrm{compute}\mathrm{the}P\left({\Pi }_{X_i}^{\mathcal{B}}=j\right)}{\underbrace{O\left(N{l}^{w-1}\right)}}+\underset{\mathrm{compute}H\left(\mathcal{B}\right)}{\underbrace{O\left(\right|\Theta \left|\right)}}=\hfill \\ \hspace{1em}O(N(w(1+l^w)+l^{w-1})+|\Theta |),\hfill \end{array}$$

which is exponential in the maximum clique size w. (The maximum clique size in a junction tree is proportional to the treewidth of the BN the junction tree is created from, which is also used in the literature to characterise computational complexity in BNs.) Interestingly, we do not need to perform belief propagation, so computing $H\left(\mathcal{B}\right)$ is more efficient than other inference tasks.

Example 4 

(Entropy of a discrete BN). For reasons of space, this example is presented as Example A3 in Appendix B.

The Kullback–Leibler divergence has a similar issue, as noted in Koller and Friedman [5] (Section 8.4.2). The best and most complete explanation of how to compute it for discrete BNs is in Moral et al. [21]. After decomposing $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$ following (8) to separate $H\left(\mathcal{B}\right)$ and $H\left(\mathcal{B},{\mathcal{B}}^{\prime }\right)$, Moral et al. [21] show that the latter takes the form

$$H\left(\mathcal{B},{\mathcal{B}}^{\prime }\right)=\sum _{i=1}^N\sum _{j\in \mathit{Val}\left({\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}\right)}\left[\sum _{k=1}^{r_i}{\pi }_{ikj}\left(\mathcal{B}\right)log{\pi }_{ik\mid j}\left({\mathcal{B}}^{\prime }\right)\right]$$

(10)

where:

• ${\pi }_{ikj}\left(\mathcal{B}\right)=P\left(X_i=k,{\Pi }_{X_i}\left({\mathcal{B}}^{\prime }\right)=j\right)$ is the probability assigned by $\mathcal{B}$ to $X_i=k$ given that the variables that are parents of $X_i$ in ${\mathcal{B}}^{\prime }$ take value j;
• ${\pi }_{ik\mid j}\left({\mathcal{B}}^{\prime }\right)=P\left(X_i=k\mid {\Pi }_{X_i}\left({\mathcal{B}}^{\prime }\right)=j\right)$ is the $(k,j)$ element of the CPT of $X_i$ in ${\mathcal{B}}^{\prime }$.

In order to compute the ${\pi }_{ikj}\left(\mathcal{B}\right)$, we need to transform $\mathcal{B}$ into its junction tree and use belief propagation to compute the joint distribution of $X_i\cup {\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}$. As a result, $H\left(\mathcal{B},{\mathcal{B}}^{\prime }\right)$ does not decompose at all: each ${\pi }_{ikj}\left(\mathcal{B}\right)$ can potentially depend on the whole BN $\mathcal{B}$.

Algorithmically, to compute $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$ we:

• Transform $\mathcal{B}$ into its junction tree.
• Compute the entropy $H\left(\mathcal{B}\right)$.
• For each node $X_i$:

  (a)

  Identify ${\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}$, the parents of $X_i$ in ${\mathcal{B}}^{\prime }$.

  (b)

  Obtain the distribution of the variables $\{X_i,{\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}\}$ from the junction tree of $\mathcal{B}$, consisting of the probabilities ${\pi }_{ikj}\left(\mathcal{B}\right)$.

  (c)

  Read the ${\pi }_{ik\mid j}\left({\mathcal{B}}^{\prime }\right)$ from the local distribution of $X_i$ in ${\mathcal{B}}^{\prime }$.

• Use the ${\pi }_{ikj}\left(\mathcal{B}\right)$ and the ${\pi }_{ik\mid j}\left({\mathcal{B}}^{\prime }\right)$ to compute (10).

The computational complexity of this procedure is as follows:

$$\begin{array}{c}\underset{\mathrm{create}\mathrm{the}\mathrm{junction}\mathrm{tree}\mathrm{of}B\mathrm{and}\mathrm{computing}H\left(\mathcal{B}\right)}{\underbrace{O(N(w(1+l^w)+l^{w-1})+|\Theta |)}}+\underset{\mathrm{produce}\mathrm{the}{\pi }_{ikj}\left(\mathcal{B}\right)}{\underbrace{O\left(N{l}^c(Nw{l}^w+|\Theta |)\right)}}+\underset{\mathrm{compute}H\left(\mathcal{B},{\mathcal{B}}^{\prime }\right)}{\underbrace{O\left(\right|\Theta \left|\right)}}=\hfill \\ \hfill O(N^{2}w{l}^{w+c}+N(w+w{l}^w+l^{w-1})+(N{l}^c+2)|\Theta |).\end{array}$$

(11)

As noted in Moral et al. [21], computing the ${\pi }_{ikj}\left(\mathcal{B}\right)$ requires a separate run of belief propagation for each configuration of the ${\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}$, for a total of ${\sum }_{i=1}^N{l}^{\left|{\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}\right|}$ times. If we assume that the DAG underlying ${\mathcal{B}}^{\prime }$ is sparse, we have that $\left|{\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}\right|⩽c$ and the overall complexity of this step becomes $O(N{l}^c·(Nw{l}^w+|\Theta |))$, N times that listed in Section 3.4. The caching scheme devised by Moral et al. [21] is very effective in limiting the use of belief propagation, but it does not alter its exponential complexity.

Example 5 

(KL between two discrete BNs). Consider the discrete BN $\mathcal{B}$ from Figure 2 top. Furthermore, consider the BN ${\mathcal{B}}^{\prime }$ from Figure 2 bottom. We constructed the global distribution of $\mathcal{B}$ in Example A1; we can similarly compose the global distribution of ${\mathcal{B}}^{\prime }$, shown below.

$X_1=a$

$X_1=b$

$X_2=c$

$X_2=d$

$X_2=c$

$X_2=d$

$X_3$

$X_3$

$X_3$

$X_3$

$X_4$

e

f

$X_4$

e

f

$X_4$

e

f

$X_4$

e

f

g

$0.013$

$0.033$

g

$0.022$

$0.054$

g

$0.016$

$0.144$

g

$0.016$

$0.139$

h

$0.072$

$0.062$

h

$0.029$

$0.025$

h

$0.013$

$0.040$

h

$0.079$

$0.243$

Since both global distributions are limited in size, we can then compute the Kullback–Leibler divergence between $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$ using (7).

$$\begin{array}{c}KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)=-0.013log0.013-0.016log0.016-0.022log0.022-0.016log0.016-\hfill \\ 0.072log0.072-0.013log0.013-0.029log0.029-0.079log0.079-\\ 0.033log0.033-0.144log0.144-0.054log0.054-0.139log0.139-\\ \hfill -0.062log0.062-0.04log0.04-0.025log0.025-0.243log0.243=0.687\end{array}$$

In the general case, when we cannot use the global distributions, we follow the approach described in Section 4.1. Firstly, we apply (8) to write

$$KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)=H\left(\mathcal{B}\right)-H\left(\mathcal{B},{\mathcal{B}}^{\prime }\right);$$

we have from Example A3 that $H\left(\mathcal{B}\right)=2.440$. As for the cross-entropy $H\left(\mathcal{B},{\mathcal{B}}^{\prime }\right)$, we apply (10):

1. 

We identify the parents of each node in ${\mathcal{B}}^{\prime }$:

$$\begin{array}{cccccccc}\hfill & {\Pi }_{X_1}^{{\mathcal{B}}^{\prime }}=\{⌀\},\hfill & \hfill & {\Pi }_{X_2}^{{\mathcal{B}}^{\prime }}=\{X_1,X_4\},\hfill & \hfill & {\Pi }_{X_3}^{{\mathcal{B}}^{\prime }}=\left\{X_1\right\},\hfill & \hfill & {\Pi }_{X_4}^{{\mathcal{B}}^{\prime }}=\left\{X_3\right\}.\hfill \end{array}$$

2. 

We construct a junction tree from $\mathcal{B}$ and we use it to compute the distributions $P\left(X_1\right)$, $P\left(X_2,X_1,X_4\right)$, $P\left(X_3,X_1\right)$ and $P\left(X_4,X_3\right)$.

$X_1$ a b $0.53$ $0.47$

$\{X_1,X_4\}$ $\{a,g\}$ $\{a,h\}$ $\{b,g\}$ $\{b,h\}$ $X_2$ c $0.070$ $0.110$ $0.053$ $0.107$ d $0.089$ $0.261$ $0.076$ $0.235$

$X_1$ a b $X_3$ e $0.289$ $0.312$ f $0.241$ $0.158$

$X_3$ e f $X_4$ g $0.120$ $0.167$ h $0.481$ $0.231$

3. 

We compute the cross-entropy terms for the individual variables in $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$:

$$\begin{array}{cc}\hfill H\left(X_1^{\mathcal{B}},X_1^{{\mathcal{B}}^{\prime }}\right)& =0.53log0.31+0.47log0.69=-0.795;\hfill \\ \hfill H\left(X_2^{\mathcal{B}},X_2^{{\mathcal{B}}^{\prime }}\right)& =0.070log0.38+0.089log0.62+0.110log0.71+0.261log0.29+\hfill \\ \hfill & 0.053log0.51+0.076log0.49+0.107log0.14+0.235log0.86\hfill \\ \hfill & =-0.807;\hfill \\ \hfill H\left(X_3^{\mathcal{B}},X_3^{{\mathcal{B}}^{\prime }}\right)& =0.289log0.44+0.241log0.56+0.312log0.18+0.158log0.82\hfill \\ \hfill & =-0.943;\hfill \\ \hfill H\left(X_4^{\mathcal{B}},X_4^{{\mathcal{B}}^{\prime }}\right)& =0.120log0.26+0.481log0.74+0.167log0.50+0.231log0.50\hfill \\ \hfill & =-0.582;\hfill \end{array}$$

which sum up to $H\left(\mathcal{B},{\mathcal{B}}^{\prime }\right)={\sum }_{i=1}^{N}H\left(X_i^{\mathcal{B}},X_i^{{\mathcal{B}}^{\prime }}\right)=-3.127$.

4. 

We compute $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)=2.440-3.127=0.687$, which matches the value we previously computed from the global distributions.

4.2. Gaussian BNs

$H\left(\mathcal{B}\right)$ decomposes along with the local distributions $X_i\mid {\Pi }_{X_i}$ in the case of GBNs: from (3), each $X_i\mid {\Pi }_{X_i}$ is a univariate normal with variance ${\sigma }_{X_i}^2\left(\mathcal{B}\right)$ and therefore

$$H\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}\right)=\frac{1}{2}log\left(2\pi {\sigma }_{X_i}^2\left(\mathcal{B}\right)\right)+\frac{1}{2}$$

(12)

which has a computational complexity of O(1) for each node, O(N) overall. Equivalently, we can start from the global distribution of $\mathcal{B}$ from Section 3.2 and consider that

$$det(\Sigma )=det\left(C_{{\mathcal{B}}^{\prime }}{C}_{{\mathcal{B}}^{\prime }}^{\mathrm{T}}\right)=det{\left(C_{\mathcal{B}}\right)}^2={\left({\prod }_{i=1}^N{C}_{\mathcal{B}}[i;i]\right)}^2={\prod }_{i=1}^N{\sigma }_{X_i}^2\left(\mathcal{B}\right)$$

(13)

because $C_{\mathcal{B}}$ is lower triangular. The (multivariate normal) entropy of $\mathbf{X}$ then becomes

$$\begin{array}{c}H\left(\mathcal{B}\right)=\frac{N}{2}+\frac{N}{2}log2\pi +\frac{1}{2}logdet(\Sigma )=\frac{N}{2}+\frac{N}{2}log2\pi +\frac{1}{2}\sum _{i=1}^{N}log{\sigma }_{X_i}^2\left(\mathcal{B}\right)\hfill \\ =\sum _{i=1}^N\frac{1}{2}+\frac{1}{2}log\left(2\pi {\sigma }_{X_i}^2\left(\mathcal{B}\right)\right)=\sum _{i=1}^{N}H\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}\right)\hfill \end{array}$$

in agreement with (12).

Example 6 

(Entropy of a GBN). For reasons of space, this example is presented as Example A4 in Appendix B.

In the literature, the Kullback–Leibler divergence between two GBNs $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$ is usually computed using the respective global distributions $N({\mathit{\mu }}_{\mathcal{B}},{\Sigma }_{\mathcal{B}})$ and $N({\mathit{\mu }}_{{\mathcal{B}}^{\prime }},{\Sigma }_{{\mathcal{B}}^{\prime }})$ [50,51,52]. The general expression is

$$KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)=\frac{1}{2}\left[tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)+{({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})}^{\mathrm{T}}{\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})-N+log\frac{det\left({\Sigma }_{{\mathcal{B}}^{\prime }}\right)}{det\left({\Sigma }_{\mathcal{B}}\right)}\right],$$

(14)

which has computational complexity

$$\begin{array}{c}\underset{\mathrm{compute}{\mathit{\mu }}_{\mathcal{B}},{\mathit{\mu }}_{{\mathcal{B}}^{\prime }}{\Sigma }_{\mathcal{B}},{\Sigma }_{{\mathcal{B}}^{\prime }}}{\underbrace{O(2N^3+2N)}}+\underset{\mathrm{invert}{\Sigma }_{{\mathcal{B}}^{\prime }}}{\underbrace{O\left(N^3\right)}}+\underset{\mathrm{multiply}{\Sigma }_{\mathcal{B}}^{-1}\mathrm{and}{\Sigma }_{{\mathcal{B}}^{\prime }}}{\underbrace{O\left(N^3\right)}}+\underset{\mathrm{trace}\mathrm{of}{\Sigma }_{\mathcal{B}}^{-1}{\Sigma }_{{\mathcal{B}}^{\prime }}}{\underbrace{O\left(N\right)}}+\hfill \\ \underset{\mathrm{compute}{({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})}^{\mathrm{T}}{\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})}{\underbrace{O(N^2+2N)}}+\underset{\mathrm{determinant}\mathrm{of}{\Sigma }_{{\mathcal{B}}^{\prime }}}{\underbrace{O\left(N^3\right)}}+\underset{\mathrm{determinant}\mathrm{of}{\Sigma }_{\mathcal{B}}}{\underbrace{O\left(N^3\right)}}=\\ \hfill O(6N^3+N^2+5N).\end{array}$$

(15)

The spectral decomposition ${\Sigma }_{{\mathcal{B}}^{\prime }}=U{\Lambda }_{{\mathcal{B}}^{\prime }}{U}^{\mathrm{T}}$ gives the eigenvalues $diag\left({\Lambda }_{{\mathcal{B}}^{\prime }}\right)=\{{\lambda }_1\left({\mathcal{B}}^{\prime }\right),\dots ,$${\lambda }_N\left({\mathcal{B}}^{\prime }\right)\}$ to compute ${\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}$ and $det\left({\Sigma }_{{\mathcal{B}}^{\prime }}\right)$ efficiently as illustrated in the example below. (Further computing the spectral decomposition of ${\Sigma }_{\mathcal{B}}$ to compute $det\left({\Sigma }_{\mathcal{B}}\right)$ from the eigenvalues $\{{\lambda }_1\left(\mathcal{B}\right),\dots ,{\lambda }_N\left(\mathcal{B}\right)\}$ does not improve complexity because it just replaces a single $O\left(N^3\right)$ operation with another one.) We thus somewhat improve the overall complexity of $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$ to $O(5N^3+N^2+6N)$.

Example 7 

(General-case KL between two GBNs). Consider the GBN $\mathcal{B}$ Figure 1 top, which we know has global distribution

$$\begin{array}{c}\hfill \left[\begin{array}{c}{X}_1\\ X_2\\ X_3\\ X_4\end{array}\right]\sim N\left(\left[\begin{array}{c}2.400\\ 1.800\\ 12.276\\ 8.848\end{array}\right],\left[\begin{array}{cccc}0.800& 0& 1.440& 1.200\\ 0& 0.600& 1.872& 1.560\\ 1.440& 1.872& 10.916& 8.347\\ 1.200& 1.560& 8.347& 6.956\end{array}\right]\right)\end{array}$$

from Example 2. Furthermore, consider the GBN ${\mathcal{B}}^{\prime }$ from Figure 1 bottom, which has global distribution

$$\begin{array}{c}\hfill \left[\begin{array}{c}{X}_1\\ X_2\\ X_3\\ X_4\end{array}\right]\sim N\left(\left[\begin{array}{c}2.400\\ 11.324\\ 6.220\\ 4.620\end{array}\right],\left[\begin{array}{cccc}0.800& 2.368& 1.040& 0.640\\ 2.368& 8.541& 3.438& 1.894\\ 1.040& 3.438& 1.652& 0.832\\ 0.640& 1.894& 0.832& 1.012\end{array}\right]\right).\end{array}$$

In order to compute $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$, we first invert ${\Sigma }_{{\mathcal{B}}^{\prime }}$ to obtain

$${\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}=\left[\begin{array}{cccc}9.945& -1.272& -2.806& -1.600\\ -1.272& 0.909& -1.091& 0\\ -2.806& -1.091& 4.642& 0\\ -1.600& 0& 0& 2.000\end{array}\right],$$

which we then multiply by ${\Sigma }_{\mathcal{B}}$ to compute the trace $tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)=57.087$. We also use ${\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}$ to compute ${({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})}^{\mathrm{T}}{\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})=408.362$. Finally, $det\left({\Sigma }_{{\mathcal{B}}^{\prime }}\right)=0.475$, $det\left({\Sigma }_{\mathcal{B}}\right)=0.132$ and therefore

$$KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)=\frac{1}{2}\left[57.087+408.362-4+log\left(\frac{0.475}{0.132}\right)\right]=230.0846.$$

(16)

As an alternative, we can compute the spectral decompositions ${\Sigma }_{\mathcal{B}}=U_{\mathcal{B}}{\Lambda }_{\mathcal{B}}{U}_{\mathcal{B}}^T$ and ${\Sigma }_{{\mathcal{B}}^{\prime }}=U_{{\mathcal{B}}^{\prime }}{\Lambda }_{{\mathcal{B}}^{\prime }}{U}_{{\mathcal{B}}^{\prime }}^T$ as an intermediate step. Multiplying the sets of eigenvalues

$$\begin{array}{cccccc}\hfill & {\Lambda }_{\mathcal{B}}=diag\left(\{18.058,0.741,0.379,0.093\}\right)\hfill & \hfill & \mathrm{and}\hfill & \hfill & {\Lambda }_{{\mathcal{B}}^{\prime }}=diag\left(\{11.106,0.574,0.236,0.087\}\right)\hfill \end{array}$$

gives the corresponding determinants; and it allows us to easily compute

$$\begin{array}{cccccc}\hfill & {\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}=U_{{\mathcal{B}}^{\prime }}{\Lambda }_{{\mathcal{B}}^{\prime }}^{-1}{U}_{{\mathcal{B}}^{\prime }}^T,\hfill & \hfill & \mathit{where}\hfill & \hfill & {\Lambda }_{{\mathcal{B}}^{\prime }}^{-1}=diag\left(\left\{\frac{1}{11.106},\frac{1}{0.574},\frac{1}{0.236},\frac{1}{0.087}\right\}\right)\hfill \end{array}$$

for use in both the quadratic form and in the trace.

However, computing $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$ from the global distributions $N({\mathit{\mu }}_{\mathcal{B}},{\Sigma }_{\mathcal{B}})$ and $N({\mathit{\mu }}_{{\mathcal{B}}^{\prime }},{\Sigma }_{{\mathcal{B}}^{\prime }})$ disregards the fact that BNs are sparse models that can be characterised more compactly by $({\mathit{\mu }}_{\mathcal{B}},C_{\mathcal{B}})$ and $({\mathit{\mu }}_{{\mathcal{B}}^{\prime }},C_{{\mathcal{B}}^{\prime }})$ as shown in Section 3.2. In particular, we can revisit several operations that are in the high-order terms of (15):

• Composing the global distribution from the local ones. We avoid computing ${\Sigma }_{\mathcal{B}}$ and ${\Sigma }_{{\mathcal{B}}^{\prime }}$, thus reducing this step to $O\left(2N\right)$ complexity.
• Computing the trace $tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)$. We can reduce the computation of the trace as follows.
  • We can replace ${\Sigma }_{\mathcal{B}}$ and ${\Sigma }_{{\mathcal{B}}^{\prime }}$ in the trace with any reordered matrix [53] (Result 8.17): we choose to use ${\tilde{\Sigma }}_{{\mathcal{B}}^{\prime }}$ and ${\tilde{\Sigma }}_{\mathcal{B}}^{\ast }$ where ${\tilde{\Sigma }}_{{\mathcal{B}}^{\prime }}$ is defined as before and ${\tilde{\Sigma }}_{\mathcal{B}}^{\ast }$ is ${\Sigma }_{\mathcal{B}}$ with the rows and columns reordered to match ${\tilde{\Sigma }}_{{\mathcal{B}}^{\prime }}$. Formally, this is equivalent to ${\tilde{\Sigma }}_{\mathcal{B}}^{\ast }=P{\tilde{\Sigma }}_{\mathcal{B}}{P}^{\mathrm{T}}$ where P is a permutation matrix that imposes the desired node ordering: since both the rows and the columns are permuted in the same way, the diagonal elements of ${\tilde{\Sigma }}_{\mathcal{B}}$ are the same as those of ${\tilde{\Sigma }}_{\mathcal{B}}^{\ast }$ and the trace is unaffected.
  • We have ${\tilde{\Sigma }}_{{\mathcal{B}}^{\prime }}=C_{{\mathcal{B}}^{\prime }}{C}_{{\mathcal{B}}^{\prime }}^{\mathrm{T}}$.
  • As for ${\tilde{\Sigma }}_{\mathcal{B}}^{\ast }$, we can write ${\tilde{\Sigma }}_{\mathcal{B}}^{\ast }=P{\tilde{\Sigma }}_{\mathcal{B}}P=\left(P{C}_{\mathcal{B}}\right){\left(P{C}_{\mathcal{B}}\right)}^{\mathrm{T}}=C_{\mathcal{B}}^{\ast }{\left(C_{\mathcal{B}}^{\ast }\right)}^{\mathrm{T}}$ where $C_{\mathcal{B}}^{\ast }=P{C}_{\mathcal{B}}$ is the lower triangular matrix $C_{\mathcal{B}}$ with the rows re-ordered to match ${\tilde{\Sigma }}_{{\mathcal{B}}^{\prime }}$. Note that $C_{\mathcal{B}}^{\ast }$ is not lower triangular unless $\mathcal{G}$ and ${\mathcal{G}}^{\prime }$ have the same partial node ordering, which implies $P={\mathrm{I}}_N$.
  Therefore

  $$\begin{array}{c}\hfill tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)=tr\left({\left(C_{{\mathcal{B}}^{\prime }}^{-1}{C}_{\mathcal{B}}^{\ast }\right)}^{\mathrm{T}}\left(C_{{\mathcal{B}}^{\prime }}^{-1}{C}_{\mathcal{B}}^{\ast }\right)\right)={\parallel C_{{\mathcal{B}}^{\prime }}^{-1}{C}_{\mathcal{B}}^{\ast }\parallel }_F^2\end{array}$$

  (17)

  where the last step rests on Seber [53] (Result 4.15). We can invert $C_{{\mathcal{B}}^{\prime }}$ in $O\left(N^2\right)$ time following Stewart [54] (Algorithm 2.3). Multiplying $C_{{\mathcal{B}}^{\prime }}^{-1}$ and $C_{\mathcal{B}}^{\ast }$ is still $O\left(N^3\right)$. The Frobenius norm ${\parallel ·\parallel }_F$ is $O\left(N^2\right)$ since it is the sum of the squared elements of $C_{{\mathcal{B}}^{\prime }}^{-1}{C}_{\mathcal{B}}^{\ast }$.
• Computing the determinants $det\left({\Sigma }_{{\mathcal{B}}^{\prime }}\right)$ and $det\left({\Sigma }_{\mathcal{B}}\right)$. From (13), each determinant can be computed in $O\left(N\right)$.
• Computing the quadratic term ${({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})}^{\mathrm{T}}{\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})$. Decomposing ${\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}$ leads to

  $${({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})}^{\mathrm{T}}{\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})={\left(C_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}^{\ast }-{\mathit{\mu }}_{\mathcal{B}}^{\ast })\right)}^{\mathrm{T}}{C}_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}^{\ast }-{\mathit{\mu }}_{\mathcal{B}}^{\ast }),$$

  (18)

  where ${\mathit{\mu }}_{{\mathcal{B}}^{\prime }}^{\ast }$ and ${\mathit{\mu }}_{\mathcal{B}}^{\ast }$ are the mean vectors re-ordered to match $C_{{\mathcal{B}}^{\prime }}^{-1}$. The computational complexity is still $O(N^2+2N)$ because $C_{{\mathcal{B}}^{\prime }}^{-1}$ is available from previous computations.

Combining (17), (13) and (18), the expression in (14) becomes

$$\begin{array}{c}KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)=\hfill \\ \hspace{1em}\frac{1}{2}\left[\parallel C_{{\mathcal{B}}^{\prime }}^{-1}{C}_{\mathcal{B}}^{\ast }{\parallel }_F^2+{\left(C_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}^{\ast }-{\mathit{\mu }}_{\mathcal{B}}^{\ast })\right)}^{\mathrm{T}}{C}_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}^{\ast }-{\mathit{\mu }}_{\mathcal{B}}^{\ast })-N+2log\frac{{\prod }_{i=1}^N{C}_{{\mathcal{B}}^{\prime }}[i;i]}{{\prod }_{i=1}^N{C}_{\mathcal{B}}[i;i]}\right].\hfill \end{array}$$

(19)

The overall complexity of (19) KL is

$$\begin{array}{c}\underset{\mathrm{compute}{\mathit{\mu }}_{\mathcal{B}},{\mathit{\mu }}_{{\mathcal{B}}^{\prime }}{C}_{\mathcal{B}},C_{{\mathcal{B}}^{\prime }}}{\underbrace{O(2N^2+2N)}}+\underset{\mathrm{compute}\parallel C_{{\mathcal{B}}^{\prime }}^{-1}{C}_{\mathcal{B}}{\parallel }_F^2}{\underbrace{O(2N^2+N^3)}}+\underset{\mathrm{compute}\mathrm{the}\mathrm{quadratic}\mathrm{form}}{\underbrace{O(N^2+2N)}}+\hfill \\ \underset{\mathrm{compute}det\left({\Sigma }_{\mathcal{B}}\right),det\left({\Sigma }_{{\mathcal{B}}^{\prime }}\right)}{\underbrace{O\left(2N\right)}}=O(N^3+5N^2+6N);\hfill \end{array}$$

(20)

while still cubic, the leading coefficient suggests that it should be about 5 times faster than the variant of (15) using the spectral decomposition.

Example 8 

(Sparse KL between two GBNs). Consider again the two GBNs from Example 7. The corresponding matrices

$$\begin{array}{cc}{C}_{\mathcal{B}}=\begin{array}{c}\\ X_1\\ X_2\\ X_4\\ X_3\end{array}\begin{array}{cccccc}& X_1& X_2& X_4& X_3& \\ (& \begin{array}{c}0.894\\ 0\\ 1.341\\ 1.610\end{array}& \begin{array}{c}0\\ 0.774\\ 2.014\\ 2.416\end{array}& \begin{array}{c}0\\ 0\\ 1.049\\ 1.258\end{array}& \begin{array}{c}0\\ 0\\ 0\\ 0.948\end{array}& )\end{array}& ,\hspace{1em}{C}_{\mathcal{B}\prime }=\begin{array}{c}\\ X_1\\ X_3\\ X_4\\ X_2\end{array}\begin{array}{cccccc}& X_1& X_3& X_4& X_2& \\ (& \begin{array}{c}0.894\\ 1.163\\ 0.715\\ 2.647\end{array}& \begin{array}{c}0\\ 0.548\\ 0\\ 0.657\end{array}& \begin{array}{c}0\\ 0\\ 0.707\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\\ 1.049\end{array}& )\end{array}\end{array}$$

readily give the determinants of ${\Sigma }_{\mathcal{B}}$ and ${\Sigma }_{{\mathcal{B}}^{\prime }}$ following (13):

$$\begin{array}{cc}\hfill det\left(C_{\mathcal{B}}\right)& ={(0.894·0.774·1.049·0.948)}^2=0.475,\hfill \\ \hfill det\left(C_{{\mathcal{B}}^{\prime }}\right)& ={(0.894·0.548·0.707·1.049)}^2=0.132.\hfill \end{array}$$

As for the Frobenius norm in (17), we first invert $C_{{\mathcal{B}}^{\prime }}$ to obtain

$$C_{{\mathcal{B}}^{\prime }}^{-1}=\begin{array}{c}\\ X_1\\ X_3\\ X_4\\ X_2\end{array}\begin{array}{cccccc}& X_1& X_3& X_4& X_2& \\ (& \begin{array}{c}1.118\\ \mathrm{-2.373}\\ \mathrm{-1.131}\\ \mathrm{-1.334}\end{array}& \begin{array}{c}0\\ 1.825\\ 0\\ \mathrm{-1.144}\end{array}& \begin{array}{c}0\\ 0\\ 1.414\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\\ 0.953\end{array}& )\end{array};$$

then we reorder the rows and columns of $C_{\mathcal{B}}$ to follow the same node ordering as $C_{{\mathcal{B}}^{\prime }}$ and compute

$${∥\left(\begin{array}{cccc}1.118& 0& 0& 0\\ -2.373& 1.825& 0& 0\\ -1.131& 0& 1.414& 0\\ -1.334& -1.144& 0& 0.953\end{array}\right)\left(\begin{array}{cccc}0.894& 0& 0& 0\\ 1.610& 0.948& 1.258& 2.416\\ 1.341& 0& 1.049& 2.014\\ 0& 0& 0& 0.774\end{array}\right)∥}_F^2=57.087$$

which, as expected, matches the value of $tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)$ we computed in Example 7. Finally, $C_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}^{\ast }-{\mathit{\mu }}_{\mathcal{B}}^{\ast })$ in (18) is

$$\left(\begin{array}{cccc}1.118& 0& 0& 0\\ -2.373& 1.825& 0& 0\\ -1.131& 0& 1.414& 0\\ -1.334& -1.144& 0& 0.953\end{array}\right)\left[\left(\begin{array}{c}2.400\\ 6.220\\ 4.620\\ 11.324\end{array}\right)-\left(\begin{array}{c}2.400\\ 12.1276\\ 8.848\\ 1.800\end{array}\right)\right]=\left(\begin{array}{c}0\\ -11.056\\ -5.459\\ 16.010\end{array}\right).$$

The quadratic form is then equal to $408.362$, which matches the value of ${({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})}^{\mathrm{T}}{\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}({\mathit{\mu }}_{{\mathcal{B}}^{\prime }}-{\mathit{\mu }}_{\mathcal{B}})$ in Example 7. As a result, the expression for $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$ is the same as in (16).

We can further reduce the complexity (20) of (19) when an approximate value of KL is suitable for our purposes. The only term with cubic complexity is $tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)={\parallel C_{{\mathcal{B}}^{\prime }}^{-1}{C}_{\mathcal{B}}^{\ast }\parallel }_F^2$: reducing it to quadratic complexity or lower will eliminate the leading term of (20), making it quadratic in complexity. One way to do this is to compute a lower and an upper bound for $tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)$, which can serve as an interval estimate, and take their geometric mean as an approximate point estimate.

A lower bound is given by Seber [53] (Result 10.39):

$$tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)⩾logdet\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)+N=-logdet\left({\Sigma }_{{\mathcal{B}}^{\prime }}\right)+logdet\left({\Sigma }_{\mathcal{B}}\right)+N,$$

(21)

which conveniently reuses the values of $det\left({\Sigma }_{\mathcal{B}}\right)$ and $det\left({\Sigma }_{{\mathcal{B}}^{\prime }}\right)$ we have from (13). For an upper bound, Seber [53] (Result 10.59) combined with Seber [53] (Result 4.15) gives

$$tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)⩽tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}\right)tr\left({\Sigma }_{\mathcal{B}}\right)=tr\left({\left(C_{{\mathcal{B}}^{\prime }}{C}_{{\mathcal{B}}^{\prime }}^{\mathrm{T}}\right)}^{-1}\right)tr\left(C_{\mathcal{B}}{C}_{\mathcal{B}}^{\mathrm{T}}\right)=\parallel C_{{\mathcal{B}}^{\prime }}^{-1}{\parallel }_F^2{\parallel C_{\mathcal{B}}\parallel }_F^2,$$

(22)

a function of $C_{\mathcal{B}}$ and $C_{{\mathcal{B}}^{\prime }}$ that can be computed in $O\left(2N^2\right)$ time. Note that, as far as the point estimate is concerned, we do not care about how wide the interval is: we only need its geometric mean to be an acceptable approximation of $tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)$.

Example 9 

(Approximate KL). From Example 7, we have that $tr\left({\Sigma }_{{\mathcal{B}}^{\prime }}^{-1}{\Sigma }_{\mathcal{B}}\right)=57.087$, $det\left({\Sigma }_{{\mathcal{B}}^{\prime }}\right)=0.475$ and $det\left({\Sigma }_{\mathcal{B}}\right)=0.132$. The lower bound in (21) is then

$$-logdet\left({\Sigma }_{{\mathcal{B}}^{\prime }}\right)+logdet\left({\Sigma }_{\mathcal{B}}\right)+4=5.281$$

and the upper bound in (22) is

$$\parallel C_{{\mathcal{B}}^{\prime }}^{-1}{\parallel }_F^2{\parallel C_{\mathcal{B}}\parallel }_F^2=17.496·19.272=337.207.$$

Their geometric mean is $42.199$, which can serve as an approximate value for $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$.

If we are comparing two GBNs whose parameters (but not necessarily network structures) have been learned from the same data, we can sometimes approximate $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$ using the local distributions $X_i\mid {\Pi }_{X_i}^{\mathcal{B}}$ and $X_i\mid {\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}$ directly. If $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$ have compatible partial orderings, we can define a common total node ordering for both such that

$$\begin{array}{cc}\hfill KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)& =KL\left(X_{\left(1\right)}\mid \{X_{\left(2\right)},\dots ,X_{\left(N\right)}\}\cdots X_N\parallel X_{\left(1\right)}\mid \{X_{\left(2\right)},\dots ,X_{\left(N\right)}\}\cdots X_N\right)\hfill \\ \hfill & =KL\left(X_{\left(1\right)}\mid {\Pi }_{X_{\left(1\right)}}^{\mathcal{B}}·\dots ·X_{\left(N\right)}\mid {\Pi }_{X_{\left(N\right)}}^{\mathcal{B}}\parallel X_{\left(1\right)}\mid {\Pi }_{X_{\left(1\right)}}^{{\mathcal{B}}^{\prime }}·\dots ·X_{\left(N\right)}\mid {\Pi }_{X_{\left(N\right)}}^{{\mathcal{B}}^{\prime }}\right).\hfill \end{array}$$

By “compatible partial orderings”, we mean two partial orderings that can be sorted into at least one shared total node ordering that is compatible with both. The product of the local distributions in the second step is obtained from the chain decomposition in the first step by considering the nodes in the conditioning other than the parents to have associated regression coefficients equal to zero. Then, following the derivations in Cavanaugh [55] for a general linear regression model, we can write the empirical approximation

$$KL\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}\parallel X_i\mid {\Pi }_{X_i}^{{\mathcal{B}}^{\prime }}\right)\approx \frac{1}{2}\left(log\frac{{\widehat{\sigma }}_{X_i}^2\left({\mathcal{B}}^{\prime }\right)}{{\widehat{\sigma }}_{X_i}^2\left(\mathcal{B}\right)}+\frac{{\widehat{\sigma }}_{X_i}^2\left(\mathcal{B}\right)}{{\widehat{\sigma }}_{X_i}^2\left({\mathcal{B}}^{\prime }\right)}-1\right)+\frac{1}{2n}\left(\frac{\parallel {\widehat{\mathbf{x}}}_i\left(\mathcal{B}\right)-{\widehat{\mathbf{x}}}_i\left({\mathcal{B}}^{\prime }\right){\parallel }_2^2}{{\widehat{\sigma }}_{X_i}^2\left({\mathcal{B}}^{\prime }\right)}\right)$$

(23)

where, following a similar notation to (4):

• ${\widehat{\mu }}_{X_i}\left(\mathcal{B}\right)$, ${\widehat{\mathit{\beta }}}_{X_i}\left(\mathcal{B}\right)$, ${\widehat{\mu }}_{X_i}\left({\mathcal{B}}^{\prime }\right)$, ${\widehat{\mathit{\beta }}}_{X_i}\left({\mathcal{B}}^{\prime }\right)$ are the estimated intercepts and regression coefficients;
• ${\widehat{\mathbf{x}}}_i\left(\mathcal{B}\right)$ and ${\widehat{\mathbf{x}}}_i\left({\mathcal{B}}^{\prime }\right)$ are the $n\times 1$ vectors

  $${\widehat{\mathbf{x}}}_i\left(\mathcal{B}\right)={\widehat{\mu }}_{X_i}\left(\mathcal{B}\right)+\mathbf{x}[•;{\prod }_{X_i}\left(\mathcal{B}\right)]{\widehat{\mathit{\beta }}}_{X_i}\left(\mathcal{B}\right),\hspace{1em}{\widehat{\mathbf{x}}}_i(\mathcal{B}\prime )={\widehat{\mu }}_{X_i}\left({\mathcal{B}}^{\prime }\right)+\mathbf{x}[•;{\prod }_{X_i}(\mathcal{B}\prime )]{\widehat{\mathit{\beta }}}_{X_i}(\mathcal{B}\prime ),$$

  the fitted values computed from the data observed for $X_i$, ${\Pi }_{X_i}\left(\mathcal{B}\right)$, ${\Pi }_{X_i}\left({\mathcal{B}}^{\prime }\right)$;
• ${\sigma }_{X_i}^2\left(\mathcal{B}\right)$ and ${\sigma }_{X_i}^2\left({\mathcal{B}}^{\prime }\right)$ are the residual variances in $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$.

We can compute the expression in (23) for each node in

$$\begin{array}{c}\underset{\mathrm{compute}{\widehat{\mathbf{x}}}_i\left(\mathcal{B}\right)\mathrm{and}{\widehat{\mathbf{x}}}_i\left({\mathcal{B}}^{\prime }\right)}{\underbrace{O\left(n\right(\left|{\Pi }_{X_i}\left(\mathcal{B}\right)\right|+\left|{\Pi }_{X_i}\left({\mathcal{B}}^{\prime }\right)\right|+2\left)\right)}}+\underset{\mathrm{compute}\mathrm{the}\mathrm{norm}\parallel {\widehat{\mathbf{x}}}_i\left(\mathcal{B}\right)-{\widehat{\mathbf{x}}}_i\left({\mathcal{B}}^{\prime }\right){\parallel }_2^2}{\underbrace{O\left(n\right)}}=\hfill \\ O\left(n\right(\left|{\Pi }_{X_i}\left(\mathcal{B}\right)\right|+\left|{\Pi }_{X_i}\left({\mathcal{B}}^{\prime }\right)\right|+\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left)\right),\hfill \end{array}$$

which is linear in the sample size if both $\mathcal{G}$ and ${\mathcal{G}}^{\prime }$ are sparse because $\left|{\Pi }_{X_i}\left(\mathcal{B}\right)\right|⩽c$, $\left|{\Pi }_{X_i}\left({\mathcal{B}}^{\prime }\right)\right|⩽c$. In this case, the overall computational complexity simplifies to $O\left(nN(2c+\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$2$}\right.)\right)$. Furthermore, as we pointed out in Scutari et al. [29], the fitted values ${\widehat{\mathbf{x}}}_i\left(\mathcal{B}\right)$, ${\widehat{\mathbf{x}}}_i\left({\mathcal{B}}^{\prime }\right)$ are computed as a by-product of parameter learning: if we consider them to be already available, the above computational complexity is reduced to just O(n) for a single node and $O\left(nN\right)$ overall. We can also replace the fitted values ${\widehat{\mathbf{x}}}_i\left(\mathcal{B}\right)$, ${\widehat{\mathbf{x}}}_i\left({\mathcal{B}}^{\prime }\right)$ in (23) with the corresponding residuals ${\widehat{\epsilon }}_i\left(\mathcal{B}\right)$, ${\widehat{\epsilon }}_i\left({\mathcal{B}}^{\prime }\right)$ because

$$\parallel {\widehat{\mathbf{x}}}_i\left(\mathcal{B}\right)-{\widehat{\mathbf{x}}}_i\left({\mathcal{B}}^{\prime }\right){\parallel }_2^2=\parallel \left(\mathbf{x}\right[•;X_i]-{\widehat{\mathbf{x}}}_i\left(\mathcal{B}\right))-(\mathbf{x}[•;X_i]-{\widehat{\mathbf{x}}}_i\left({\mathcal{B}}^{\prime }\right)){\parallel }_2^2=\parallel {\widehat{\mathbf{\epsilon }}}_i\left(\mathcal{B}\right)-{\widehat{\mathbf{\epsilon }}}_i\left({\mathcal{B}}^{\prime }\right){\parallel }_2^2$$

if the latter are available but the former are not.

Example 10 

(KL between GBNs with parameters estimated from data). For reasons of space, this example is presented as Example A5 in Appendix B.

4.3. Conditional Gaussian BNs

The entropy $H\left(\mathcal{B}\right)$ decomposes into a separate $H\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}\right)$ for each node, of the form (9) for discrete nodes and (12) for continuous nodes with no discrete parents. For continuous nodes with both discrete and continuous parents,

$$H\left(X_i\mid {\Pi }_{X_i}^{\mathcal{B}}\right)=\frac{1}{2}\sum _{{\delta }_{X_i}\in \mathit{Val}\left({\Delta }_{X_i}\right)}{\pi }_{{\delta }_{X_i}}log\left(2\pi {\sigma }_{X_i,{\delta }_{X_i}}^2\left(\mathcal{B}\right)\right)+\frac{1}{2},$$

(24)

where ${\pi }_{{\delta }_{X_i}}$ represents the probability associated with the configuration ${\delta }_{X_i}$ of the discrete parents ${\Delta }_{X_i}$. This last expression can be computed in $\left|\mathit{Val}\left({\Delta }_{X_i}\right)\right|$ time for each node. Overall, the complexity of computing $H\left(\mathcal{B}\right)$ is

$$O(\sum _{X_i\in {\mathbf{X}}_D}\left|{\Theta }_{X_i}\right|+\sum _{X_i\in {\mathbf{X}}_G}max\left\{1,\left|\mathit{Val}\left({\Delta }_{X_i}\right)\right|\right\}).$$

where the max accounts for the fact that $\left|\mathit{Val}\left({\Delta }_{X_i}\right)\right|=0$ when ${\Delta }_{X_i}=⌀$ but the computational complexity is O(1) for such nodes.

Example 11 

(Entropy of a CLGBN). For reasons of space, this example is presented as Example A6 in Appendix B.

As for $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$, we could not find any literature illustrating how to compute it. The partition of the nodes in (5) implies that

$$KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)=\underset{\mathrm{discrete}\mathrm{nodes}}{\underbrace{KL\left({\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)}}+\underset{\mathrm{continuous}\mathrm{nodes}}{\underbrace{KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)}}.$$

(25)

We can compute the first term following Section 4.1: ${\mathbf{X}}_D^{\mathcal{B}}$ and ${\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}$ form two discrete BNs whose DAGs are the spanning subgraphs of $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$ and whose local distributions are the corresponding ones in $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$, respectively. The second term decomposes into

$$KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)=\sum _{{\mathbf{x}}_D\in \mathit{Val}\left({\mathbf{X}}_D\right)}P\left({\mathbf{X}}_D^{\mathcal{B}}={\mathbf{x}}_D\right)KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}={\mathbf{x}}_D\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}={\mathbf{x}}_D\right)$$

(26)

similarly to (10) and (24). We can compute it using the multivariate normal distributions associated with the ${\mathbf{X}}_D^{\mathcal{B}}={\mathbf{x}}_D$ and the ${\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}={\mathbf{x}}_D$ in the global distributions of $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$.

Example 12 

(General-case KL between two CLGBNs). Consider the CLGBNs $\mathcal{B}$ from Figure 3 top, which we already used in Examples 3 and 11, and ${\mathcal{B}}^{\prime }$ from Figure 3 bottom. The variables ${\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}$ identify the following mixture components in the global distribution of ${\mathcal{B}}^{\prime }$:

$$\begin{array}{cc}\hfill \{a,c,e\},\{b,c,e\},\{a,d,e\},\{b,d,e\}& \mapsto \left\{e\right\},\hfill \\ \hfill \{a,c,f\},\{b,c,f\},\{a,d,f\},\{b,d,f\}& \mapsto \left\{f\right\}.\hfill \end{array}$$

Therefore, ${\mathcal{B}}^{\prime }$ only encodes two different multivariate normal distributions.

Firstly, we construct two discrete BNs using the subgraphs spanning ${\mathbf{X}}_D^{\mathcal{B}}={\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}=\{X_1,X_2,X_3\}$ in $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$, which have arcs $\{X_1\to X_2\}$ and $\{X_1\to X_2,X_2\to X_3\}$, respectively. The CPTs for $X_1$, $X_2$ and $X_3$ are the same as in $\mathcal{B}$ and in ${\mathcal{B}}^{\prime }$. We then compute $KL\left({\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)=0.577$ following Example 5.

Secondly, we construct the multivariate normal distributions associated with the components of ${\mathcal{B}}^{\prime }$ following Example 3 (in which we computed those of $\mathcal{B}$). For $\left\{e\right\}$, we have

$$\begin{array}{c}\hfill \left[\begin{array}{c}{X}_4\\ X_5\\ X_6\end{array}\right]\sim N\left(\left[\begin{array}{c}0.300\\ 1.400\\ 1.140\end{array}\right],\left[\begin{array}{ccc}0.160& 0.000& 0.032\\ 0.000& 1.690& 1.183\\ 0.032& 1.183& 2.274\end{array}\right]\right);\end{array}$$

for $\left\{f\right\}$, we have

$$\begin{array}{c}\hfill \left[\begin{array}{c}{X}_4\\ X_5\\ X_6\end{array}\right]\sim N\left(\left[\begin{array}{c}1.000\\ 0.500\\ 0.650\end{array}\right],\left[\begin{array}{ccc}0.090& 0.000& 0.018\\ 0.000& 2.250& 1.575\\ 0.018& 1.575& 2.546\end{array}\right]\right).\end{array}$$

Then,

$$\begin{array}{cc}\hfill & KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)\hfill \\ \hfill & =\sum _{x_1\in \{a,b\}}\sum _{x_2\in \{c,d\}}\sum _{x_3\in \{e,f\}}P\left({\mathbf{X}}_D^{\mathcal{B}}=\{x_1,x_2,x_3\}\right)·\hfill \\ \hfill & KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}=\{x_1,x_2,x_3\}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}=\{x_1,x_2,x_3\}\right)\hfill \\ \hfill & =\underset{\{a,c,e\}}{\underbrace{0.040\times 1.721}}+\underset{\{b,c,e\}}{\underbrace{0.036\times 1.721}}+\underset{\{a,d,e\}}{\underbrace{0.040\times 2.504}}+\underset{\{b,d,e\}}{\underbrace{0.084\times 2.504}}+\hfill \\ \hfill & \underset{\{a,c,f\}}{\underbrace{0.16\times 4.303}}+\underset{\{b,c,f\}}{\underbrace{0.144\times 4.303}}+\underset{\{a,d,f\}}{\underbrace{0.16\times 6.31}}+\underset{\{b,d,f\}}{\underbrace{0.336\times 6.31}}\hfill \\ \hfill & =4.879\hfill \end{array}$$

and $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)=KL\left({\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)+KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)=0.577+4.879=5.456$.

The computational complexity of this basic approach to computing $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$ is

$$\begin{array}{c}\underset{\mathrm{compute}KL\left({\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)}{\underbrace{O(Mw{l}^{w+c}+M(w+w{l}^w+l^{w-1})+(M{l}^c+2)\left|{\Theta }_{{\mathbf{X}}_D}\right|)}}+\hfill \\ \underset{\mathrm{compute}\mathrm{all}\mathrm{the}KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}={\mathbf{x}}_D\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}={\mathbf{x}}_D\right)}{\underbrace{O(l^M·\left(6{(N-M)}^3+{(N-M)}^2+5(N-M)\right))}},\hfill \end{array}$$

(27)

which we obtain by adapting (11) and (15) to follow the notation $\left|{\mathbf{X}}_D\right|=M$ and $\left|{\mathbf{X}}_G\right|=N-M$ we established in Section 3.3. The first term implicitly covers the cost of computing the $P\left({\mathbf{X}}_D^{\mathcal{B}}={\mathbf{x}}_D\right)$, which relies on exact inference like the computation of $KL\left({\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)$. The second term is exponential in M, which would lead us to conclude that it is computationally unfeasible to compute $KL\left(\mathcal{B}\parallel {\mathcal{B}}^{\prime }\right)$ whenever we have more than a few discrete variables in $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$. Certainly, this would agree with Hershey and Olsen [22], who reviewed various scalable approximations of the KL divergence between two Gaussian mixtures.

However, we would again disregard the fact that BNs are sparse models. Two properties of CLGBNs that are apparent from Examples 3 and 12 allow us to compute (26) efficiently:

• We can reduce ${\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}$ to ${\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{\Delta }}^{\mathcal{B}}$ where ${\mathbf{\Delta }}^{\mathcal{B}}={\bigcup }_{X_i\in {\mathbf{X}}_G}{\Delta }_{X_i}^{\mathcal{B}}\subseteq {\mathbf{X}}_D^{\mathcal{B}}$. In other words, the continuous nodes are conditionally independent on the discrete nodes that are not their parents (${\mathbf{X}}_D^{\mathcal{B}}\setminus {\mathbf{\Delta }}^{\mathcal{B}}$) given their parents (${\mathbf{\Delta }}^{\mathcal{B}}$). The same is true for ${\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}$. The number of distinct terms in the summation in (26) is then given by $\left|\mathit{Val}({\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }})\right|$ which will be smaller than $\left|\mathit{Val}\left({\mathbf{X}}_D^{\mathcal{B}}\right)\right|$ in sparse networks.
• The conditional distributions ${\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}=\mathit{\delta }$ and ${\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}=\mathit{\delta }$ are multivariate normals (not mixtures). They are also faithful to the subgraphs spanning the continuous nodes ${\mathbf{X}}_G$, and we can represent them as GBNs whose parameters can be extracted directly from $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$. Therefore, we can use the results from Section 4.2 to compute their Kullback–Leibler divergences efficiently.

As a result, (26) simplifies to

$$\begin{array}{c}KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)=\hfill \\ \hfill \sum _{\mathit{\delta }\in \mathit{Val}({\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }})}P\left(\{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\mathit{\delta }\right)KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\mathit{\delta }\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\mathit{\delta }\right).\end{array}$$

where $P\left(\{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\mathit{\delta }\right)$ is the probability that the nodes ${\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}$ take value $\mathit{\delta }$ as computed in $\mathcal{B}$. In turn, (27) reduces to

$$\begin{array}{c}\underset{\mathrm{compute}KL\left({\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)}{\underbrace{O(Mw{l}^{w+c}+M(w+w{l}^w+l^{w-1})+(M{l}^c+2)\left|{\Theta }_{{\mathbf{X}}_D}\right|)}}+\hfill \\ \underset{\mathrm{compute}\mathrm{all}\mathrm{the}KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\mathit{\delta }\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\mathit{\delta }\right)}{\underbrace{O(l^{\left|\mathit{Val}\right({\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\left)\right|}·\left({(N-M)}^3+5{(N-M)}^2+6(N-M)\right))}}.\hfill \end{array}$$

because we can replace $l^M$ with $l^{\left|\mathit{Val}\right({\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\left)\right|}$, which is an upper bound to the unique components in the mixture, and because we replace the complexity in (15) with that (20). We can also further reduce the second term to quadratic complexity as we discussed in Section 4.2. The remaining drivers of the computational complexity are:

• the maximum clique size w in the subgraph spanning ${\mathbf{X}}_D^{\mathcal{B}}$;
• the number of arcs from discrete nodes to continuous nodes in both $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$ and the overlap between ${\mathbf{\Delta }}^{\mathcal{B}}$ and ${\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}$.

Example 13 

(Sparse KL between two CLGBNs). Consider again the CLGBNs $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$ from Example 12. The node sets ${\mathbf{\Delta }}^{\mathcal{B}}=\{X_2,X_3\}$ and ${\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}=\left\{X_3\right\}$ identify four KL divergences to compute: $\mathit{Val}({\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }})=\left\{\{c,e\},\{c,f\},\{d,e\},\{d,f\}\right\}$.

$$\begin{array}{c}KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)=\hfill \\ P\left(\{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{c,e\}\right)KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{c,e\}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{c,e\}\right)+\\ P\left(\{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{c,f\}\right)KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{c,f\}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{c,f\}\right)+\\ P\left(\{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{d,e\}\right)KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{d,e\}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{d,e\}\right)+\\ \hfill P\left(\{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{d,f\}\right)KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{d,f\}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}=\{d,f\}\right)\end{array}$$

All the BNs in the Kullback–Leibler divergences are GBNs whose structure and local distributions can be read from $\mathcal{B}$ and ${\mathcal{B}}^{\prime }$. The four GBNs associated with ${\mathbf{X}}_G^{\mathcal{B}}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}$ have nodes ${\mathbf{X}}_G^{\mathcal{B}}=\{X_4,X_5,X_6\}$, arcs $\{X_5\to X_4,X_4\to X_6\}$ and the local distributions listed in Figure 3. The corresponding GBNs associated with ${\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid \{{\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\}$ are, in fact, only two distinct GBNs associated with $\left\{e\right\}$ and $\left\{f\right\}$. They have arcs $\{X_4\to X_6,X_5\to X_6\}$ and local distributions: for $\left\{e\right\}$,

$$\begin{array}{cccc}\hfill & X_4=0.3+{\epsilon }_{X_4},\hfill & \hfill & {\epsilon }_{X_4}\sim N(0,0.16),\hfill \\ \hfill & X_5=1.4+{\epsilon }_{X_5},\hfill & \hfill & {\epsilon }_{X_5}\sim N(0,1.69),\hfill \\ \hfill & X_6=0.1+0.2X_4+0.7X_5+{\epsilon }_{X_6},\hfill & \hfill & {\epsilon }_{X_6}\sim N(0,1.44);\hfill \end{array}$$

for $\left\{f\right\}$,

$$\begin{array}{cccc}\hfill & X_4=1.0+{\epsilon }_{X_4},\hfill & \hfill & {\epsilon }_{X_4}\sim N(0,0.09),\hfill \\ \hfill & X_5=0.5+{\epsilon }_{X_5},\hfill & \hfill & {\epsilon }_{X_5}\sim N(0,2.25),\hfill \\ \hfill & X_6=0.1+0.2X_4+0.7X_5+{\epsilon }_{X_6},\hfill & \hfill & {\epsilon }_{X_6}\sim N(0,1.44).\hfill \end{array}$$

Plugging in the numbers,

$$\begin{array}{c}KL\left({\mathbf{X}}_G^{\mathcal{B}}\mid {\mathbf{X}}_D^{\mathcal{B}}\parallel {\mathbf{X}}_G^{{\mathcal{B}}^{\prime }}\mid {\mathbf{X}}_D^{{\mathcal{B}}^{\prime }}\right)=\underset{\{c,e\}}{\underbrace{0.076\times 1.721}}+\underset{\{c,f\}}{\underbrace{0.304\times 4.303}}+\hfill \\ \underset{\{d,e\}}{\underbrace{0.124\times 2.504}}+\underset{\{d,f\}}{\underbrace{0.496\times 6.310}}=4.879\hfill \end{array}$$

which matches the value we computed in Example 12.

5. Conclusions

We started this paper by reviewing the three most common distributional assumptions for BNs: discrete BNs, Gaussian BNs (GBNs) and conditional linear Gaussian BNs (CLGBNs). Firstly, we reviewed the link between the respective global and local distributions, and we formalised the computational complexity of decomposing the former into the latter (and vice versa).

We then leveraged these results to study the complexity of computing Shannon’s entropy. We can, of course, compute the entropy of a BN from its global distribution using standard results from the literature. (In the case of discrete BNs and CLGBNS, only for small networks because $|\Theta |$ grows combinatorially.) However, this is not computationally efficient because we incur the cost of composing the global distribution. While the entropy does not decompose along with the local distributions for either discrete BNs or CLGBNS, we show that it is nevertheless efficient to compute it from them.

Computing the Kullback–Leibler divergence between two BNs following the little material found in the literature is more demanding. The discrete case has been thoroughly investigated by Moral et al. [21]. However, the literature typically relies on composing the global distributions for GBNs and CGBNs. Using the local distributions, thus leveraging the intrinsic sparsity of BNs, we showed how to compute the Kullback–Leibler divergence exactly with greater efficiency. For GBNs, we showed how to compute the Kullback–Leibler divergence approximately with quadratic complexity (instead of cubic). If the two GBNs have compatible node orderings and their parameters are estimated from the same data, we can also approximate their Kullback–Leibler divergence with complexity that scales with the number of parents of each node. All these results are summarised in

Table A1. Summary of all the computational complexity results in this paper, including the type of BN and the page where they have been derived.

Composing and decomposing the global distributions
$O(N{l}^N+l{\sum }_{i=1}^N{l}^{|{\Pi }_{X_i}|})$	discrete BNs	Section 3.1
$O\left(N^3\right)$	GBNs	Section 3.2
$O(M{l}^M+{(N-M)}^3{l}^{{\mathbf{\Delta }}^{}})$	CLGBNs	Section 3.3
Computing Shannon’s entropy
$O(N(w(1+l^w)+l^{w-1})+|\Theta |)$	discrete BNs	Section 4.1
O(N)	GBNs	Section 4.2
$O({\sum }_{X_i\in {\mathbf{X}}_D}\left|{\Theta }_{X_i}\right|+{\sum }_{X_i\in {\mathbf{X}}_G}max\left\{1,\left|\mathit{Val}\left({\Delta }_{X_i}\right)\right|\right\})$	CLGBNs	Section 4.3
Computing the Kullback–Leibler divergence
$O(N^{2}w{l}^{w+c}+N(w+w{l}^w+l^{w-1})+(N{l}^c+2)|\Theta |)$	discrete BNs	Section 4.1
$O(6N^3+N^2+5N)$	GBNs	Section 4.2
$O(Mw{l}^{w+c}+M(w+w{l}^w+l^{w-1})+(M{l}^c+2)\left|{\Theta }_{{\mathbf{X}}_D}\right|)+$
$O\left(l^M·\left(6{(N-M)}^3+{(N-M)}^2+5(N-M)\right)\right)$	CLGBNs	Section 4.3
Sparse Kullback–Leibler divergence
$O(N^3+5N^2+6N)$	GBNs	Section 4.2
$O(Mw{l}^{w+c}+M(w+w{l}^w+l^{w-1})+(M{l}^c+2)\left|{\Theta }_{{\mathbf{X}}_D}\right|)+$
$O\left(l^{\left|\mathit{Val}\right({\mathbf{\Delta }}^{\mathcal{B}}\cup {\mathbf{\Delta }}^{{\mathcal{B}}^{\prime }}\left)\right|}·\left({(N-M)}^3+5{(N-M)}^2+6(N-M)\right)\right)$	CLGBNs	Section 4.3
Approximate Kullback–Leibler divergence
$O(7N^2+6N)$	GBNs	Section 4.2
Efficient empirical Kullback–Leibler divergence
$O\left(nN\right(2c+\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left)\right)$	GBNs	Section 4.2

in Appendix A.

Finally, we provided step-by-step numeric examples of how to compute Shannon’s entropy and the Kullback–Leibler divergence for discrete BNs, GBNs and CLGBNs. (See also Appendix B). Considering this is a highly technical topic, and no such examples are available anywhere in the literature, we feel that they are helpful in demystifying this topic and in integrating BNs into many general machine learning approaches.