# On Neural Networks as Infinite Tree-Structured Probabilistic Graphical
Models

Boyao Li111Boyao Li and Alexander J. Thomson contributed equally to this work.
Alexandar J. Thomson111Boyao Li and Alexander J. Thomson contributed equally
to this work. Matthew M. Engelhard222Matthew M. Engelhard and David Page
contributed equally in supervising this work. David Page222Matthew M.
Engelhard and David Page contributed equally in supervising this work.
Department of Biostatistics and Bioinformatics, Duke University
<EMAIL_ADDRESS>
###### 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.

##  2 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).***Note that
when viewing the DNN as a function approximator, these conditional probability
distributions are not explicitly defined. 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:

| $B$ | $\neg B$
---|---|---
$A$ | $e^{w}$ | 1
$\neg A$ | $1$ | 1

######  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 $\cal{D}$ defined by the BN
or MN to be the correct meaning of the DNN, the DNN must be using an
approximation $\cal{D^{\prime}}$ to $\cal{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
$\cal{D^{\prime}}$ to $\cal{D}$ is precise, we prefer in the PGM literature to
characterize approximate distributions such as $\cal{D^{\prime}}$ with an
alternative PGM that precisely corresponds to $\cal{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. 1.

Consider the observed nodes in the Bayesian network that correspond to the
input of the neural network and their outgoing edges.

  2. 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. 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. 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.

(a) The neural network’s graphical structure before applying the first step of
this PGM construction.

(b) The neural network’s graphical structure after applying the first step of
this PGM construction.

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.

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\to\infty$,
$P(H=1|\vec{x})\to\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}$, $M$
observed parents $g_{1},...,g_{M}$, and a single child node $C$. 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, and $\gamma$, for the child node. 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:

$h_{i}$ | $\neg h_{i}$
---|---
$e^{p_{i}}$ | 1

Since $g_{j}$ are observed, their values are found in $\vec{x}$. The potential
for the ’child’ node can be represented as:

$C$ | $\neg C$
---|---
$c$ | 1

where $c$ is a positive real number that can be arbitrarily large.

| $H$ | $\neg H$
---|---|---
$h_{i}$ | $e^{\theta_{i}}$ | 1
$\neg h_{i}$ | $1$ | 1
| $H$ | $\neg H$
---|---|---
$g_{j}$ | $e^{w_{j}}$ | 1
$\neg g_{j}$ | $1$ | 1
| $H$ | $\neg H$
---|---|---
$C$ | $e^{\gamma}$ | 1
$\neg C$ | $1$ | 1

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},...,g_{M}^{L}$ with weights
$\theta_{1}/L,...,\theta_{N}/L$ and $w_{1}/L,...,w_{M}/L$ respectively. Note,
as per the construction, additional copies of $C$ are not connected to its
parent $H$, however, its weight becomes $\gamma/L$. The potential on $H$ and
these copied nodes is then:

| $H$ | $\neg H$
---|---|---
$h_{i}^{k}$ | $e^{\theta_{i}/L}$ | 1
$\neg h_{i}^{k}$ | $1$ | 1
| $H$ | $\neg H$
---|---|---
$g_{j}^{k}$ | $e^{w_{j}/L}$ | 1
$\neg g_{j}^{k}$ | $1$ | 1

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,

$\displaystyle
P(H,C,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})$
$\displaystyle=\frac{1}{Z}\times e^{(\gamma/L)H\times
C}\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}}$ $\displaystyle=\frac{1}{Z}\times e^{(\gamma/L)H\times C}\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}$:

$\displaystyle\sum_{h_{\alpha}^{\beta},\neg
h_{\alpha}^{\beta}}P(H,C,h_{1}^{1},...,h_{1}^{L},...,h_{N}^{1},...,h_{N}^{L}|\vec{x})$
$\displaystyle=\frac{1}{Z}\times e^{(\gamma/L)H\times C}\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}}$ $\displaystyle=\frac{1}{Z}\times e^{(\gamma/L)H\times C}\times
e^{\sum_{j=1}^{M}w_{j}g_{j}\times H}\times$
$\displaystyle\left[e^{p_{\alpha}}e^{(\theta_{\alpha}/L)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.$ $\displaystyle\hskip
56.9055pt\left.+\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]$ $\displaystyle=\frac{1}{Z}\times e^{(\gamma/L)H\times
C}\times e^{\sum_{j=1}^{M}w_{j}g_{j}\times
H}\times(e^{p_{\alpha}}e^{(\theta_{\alpha}/L)H}+1)$ $\displaystyle\hskip
56.9055pt\times\left[\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]_{\textstyle.}$

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

$\displaystyle\frac{1}{Z}\times e^{(\gamma/L)H\times C}\times
e^{\sum_{j=1}^{M}w_{j}g_{j}\times
H}\times(e^{p_{\alpha}}e^{(\theta_{\alpha}/L)H}+1)^{L}$ $\displaystyle\hskip
42.67912pt\times\left[\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]_{\textstyle.}$

Summing out the $L$ copies of each latent parent would then yield:

$\frac{1}{Z}\times e^{(\gamma/L)H\times C}\times
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}\;,$

We then sum out $C$:

$\frac{1}{Z}\times(ce^{(\gamma/L)H}+1)\times 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:

$\displaystyle P(H=1|\vec{x})$ $\displaystyle=\left[(ce^{(\gamma/L)\times
1}+1)\times 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}\right]\bigg{/}$ $\displaystyle\hskip
11.38109pt\left[(ce^{(\gamma/L)\times 1}+1)\times
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}\right.$
$\displaystyle\hskip 17.07164pt\left.+(ce^{(\gamma/L)\times 0}+1)\times
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}\right]$
$\displaystyle=\frac{(ce^{(\gamma/L)}+1)\times
e^{\sum_{j=1}^{M}w_{j}g_{j}}\times\prod_{i}^{N}(e^{p_{i}}e^{(\theta_{i}/L)}+1)^{L}}{(ce^{(\gamma/L)}+1)\times
e^{\sum_{j=1}^{M}w_{j}g_{j}}\times\prod_{i}^{N}(e^{p_{i}}e^{(\theta_{i}/L)}+1)^{L}+(c+1)\prod_{i}^{N}(e^{p_{i}}+1)^{L}}$
$\displaystyle=\left[\frac{1}{1+\frac{(c+1)\prod_{i}^{N}(e^{p_{i}}+1)^{L}}{(ce^{(\gamma/L)}+1)\times
e^{\sum_{j=1}^{M}w_{j}g_{j}}\times\prod_{i}^{N}(e^{p_{i}}e^{(\theta_{i}/L)}+1)^{L}}}\right]_{\textstyle.}$

We then consider:

$\displaystyle\lim_{L\to\infty}\frac{(c+1)\prod_{i}^{N}(e^{p_{i}}+1)^{L}}{(ce^{(\gamma/L)}+1)\times
e^{\sum_{j=1}^{M}w_{j}g_{j}}\times\prod_{i}^{N}(e^{p_{i}}e^{(\theta_{i}/L)}+1)^{L}}=$
$\displaystyle\lim_{L\to\infty}\exp\left(-\sum_{j=1}^{M}w_{j}g_{j}+\sum_{i=1}^{N}L\times\log(e^{p_{i}}+1)\right.$
$\displaystyle\hskip
51.21495pt\left.-\sum_{i=1}^{N}L\times\log(e^{p_{i}}e^{(\theta_{i}/L)}+1)-\log(ce^{(\gamma/L)}+1)+\log(c+1)\right)_{\textstyle,}$

and the logarithm of this limit is,

$\displaystyle\lim_{L\to\infty}\left[-\sum_{j=1}^{M}w_{j}g_{j}+\sum_{i=1}^{N}L\times\log(e^{p_{i}}+1)\right.$
$\displaystyle\hskip
28.45274pt\left.-\sum_{i=1}^{N}L\times\log(e^{p_{i}}e^{(\theta_{i}/L)}+1)-\log(ce^{(\gamma/L)}+1)+\log(c+1)\right]$
$\displaystyle=-\sum_{j=1}^{M}w_{j}g_{j}\hskip
5.69054pt+\lim_{L\to\infty}\left[-\log(ce^{(\gamma/L)}+1)+\log(c+1)\right]+$
$\displaystyle\hskip
28.45274pt\lim_{L\to\infty}\frac{\sum_{i=1}^{N}\left(\log(e^{p_{i}}+1)-\log(e^{p_{i}}e^{(\theta_{i}/L)}+1)\right)}{1/L}.$

Note that $\lim_{L\to\infty}\left[-\log(ce^{(\gamma/L)}+1)+\log(c+1)\right]=0$
and the second limit in the previous expression clearly has the indeterminate
form of $\frac{0}{0}$. Let $G=-\sum_{j=1}^{M}w_{j}g_{j}$ and consider, then
the following change of variables, $S=1/L$, and subsequent use of l’Hôspital’s
rule.

$\displaystyle G+\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}$
$\displaystyle=G+\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}$ $\displaystyle=G+\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}$ $\displaystyle=G+\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}$
$\displaystyle=G-\sum_{i}^{N}\frac{e^{p_{i}}}{e^{p_{i}}+1}\times\theta_{i}$
$\displaystyle=-\sum_{j=1}^{M}w_{j}g_{j}-\sum_{i}^{N}\sigma(p_{i})\theta_{i}.$

Therefore,

$\displaystyle\lim_{L\to\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}}$
$\displaystyle=\exp(-\sum_{j=1}^{M}w_{j}g_{j}-\sum_{i}^{N}\sigma(p_{i})\theta_{i}),$

and,

$\displaystyle\lim_{L\to\infty}P(H=1|\vec{x})$
$\displaystyle=\frac{1}{1+e^{-\sum_{j=1}^{M}w_{j}g_{j}-\sum_{i}^{N}\sigma(p_{i})\theta_{i}}}$
$\displaystyle=\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
$\bm{x}=\bm{h}_{0}$, a sequence of layers of hidden variables
$\bm{h}_{1},...,\bm{h}_{L}$, and output variables $\bm{y}$. Each pair of
consecutive layers forms a bipartite subgraph of the network as a whole, and
the variables $\bm{h}_{i}=(h_{i1},...,h_{iM_{i}})$ follow a multivariate
continuous Bernoulli distribution with parameters
$\bm{\lambda}_{i}=(\lambda_{i1},...,\lambda_{iM_{i}})$ that depend on
variables in the previous layer $\bm{h}_{i-1}$ as follows:

$h_{ij}\sim\mathcal{CB}(\lambda_{ij})\text{, where
}\bm{\lambda_{i}}=\sigma(\bm{W}_{i-1}\bm{h}_{i-1}+\bm{b}_{i-1}).$ (1)

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

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

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

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

Let $\bm{h}^{(0)},\bm{h}^{(1)},\bm{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:

$\displaystyle\frac{d\rho_{i}}{dt}=\frac{\partial{H}}{\partial{\mu_{i}}}$
$\displaystyle\frac{d\mu_{i}}{dt}=-\frac{\partial{H}}{\partial{\rho_{i}}}$ (3)

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

$\displaystyle
H=H(\bm{\rho},\bm{\mu})=U(\bm{\rho})+\frac{1}{2}\bm{\mu}^{T}M^{-1}\bm{\mu}$
(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
$\bm{\rho}=\bm{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(\bm{h})=-\log p(\bm{y},\bm{h}|\bm{x};\bm{\theta})=-\log
p(\bm{y}|\bm{h}_{L};\bm{\theta}_{L})-\sum_{i=1}^{L}\sum_{j=1}^{M_{i}}\log
p_{\mathcal{CB}}(h_{ij}|\lambda_{ij}(\bm{h}_{i-1};\bm{\theta}_{i-1})).$ (5)

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

The initial state of the chain $\bm{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(\bm{W}_{i-1}^{(0)}\bm{h}_{i-1}^{(0)}+\bm{b}_{i-1}^{(0)})_{j})$
for $i\in\\{1,...L\\}$, where $\bm{h}_{0}=\bm{x}$ are the input variables, and
the values of $\bm{W}_{i}^{(0)}$ and $\bm{b}_{i}^{(0)}$ are manually set or
drawn from a standard normal or uniform distribution. We update $\bm{h}$
through a number of burn-in steps before beginning to update our parameters to
ensure that $\bm{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}(\bm{\theta}^{(n)})=-\log
p(\bm{y},\bm{h}|\bm{x};\bm{\theta}^{(n)}).$ (6)

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

$\displaystyle\bm{W}_{i}^{(n+1)}=\bm{W}_{i}^{(n)}-\eta\frac{\partial\mathcal{L}}{\partial\bm{W}_{i}^{(n)}}$
$\displaystyle\bm{b}_{i}^{(n+1)}=\bm{b}_{i}^{(n)}-\eta\frac{\partial\mathcal{L}}{\partial\bm{b}_{i}^{(n)}}$
(7)

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

###  5.2 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.

Network | Algorithm | Accuracy | Training Steps (Epochs) | Training time (s)
---|---|---|---|---
A | SGD (BP) | 100% | 1693 | 19.8
CD-HMC | 100% | 263 | 14.0
B | SGD (BP) | 100% |  4000 |  40
CD-HMC | 100% | 445 | 22.8
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).

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.

(a) Make Moons, 32 hidden nodes

(b) Make Moons, 128 hidden nodes

(c) MNIST, 32 hidden nodes

(d) MNIST, 128 hidden nodes

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). Network | Algorithm | Make Moons | MNIST
---|---|---|---
Test Accuracy | Test Loss | Test Accuracy | Test Loss
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
128 nodes | SGD (BP) | 80.5% | 0.3751 | 100% | 0.0017
CD-HMC | 80.5% | 0.4822 | 99.75% | 0.1815
CD-Gibbs | 81% | 0.6408 | 99.75% | 0.4733
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.

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 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 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.

## References

  * Buhrmester et al. [2021] Vanessa Buhrmester, David Münch, and Michael Arens.  Analysis of explainers of black box deep neural networks for computer vision: A survey.  _Machine Learning and Knowledge Extraction_ , 3(4):966–989, 2021.
  * Kingma and Welling [2014] Diederik P. Kingma and Max Welling.  Auto-Encoding Variational Bayes.  In _2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings_ , 2014.
  * Choe et al. [2017] Yo Joong Choe, Jaehyeok Shin, and Neil Spencer.  Probabilistic interpretations of recurrent neural networks.  _Probabilistic Graphical Models_ , 2017.
  * Garriga-Alonso et al. [2018] Adrià Garriga-Alonso, Carl Edward Rasmussen, and Laurence Aitchison.  Deep convolutional networks as shallow gaussian processes.  _arXiv preprint arXiv:1808.05587_ , 2018.
  * Amari [1995] Shun-ichi Amari.  Information geometry of the em and em algorithms for neural networks.  _Neural networks_ , 8(9):1379–1408, 1995.
  * Song et al. [2016] Zhao Song, Ricardo Henao, David Carlson, and Lawrence Carin.  Learning sigmoid belief networks via monte carlo expectation maximization.  In _Artificial Intelligence and Statistics_ , pages 1347–1355. PMLR, 2016.
  * Dutordoir et al. [2021] Vincent Dutordoir, James Hensman, Mark van der Wilk, Carl Henrik Ek, Zoubin Ghahramani, and Nicolas Durrande.  Deep neural networks as point estimates for deep gaussian processes.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Sun et al. [2020] Shengyang Sun, Jiaxin Shi, and Roger Baker Grosse.  Neural networks as inter-domain inducing points.  In _Third Symposium on Advances in Approximate Bayesian Inference_ , 2020.
  * Damianou and Lawrence [2013] Andreas Damianou and Neil D Lawrence.  Deep gaussian processes.  In _Artificial intelligence and statistics_ , pages 207–215. PMLR, 2013.
  * Nair and Hinton [2010] Vinod Nair and Geoffrey E Hinton.  Rectified linear units improve restricted boltzmann machines.  In _Proceedings of the 27th international conference on machine learning (ICML-10)_ , pages 807–814, 2010.
  * Loaiza-Ganem and Cunningham [2019] Gabriel Loaiza-Ganem and John P. Cunningham.  The continuous bernoulli: fixing a pervasive error in variational autoencoders, 2019.
  * Hinton et al. [2006] Geoffrey E Hinton, Simon Osindero, and Yee-Whye Teh.  A fast learning algorithm for deep belief nets.  _Neural computation_ , 18(7):1527–1554, 2006\.
  * Frankle and Carbin [2019] Jonathan Frankle and Michael Carbin.  The lottery ticket hypothesis: Finding sparse, trainable neural networks, 2019.

##  Appendix A Proof of Theorem 1: Bayesian Belief Net and Markov Net
Equivalence

######  Theorem.

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.

###### Proof.

According to $M$ the probability of a setting $\vec{V}$ of its variables is

$\frac{1}{Z}\Pi_{i}\phi_{i}(\vec{V})$

where $\phi_{i}$ are the potentials in $M$, and $Z$ is the partition function,
defined as

$Z=\Sigma_{\vec{V}}\Pi_{i}\phi_{i}(\vec{V})$

We use $DOM(\phi)$ to designate the variables in a potential $\phi$. Because
the nodes and structures of $M$ and $N$ agree, we will refer to the parents,
children, ancestors, and descendants of any node in $M$ to designate the
corresponding nodes in $N$. Likewise we will refer to the input and output
variables of $M$ as those nodes of $N$ that have no parents and no children,
respectively. Because $M$ has treewidth 1, each node of $M$ d-separates its
set of ancestors from its set of descendants and indeed from all other nodes
in $M$. As a result, it is known that the partition function can be computed
efficiently in treewidth-1 Markov networks, for example by the following
recursive procedure $f$ defined below. Let $V_{0}$ be the empty set of
variables, and let $V_{1}$ be the input variables of $M$. Let $Ch(V)$ denote
the children of any set $V$ of variables in $M$, and similaly let $Pa(V)$
denote the parents of $V$. For convenience, when $V$ is a singleton we drop
the set notation and let $V$ denote the variable itself. For all natural
numbers $i\geq 0$:

$f(V_{i})=\Pi_{N\in
Ch(V_{i})}\Sigma_{N=0,1}\Pi_{\phi_{j}:DOM(\phi_{j})\subseteq
V_{i},DOM(\phi_{j})\not\subseteq V_{i-1}}\phi_{j}(V_{i})$ $f(V_{m+1})=1$

where $V_{m}$ is not the full set of variable in $M$ but $V_{m+1}$ is the full
set. Then $Z=f(V_{1})$.

For each variable $v\in\vec{V}$, we can multiply the potentials on the edges
between $v$ and its parents, to get a single potential $\phi_{\\{v,Pa(v)\\}}$
over $\\{v,Pa(v)\\}$. For a given setting of the parents of $v$ in $\vec{V}$,
let $\phi_{v|Pa(v)}$ denote the result of conditioning on this setting of the
parents, and let $\phi_{v,\neg v|Pa(v)}$ denote the result of summing out
variable $v$. Using these product potentials of $M$, and given the method
above for computing $Z$ for a tree-structured Markov network, we can define
the probability of a particular setting $\vec{V}$ as

$P(\vec{V})=\Pi_{v\in\vec{V}}\frac{\phi_{v|Pa(v)}}{\phi_{v,\neg v|Pa(v)}}$

These terms are exactly the terms of the logistic conditional probabilities of
the Bayesian belief network $N$:

$P(\vec{V})=\Pi_{v\in\vec{V}}P(v|Pa(v))$

∎

Note that in general when converting a Bayes net structure to a Markov net
structure, to empower the Markov net to represent any probability distribution
representable by the Bayes net we have to moralize. A corollary of the above
theorem is that in the special case where the Bayes net uses only sigmoid
activations, and its underlying undirected graph is tree-structured,
moralization is not required.

##  Appendix B HMC Trajectories

Suppose the current chain state is $\bm{h}^{(n)}=\bm{\rho}_{n}(0)$. We then
draw a momentum $\bm{\mu}_{n}(0)\sim\mathcal{N}(0,M)$. The HMC trajectories
imply that after $\Delta t$, we have:

$\begin{split}\bm{\mu}_{n}(t+\frac{\Delta t}{2})&=\bm{\mu}_{n}(t)-\frac{\Delta
t}{2}\nabla U(\bm{\rho})\bigg{|}_{\bm{\rho}=\bm{\rho}_{n}(t)}\\\
\bm{\rho}_{n}(t+\Delta t)&=\bm{\rho}_{n}(t)+\Delta
tM^{-1}\bm{\mu}_{n}(t+\frac{\Delta t}{2})\\\ \bm{\mu}_{n}(t+\Delta
t)&=\bm{\mu}_{n}(t+\frac{\Delta t}{2})-\frac{\Delta t}{2}\nabla
U(\bm{\rho})\bigg{|}_{\bm{\rho}=\bm{\rho}_{n}(t+\Delta t)}.\end{split}$ (B.1)

We may then apply these equations to $\bm{\rho}_{n}(0)$ and $\bm{\mu}_{n}(0)$
$L$ times to get $\bm{\rho}_{n}(L\Delta t)$ and $\bm{\mu}_{n}(L\Delta t)$.
Thus, the transition from $\bm{h}_{(n)}=\bm{\rho}_{n}$ to the next state
$\bm{h}^{(n+1)}$ is given by:

$\bm{h}^{(n+1)}\Big{|}(\bm{h}^{(n)}=\bm{\rho}_{n}(0))=\begin{cases}\bm{\rho}_{n}(L\Delta
t)&\text{with probability $\alpha(\bm{\rho}_{n}(0),\bm{\rho}_{n}(L\Delta
t))$}\\\ \bm{\rho}_{n}(0)&\text{otherwise}\end{cases}$ (B.2)

where

$\alpha(\bm{\rho}_{n}(0),\bm{\rho}_{n}(L\Delta
t))=\min(1,\exp(H(\bm{\rho}_{n}(0),\bm{\mu}_{n}(0))-H(\bm{\rho}_{n}(L\Delta
t),\bm{\mu}_{n}(L\Delta t)))).$ (B.3)

##  Appendix C Algorithm 1

Algorithm 1 CD-HMC Training for the Continuous Bernoulli Belief Network

1:Initialized $\bm{h}^{(0)}$,
$\bm{W}^{(0)}=\\{\bm{W}_{i}^{(0)},i=1,2,...,L\\}$,
$\bm{b}^{(0)}=\\{\bm{b_{i}}^{(0)},i=1,2,...,L\\}$

2:$\bm{W}^{(n)},\bm{b}^{(n)}$ when loss converges.

3:

4:procedure Burn-in($N$)

5: for $i\leftarrow 0$ to $N-1$ do

6: $\bm{h}^{(i+1)}\leftarrow$
Sampling$(\bm{h}^{(i)},\bm{W}^{(0)}$,$\bm{b}^{(0)})$

7: end for

8:end procedure

9:procedure Training($M$)

10: for $i\leftarrow 0$ to $M-1$ do

11: $\bm{W}^{(i+1)},\bm{b}^{(i+1)}\leftarrow$ Weight-
updating$(\bm{h}^{(N+iK)},\bm{W}^{(i)},\bm{b}^{(i)})$

12: for $j\leftarrow 0$ to $K-1$ do

13: $\bm{h}^{(N+ik+j+1)}\leftarrow$
Sampling$(\bm{h}^{(N+iK+j)},\bm{W}^{(i+1)},\bm{b}^{(i+1)})$

14: end for

15: end for

16:end procedure