On Neural Networks as Infinite Tree-Structured Probabilistic Graphical
  Models
Deep neural networks (DNNs) lack the precise semantics and definitive
probabilistic interpretation of probabilistic graphical models (PGMs). In this
paper, we propose an innovative solution by constructing infinite
tree-structured PGMs that correspond exactly to neural networks. Our research
reveals that DNNs, during forward propagation, indeed perform approximations of
PGM inference that are precise in this alternative PGM structure. Not only does
our research complement existing studies that describe neural networks as
kernel machines or infinite-sized Gaussian processes, it also elucidates a more
direct approximation that DNNs make to exact inference in PGMs. Potential
benefits include improved pedagogy and interpretation of DNNs, and algorithms
that can merge the strengths of PGMs and DNNs.
Boyao Li, Alexandar J. Thomson, Matthew M. Engelhard, David Page
# On Neural Networks as Infinite Tree-Structured Probabilistic Graphical Models

###### Abstract

Deep neural networks (DNNs) lack the precise semantics and definitive probabilistic interpretation of probabilistic graphical models (PGMs). In this paper, we propose an innovative solution by constructing infinite tree-structured PGMs that correspond exactly to neural networks. Our research reveals that DNNs, during forward propagation, indeed perform approximations of PGM inference that are precise in this alternative PGM structure. Not only does our research complement existing studies that describe neural networks as kernel machines or infinite-sized Gaussian processes, it also elucidates a more direct approximation that DNNs make to exact inference in PGMs. Potential benefits include improved pedagogy and interpretation of DNNs, and algorithms that can merge the strengths of PGMs and DNNs.

## 1 Introduction

Deep neural networks (DNNs), including large language models, offer state-of-the-art prediction performance, but they are difficult to interpret due to their complex multilayer structure, large number of latent variables, and the presence of nonlinear activation functions (Buhrmester et al., 2021). To gain a precise statistical interpretation for DNNs, much progress has been made in linking them to probabilistic graphical models (PGMs). Variational autoencoders (VAEs) (Kingma and Welling, 2014) are an early example; more recent examples relate recurrent neural networks (RNNs) with hidden Markov models (HMMs) (Choe et al., 2017) and convolutional neural networks (CNNs) with Gaussian processes (GPs) (Garriga-Alonso et al., 2018). When such a connection is possible, benefits include:
* Clear statistical semantics for a trained DNN model beyond providing the conditional distribution over output variables given input variables. Instead, PGMs provide a joint distribution over all variables including latent variables.
* Ability to make inferences about how evidence of some nodes influences probabilities at others, including how later nodes influence earlier ones, as in Bayes nets or Markov nets.
* Ability to understand weight initializations of DNNs as representing prior distributions and trained DNNs as representing posterior distributions, or ensembles of models.
* Proposal of new algorithms by importing algorithmic approaches from PGMs into DNNs.

In this paper, we establish a correspondence between DNNs and PGMs. Given an arbitrary DNN, we first construct an infinite-width tree-structured PGM. We then demonstrate that during training, the DNN executes approximations of precise inference in the PGM during the forward propagation. We prove our result in the case of sigmoid activations and then indicate how the proof can be expanded to other activation functions, provided that some form of normalization is employed. These findings provide immediate benefits such as those listed above.

This work stands apart from most theoretical analyses of DNNs, which typically view DNNs purely as _function approximators_ and prove theorems about the quality of function approximation. Here we instead show that DNNs may be viewed as statistical models, specifically PGMs. This work is also different from the field of _Bayesian neural networks_, where the goal is to seek and model a probability distribution over neural network parameters. In our work, the neural network itself defines a joint probability distribution over its variables (nodes). Our work therefore is synergistic with Bayesian neural networks but more closely related to older work to learn stochastic neural networks via expectation maximization (EM) (Amari, 1995) or approximate EM (Song et al., 2016).

Although the approach is different, our motivation is similar to that of Dutordoir et al. (2021) and Sun et al. (2020) in their work to link DNNs to deep Gaussian processes (GPs) (Damianou and Lawrence, 2013). By identifying the forward pass of a DNN with the mean of a deep GP layer, they aim to augment DNNs with advantages of GPs, notably the ability to quantify uncertainty over both output and latent nodes. What distinguishes our work is that we make the DNN-PGM approximation explicit and include _all_ sigmoid DNNs, not just unsupervised belief networks or other specific cases.
Background: Comparison to Bayesian Networks and Markov Networks

Syntactically a Bayesian network (BN) is a directed acyclic graph, like a neural network, whose nodes are random variables. Semantically, a BN represents a full joint probability distribution over its variables as \(P(\vec{v})=\prod_{i}P(v_{i}|pa(v_{i}))\), where \(\vec{v}\) is a complete setting of the variables, and \(pa(v_{i})\) denotes the parents of variable \(v_{i}\). If the conditional probability distributions (CPDs) \(P(v_{i}|pa(v_{i}))\) are all logistic regression models, we refer to the network as a sigmoid BN.

It is well known that given sigmoid activation and a cross-entropy error, training a single neuron by gradient descent is identical to training a logistic regression model. Hence, a neural network under such conditions can be viewed as a "stacked logistic regression model", and also as a Bayesian network with logistic regression CPDs at the nodes. Technically, the sigmoid BN has a distribution over the input variables (variables without parents), whereas the neural network does not, and all nodes are treated as random variables. These distributions are easily added, and distributions of the input variables can be viewed as represented by the joint sample over them in our training set.

A Markov network (MN) syntactically is an undirected graph with potentials \(\phi_{i}\) on its cliques, where each potential gives the relative probabilities of the various settings for its variables (the variables in the clique). Semantically, it defines the full joint distribution on the variables as \(P(\vec{v})=\frac{1}{Z}\prod_{i}\phi_{i}(\vec{v})\) where the partition function \(Z\) is defined as \(\sum_{\vec{v}}\prod_{i}\phi_{i}(\vec{v})\). It is common to use a loglinear form of the same MN, which can be obtained by treating a setting of the variables in a clique as a binary feature \(f_{i}\), and the natural log of the corresponding entry for that setting in the potential for that clique as a weight \(w_{i}\) on that feature; the equivalent definition of the full joint is then \(P(\vec{v})=\frac{1}{Z}e^{\sum_{i}w_{i}f_{i}(\vec{v})}\). For training and prediction at this point the original graph itself is superfluous.

The potentials of an MN may be on subsets of cliques; in that case we simply multiply all potentials on subsets of a clique to derive the potential on the clique itself. If the MN can be expressed entirely as potentials on edges or individual nodes, we call it a "pairwise" MN. An MN whose variables are all binary is a binary MN.

A DNN of any architecture is, like a Bayesian network, a directed acyclic graph. A sigmoid activation can be understood as a logistic model, thus giving a conditional probability distribution for a binary variable given its parents. Thus, there is a natural interpretation of a DNN with sigmoid activations as a Bayesian network (e.g., Bayesian belief network).* As reviewed in theorem 1, this Bayes net in turn is equivalent to (represents the same probability distribution) as a Markov network where every edge of weight \(w\) from variable \(A\) to variable \(B\) has a potential of the following form:
\begin{tabular}{|c|c|c|} \hline  & \(B\) & \(\neg B\) \\ \hline \(A\) & \(e^{w}\) & \(1\) \\ \hline \(\neg A\) & \(1\) & \(1\) \\ \hline \end{tabular}

**Theorem 1**.: _Let \(N\) be a Bayesian belief network whose underlying undirected graph has treewidth 1, and let \(w_{AB}\) denote the coefficient of variable \(A\) in the logistic CPD for its child \(B\). Let \(M\) be a binary pairwise Markov random field with the same nodes and edges (now undirected) as \(N\). Let \(M\)'s potentials all have the value \(e^{w_{AB}}\) if the nodes \(A\) and \(B\) on either side of edge \(AB\) are true, and the value \(1\) otherwise. \(M\) and \(N\) represent the same joint probability distribution over their nodes._

We don't claim theorem 1 is new, but we provide a proof in Appendix A because it captures several components of common knowledge to which we couldn't find a single reference.

For space reasons, we assume the reader is already familiar with the Variable Elimination (VE) algorithm for computing the probability distribution over any query variable(s) given evidence (known values) at other variables in the network. This algorithm is identical for Bayes nets and Markov nets. It repeatedly multiplies together all the potentials (in a Bayes net, conditional probability distributions) involving the variable to be eliminated, and then sums that variable out of the resulting table, until only the query variable(s) remain. Normalization of the resulting table yields the final answer. VE is an exact inference algorithm, meaning its answers are exactly correct.

## 3 The Construction of Tree-structured PGMs

Although both a binary pairwise Markov network (MN) and a Bayesian network (BN) share the same sigmoid functional structure as a DNN with sigmoid activations, it can be shown that the DNN does not in general define the same probability for the output variables given the input variables: forward propagation in the DNN is very fast but yields a different result than VE in the MN or BN, which can be much slower because the inference task is NP-complete. Therefore, if we take the distribution \(\mathcal{D}\) defined by the BN or MN to be the correct meaning of the DNN, the DNN must be using an approximation \(\mathcal{D}^{\prime}\) to \(\mathcal{D}\). Procedurally, the approximation can be shown to be exactly the following: the DNN repeatedly treats the _expectation_ of a variable \(V\), given the values of \(V\)'s parents, as if it were the actual _value_ of \(V\). Thus previously binary variables in the Bayesian network view and binary features in the Markov network view become continuous. While this procedural characterization of the approximation of \(\mathcal{D}^{\prime}\) to \(\mathcal{D}\) is precise, we prefer in the PGM literature to characterize approximate distributions such as \(\mathcal{D}^{\prime}\) with an alternative PGM that precisely corresponds to \(\mathcal{D}^{\prime}\); for example, in some variational methods we may remove edges from a PGM to obtain a simpler PGM in which inference is more efficient. Treewidth-1 (tree-structured or forest-structured) PGMs are among the most desirable because in those exact inference by
VE or other algorithms becomes efficient. We seek to so characterize the DNN approximation here.

To begin, we consider the Bayesian network view of the DNN. Our first step in this construction is to copy the shared parents in the network into separate nodes whose values are not tied. The algorithm for this step is as follows:

1. Consider the observed nodes in the Bayesian network that correspond to the input of the neural network and their outgoing edges.
2. At each node, for each outgoing edge, create a copy of the current node that is only connected to one of the original node's children with that edge. Since these nodes are observed at this step, these copies do all share the same values. The weights on these edges remain the same.
3. Consider then the children of these nodes. Again, for each outgoing edge, make a copy of this node that is only connected to one child with that edge. In this step, for each copied node, we then also copy the entire subgraph formed by all ancestor nodes of the current node. Note that while weights across copies are tied, the values of the copies of any node are not tied. However, since we also copy the subtree of all input and intermediary hidden nodes relevant in the calculation of this node for each copy, the probability of any of these copied nodes being true remains the same across copies.
4. We repeat this process until we have separate trees for each output node in the original deep neural network graph.

This process ultimately creates a graph whose undirected structure is a tree or forest. In the directed structure, trees converge at the output nodes. The probability of any copy of a latent node given the observed input is the same across all the copies, but when sampling, their values may not be the same.

The preceding step alone is still not sufficient to accurately express the deep neural network as a PGM. Recall that in the Markov network view, we have seen that the neural network makes a mean-field approximation where it uses the expected value of a node in place of its actual value. The following additional step in the construction yields this same behavior. This next step of the construction creates \(L\) copies of every non-output node in the network while also copying the entire subtrees of each of these nodes, as was done in step 1. The weight of a copied edges is then set to its original value divided by \(L\). As \(L\) approaches infinity, we show that the gradient in this PGM construction matches the gradient in the neural network exactly.

This second step in the construction can be thought of intuitively by considering the behavior of sampling in the Bayesian network view. Since we make \(L\) copies of each node while also copying the subgraph of its ancestors, these copied nodes all share the same probabilities. As \(L\) grows large, even if we sampled every copied node only once, we would
expect the average value across these \(L\) copies to match the probability of an individual copied node being true. Given that we set the new weights between these copies and their parents as the original weights divided by \(L\), the sum of products (new weights times parent values) yields the average parent value multiplied by the original weight. As \(L\) goes to infinity, we remove sampling bias and the result exactly matches the value of the sigmoid activation function of the neural network, where this expectation in the PGM view is passed repeatedly to the subsequent neurons. The formal proof of this result, based on variable elimination, is found below. There, we show the following:

**Theorem 2**.: _In the PGM construction, as \(L\rightarrow\infty\), \(P(H=1|\vec{x})\rightarrow\sigma(\sum_{j=1}^{M}w_{j}g_{j}+\sum_{i}^{N}\theta_ {i}\sigma(p_{i}))\), for an arbitrary latent node \(H\) in the DNN that has observed parents \(g_{1},...,g_{M}\) and latent parents \(h_{1},...,h_{N}\) that are true with probabilities \(\sigma(p_{1}),...,\sigma(p_{N})\). \(w_{1},...,w_{M}\) and \(\theta_{1},...,\theta_{N}\) are the weights on edges between these nodes and \(H\)._

In order to prove that as \(L\) goes to infinity, this PGM construction does indeed match the neural network's forward propagation, we consider an arbitrary latent node \(H\) with \(N\) unobserved parents \(h_{1},...,h_{N}\) and \(M\) observed parents \(g_{1},...,g_{M}\). The edges between these parents and \(H\) then each have a weight \(\theta_{i}\), \(1\leq i\leq N\), for the unobserved nodes, and \(w_{j}\), \(1\leq j\leq M\), for the observed nodes. The network as a whole has observed evidence \(\vec{x}\). For the rest of this problem we use a Markov network view of the neural network. The potentials for these nodes in the network are as follows:

Figure 1: The first step of the PGM construction where shared latent parents are separated into copies along with the subtree of their ancestors. Copies of nodes H1 and H2 are made in this example.
\begin{tabular}{|c|c|} \hline \(h_{i}\) & \(\neg h_{i}\) \\ \hline \(e^{p_{i}}\) & 1 \\ \hline \end{tabular} Since \(g_{j}\) are observed, their values are found in \(\vec{x}\).

\begin{tabular}{|c|c|c|} \hline  & \(H\) & \(\neg H\) \\ \hline \(h_{i}\) & \(e^{\theta_{i}}\) & 1 \\ \hline \(\neg h_{i}\) & 1 & 1 \\ \hline \end{tabular}  & 
\begin{tabular}{|c|c|c|} \hline  & \(H\) & \(\neg H\) \\ \hline \(g_{j}\) & \(e^{w_{j}}\) & 1 \\ \hline \(\neg g_{j}\) & 1 & 1 \\ \hline \end{tabular} Suppose, then, using the second step of our construction, we make \(L\) copies of all the nodes that were parents of \(H\) in the Bayesian network view of the DNN, \(h_{1}^{1},...,h_{1}^{L},...,h_{N}^{1},...,h_{N}^{L}\) and \(g_{1}^{1},...,g_{1}^{L},...,g_{M}^{1}\) with weights \(\theta_{1}/L,...,\theta_{N}/L\) and \(w_{1}/L,...,w_{M}/L\) respectively. The potential on \(H\) and these copied nodes is then:

\begin{tabular}{|c|c|c|} \hline  & \(H\) & \(\neg H\) \\ \hline \(h_{i}^{k}\) & \(e^{\theta_{i}/L}\) & 1 \\ \hline \(\neg h_{i}^{k}\) & 1 & 1 \\ \hline \end{tabular}  & 
\begin{tabular}{|c|c|c|} \hline  & \(H\) & \(\neg H\) \\ \hline \(g_{j}^{k}\) & \(e^{w_{j}/L}\) & 1 \\ \hline \(\neg g_{j}^{k}\) & 1 & 1 \\ \hline \end{tabular} where \(1\leq i\leq N\), \(1\leq j\leq M\), and \(1\leq k\leq L\). The potentials for each of the copied nodes are the same as the nodes they were originally copied from. We then have that,

\[P(H,h_{1}^{1},...,h_{1}^{L},...,h_{N}^{1},...,h_{N}^{L},g_{1},...,g_{1}^{L},...,g_{M}^{1},...,g_{M}^{L}|\vec{x})\] \[=\frac{1}{Z}\times\prod_{j=1}^{M}\prod_{k=1}^{L}e^{(w_{j}/L)H \times g_{j}^{k}}\times\prod_{i=1}^{N}\prod_{k=1}^{L}e^{(\theta_{i}/L)H\times h _{i}^{k}}\] \[=\frac{1}{Z}\times e^{\sum_{j=1}^{M}w_{j}g_{j}\times H}\times e^{ H(\frac{\theta_{1}}{L}\sum_{k=1}^{L}h_{1}^{k}+...+\frac{\theta_{N}}{L}\sum_{k=1}^{L}h_{ N}^{k})}\.\]

Summing out an arbitrary, copied latent node, \(h_{\alpha}^{\beta}\):
\[\sum_{h_{\alpha}^{\beta},\neg h_{\alpha}^{\beta}}P(H,h_{1}^{1},...,h_{1}^{L},...,h_{N}^{1},...,h_{N}^{L}|\vec{x})\] \[=\frac{1}{Z}\times e^{\sum_{j=1}^{M}w_{j}g_{j}\times H}\times\sum_ {h_{\alpha}^{\beta},\neg h_{\alpha}^{\beta}}\prod_{i=1}^{N}\prod_{k=1}^{L}e^{( \theta_{i}/L)H\times h_{i}^{k}}\] \[=\left(\frac{1}{Z}\times e^{\sum_{j=1}^{M}w_{j}g_{j}\times H} \times e^{p_{\alpha}}e^{(\theta_{\alpha}/L)H}\prod_{\begin{subarray}{c}i=1,.., N\\ (i,k)\neq(\alpha,\beta)\end{subarray}}\prod_{k=1,...L}e^{(\theta_{i}/L)H\times h_{i}^{k}}\right.\] \[\left.+\frac{1}{Z}\times e^{\sum_{j=1}^{M}w_{j}g_{j}\times H} \times\prod_{\begin{subarray}{c}i=1,..,N\\ (i,k)\neq(\alpha,\beta)\end{subarray}}\prod_{k=1,...L}e^{(\theta_{i}/L)H \times h_{i}^{k}}\right)\] \[=\left(\frac{1}{Z}\times e^{\sum_{j=1}^{M}w_{j}g_{j}\times H} \times(e^{p_{\alpha}}e^{(\theta_{\alpha}/L)H}+1)\times\prod_{\begin{subarray} {c}i=1,..,N\\ (i,k)\neq(\alpha,\beta)\end{subarray}}\prod_{k=1,...L}e^{(\theta_{i}/L)H\times h _{i}^{k}}\right)\.\]

Summing out all \(L\) copies of \(h_{\alpha}\):

\[\left(\frac{1}{Z}\times e^{\sum_{j=1}^{M}w_{j}g_{j}\times H}\times(e^{p_{ \alpha}}e^{(\theta_{\alpha}/L)H}+1)^{L}\times\prod_{\begin{subarray}{c}i=1,.., N\\ i\neq\alpha\end{subarray}}\prod_{k=1,...L}e^{(\theta_{i}/L)H\times h_{i}^{k}} \right)\.\]

Summing out the \(L\) copies of each latent parent would then yield:

\[\frac{1}{Z}e^{\sum_{j=1}^{M}w_{j}g_{j}\times H}\times\prod_{i}^{N}(e^{p_{i}}e^ {(\theta_{i}/L)H}+1)^{L}\,\]

which, in turn, gives us:
\[P(H=1|\vec{x}) =\frac{e^{\sum_{j=1}^{M}w_{j}g_{j}\times 1}\times\prod_{i}^{N}(e^{p_ {i}}e^{(\theta_{i}/L)\times 1}+1)^{L}}{e^{\sum_{j=1}^{M}w_{j}g_{j}\times 1}\times \prod_{i}^{N}(e^{p_{i}}e^{(\theta_{i}/L)\times 1}+1)^{L}+e^{\sum_{j=1}^{M}w_{j}g_{j} \times 0}\times\prod_{i}^{N}(e^{p_{i}}e^{(\theta_{i}/L)\times 0}+1)^{L}}\] \[=\frac{e^{\sum_{j=1}^{M}w_{j}g_{j}}\times\prod_{i}^{N}(e^{p_{i}}e^ {(\theta_{i}/L)}+1)^{L}}{e^{\sum_{j=1}^{M}w_{j}g_{j}}\times\prod_{i}^{N}(e^{p_ {i}}e^{(\theta_{i}/L)}+1)^{L}+\prod_{i}^{N}(e^{p_{i}}+1)^{L}}\] \[=\left(\frac{1}{1+\frac{\prod_{i}^{N}(e^{p_{i}}+1)^{L}}{e^{\sum_ {j=1}^{M}w_{j}g_{j}}\times\prod_{i}^{N}(e^{p_{i}}e^{(\theta_{i}/L)}+1)^{L}}} \right)\;.\]

We then consider:

\[\lim_{L\rightarrow\infty}\frac{\prod_{i}^{N}(e^{p_{i}}+1)^{L}}{e ^{\sum_{j=1}^{M}w_{j}g_{j}}\times\prod_{i}^{N}(e^{p_{i}}e^{(\theta_{i}/L)}+1) ^{L}}\] \[=\lim_{L\rightarrow\infty}e^{-\sum_{j=1}^{M}w_{j}g_{j}+\sum_{i=1 }^{N}L\times\log(e^{p_{i}}+1)-\sum_{i=1}^{N}L\times\log(e^{p_{i}}e^{(\theta_{ i}/L)}+1)}\;,\]

\[\lim_{L\rightarrow\infty}-\sum_{j=1}^{M}w_{j}g_{j}+\sum_{i=1}^{N} L\times\log(e^{p_{i}}+1)-\sum_{i=1}^{N}L\times\log(e^{p_{i}}e^{(\theta_{i}/L)}+1)\] \[=-\sum_{j=1}^{M}w_{j}g_{j}+\lim_{L\rightarrow\infty}\frac{\sum_ {i=1}^{N}\left(\log(e^{p_{i}}+1)-\log(e^{p_{i}}e^{(\theta_{i}/L)}+1)\right)}{ 1/L}\]

This limit clearly has the indeterminate form of \(\frac{0}{0}\). Consider, then the following change of variables, \(S=1/L\), and subsequent use of l'Hospital's rule.
\[-\sum_{j=1}^{M}w_{j}g_{j}+\lim_{S\to 0^{+}}\frac{\sum_{i=1}^{N} \left(\log(e^{p_{i}}+1)-\log(e^{p_{i}}e^{(\theta_{i}S)}+1)\right)}{S}\] \[=-\sum_{j=1}^{M}w_{j}g_{j}+\lim_{S\to 0^{+}}\frac{\frac{\partial}{ \partial S}\sum_{i=1}^{N}\left(\log(e^{p_{i}}+1)-\log(e^{p_{i}}e^{(\theta_{i}S) }+1)\right)}{\frac{\partial}{\partial S}S}\] \[=-\sum_{j=1}^{M}w_{j}g_{j}+\lim_{S\to 0^{+}}\frac{\sum_{i=1}^{N} \frac{-1}{e^{p_{i}}e^{(\theta_{i}S)}+1}\times e^{p_{i}}e^{\theta_{i}S}\times \theta_{i}}{1}\] \[=-\sum_{j=1}^{M}w_{j}g_{j}+\lim_{S\to 0^{+}}\sum_{i=1}^{N} \frac{-e^{p_{i}}e^{\theta_{i}S}\times\theta_{i}}{e^{p_{i}}e^{\theta_{i}S}+1}\] \[=-\sum_{j=1}^{M}w_{j}g_{j}-\sum_{i}^{N}\frac{e^{p_{i}}}{e^{p_{i} }+1}\times\theta_{i}=-\sum_{j=1}^{M}w_{j}g_{j}-\sum_{i}^{N}\sigma(p_{i})\theta _{i}\;.\]

Therefore,

\[\lim_{L\rightarrow\infty}\frac{\prod_{i}^{N}(e^{p_{i}}+1)^{L}}{e^{\sum_{j=1}^ {M}w_{j}g_{j}}\times\prod_{i}^{N}(e^{p_{i}}e^{(\theta_{i}/L)}+1)^{L}}=e^{- \sum_{j=1}^{M}w_{j}g_{j}-\sum_{i}^{N}\sigma(p_{i})\theta_{i}}\;,\]

and,

\[\lim_{L\rightarrow\infty}P(H=1|\vec{x})=\frac{1}{1+e^{-\sum_{j=1}^{M}w_{j}g_{ j}-\sum_{i}^{N}\sigma(p_{i})\theta_{i}}}=\sigma(\sum_{j=1}^{M}w_{j}g_{j}+ \sum_{i}^{N}\sigma(p_{i})\theta_{i})\;.\]

This is exactly the result of the deep neural network. Suppose then that \(z\) is a hidden node whose parents in the Bayesian network view are all observed. By our PGM construction, we have that the potential for \(z\) true is \(e^{\sum_{x\in\vec{x}}w_{zx}x}\) where \(w_{zx}\) is the weight between nodes \(z\) and \(x\), and, for \(z\) false, 1. This clearly matches the deep neural network's sigmoid activation in this 'first layer'. Consider, then, the nodes whose parents in the Bayesian network view are either one of these first layer hidden nodes, or an observed node. By our PGM construction, we have shown that so long as nodes in the previous layer are either observed or have sigmoid conditional probabilities, as is the case here, the conditional probability of any nodes that immediately follow will also have a sigmoid conditional probability. Repeating this argument up to the output nodes gives us that the conditional probability in this PGM construction and the activation values of the DNN match for any layer in the DNN.

Consider then the view of the DNN where each layer is defined such that the values of the activation functions for all neurons in a given layer can be calculated using neurons
of the preceding layers. The DNN is structured then such that the first layer depends only on observed evidence and all layers can be calculated sequentially from that starting point. We have already established that nodes with only observed parents have this sigmoid conditional probability. Given the structure of the DNN and Theorem 2, we then have that the corresponding layers in our PGM construction of the DNN can be similarly computed sequentially from that first layer and have conditional probabilities that exactly match the DNN's activations.

## 4 Implications and Extensions

We are not claiming that one should actually carry out the PGM construction used in the preceding section, since that PGM is infinite. Rather, its contribution is to let us understand precisely the approximation that SGD in a DNN is making; although a DNN itself can be understood as a BN or MN, SGD is not using that BN or MN but rather the infinite tree-structured one. While that PGM is infinite, is it built using the original as a template in a straightforward fashion, and hence is easy to understand. Beyond this contribution to comprehensibility and pedagogy, are there other applications?

One application is an ability to use standard PGM algorithms such as Markov chain Monte Carlo (MCMC) to sample latent variables given observed values of input and output variables, such as for producing confidence intervals or understanding relationships among variables. One could already do so using Gibbs sampling in the BN or MN directly represented by the DNN itself (which we will call the "direct PGM"), but then one wouldn't be using the BN or MN that SGD in the DNN actually used during training. For that, our result has shown that one instead needs to use Gibbs sampling in the infinite tree-structured PGM, which is impractical. Nevertheless, for any variable \(V\) in the original DNN, on each iteration a Gibbs sampler takes infinitely many samples of \(V\) given infinitely many samples of each of the members of \(V\)'s Markov blanket in the original DNN. By treating the variables of the original DNN as continuous, with their values approximating their sampled probabilities in the Gibbs sampler, we can instead apply Hamiltonian Monte Carlo or other MCMC methods for continuous variables in the much smaller DNN structure. We explore this approach empirically rather than theoretically in the next section. Another, related application of our result is that one could further fine-tune the trained DNN using other PGM algorithms, such as contrastive divergence. We also explore this use in the next section.

One might object that most results in this paper use sigmoid activation functions. Nair and Hinton showed that rectified linear units (ReLU) can be thought of as a combination of infinitely many sigmoid units with varying biases (Nair and Hinton, 2010). Hence our result in the previous section can be extended to ReLU activations by the same argument. More generally, with any non-negative activation function that can yield values greater than one, while our BN argument no longer holds, the MN version of the argument can be extended. An MN already requires normalization to represent a probability distribution. While Batch
Normalization and Layer Normalization typically are motivated procedurally, to keep nodes from "saturating," and consequently to keep gradients from "exploding" or "vanishing," as the names suggest, they also can be used to bring variables into the range \([0,1]\) and hence to being considered as probabilities. Consider an idealized variant of these that begins by normalizing all the values coming from a node \(h\) of a neural network, over a given minibatch, to sum to \(1.0\); the argument can be extended to a set of \(h\) and all its siblings in a layer (or other portion of the network structure) assumed to share their properties. It is easily shown that if the parents of any node \(h\) in the neural network provide to \(h\) approximate probabilities that those parent variables are true in the distribution defined by the Markov network given the inputs, then \(h\) in turn provides to its children an approximate probability that \(h\) is true in the distribution defined by the Markov network given the inputs. Use of Batch or Layer Normalization is only approximate and hence adds an additional source of approximation to the result of the preceding section. Detailed consideration of other activation functions is left for further work; in the next section we return to the sigmoid case.

## 5 Alternative Training Algorithms: The Sigmoid Case

To illustrate the potential utility of the infinite tree-structured PGM view of a DNN, in this section we pursue one of its potential implications in depth. We have already noted we can view forward propagation in an all-sigmoid DNN as exact inference in a tree-structured BN, such that the CPD of each hidden variable is a logistic regression. In other words, each hidden node is a Bernoulli random variable, with parameter \(\lambda\) being a sigmoid activation (i.e. logistic function) applied to a linear function of the parent nodes. This view suggests alternative learning algorithms such as contrastive divergence (CD) that use sampling methods for high-dimensional binary random variables such as Gibbs sampling. Doing so has a natural advantage over SGD, which samples the values of the hidden variables using only the evidence at the input values. Instead, Gibbs sampling uses _all_ the available evidence, both at input and output variables. MCMC has now advanced beyond Gibbs sampling with methods such as Hamiltonian Monte Carlo (HMC), but HMC will sample values in \(\{0,1\}\) rather than in \([0,1]\).

To use HMC as proposed, we define hidden variables using the recently-developed _continuous_ Bernoulli distribution (Loaiza-Ganem and Cunningham, 2019), where the single parameter \(\lambda\) of the distribution is defined as in the Bernoulli case. Whereas the Bernoulli density is unnormalized when viewed (improperly) as a density over \((0,1)\), the continuous Bernoulli distribution is a proper density. Somewhat counter-intuitively, the expected value of this distribution is _not_ equal to the parameter \(\lambda\), which has important implications for inference. This option leads to learning algorithms that are able to take advantage of sampling methods effective for high-dimensional continuous random variables, including HMC.

With respect to this BN, whose variables correspond exactly with those of the DNN, we can see that our previous CD-1 algorithm, which is standard SGD for DNNs, samples
settings of the input variables, computes expectations on all the remaining variables (latent and output variables), and adjusts the weights (CPDs) toward maximizing the probability conditional log likelihood, or minimizing cross-entropy. If instead of one gradient step we continued to convergence, the resulting algorithm _almost_ would be the well-known Expectation Maximization (EM) algorithm for training any BN's parameters from a data set with missing or hidden variables, given the BN structure. We say _almost_ because this analysis reveals one important "error" or shortcoming in the SGD algorithm: in the E step where we compute the expected values of the latent variables, it _entirely ignores_ the values of the output variables. In other words, we can view SGD as an approximation to EM in which evidence from output variables is ignored during the E step, and therefore gradients must be backpropagated across _all_ layers to account for this evidence when updating the weights in the M step. This strategy is effective in later layers, which are closer to the output evidence, but highly ineffective in earlier layers. This limitation has been recognized as the vanishing gradient problem and addressed through ad-hoc initialization and pre-training strategies as well as Nesterov momentum. Nesterov momentum can be seen as "looking ahead" one step of SGD when computing the gradient, which partially incorporates evidence at output variables into the expectations at the hidden nodes.

More precisely and formally correcting this shortcoming is not easy: computing the correct expectation is NP-complete, and the most obvious algorithm always requires time exponential in the number of latent variables. In such situations in PGMs, researchers have replaced the E step with MCMC sampling (e.g., Gibbs) (Hinton et al., 2006; Song et al., 2016). However, running these MCMC chains to convergence is impractical, therefore it is common in practice to take a small number \(k\) of steps in the chain between gradient steps, which again gives rise to the CD-1 or CD-\(k\) algorithm. However, a different training algorithm for DNNs, motivated by the natural correspondence between DNNs and BNs described above - and which correctly accounts for evidence from the output variables through proper EM updates - converges to the correct answer in fewer training epochs compared to SGD.

The Continuous Bernoulli Bayes net (CBBN) is similar to the sigmoid BN (_i.e._, "stacked logistic regression"), but with logistic regression CPDs replaced by their continuous Bernoulli analogues. Equivalently, it is a feedforward, stochastic neural network in which hidden variables are continuous Bernoulli distributed. Consider a Bayesian network composed of input variables \(\mathbf{x}=\mathbf{h}_{0}\), a sequence of layers of hidden variables \(\mathbf{h}_{1},...,\mathbf{h}_{L}\), and output variables \(\mathbf{y}\). Each pair of consecutive layers forms a bipartite subgraph of the network as a whole, and the variables \(\mathbf{h}_{i}=(h_{i1},...,h_{iM_{i}})\) follow a multivariate continuous Bernoulli distribution with parameters \(\mathbf{\lambda}_{i}=(\lambda_{i1},...,\lambda_{iM_{i}})\) that depend on variables in the previous layer \(\mathbf{h}_{i-1}\) as follows:

\[h_{ij}\sim\mathcal{CB}(\lambda_{ij}),\,\text{where}\,\,\mathbf{\lambda}_{i}= \sigma(\mathbf{W}_{i-1}\mathbf{h}_{i-1}+\mathbf{b}_{i-1}). \tag{1}\]

\(\sigma:\mathbb{R}\rightarrow(0,1)\) is a non-linearity - here the logistic function - that is applied element-wise, and \(\mathbf{\theta}_{i}=(\mathbf{W}_{i},\mathbf{b}_{i})\) are parameters to be learned. For a complete setting of the variables
\(\{\mathbf{x},\mathbf{h},\mathbf{y}\}\), where \(\mathbf{h}=\{\mathbf{h}_{1},...,\mathbf{h}_{L}\}\), and parameters \(\mathbf{\theta}=\{\mathbf{\theta}_{i}\}_{i=0}^{L}\), the likelihood \(p(\mathbf{y},\mathbf{h}|\mathbf{x};\mathbf{\theta})\) may be decomposed as:

\[p(\mathbf{y},\mathbf{h}|\mathbf{x};\mathbf{\theta})=p(\mathbf{y}|\mathbf{h}_{L};\mathbf{ \theta}_{L})\cdot\prod_{i=1}^{L}\prod_{j=1}^{M_{i}}p_{\mathcal{CB}}(h_{ij}| \lambda_{ij}(\mathbf{h}_{i-1};\mathbf{\theta}_{i-1})), \tag{2}\]

where \(p_{\mathcal{CB}}(\cdot|\cdot)\) denotes the continuous Bernoulli density, and a specific form for \(p(\mathbf{y}|\mathbf{h}_{L},\mathbf{\theta}_{L})\) has been omitted to allow variability in the output variables. In our experiments, \(\mathbf{y}\) is a Bernoulli or categorical random variable parameterized via the logistic or softmax function, respectively.

### Learning via Contrastive Divergence with Hamiltonian Monte Carlo Sampling

Let \(\mathbf{h}^{(0)},\mathbf{h}^{(1)},\mathbf{h}^{(2)},...\) denote a chain of MCMC samples of the complete setting of hidden variables in our CBBN. As previously noted, we allow hidden variables \(h_{ij}\in(0,1)\) for \(i\in\{1,...,L\}\) and \(j\in\{1,...,M_{i}\}\), and use Hamiltonian Monte Carlo (HMC) to generate the next state due to its fast convergence. Since HMC samples are unbounded, we sample the _logit_ associated with \(h_{ij}\in(0,1)\), i.e. \(\sigma^{-1}(h_{ij})\in(-\infty,\infty)\), rather than sampling the \(h_{ij}\) directly.

The HMC trajectories are defined by Hamilton's Equations:

\[\frac{d\rho_{i}}{dt}=\frac{\partial H}{\partial\mu_{i}} \frac{d\mu_{i}}{dt}=-\frac{\partial H}{\partial\rho_{i}} \tag{3}\]

where \(\rho_{i},\mu_{i}\) are the \(i\)th component of the position and momentum vector. The Hamiltonian \(H\) is

\[H=H(\mathbf{\rho},\mathbf{\mu})=U(\mathbf{\rho})+\frac{1}{2}\mathbf{\mu}^{T}M^{- 1}\mathbf{\mu} \tag{4}\]

where \(M\) is a positive definite and symmetric mass matrix, and \(M^{-1}\) could represent a diagonal estimate of the covariance. Defining the position \(\mathbf{\rho}=\mathbf{h}\), the complete set of hidden variables of our CBBN, we have that the potential energy \(U\) is the negative log-likelihood associated with equation (2):

\[U(\mathbf{h})=-\log p(\mathbf{y},\mathbf{h}|\mathbf{x};\mathbf{\theta})=-\log p(\mathbf{ y}|\mathbf{h}_{L};\mathbf{\theta}_{L})-\sum_{i=1}^{L}\sum_{j=1}^{M_{i}}\log p_{ \mathcal{CB}}(h_{ij}|\lambda_{ij}(\mathbf{h}_{i-1};\mathbf{\theta}_{i-1})). \tag{5}\]

We set the leap frog size \(L>0\), step size \(\Delta t>0\). A description of the HMC trajectories (_i.e._, evolution of \(\mathbf{h}\)) is provided in the supplementary material.

The initial state of the chain \(\mathbf{h}^{(0)}\) is drawn with a simple forward pass through the network, ignoring the output variables; in other words, we have \(h_{ij}^{(0)}\sim\mathcal{CB}(\sigma(\mathbf{W}_{i-1}^{(0)}\mathbf{h}_{i-1}^{(0)}+\mathbf{ b}_{i-1}^{(0)})_{j})\) for
\(i\in\{1,...L\}\), where \(\mathbf{h}_{0}=\mathbf{x}\) are the input variables, and the values of \(\mathbf{W}_{i}^{(0)}\) and \(\mathbf{b}_{i}^{(0)}\) are manually set or drawn from a standard normal or uniform distribution. We update \(\mathbf{h}\) through a number of burn-in steps before beginning to update our parameters to ensure that \(\mathbf{h}\) is first consistent with evidence from the output variables. After \(k\) steps, corresponding to CD-\(k\), we define the loss based on equation (2):

\[\mathcal{L}(\mathbf{\theta}^{(n)})=-\log p(\mathbf{y},\mathbf{h}|\mathbf{x};\mathbf{ \theta}^{(n)}). \tag{6}\]

We then apply the following gradients to update the parameters \(\{\mathbf{W}_{i}^{(n)}\}_{i=0}^{L}\) and \(\{\mathbf{b}_{i}^{(n)}\}_{i=0}^{L}\):

\[\mathbf{W}_{i}^{(n+1)}=\mathbf{W}_{i}^{(n)}-\eta\frac{\partial\mathcal{ L}}{\partial\mathbf{W}_{i}^{(n)}} \mathbf{b}_{i}^{(n+1)}=\mathbf{b}_{i}^{(n)}-\eta\frac{\partial\mathcal{L}}{ \partial\mathbf{b}_{i}^{(n)}} \tag{7}\]

where \(\eta\) is the learning rate. Algorithm 1 (see supplementary material) summarizes this procedure.

### Experimental Results

The preceding algorithm shows the potential of the DNN-as-PGM view to generate new algorithms and approaches, but does it work? Taking the view of neural net as BN, we start with the hard problem of learning the exclusive-or function, but with the "right" prior (up to magnitude) on the BN parameters. This algorithm is also CD-k - here we use CD-1 - but under this alternative correspondence between PGM and neural network. We therefore call it CD-HMC to distinguish it from the earlier CD-k algorithm that is identical to SGD. As shown in Table 1, CD-HMC converges in half as many training epochs as SGD using the same minibatch size, and each epoch takes little longer than for SGD. But what if we did not know the correct prior?

Using a hidden layer of thirty variables makes it highly probable that there exists a pair latent variables that, together with the inputs and output, have their weights randomly initialized to the correct prior. The empirical results below bear this out: results are similar to experiments using the correct prior (Table 1). If more than one combination of two such latent variables exist, the resulting trained model becomes an ensemble of accurate posteriors. The argument scales to more complex target functions by using more hidden layers in the neural network. These empirical results support the following simple view of why overparameterization works: more parameters, from more latent units and layers, provides a higher probability of having the correct prior embedded somewhere. And if it is embedded more than once, all the better since the final logistic regression layer simply learns a weighted vote among the various possible ensemble components. Further empirical support for this view can be found elsewhere Frankle and Carbin (2019), but for the first time here we view initialization as encoding multiple priors. This in turn suggests viewing the trained DNN as encoding a posterior distribution, from which we can make inferences about uncertainty.
In addition to the experiments on learning the exclusive-or function, we also explore our method on two other datasets. One is a synthetic dataset generated using the _Make Moons_ from _sklearn_ (1k data, 30% noise). The other is MNIST, where we randomly choose 2k images of the digits 0 and 1. Using networks with one hidden layer of 32 or 128 hidden variables and sigmoid activation, we test CD-HMC on both datasets and compare it to SGD and CD-Gibbs, in which we return to logistic regression CPDs and use Gibbs sampling. The training and test datasets are split (80:20 ratio), and each model is trained for 400 epochs with weights updated by gradient descent (learning rate=0.01). For CD-HMC and CD-Gibbs, we draw 500 "burn-in" samples for each data point before the first weight update. The results below illustrate CD-HMC has similar accuracy to SGD and CD-Gibbs on the test set and converges in fewer epochs on _Make Moons_ dataset (Figure 2). Table 2 also shows that CD-HMC has a higher test loss than SGD under all the settings for networks and datasets. This is likely due to variability in sampling in contrast to SGD, which is deterministic.

As the size of the dataset increases, the training time of CD-HMC and CD-Gibbs becomes considerably longer than SGD due to the high cost of drawing high-dimensional samples. However, HMC samples all hidden nodes together, whereas CD-Gibbs samples them one by one, and as a result, CD-HMC is much faster.

As these results suggest, and consistent with prior work on learning in sigmoid BNs, it is difficult (or perhaps impossible) to match the computational efficiency of SGD while

\begin{table}
\begin{tabular}{c|c|c|c|c} Network & Algorithm & Accuracy & Training Steps (Epochs) & Training time (s) \\ \hline \hline \multirow{2}{*}{A} & SGD (BP) & 100\% & 1693 & 19.8 \\  & CD-HMC & 100\% & 263 & 14.0 \\ \hline \multirow{2}{*}{B} & SGD (BP) & 100\% & 4000 & 40 \\  & CD-HMC & 100\% & 445 & 22.8 \\ \end{tabular}
\end{table}
Table 1: Comparing SGD (CD-1) and new algorithm CD-HMC on learning exclusive-or using 30 hidden variables and a random initialization, given correct initialization (A), or all possible priors (B).

\begin{table}
\begin{tabular}{c|c|c|c|c|c} Network & Algorithm & \multicolumn{2}{c|}{_Make Moons_} & \multicolumn{2}{c}{MNIST} \\ \cline{3-6}  & & Test Accuracy & Test Loss & Test Accuracy & Test Loss \\ \hline \hline \multirow{3}{*}{32 nodes} & SGD (BP) & 80\% & 0.3753 & 100\% & 0.0021 \\  & CD-HMC & 80.5\% & 0.5131 & 99.25\% & 0.2606 \\  & CD-Gibbs & 80.5\% & 0.6425 & 100\% & 0.4663 \\ \hline \multirow{3}{*}{128 nodes} & SGD (BP) & 80.5\% & 0.3751 & 100\% & 0.0017 \\  & CD-HMC & 80.5\% & 0.4822 & 99.75\% & 0.1815 \\ \cline{1-1}  & CD-Gibbs & 81\% & 0.6408 & 99.75\% & 0.4733 \\ \end{tabular}
\end{table}
Table 2: SGD, CD-HMC and CD-Gibbs performance on synthetic data (_Make Moons_, 1k samples, noise=0.3) and MNIST (2k images, \(\{0,1\}\) only) with 32 or 128 hidden units and random initialization.
maintaining the view of DNNs as PGMs. However, a hybrid approach - for example, using SGD (CD-1) as pre-training and then applying CD-HMC - preserves the benefits of understanding the DNN as a PGM, or even as a distribution (ensemble) of PGMs, while also allowing exact inference and quantifying uncertainty.

## 6 Limitations and Future Work

Limitations of the present work and directions for future work include establishing formal results about how closely batch- and layer-normalization approximate Markov network normalization when using non-sigmoid activations, establishing theoretical results relating HMC in the neural network to Gibbs sampling in the large treewidth-1 Markov network, and obtaining empirical results for HMC with non-sigmoid activations. Also of great interest is comparing HMC and other PGM algorithms to Shapley values, Integrated Gradients, and other approaches for assessing the relationship of some latent variables to each other

Figure 2: Test accuracy of SGD, CD-HMC and CD-Gibbs on synthetic data (_Make Moons_, 1k, noise=0.3) and MNIST (2k images, \(\{0,1\}\) only).
or to inputs and/or outputs in a neural network. Finally, the large treewidth-1 PGM is a substantial approximation to the direct PGM of a DNN. After training the DNN and hence the large treewidth-1 model, can we fine-tune with a less-approximate approach, perhaps based on loopy belief propagation or other approximate algorithms often used in PGMs?

## Acknowledgements

The authors would like to thank Sayan Mukherjee, Samuel I. Berchuck, Youngsoo Baek, Andrew Allen, William H. Majoros, David Carlson, Juan Restrepo and Houssam Nassif for their helpful discussion about the theoretical work. We are also grateful to Mengyue Han and Jinyi Zhou for their technical support.

This project is in part supported by Impact of Genomic Variation on Function (IGVF) Consortium of the National Institutes of Health via grant U01HG011967.
