# Explicitly Bayesian Regularizations in Deep Learning

Xinjie Lan, Kenneth E. Barner
Department of Electrical and Computer Engineering
University of Delaware
Newark, DE 19713
<EMAIL_ADDRESS>

###### Abstract

Generalization is essential for deep learning. In contrast to previous works
claiming that Deep Neural Networks (DNNs) have an implicit regularization
implemented by the stochastic gradient descent, we demonstrate explicitly
Bayesian regularizations in a specific category of DNNs, i.e., Convolutional
Neural Networks (CNNs). First, we introduce a novel probabilistic
representation for the hidden layers of CNNs and demonstrate that CNNs
correspond to Bayesian networks with the serial connection. Furthermore, we
show that the hidden layers close to the input formulate prior distributions,
thus CNNs have explicitly Bayesian regularizations based on the Bayesian
regularization theory. In addition, we clarify two recently observed empirical
phenomena that are inconsistent with traditional theories of generalization.
Finally, we validate the proposed theory on a synthetic dataset.

##  1 Introduction

Generalization is a fundamental problem of deep learning theory. Since
traditional generalization theories cannot clarify over-parameterized Deep
Neural Networks (DNNs), especially in some cases, e.g., the random label
classification [35, 20, 21], most recent works ascribe the generalization
performance of DNNs to the Stochastic Gradient Descent (SGD) optimization
implementing an implicit regularization, i.e., a bias towards models with low
complexity [11, 3, 7, 30, 1].

Nevertheless, the implicit regularization explanations have lots of
controversies. For example, Keskar et al. claim that the sharp minima learned
by SGD lead to poor generalization [11], but Dinh et al. demonstrate that
sharp minima can generalize DNNs well [3]. Sourdry et al. argue that SGD
learns the minimum norm as the implicit regularization [7, 30], but Arora et
al. doubt that the implicit regularization cannot be captured by the
mathematical norm [1, 20].

In the context of regularization theories, Bayesian theory indicates that an
explicit regularization can be formulated as a prior distribution [31].
However, it is scarce that explaining the regularization of DNNs from the
Bayesian viewpoint, because the extremely complicated architecture of DNNs and
the difficulty of probabilistic modeling high-dimensional activations impede
establishing a comprehensive probabilistic explanation for the architecture of
DNNs.

Numerous efforts have been devoted to explain the architecture of DNNs from a
probabilistic perspective [29, 19, 25, 32, 26, 16, 18, 4, 24]. However,
existing probabilistic explanations have several limitations. The Boltzmann
explanations only clarify the restricted Boltzmann machine [29, 19]. The
mixture model explanations require too strict assumptions [25, 32, 26]. The
gaussian process explanations cannot clarify the hierarchical representation,
an essential of deep learning [18]. In summary, existing explanations cannot
establish a solid probabilistic foundation for investigating the
regularization of DNNs from the Bayesian perspective.

To establish a solid probabilistic foundation for the study of Bayesian
regularization in deep learning, we introduce a novel probabilistic
representation for the hidden layers of Convolutional Neural Networks (CNNs)
[13] by extending the Boltzmann explanation [19]. More specifically, a fully
connected hidden layer formulate a multivariate discrete Boltzmann
distribution [6], and its energy function is equal to the negative activations
of the layer. In addition, a convolutional layer formulate a specific
Boltzmann explanation, the Markov Random Fields (MRFs) model [5, 33], and its
potential functions are determined by the linear convolutional filters.

Based on the proposed probabilistic representation, we demonstrate that CNNs
formulate a directed acyclic graph, i.e., a Bayesian network [23, 12], and SGD
learns the causality between the given training dataset and training labels.
In the context of Bayesian networks, the causal links are described by
conditional probabilities, thus the hidden layers close to the training
dataset formulate prior distributions to describe the cause and the output
layer formulates a probability to describe the expected effect. As a result,
we conclude that CNNs have explicitly Bayesian regularizations based on the
correspondence between prior distributions and regularizations.

We clarify two recently observed empirical phenomena that are inconsistent
with traditional theories of generalization. First, increasing the number of
hidden units can lead to a decrease in generalization error [21]. Bayesian
theory indicates that CNNs can use more hidden units to generate a prior
distribution including more prior knowledge, thereby achieving the lower
generalization error. Second, CNNs with good generalization performance on
real labels achieve zero training error but very high generalization error on
random labels [35]. A Bayesian network can be simplified as two components,
i.e., a prior distribution and a likelihood distribution, and regularization
merely corresponds to the prior distribution. Given the same CNN architecture
and the training dataset, CNNs should learn the same prior distributions to
formulate the same regularization, different generalization errors result from
the likelihood distribution but not related to the prior distribution, i.e.,
regularization.

##  2 Background

###  2.1 Boltzmann distribution

In mathematics, the Boltzmann distribution (a.k.a, the Gibbs distribution or
the renormalization group) formulates the probability of a state, e.g.,
$p(\boldsymbol{x};\boldsymbol{\theta})$, as a function of the state’s energy
[6, 19, 34, 15].

${\textstyle
p(\boldsymbol{x};\boldsymbol{\theta})=\frac{1}{Z(\boldsymbol{\theta})}\text{exp}[-E(\boldsymbol{x;\theta})]\text{,}}$
(1)

where $E(\boldsymbol{x;\theta})$ is the energy function,
$Z(\boldsymbol{\theta})=\sum_{\boldsymbol{x}}{\text{exp}[-E(\boldsymbol{x;\theta})]}$
is the partition function, and $\boldsymbol{\theta}$ denote all parameters. A
classical example of the Boltzmann distribution in machine learning is the
Restricted Boltzmann Machine (RBM) [29].

The Boltzmann distribution can be reformulated as various probabilistic models
through redefining the energy function $E(\boldsymbol{x;\theta})$. In the
context of probabilistic graphical model, a fundamental Boltzmann distribution
is the Markov Random Fields (MRFs) and its energy function is equivalent to
the summation of multiple potential functions, i.e., ${\scriptstyle
E(\boldsymbol{x;\theta})=-\sum_{k}f_{k}(\boldsymbol{x};\boldsymbol{\theta}_{k})}$
[5, 33]. Moreover, assuming a single potential function defines an expert,
i.e.,
${\scriptstyle\boldsymbol{F}_{k}=\frac{1}{Z(\boldsymbol{\theta}_{k})}\text{exp}[-f_{k}(\boldsymbol{x};\boldsymbol{\theta}_{k})]}$,
the Boltzmann distribution becomes a Product of Experts (PoE) model , i.e.,
${\scriptstyle
p(\boldsymbol{x};\boldsymbol{\theta})=\frac{1}{Z(\boldsymbol{\theta})}\prod_{k}\boldsymbol{F}_{k}}$,
where ${\scriptstyle
Z(\boldsymbol{\theta})=\prod_{k}Z(\boldsymbol{\theta}_{k})}$ [9, 10].

###  2.2 Bayesian networks

A Bayesian network represents the causal probabilistic relationship among a
set of random variables based on their conditional distributions. It consists
of two major parts: a directed acyclic graph and a set of conditional
probability distributions [23, 12]. For example,
$\boldsymbol{A}\rightarrow\boldsymbol{B}$ forms a directed acyclic graph to
show the possible causality between two random variables $\boldsymbol{A}$ and
$\boldsymbol{B}$, which can be formulated as
$P(\boldsymbol{A},\boldsymbol{B})=P(\boldsymbol{B|A})P(\boldsymbol{A})$, where
$P(\boldsymbol{A})$ formulates a prior distribution to expresses the
uncertainty of the cause and $P(\boldsymbol{B|A})$ formulates the uncertainty
of the effect $\boldsymbol{B}$ given the cause $\boldsymbol{A}$. In
particular, since the prior distribution $P(\boldsymbol{A})$ expresses the
prior belief of the cause that lead to the effect, $P(\boldsymbol{A})$
corresponds to the regularization in terms of the Bayesian regularization
theory.

In this paper, we show that the entire architecture of CNNs forms a Bayesian
network and the hidden layers close to the input formulate prior distributions
to express the prior belief leading to the expected effect (i.e., labels),
thus these hidden layers form explicitly Bayesian regularizations.

##  3 Main results

###  3.1 Preliminaries

We assume that $\boldsymbol{X}$ and $\boldsymbol{Y}$ are two random variables
and
$P_{\boldsymbol{\theta}}(\boldsymbol{X},\boldsymbol{Y})=P(\boldsymbol{Y|X})P(\boldsymbol{X})$
is an unknown joint distribution between $\boldsymbol{X}$ and
$\boldsymbol{Y}$, where $P(\boldsymbol{X})$ describes the prior knowledge of
$\boldsymbol{X}$, $P(\boldsymbol{Y|X})$ describes the statistical connection
between $\boldsymbol{X}$ and $\boldsymbol{Y}$, and $\boldsymbol{\theta}$
indicate the parameters of
$P_{\boldsymbol{\theta}}(\boldsymbol{X},\boldsymbol{Y})$.

A dataset
$\boldsymbol{\mathcal{D}}=\\{(\boldsymbol{x}^{j},\boldsymbol{y}^{j})|\boldsymbol{x}^{j}\in{\mathbb{R}}^{M},\boldsymbol{y}^{j}\in{\mathbb{R}}^{L}\\}_{j=1}^{J}$
is composed of i.i.d. samples generated from
$P_{\theta}(\boldsymbol{X},\boldsymbol{Y})$. In addition, a CNN with $I$
hidden layers is denoted as
$\text{CNN}=\\{\boldsymbol{x;f_{1};...;f_{I};f_{Y}}\\}$ and trained by
$\boldsymbol{\mathcal{D}}$, where
$(\boldsymbol{x},\boldsymbol{y})\in\boldsymbol{\mathcal{D}}$ are the input of
the CNN and the corresponding training label, and $\boldsymbol{f_{Y}}$ is an
estimation of the true distribution $P(\boldsymbol{Y|X})$. As a result,
$\boldsymbol{x}\sim P(\boldsymbol{X})$ and $\boldsymbol{f_{Y}}\approx
P(\boldsymbol{Y|X})$.

###  3.2 The probabilistic representation for the hidden layers of CNNs

In this section, we extend the Boltzmann explanations to CNNs and introduce a
novel probabilistic representation for the hidden layers of CNNs. The detailed
derivation is included in Appendix 6.1. Since most CNNs only contain two
categories of hidden layers, namely fully connected layers and convolutional
layers with non-linear operators, e.g., ReLU and pooling, the proposed
probabilistic representation can be summarized as the following propositions.

Proposition 1: A fully connected layer formulates a multivariate discrete
Boltzmann distribution and the energy function of which is equivalent to the
negative activations of the fully connected layer.

A fully connected layer $\boldsymbol{f_{i}}=\\{f_{in}\\}_{n=1}^{N}$ with $N$
neurons formulates a multivariate discrete Boltzmann distribution to describe
the probability of $N$ features defined by the neurons occurring in the input.
Their energy functions are equal to the negative of their respective neurons,
i.e., $E_{f_{in}}=-f_{in}(\boldsymbol{f_{j}})$.

$P(\boldsymbol{F_{i}})=\\{P(f_{in})=\frac{1}{Z_{\boldsymbol{F_{i}}}}\text{exp}[f_{in}(\boldsymbol{f_{j}})]\\}_{n=1}^{N}$
(2)

where the previous layer $\boldsymbol{f_{j}}=\\{f_{jm}\\}_{m=1}^{M}$ is the
input, $f_{in}(\boldsymbol{f_{j}})=\sigma[g_{in}(\boldsymbol{f_{j}})]$ is the
activation of the $n$th neuron,
${g_{in}(\boldsymbol{f_{j}})=\sum_{m=1}^{M}\alpha_{mn}\cdot f_{jm}+b_{in}}$ is
a linear filter with $\alpha_{mn}$ being the weights and $b_{in}$ being the
bias, $\sigma(\cdot)$ is an activation function, and
${Z_{\boldsymbol{F}}=\sum_{n=1}^{N}\text{exp}[f_{in}(\boldsymbol{f_{j}})]}$ is
the partition function.

Proposition 2: A convolutional layer formulates a specific Boltzmann
distribution, i.e., the MRF model, and its energy function is equal to the
negative summation of all the convolutional channels.

A convolutional layer with $N$ convolutional channels, i.e.,
$\boldsymbol{f_{i}}=\\{\boldsymbol{f^{n}_{i}}=\sigma(\varphi_{i}^{n}(\boldsymbol{f_{j}}))\\}_{n=1}^{N}$,
formulates a specific Boltzmann distribution, namely a MRF model.

$P(\boldsymbol{F_{i}})=\frac{1}{Z_{\boldsymbol{F_{i}}}}\text{exp}[\sum_{n=1}^{N}\varphi_{i}^{n}(\boldsymbol{f_{j}})]$
(3)

where $\boldsymbol{f_{j}}=\\{\boldsymbol{f_{j}^{q}}\\}_{q=1}^{Q}$ is the input
with $Q$ channels,
$\varphi_{i}^{n}(\boldsymbol{f_{j}})=\sum_{q=1}^{Q}\boldsymbol{S^{(n,q)}}\circ\boldsymbol{f_{j}^{q}}+b_{i}^{n}\cdot\boldsymbol{1})$
is the $n$th linear convolutional channel corresponding to the potential
function of a MRF model, in which $\boldsymbol{S^{(n,q)}}$ is the $n$th
convolutional filter applying to the $q$th channel, $b_{i}^{n}$ is the bias,
$\sigma(\cdot)$ is an non-linear operator, and the partition function
$Z_{\boldsymbol{F_{i}}}=\int_{\boldsymbol{f_{j}}}\text{exp}[\sum_{n=1}^{N}\varphi_{i}^{n}(\boldsymbol{f_{j}})]d{\boldsymbol{f_{j}}}$.
In particular, based on the equivalence between SGD and the first order
approximation [2], Equation 3 indicates that the statistical property of a
convolutional layer is determined by the linear convolutional operations
$\\{\varphi_{i}^{n}(\boldsymbol{f_{j}})\\}_{n=1}^{N}$.

In summary, the Boltzmann distribution establishes a novel connection between
the hidden layers of CNNs and probabilistic models through defining the energy
function of a Boltzmann distribution as a function of the activations in a
hidden layer. It is important to note that the energy function is a sufficient
statistics of a Boltzmann distribution, i.e., the later is entirely dependent
on the former [22]. As a result, we can specify a unique Boltzmann
distribution given a hidden layer of CNNs.

Moreover, the architecture of hidden layers determines the statistical
property of the corresponding Boltzmann distribution. For example, a fully
connected layer assumes all the input data being connected to each other, but
a convolutional layer assumes local connectivity, i.e., the input data only
depends on its neighbors. Therefore, their corresponding Boltzmann
distributions have different statistical properties: Equation 2 uses the
linear filter $g_{in}(\boldsymbol{f_{j}})$ to describe the fully connected
features, and Equation 3 uses the potential functions of MRF to model the
local connected features.

###  3.3 Explicitly Bayesian regularizations in CNNs

Proposition 3: The entire architecture of CNNs can be explained as a Bayesian
network.

Given a $\text{CNN}=\\{\boldsymbol{x;f_{1};...;f_{I};f_{Y}}\\}$, the output of
the hidden layer $\boldsymbol{f_{i}}$ is the input of the next hidden layer
$\boldsymbol{f_{i+1}}$, thus the CNN forms a serial connection in terms of
Bayesian networks.

$\boldsymbol{f_{1}\rightarrow\cdots\rightarrow f_{I}\rightarrow f_{Y}}$ (4)

Since the proposed probabilistic representation specifies the distribution of
the hidden layers of CNNs, the entire architecture of the CNN corresponds to a
Bayesian network as

${\textstyle
P(\boldsymbol{F_{Y};F_{I};\cdots;F_{1}|X})=P(\boldsymbol{F_{Y}|F_{I}})\cdot...\cdot
P(\boldsymbol{F_{i+1}|F_{i}})\cdot...\cdot P(\boldsymbol{F_{1}|X})}$ (5)

where $P(\boldsymbol{F_{1}|X})$ formulates a prior Boltzmann distribution to
describe the cause, the $P(\boldsymbol{F_{Y}|F_{I}})$ formulates a Boltzmann
distribution to express the expected effect, and the intermediate hidden
layers formulate multiple conditional distributions to transfer the influence
of the cause to the effect.

Previous works demonstrate that the output of CNNs estimates a Bayesian
posterior distribution $P(\boldsymbol{F_{Y}|X})$ [27, 36]. In the context of
Bayesian networks, $P(\boldsymbol{F_{Y}|X})$ can be explained as a variable
inference given the Bayesian networks corresponding to CNNs.

$P(\boldsymbol{F_{Y}|X})=\int_{\boldsymbol{F_{I}}}\cdots\int_{\boldsymbol{F_{1}}}P(\boldsymbol{F_{Y};F_{I}..;F_{1}|X})d{\boldsymbol{F_{I}}}\cdots
d{\boldsymbol{F_{1}}}$ (6)

Since CNNs form a Bayesian network with the serial connection and the
distribution of each layer is a Boltzmann distribution, we prove that
$P(\boldsymbol{F_{Y}|X})$ is also a Boltzmann distribution and the energy
function of which is also equivalent to the functionality of the entire CNN in
Appendix 6.2.

Given the training dataset $\boldsymbol{\mathcal{D}}$ and a predefined loss
function, e.g., the cross entropy, SGD can be explained as a parameter
learning procedure for $P(\boldsymbol{F_{Y}|X})$. Notably, since CNNs form a
serial connection indicated by Equation 4, both the Bayesian posterior
distribution $P(\boldsymbol{F_{Y}|X})$ and the entire Bayesian network
$P(\boldsymbol{F_{Y};F_{I}..;F_{1}|X})$ correspond to the same architecture of
CNNs, thus SGD can learn the parameters of the entire Bayesian network by
merely optimizing $P(\boldsymbol{F_{Y}|X})$.

Proposition 4: CNNs have explicitly Bayesian regularizations

In the Bayesian approach, regularization is achieved by specifying a prior
distribution over the parameters and subsequently averaging over the posterior
distribution [31]. Therefore, the Bayesian networks explanation provides a
novel perspective to clarify the generalization of CNNs.

Above all, CNNs have explicitly Bayesian regularizations. Since we explain the
entire architecture of CNNs as a Bayesian network with the serial connection
indicated by Equation 4, the hidden layers close to the input should formulate
prior distributions, thereby representing explicit regularizations. In
contrast to previous works ascribing the generalization of CNNs to SGD
implementing implicit regularization, Proposition 3 shows that SGD merely
serves as a parameter learning for $P(\boldsymbol{F_{Y}|X})$ and explicit
Bayesian regularizations are defined by the hidden layers corresponding to
prior distributions.

Moreover, the Bayesian regularization explanation can clarify several recently
observed empirical phenomena that are inconsistent with traditional theories
of generalization [20]. First, increasing the number of hidden units of CNNs,
i.e., increasing the number of parameters, can lead to a decrease in
generalization error [21]. Bayesian theory indicates that more complex prior
distributions can include more prior belief of the training dataset. As a
result, CNNs can use many hidden units to generate much powerful prior
distributions for regularizing the entire network, thereby guaranteeing the
generalization performance though they are over-parametrized.

Second, CNNs with good generalization performance on real labels achieve zero
training error but very high generalization error on random labels [35].
Bayesian theory enables us to utilize the ‘separate and conquer’ rule to
understand the special case. In terms of Bayesian theory, a Bayesian network
can be simplified as two components, namely a prior distribution and a
likelihood distribution. Given the same CNN architecture and the same training
dataset with different labels, e.g., real labels and random labels, the CNN
should formulate the same regularizations, because regularization merely
corresponds the prior distribution. Different generalization errors for
different labels result from the learned likelihood distribution but not
related to the prior distribution.

##  4 Simulations

In this section, we use a synthetic dataset to validate explicitly Bayesian
regularizations in CNNs, and demonstrate that the proposed Bayesian
regularization can explain two noteworthy empirical phenomena that are
inconsistent with traditional theories of generalization.

###  4.1 Setup

Figure 1:  The first row shows five synthetic images of handwritten digits,
the second row shows their respective empirical distributions, and the red
curve indicates the true distribution $\mathcal{N}(0,1024)$.

To validate explicitly Bayesian regularizations in CNNs, we need to firstly
validate the probabilistic representation for the hidden layers of CNNs.
However, since the distributions of benchmark datasets are unknown, it is
impossible to directly use them to validate the proposed probabilistic
representation. Alternatively, we generate a synthetic dataset obeying the
Gaussian distribution $\mathcal{N}(0,1024)$ based on the NIST dataset of
handwritten digits [14]. It consists of 20,000 $32\times 32$ grayscale images
in 10 classes (digits from 0 to 9), and each class has 1,000 training images
and 1,000 testing images. The method for generating the synthetic dataset is
included in Appendix 6.3.

Most image datasets satisfy two fundamental assumptions: stationary and Markov
[17]. The former means that the distribution is independent on the location,
i.e., the distribution of pixels at different locations are identically
distributed. The later indicates that the distribution of a single pixel is
independent on the other pixels given its neighbors. The two assumptions allow
us to use an empirical distribution to simulate the true distribution of the
synthetic dataset. Figure 1 shows five synthetic images and their perspective
empirical distributions. It shows that their empirical distributions precisely
simulate the true distribution of synthetic images, i.e.,
$\mathcal{N}(0,1024)$.

We design three simple by comprehensive CNNs for classifying the synthetic
dataset. Every CNN has two convolutional hidden layers, two max pooling
layers, and a fully connected hidden layer. The only difference within these
CNNs is that the first convolutional layer $\boldsymbol{f_{1}}$ of different
CNNs has different number of convolutional channels.

Table 1: The architectures of CNNs for classifying the synthetic dataset

R.V. | Layer | Description | CNN1 | CNN2 | CNN3
---|---|---|---|---|---
$\boldsymbol{X}$ | $\boldsymbol{x}$ | Input | $32\times 32\times 1$ | $32\times 32\times 1$ | $32\times 32\times 1$
$\boldsymbol{F_{1}}$ | $\boldsymbol{f_{1}}$ |  Conv ($3\times 3$) +ReLU | $30\times 30\times\boldsymbol{4}$ | $30\times 30\times\boldsymbol{12}$ | $30\times 30\times\boldsymbol{20}$
$\boldsymbol{f_{2}}$ | Maxpool | $15\times 15\times{4}$ | $15\times 15\times 12$ | $15\times 15\times 20$
$\boldsymbol{F_{2}}$ | $\boldsymbol{f_{3}}$ |  Conv ($5\times 5$) + ReLU | $11\times 11\times 20$ | $11\times 11\times 20$ | $11\times 11\times 20$
$\boldsymbol{f_{4}}$ | Maxpool | $5\times 5\times{20}$ | $5\times 5\times 20$ | $5\times 5\times 20$
$\boldsymbol{F_{3}}$ | $\boldsymbol{f_{5}}$ | Fully connected | $1\times 1\times 20$ | $1\times 1\times 20$ | $1\times 1\times 20$
$\boldsymbol{F_{Y}}$ | $\boldsymbol{f_{Y}}$ | Output(softmax) | $1\times 1\times 10$ | $1\times 1\times 10$ | $1\times 1\times 10$

  * •

R.V. is the random variable of the hidden layer(s).

Based on the proposed probabilistic representation, we identify the
distribution of each hidden layer in the above
$\text{CNNs}=\\{\boldsymbol{x;f_{1};...;f_{I};f_{Y}}\\}$. First,
$\boldsymbol{f_{1}}$ formulates a specific Boltzmann distribution, namely, a
MRF model. Since the max pooling layer $\boldsymbol{f_{2}}$ can be viewed as
an non-linear differentiable operation as ReLU, the distribution of
$\boldsymbol{f_{2}}$ can be approximated by $\boldsymbol{f_{1}}$ based on the
equivalence between SGD and the first order approximation [2]. Therefore, we
can use a single random variable $\boldsymbol{F_{1}}$ to represent the
statistical properties of both $\boldsymbol{f_{1}}$ and $\boldsymbol{f_{2}}$.
Similarly, we use $\boldsymbol{F_{2}}$ to represent the statistical properties
of both $\boldsymbol{f_{3}}$ and $\boldsymbol{f_{4}}$. As a result,
$P(\boldsymbol{F_{1}})$ and $P(\boldsymbol{F_{2}|F_{1}})$ can be formulated as

$\begin{split}P(\boldsymbol{F_{1}})=\frac{1}{Z_{\boldsymbol{F_{1}}}}\text{exp}[\sum_{n=1}^{N}\varphi_{1}^{n}(\boldsymbol{x})]\text{
and
}P(\boldsymbol{F_{2}|F_{1}})=\frac{1}{Z_{\boldsymbol{F_{2}}}}\text{exp}[\sum_{k=1}^{K}\varphi_{3}^{k}(\boldsymbol{f_{2}})]\end{split}$
(7)

where $\boldsymbol{x}$ only has a single channel since it denotes a gray scale
image, and
$\varphi_{1}^{n}(\boldsymbol{x})=\boldsymbol{S_{1}^{n}}\circ\boldsymbol{x}+b_{1}^{n}\cdot\boldsymbol{1}$
is the $n$th linear convolutional channel in $\boldsymbol{f_{1}}$. Similarly,
$\boldsymbol{f_{2}}=\\{\boldsymbol{f^{n}_{2}}\\}_{n=1}^{N}$ has $N$ channels,
and
$\varphi_{3}^{k}(\boldsymbol{f_{2}})=\sum_{n=1}^{N}\boldsymbol{S_{3}^{(k,n)}}\circ\boldsymbol{f_{2}^{k}}+b_{2}^{k}\cdot\boldsymbol{1})$
is the $k$th linear convolutional channel in $\boldsymbol{f_{3}}$. The two
partition functions are
$Z_{\boldsymbol{F_{1}}}=\int_{\boldsymbol{x}}\text{exp}[\sum_{n=1}^{N}\varphi_{1}^{n}(\boldsymbol{x})]d{\boldsymbol{x}}$
and
$Z_{\boldsymbol{F_{2}}}=\int_{\boldsymbol{f_{2}}}\text{exp}[\sum_{k=1}^{K}\varphi_{3}^{k}(\boldsymbol{f_{2}})]d{\boldsymbol{f_{2}}}$.

Although we derive $P(\boldsymbol{F_{1}})$ and $P(\boldsymbol{F_{2}|F_{1}})$,
it is still hard to calculate $P(\boldsymbol{F_{1}})$ and
$P(\boldsymbol{F_{2}|F_{1}})$ because $Z_{\boldsymbol{F_{1}}}$ and
$Z_{\boldsymbol{F_{2}}}$ are intractable for the high dimensional datasets
$\boldsymbol{x}$ and $\boldsymbol{f_{2}}$. Alternatively, since MRF models
satisfy the two assumptions (i.e., stationary and Markov) and the energy
function is a sufficient statistics of a Boltzmann distribution, we can use
the empirical distributions of their respective energy functions
(${\scriptstyle\boldsymbol{E}_{\boldsymbol{F_{1}}}=-\sum_{n=1}^{N}\varphi_{1}^{n}(\boldsymbol{x})}$
and
${\scriptstyle\boldsymbol{E}_{\boldsymbol{F_{2}}}=-\sum_{k=1}^{K}\varphi_{2}^{k}(\boldsymbol{f_{2}})}$)
to estimate $P(\boldsymbol{F_{1}})$ and $P(\boldsymbol{F_{2}|F_{1}})$.
Specifically, the two empirical distributions can be expressed as

$\begin{split}\hat{P}[\boldsymbol{F_{1}}(B_{i})]=\frac{1}{C_{\boldsymbol{E}_{\boldsymbol{F_{1}}}}}[\sum_{nn=1}^{C_{\boldsymbol{E}_{\boldsymbol{F_{1}}}}}B_{i}\leq\boldsymbol{1}[\boldsymbol{E}_{\boldsymbol{F_{1}}}(nn)]\leq
B_{i+1}]\\\
\hat{P}[\boldsymbol{F_{2}|F_{1}}(B_{i})]=\frac{1}{C_{\boldsymbol{E}_{\boldsymbol{F_{2}}}}}[\sum_{mm=1}^{C_{\boldsymbol{E}_{\boldsymbol{F_{2}}}}}B_{i}\leq\boldsymbol{1}[\boldsymbol{E}_{\boldsymbol{F_{2}}}(mm)]\leq
B_{i+1}]\end{split}$ (8)

where $C_{\boldsymbol{E}_{\boldsymbol{F_{1}}}}$ and
$C_{\boldsymbol{E}_{\boldsymbol{F_{2}}}}$ are the number of elements in
$\boldsymbol{E}_{\boldsymbol{F_{1}}}$ and
$\boldsymbol{E}_{\boldsymbol{F_{2}}}$, respectively. $B_{i}$ and $B_{i+1}$ are
the pre-defined edges, and $\boldsymbol{1}(\cdot)$ is the indicator function.

Based on the Boltzmann explanation for fully connected layers (Equation 2), we
can formulate the distribution of $\boldsymbol{f_{5}}$ and
$\boldsymbol{f_{Y}}$ as follows.

$\begin{split}P(\boldsymbol{F_{3}|F_{2}})=\\{P(f_{5m})&=\frac{1}{Z_{\boldsymbol{F_{5}}}}\text{exp}[f_{5m}(\boldsymbol{f_{4}})]\\}_{m=1}^{20}\\\
P(\boldsymbol{F_{Y}|F_{3}})=\\{P(f_{yl})&=\frac{1}{Z_{\boldsymbol{F_{Y}}}}\text{exp}[f_{yl}(\boldsymbol{f_{5}})]\\}_{l=1}^{10}\end{split}$
(9)

where the flattening $\boldsymbol{f_{4}}$ is the input of
$P(\boldsymbol{F_{3}})$. In addition, the two partition functions are defined
as
$Z_{\boldsymbol{F_{5}}}=\sum_{m=1}^{20}\text{exp}[f_{5m}(\boldsymbol{f_{4}})]$
and
$Z_{\boldsymbol{F_{Y}}}=\sum_{l=1}^{10}\text{exp}[f_{yl}(\boldsymbol{f_{5}})]$.

Based on the Bayesian networks explanation for CNNs, the above CNNs form a
Bayesian network
$P(\boldsymbol{F_{Y};F_{3};F_{2};F_{1}})=P(\boldsymbol{F_{Y}|F_{3}})\cdot
P(\boldsymbol{F_{3}|F_{2}})\cdot P(\boldsymbol{F_{2}|F_{1}})\cdot
P(\boldsymbol{F_{1}})$, where $P(\boldsymbol{F_{1}})$ formulates a prior
distribution of the training dataset to describe the cause, i.e.,
$P(\boldsymbol{F_{1}})\approx P(\boldsymbol{X})$,
$P(\boldsymbol{F_{Y}|F_{3}})$ formulates the uncertainty of the expected
effect, and the intermediate layers formulate $P(\boldsymbol{F_{2}|F_{1}})$
and $P(\boldsymbol{F_{3}|F_{2}})$ to transfer the influence of the prior
belief to the expected effect.

###  4.2 Validating explicitly Bayesian regularizations in CNNs

Since $\boldsymbol{x}$ is a sample generated from $P(\boldsymbol{X})$, i.e.,
$\boldsymbol{x}\sim P(\boldsymbol{X})=\mathcal{N}(0,1024)$, we can calculate
the distance between $P(\boldsymbol{X})$ and $\hat{P}(\boldsymbol{F_{1}})$,
i.e., $\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{F_{1}})]$, to examine
whether $\boldsymbol{f_{1}}$ formulates a prior distribution of the input for
validating if explicitly Bayesian regularizations exist in CNNs.

We choose CNN2 from Table 1 to classify the synthetic dataset. After CNN2 is
well trained, i.e., the training error of CNN2 becomes zero, we randomly
choose a synthetic image as the input of CNN2 from the testing dataset, and
derive $\hat{P}(\boldsymbol{F_{1}})$, $\hat{P}(\boldsymbol{F_{2}|F_{1}})$,
$P(\boldsymbol{F_{3}|F_{2}})$, and $P(\boldsymbol{F_{Y}|F_{3}})$. Finally, all
the distributions are shown in Figure 2. We can find that
$\hat{P}(\boldsymbol{F_{1}})$ is indeed close to $P(\boldsymbol{X})$, because
the KL divergence is samll, i.e.,
$\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{F_{1}})]=0.49$.

Figure 2:  (A) the synthetic image $\boldsymbol{x}$ is the input of CNN2, (B)
$\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{x})]=0.53$, where the red
curve indicates the truly prior distribution
$P(\boldsymbol{X})=\mathcal{N}(0,1024)$ and $\hat{P}(\boldsymbol{x})$ is the
empirical distribution of $\boldsymbol{x}$, (C)
$\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{F_{1}})]=0.49$, where
$\hat{P}(\boldsymbol{F_{1}})$ is the empirical distribution of
$\boldsymbol{f_{1}}$ and $\boldsymbol{f_{2}}$, (D)
$\hat{P}(\boldsymbol{F_{2}|F_{1}})$ is the empirical distribution of
$\boldsymbol{f_{3}}$ and $\boldsymbol{f_{4}}$, (E) the distribution
${P}(\boldsymbol{F_{3}|F_{2}})$, and (F) the distribution
${P}(\boldsymbol{F_{Y}|F_{3}})$.

This experiment validates explicitly Bayesian regularizations in CNNs. Since
we can theoretically prove CNN2 formulating a Bayesian network
$P(\boldsymbol{F_{Y};F_{3};F_{2};F_{1}})$ and empirically show
$\hat{P}(\boldsymbol{F_{1}})$ modeling the prior distribution
$P(\boldsymbol{X})$, CNN2 has an explicitly Bayesian regularization defined by
the first convolutional layer $\boldsymbol{f_{1}}$. In addition, the
experiment validates the proposed probabilistic representation for the hidden
layers of CNNs, especially the convolutional layer.

###  4.3 The effect of increasing the number of hidden units on the
generalization performance

Based on the proposed Bayesian networks explanation and the experiment
presented in Section 4.2, we verify that $\boldsymbol{f_{1}}$ formulates a
prior distribution as an explicit Bayesian regularization. Since the only
difference within the three CNNs is that the first convolutional layer
$\boldsymbol{f_{1}}$ of different CNNs has different number of convolutional
channels, we can use the three CNNs to examine the effect of increasing hidden
units on the generalization performance of CNNs. We use
$\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{F_{1}})]$ to quantify the
effect of increasing the number of hidden units on the prior distribution
$\hat{P}(\boldsymbol{F_{1}})$, i.e., the Bayesian regularization. We average
$\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{F_{1}})]$ overall all
testing images and use the testing error to quantify the generalization
performance over 30 training epochs.

Figure 3 shows that increasing the number of hidden units of CNNs has a
negative correlation with the generalization error, which validates the
Bayesian regularization explanation (Proposition 4). Specifically, Figure 3
(Left) shows that $\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{F_{1}})]$
is decreasing as we increase the number of hidden units. It means that
$\boldsymbol{f_{1}}$ formulates a better Bayesian regularization as increasing
the number of hidden unit. As a result, the corresponding CNNs should achieve
the lower generalization error, which is demonstrated by Figure 3 (Right).

Figure 3:  (Left) The average
$\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{F_{1}})]$ overall all
testing images. (Right) The generalization performance of the three CNNs
represented by the testing error. The number behind the three CNNs in the
legend denotes the number of hidden units of $\boldsymbol{f_{1}}$.

###  4.4 The generalization performance of CNNs on random labels

Figure 4:  The distribution of the hidden layers in CNN3 given a synthetic
testing image with the random label. The notation is the same as Figure 2,
e.g., $\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{x})]=0.66$, and
$\text{KL}[P(\boldsymbol{X})||\hat{P}(\boldsymbol{F_{1}})]=0.47$.

Similar to the notable experiment presented in [35], we use CNN3 to classify
the synthetic dataset with random labels. Figure 5 shows CNN3 achieving zero
training error but very high testing error. We visualize the distribution of
hidden layers in CNN3 given a testing image with a random label, and Figure
4(C) shows that CNN3 still can learn an accurate prior distribution
$\hat{P}(\boldsymbol{F_{1}})$. In this sense, CNN3 still achieves good
generalization performance on random labels.

We conjecture that CNNs cannot classify random labels. In terms of Bayesian
networks, CNNs describe the causality between the given dataset and the
corresponding labels, which can be simplified as the joint distribution
$P(\boldsymbol{X},\boldsymbol{Y})=P(\boldsymbol{Y|X})P(\boldsymbol{X})$, where
$P(\boldsymbol{X})$ corresponds to the prior belief, and $P(\boldsymbol{Y|X})$
describes the statistical causality. Leveraging their powerful representation
ability, CNNs definitely can learn a $P(\boldsymbol{Y|X})$ from the training
dataset and corresponding random labels, thus the training error is zero.
However, the causality represented by the learned $P(\boldsymbol{Y|X})$ is
obviously not consistent with the testing labels, because labels are random,
thus the testing error is always around $90\%$. In summary, the high testing
error on random labels is because $P(\boldsymbol{Y|X})$ cannot model two
independent variables, but not related to the prior distribution, i.e.,
Bayesian regularization.

Figure 5:  The training and testing errors of CNN3 on random labels over 200
training epochs.

##  5 Conclusion

In this paper, we demonstrate CNNs have explicitly Bayesian regularizations.
First, we introduce a novel probabilistic representation for the hidden layers
of CNNs and demonstrate that CNNs correspond to a Bayesian network with the
serial connection. We conclude that CNNs have explicitly Bayesian
regularizations because the hidden layers close to the input formulate prior
distributions. Moreover, based on the proposed Bayesian regularization, we
clarify two recently observed empirical phenomena that are inconsistent with
traditional theories of generalization. In the future, we hope to extend the
proposed Bayesian regularization theory to general DNNs, e.g., the residual
networks [8], and to design new regularization algorithms through pre-training
the hidden layers corresponding to prior distributions for improving the
generalization performance of DNNs.

## References

  * [1] Sanjeev Arora, Nadav Cohen, Wei Hu, and Yuping Luo.  Implicit regularization in deep matrix factorization.  arXiv preprint arXiv:1905.13655, 2019.
  * [2] Roberto Battiti.  First-and second-order methods for learning: between steepest descent and newton’s method.  Neural computation, 4(2):141–166, 1992.
  * [3] Laurent Dinh, Razvan Pascanu, Samy Bengio, and Yoshua Bengio.  Sharp minima can generalize for deep nets.  In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 1019–1028. JMLR. org, 2017.
  * [4] Adrià Garriga-Alonso, Carl Edward Rasmussen, and Laurence Aitchison.  Deep convolutional networks as shallow gaussian processes.  arXiv preprint arXiv:1808.05587, 2018.
  * [5] S. Geman and D. Geman.  Stochastic relaxation, gibbs distributions, and the bayesian restoration of images.  IEEE Transactions. on Pattern Analysis and Machine Intelligence, pages 721–741, June 1984.
  * [6] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.
  * [7] Suriya Gunasekar, Jason D Lee, Daniel Soudry, and Nati Srebro.  Implicit bias of gradient descent on linear convolutional networks.  In Advances in Neural Information Processing Systems, pages 9461–9471, 2018.
  * [8] K He, X Zhang, S Ren, and J Sun.  Deep residual learning for image recognition.  In CVPR, pages 770–778, 2016.
  * [9] Geoffrey E Hinton.  Products of experts.  9th International Conference on Artificial Neural Networks, 1999\.
  * [10] Geoffrey E. Hinton.  Training products of experts by minimizing contrastive divergence.  Neural Computation, 14:1771–1800, 2002.
  * [11] Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Peter Tang.  On large-batch training for deep learning: Generalization gap and sharp minima.  arXiv preprint arXiv:1609.04836, 2016.
  * [12] Timo Koski and John Noble.  Bayesian networks: an introduction, volume 924.  John Wiley & Sons, 2011.
  * [13] Krizhevsky, Sutskever A., and G. E. Hinton.  Imagenet classification with deep convolutional neural networks.  In NeurIPS, volume 25, pages 1090–1098, 2012.
  * [14] Xinjie Lan.  A synthetic dataset for deep learning.  arXiv preprint arXiv:1906.11905, 2019.
  * [15] Yann LeCun, Sumit Chopra, Raia Hadsell, Marc’Aurelio Ranzato, and Fu Jie Huang.  A tutorial on energy-based learning.  MIT Press, 2006.
  * [16] Jaehoon Lee, Yasaman Bahri, Roman Novak, Samuel S. Schoenholz, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Deep neural networks as gaussian processes.  In ICLR, 2018.
  * [17] Siwei Lyu and Eero. P. Simoncelli.  Statistical modeling of images with fields of gaussian scale mixtures.  In NeurIPS, 2007.
  * [18] Alexander Matthews, Mark Rowland, Jiri Hron, Richard E. Turner, and Zoubin Ghahramani.  Gaussian process behaviour in wide deep neural networks.  In ICLR, 2018.
  * [19] Pankaj Mehta and David J. Schwab.  An exact mapping between the variational renormalization group and deep learning.  arXiv preprint arXiv:1410.3831, 2014.
  * [20] Behnam Neyshabur, Srinadh Bhojanapalli, David McAllester, and Nathan Srebro.  Exploring generalization in deep learning.  In NeurIPS, 2017.
  * [21] Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro.  In search of the real inductive bias: On the role of implicit regularization in deep learning.  In ICLR, 2015.
  * [22] Frank Nielsen and Vincent Garcia.  Statistical exponential families: A digest with flash cards.  arXiv preprint arXiv:0911.4863, 2009.
  * [23] Thomas Dyhre Nielsen and Finn Verner Jensen.  Bayesian networks and decision graphs.  Springer Science & Business Media, 2009.
  * [24] Roman Novak, Lechao Xiao, Jaehoon Lee, Yasaman Bahri, Greg Yang, Jiri Hron, Daniel A. Abolafia, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Bayesian deep convolutional networks with many channels are gaussian processes.  In ICLR, 2018.
  * [25] Aaron Oord and B Schrauwen.  Factoring variations in natural images with deep gaussian mixture models.  In NeurIPS, 2014.
  * [26] Ankit Patel, Minh Nguyen, and Richard Baraniuk.  A probabilistic framework for deep learning.  In NeurIPS, 2016.
  * [27] M. D. Richard and R. Lippmann.  Neural network classifiers estimate bayesian a posteriori probabilities.  Neural Compution, 3:461–483, 1991.
  * [28] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning representations by back-propagating errors.  Nature, 323:533–536, October 1986.
  * [29] Ruslan Salakhutdinov and Geoffrey Hinton.  Deep boltzmann machines.  In AISTATS 2009, pages 448–455, 2009.
  * [30] Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro.  The implicit bias of gradient descent on separable data.  The Journal of Machine Learning Research, 19(1):2822–2878, 2018\.
  * [31] Harald Steck and Tommi S. Jaakkola.  On the dirichlet prior and bayesian regularization.  In NeurIPS, 2003.
  * [32] Yichuan Tang, Ruslan Salakhutdinov, and Geoffrey Hinton.  Deep mixtures of factor analysers.  arXiv preprint arXiv:1206.4635, 2015.
  * [33] Chaohui Wang, Nikos Komodakis, and Nikos Paragios.  Markov random field modeling, inference & learning in computer vision & image understanding: A survey.  Computer Vision and Image Understanding, 117(11):1610–1627, 2013\.
  * [34] Sho Yaida.  Non-gaussian processes and neural networks at finite widths.  arXiv preprint arXiv:1910.00019, 2019.
  * [35] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning requires rethinking generalization.  In ICLR, 2016.
  * [36] G. Zhang.  Neural networks for classification: a survey.  IEEE Transactions on Systems, Man, and Cybernetics, 30:451–462, 2000.

##  6 Appendix

###  6.1 Proof: the proposed probabilisitic representation for the hidden
layers of CNNs

Figure 6:  The $\text{CNN}=\\{\boldsymbol{x;f_{1};f_{2};f_{3};f_{Y}}\\}$ has
two convolutional hidden layers $\boldsymbol{f_{1}}$ and $\boldsymbol{f_{2}}$
and a fully connected hidden layer $\boldsymbol{f_{3}}$. The flattening
$\boldsymbol{f_{2}}$ is the input of $\boldsymbol{f_{3}}$.

Based on the CNN show in Figure 6, we prove the proposed probabilistic
representation for the hidden layers of CNNs in Section 3.2. In the above CNN,
the input $\boldsymbol{x}$ has $M$ channels
$\boldsymbol{x}=\\{\boldsymbol{x}^{m}\\}_{m=1}^{M}$ (e.g., $M=3$ if
$\boldsymbol{x}$ are color images), $\boldsymbol{f_{1}}$ has $N$ convolutional
channels, i.e., $\boldsymbol{f_{1}}=\\{\boldsymbol{f_{1}^{n}}\\}_{n=1}^{N}$,
and $\boldsymbol{f_{2}}$ has $Q$ convolutional channels, i.e.,
$\boldsymbol{f_{2}}=\\{\boldsymbol{f_{2}^{q}}\\}_{q=1}^{Q}$. In addition,
$\boldsymbol{f_{3}}$ has $K$ neurons, i.e.,
$\boldsymbol{f_{3}}=\\{\boldsymbol{f_{3k}}\\}_{k=1}^{K}$.

Assuming the output layer $\boldsymbol{f_{Y}}$ is the softmax with $L$ nodes,
its distribution can be formulated as

$P(\boldsymbol{F_{Y}})=\\{P(f_{yl})=\frac{1}{Z_{\boldsymbol{F_{Y}}}}\text{exp}(f_{yl})\\}_{l=1}^{L}$
(10)

where $Z_{\boldsymbol{F_{Y}}}=\sum_{l=1}^{L}\text{exp}(f_{yl})$ is the
partition function. Since $f_{yl}=\sum_{k=1}^{K}\gamma_{kl}\cdot
f_{3k}+b_{yl}$, we have

$P(\boldsymbol{F_{Y}})=\\{P(f_{yl})=\frac{1}{Z_{\boldsymbol{Y}}}\text{exp}(\sum_{k=1}^{K}\gamma_{kl}\cdot
f_{3k}+b_{yl})\\}_{l=1}^{L}$ (11)

Based on the properties of the exponential function, i.e.,
${\textstyle\text{exp}(a+b)=\text{exp}(a)\cdot\text{exp}(b)}$ and
${\textstyle\text{exp}(a\cdot b)=[\text{exp}(b)]^{a}}$, we can reformulate
$P(\boldsymbol{F_{Y}})$ as

$P(\boldsymbol{F_{Y}})=\\{P(f_{yl})=\frac{1}{Z^{\prime}_{\boldsymbol{Y}}}\prod_{k=1}^{K}[\text{exp}(f_{3k})]^{\gamma_{kl}}\\}_{l=1}^{L}$
(12)

where $Z^{\prime}_{\boldsymbol{Y}}=Z_{\boldsymbol{Y}}/\text{exp}(b_{yl})$.
Since $\\{\text{exp}(f_{3k})\\}_{k=1}^{K}$ are scalar, we can introduce a new
partition function ${\textstyle
Z_{\boldsymbol{F_{3}}}=\sum_{k=1}^{K}\text{exp}(f_{3k})}$ such that
$\\{\frac{1}{Z_{\boldsymbol{F_{3}}}}\text{exp}(f_{3k})\\}_{k=1}^{K}$ becomes a
probability measure. As a result, we can further reformulate
$P(\boldsymbol{F_{Y}})=\\{P(f_{yl})\\}_{l=1}^{L}$ as a Product of Expert (PoE)
model as

$P(\boldsymbol{F_{Y}})=\\{P(f_{yl})=\frac{1}{Z^{\prime\prime}_{\boldsymbol{Y}}}\prod_{k=1}^{K}[\frac{1}{Z_{\boldsymbol{F_{3}}}}\text{exp}(f_{3k})]^{\gamma_{kl}}\\}_{l=1}^{L}$
(13)

where ${\textstyle
Z^{\prime\prime}_{\boldsymbol{Y}}=Z_{\boldsymbol{Y}}/[\text{exp}(b_{yl})\cdot\prod_{k=1}^{K}[Z_{\boldsymbol{F_{3}}}]^{\gamma_{kl}}}]$
and each expert is defined as
$\frac{1}{Z_{\boldsymbol{F_{3}}}}\text{exp}(f_{3k})$.

It is noteworthy that all the experts
$\\{\frac{1}{Z_{\boldsymbol{F_{3}}}}\text{exp}(f_{3k})\\}_{k=1}^{K}$ form a
probability measure and establish an exact one-to-one correspondence to all
the neurons in the third hidden layer $\boldsymbol{f_{3}}$, i.e.,
$\\{f_{3l}\\}_{k=1}^{K}$. Therefore, the distribution of $\boldsymbol{f_{3}}$
can be expressed as

$P(\boldsymbol{F_{3}})=\\{P(f_{3k})=\frac{1}{Z_{\boldsymbol{F_{3}}}}\text{exp}(f_{3k})\\}_{k=1}^{K}$
(14)

Based on the definition of Boltzmann distribution, Equation 10 and 14 show
that $\boldsymbol{f_{Y}}$ and $\boldsymbol{f_{3}}$ formulate two multivariate
Boltzmann distributions and their energy functions are equivalent to the
negative of the output nodes $\\{f_{yl}\\}_{l=1}^{L}$ and the neurons
$\\{f_{3k}\\}_{k=1}^{K}$, respectively.

In summary, the above derivation proves Proposition 1, i.e., a fully connected
layer formulates a multivariate discrete Boltzmann distribution and the energy
function of which is equivalent to the negative activations of the fully
connected layer. The subsequent derivation proves Proposition 2.

Since $\\{f_{3k}=\sum_{j=1}^{J}\beta_{jk}\cdot f_{2j}+b_{3k}\\}_{k=1}^{K}$, we
can derive

$P(\boldsymbol{F_{3}})=\\{P(f_{3k})=\frac{1}{Z_{\boldsymbol{F_{3}}}}\text{exp}(\sum_{j=1}^{J}\beta_{jk}\cdot
f_{2j}+b_{3k})\\}_{k=1}^{K}$ (15)

where $\\{f_{2j}\\}_{j=1}^{J}$ is the flattened output of
$\boldsymbol{f_{2}}$, $\beta_{jk}$ is the weight of the edge between $f_{2j}$
and $f_{3k}$, and $b_{3k}$ denotes the bias. Since $\boldsymbol{f_{2}}$ is a
convolutional layer with $Q$ convolutional channels, it can be reformulated as
$\\{f_{2j}\\}_{j=1}^{J}=\\{\boldsymbol{f^{q}_{2}}\\}_{q=1}^{Q}$. As a result,
$P(\boldsymbol{F_{3}})$ can be reformulated as

$P(\boldsymbol{F_{3}})=\\{P(f_{3k})=\frac{1}{Z_{\boldsymbol{F_{3}}}}\text{exp}(\sum_{q=1}^{Q}\boldsymbol{\beta}^{T}_{qk}\cdot\boldsymbol{f^{q}_{2}}+b_{3k})\\}_{k=1}^{K}$
(16)

where $\boldsymbol{\beta}_{qk}$ is a subset of the parameters
$\\{\beta_{jk}\\}_{j=1}^{J}$ such that
$\sum_{q=1}^{Q}\boldsymbol{\beta}^{T}_{qk}\cdot\boldsymbol{f}^{q}_{2}=\sum_{j=1}^{J}\beta_{jk}\cdot
f_{2j}$, and $\boldsymbol{\beta}^{T}_{qk}$ is the transpose of
$\boldsymbol{\beta}_{qk}$. Therefore, $P(\boldsymbol{F_{3}})$ can be
reformulated as

$P(\boldsymbol{F_{3}})=\\{P(f_{3k})=\frac{1}{Z^{\prime}_{\boldsymbol{F_{3}}}}\prod_{q=1}^{Q}\text{exp}(\boldsymbol{\beta}^{T}_{qk}\cdot\boldsymbol{f^{q}_{2}})\\}_{k=1}^{K}$
(17)

where
$Z^{\prime}_{\boldsymbol{F_{3}}}=Z_{\boldsymbol{F_{3}}}/[\text{exp}(b_{3k})]$.
Recall the element-wise matrix power, e.g.,
$\text{exp}(\boldsymbol{a})^{\boldsymbol{b}}$, where
$\boldsymbol{a}=[a_{1};a_{2};a_{3}]$ and $\boldsymbol{b}=[b_{1};b_{2};b_{3}]$,
we can derive that
$\text{exp}(\boldsymbol{a})=[\text{exp}(a_{1});\text{exp}(a_{2});\text{exp}(a_{3})]$
and
$\text{exp}(\boldsymbol{a})^{\boldsymbol{b}}=[\text{exp}(a_{1})^{b_{1}};\text{exp}(a_{2})^{b_{2}};\text{exp}(a_{3})^{b_{3}}]=[\text{exp}(a_{1}b_{1});\text{exp}(a_{2}b_{2});\text{exp}(a_{3}b_{3})]$.
As a result,
$\text{exp}(\boldsymbol{ab})=\text{exp}(a_{1}b_{1}+a_{2}b_{2}+a_{3}b_{3})=\prod_{|a|}\text{exp}(\boldsymbol{a})^{\boldsymbol{b}}$,
where $|\boldsymbol{a}|$ is the number of elements in $\boldsymbol{a}$.

Based on the element-wise matrix power, we can reformulate
$P(\boldsymbol{F_{3}})$ as

$P(\boldsymbol{F_{3}})=\\{P(f_{3k})\approx\frac{1}{Z^{\prime}_{\boldsymbol{F_{3}}}}\prod_{q=1}^{Q}\text{exp}(\boldsymbol{\beta}^{T}_{qk}\cdot\boldsymbol{f^{q}_{2}})=\frac{1}{Z^{\prime}_{\boldsymbol{F_{3}}}}\prod_{q=1}^{Q}\prod_{|\boldsymbol{f^{q}_{2}}|}[\text{exp}(\boldsymbol{f^{q}_{2}})]^{\boldsymbol{\beta}_{qk}}\\}_{k=1}^{K}$
(18)

Moreover, we introduce a new partition function
$Z_{\boldsymbol{F^{q}_{2}}}=\sum_{\boldsymbol{f^{q}_{2}}}\text{exp}(\boldsymbol{f^{q}_{2}})$
to guarantee that
$\frac{1}{Z_{\boldsymbol{F^{q}_{2}}}}\text{exp}(\boldsymbol{f^{q}_{2}})$ is a
probability measure. As a result, $P(\boldsymbol{F_{3}})$ can be reformulated
as

$P(\boldsymbol{F_{3}})=\\{P(f_{3k})=\frac{1}{Z^{\prime\prime}_{\boldsymbol{F_{3}}}}\prod_{q=1}^{Q}\prod_{|\boldsymbol{f^{q}_{2}}|}[\frac{1}{Z_{\boldsymbol{F^{q}_{2}}}}\text{exp}(\boldsymbol{f^{q}_{2}})]^{\boldsymbol{\beta}_{qk}}\\}_{k=1}^{K}$
(19)

where
$Z^{\prime\prime}_{\boldsymbol{F_{3}}}=Z^{\prime}_{\boldsymbol{F_{3}}}/\prod_{q=1}^{Q}[Z_{\boldsymbol{F^{q}_{2}}}]^{\boldsymbol{\beta}_{qk}\cdot|\boldsymbol{f^{q}_{2}}|}$.
Overall, $P(\boldsymbol{F_{3}})$ can be reformulated as a PoE model, in which
each expert is defined by every convolutional channel in $\boldsymbol{f_{2}}$

$P(\boldsymbol{f}^{q}_{2})=\frac{1}{Z_{\boldsymbol{F^{q}_{2}}}}\text{exp}(\boldsymbol{f^{q}_{2}})$
(20)

Since $\boldsymbol{f_{2}}$ has $K$ convolutional channels, i.e.,
$\boldsymbol{f_{2}}=\\{\boldsymbol{f^{q}_{2}}\\}_{q=1}^{Q}$, and each channel
can be formulated as the summation of all the convolutional channels in
$\boldsymbol{f_{1}}$, i.e.,

$\boldsymbol{f_{2}}=\\{\boldsymbol{f^{q}_{2}}=\sigma_{2}(\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2}\cdot\boldsymbol{1})\\}_{q=1}^{Q}$
(21)

where $\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}}$ is the output of
the $q$th convolutional filter applying into the $n$th channel of
$\boldsymbol{f_{1}}$, $b_{2}^{q}$ is the bias, and $\sigma_{2}(\cdot)$ is the
activation function. Therefore, we have

$P(\boldsymbol{f}^{q}_{2})=\frac{1}{Z_{\boldsymbol{F^{q}_{2}}}}\text{exp}[\sigma_{2}(\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2}\cdot\boldsymbol{1})]$
(22)

Based on the equivalence between the gradient descent learning and the first
order approximation (Appendix 6.4), Equation 22 can be approximated as

$P(\boldsymbol{f}^{q}_{2})\approx\frac{1}{Z_{\boldsymbol{F^{q}_{2}}}}\text{exp}(\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2}\cdot\boldsymbol{1})$
(23)

Let us define a linear function as
$\varphi^{q}_{2}(\boldsymbol{f_{1}})=\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2}\cdot\boldsymbol{1}$,
thus the distribution of the $q$th convolutional channel in
$\boldsymbol{f_{2}}$, i.e., $P(\boldsymbol{f}^{q}_{2})$, can be expressed as

$P(\boldsymbol{f}^{q}_{2})\approx\frac{1}{Z_{\boldsymbol{F^{q}_{2}}}}\text{exp}[\varphi^{q}_{2}(\boldsymbol{f_{1}})]$
(24)

If we regard all the linear filters
$\\{\varphi^{q}_{2}(\boldsymbol{f_{1}})\\}_{k=1}^{K}$ as potential functions
for modeling signal structures of $\boldsymbol{f_{1}}$, we can formulate
$P(\boldsymbol{F_{2}})$ as a specific Boltzmann distribution, i.e., the Markov
Random Fields (MRFs), which can be expressed as

$P(\boldsymbol{F_{2}})=\frac{1}{Z_{\boldsymbol{F_{2}}}}\prod_{q=1}^{Q}P(\boldsymbol{f}^{q}_{2})=\frac{1}{Z_{\boldsymbol{F_{2}}}}\text{exp}[\sum_{q=1}^{Q}\varphi^{q}_{2}(\boldsymbol{f_{1}})]$
(25)

where
$Z_{\boldsymbol{F_{2}}}=\int_{\boldsymbol{f_{1}}}\text{exp}[\sum_{q=1}^{Q}\varphi^{k}_{2}(\boldsymbol{f_{1}})]d\boldsymbol{f_{1}}$
is the partition function for $P(\boldsymbol{F_{2}})$. In summary, we prove
that the convolutional layer $\boldsymbol{f_{2}}$ formulates a MRF model to
describe the statistical property of $\boldsymbol{f_{1}}$.

Figure 7:  The equivalence between a convolutional operation and a matrix
multiplication.

Subsequently, we prove that the convolutional layer $\boldsymbol{f_{1}}$ also
formulates a MRF model to describe the statistical property of the input
$\boldsymbol{x}$. Based on Equation 23, we have

$P(\boldsymbol{f}^{q}_{2})\approx\frac{1}{Z^{\prime}_{\boldsymbol{F^{q}_{2}}}}\text{exp}(\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}})$
(26)

where
$Z^{\prime}_{\boldsymbol{F^{k}_{2}}}=Z_{\boldsymbol{F^{q}_{2}}}/\text{exp}(b^{q}_{2})$.
Based on the equivalence between a convolutional channel and a matrix
multiplication proven in Appendix 6.5, we can regard a convolutional operation
as a matrix multiplication, which is visualized in Figure 7. As a result, we
have

$P(\boldsymbol{f}^{q}_{2})\approx\frac{1}{Z^{\prime}_{\boldsymbol{F^{q}_{2}}}}\text{exp}(\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}})=\frac{1}{Z^{\prime}_{\boldsymbol{F^{q}_{2}}}}\text{exp}(\sum_{n=1}^{N}\boldsymbol{W^{(q,n)}_{2}}\cdot\boldsymbol{f^{n}_{1}})$
(27)

where $\boldsymbol{W^{(q,n)}_{2}}$ is a matrix corresponding to the
convolution filter $\boldsymbol{S^{(q,n)}_{2}}$.

Based on the element-wise matrix power, $P(\boldsymbol{f}^{q}_{2})$ can be
reformulated as

$P(\boldsymbol{f}^{q}_{2})\approx\frac{1}{Z^{\prime}_{\boldsymbol{F^{q}_{2}}}}\prod_{n=1}^{N}\text{exp}(\boldsymbol{W^{(q,n)}_{2}}\cdot\boldsymbol{f^{n}_{1}})=\frac{1}{Z^{\prime}_{\boldsymbol{F^{q}_{2}}}}\prod_{n=1}^{N}\text{exp}(\boldsymbol{f^{n}_{1}})^{\boldsymbol{W^{(q,n)}_{2}}}$
(28)

Similarly, we introduce a new partition function
$Z_{\boldsymbol{F^{n}_{1}}}=\int_{\boldsymbol{f^{n}_{1}}}\text{exp}(\boldsymbol{f^{n}_{1}})d\boldsymbol{f^{n}_{1}}$
to guarantee that
$\frac{1}{Z_{\boldsymbol{F^{n}_{1}}}}\text{exp}(\boldsymbol{f^{n}_{1}})$ is a
probability measure. As a result, $P(\boldsymbol{f^{q}_{2}})$ can be
reformulated as

$P(\boldsymbol{f^{q}_{2}})=P(f_{yl})\approx\frac{1}{Z^{\prime\prime}_{\boldsymbol{F^{q}_{2}}}}\prod_{n=1}^{N}[\frac{1}{Z_{\boldsymbol{F^{n}_{1}}}}\text{exp}(\boldsymbol{f^{n}_{1}})]^{\boldsymbol{W^{(q,n)}_{2}}}$
(29)

where
$Z^{\prime\prime}_{\boldsymbol{F^{q}_{2}}}=Z^{\prime}_{\boldsymbol{F^{q}_{2}}}/\prod_{n=1}^{N}[Z_{\boldsymbol{F^{n}_{1}}}]^{\boldsymbol{W^{(q,n)}_{2}}}$.
Therefore, $P(\boldsymbol{f^{q}_{2}})$ can be reformulated as a PoE model, in
which all the experts are defined as
$P(\boldsymbol{f^{n}_{1}})=\frac{1}{Z_{\boldsymbol{F^{n}_{1}}}}\text{exp}(\boldsymbol{f^{n}_{1}})$.

Since $\boldsymbol{f_{1}}$ has $N$ convolutional channels

$\boldsymbol{f_{1}}=\\{\boldsymbol{f^{n}_{1}}=\sigma_{1}(\sum_{m=1}^{M}\boldsymbol{S^{(n,m)}_{1}}\circ\boldsymbol{x^{m}}+b^{n}_{1}\cdot\boldsymbol{1})\\}_{n=1}^{N}$
(30)

Based on the same derivation as Equations 23, 24, 25, $P(\boldsymbol{F_{1}})$
also can be expressed as a MRF model

$P(\boldsymbol{F_{1}})=\frac{1}{Z_{\boldsymbol{F_{1}}}}\prod_{n=1}^{N}P(\boldsymbol{f}^{n}_{1})=\frac{1}{Z_{\boldsymbol{F_{1}}}}\text{exp}[\sum_{n=1}^{N}\varphi^{n}_{1}(\boldsymbol{x})]$
(31)

where
$\varphi^{n}_{1}(\boldsymbol{x})=\sum_{m=1}^{M}\boldsymbol{S^{(n,m)}_{1}}\circ\boldsymbol{x^{m}}+b^{n}_{1}\cdot\boldsymbol{1})$
and
$Z_{\boldsymbol{F_{1}}}=\int_{\boldsymbol{x}}\text{exp}[\sum_{n=1}^{N}\varphi^{n}_{1}(\boldsymbol{x})]d\boldsymbol{x}$.
Overall, the distribution of a convolutional layer can be formulated as a MRF
model.

###  6.2 Proof: the variable inference of the Bayesian networks corresponding
to the CNNs

Given a $\text{CNN}=\\{\boldsymbol{x;f_{1};...;f_{I};f_{Y}}\\}$, we
demonstrate the entire architecture of the CNN corresponds to a Bayesian
network as

${\textstyle
P(\boldsymbol{F_{Y};F_{I};\cdots;F_{1}|X})=P(\boldsymbol{F_{Y}|F_{I}})\cdot...\cdot
P(\boldsymbol{F_{i+1}|F_{i}})\cdot...\cdot P(\boldsymbol{F_{1}|X})}$ (32)

Here we prove that the posterior distribution $P(\boldsymbol{F_{Y}|X})$ is a
Boltzmann distribution, and the energy function of which is also equivalent to
the functionality of the entire CNN.

Based on the variable elimination, $P(\boldsymbol{F_{Y}|X})$ can be formulated
as

$\begin{split}P(\boldsymbol{F_{Y}|X})&=\int_{\boldsymbol{F_{I}}}\cdots\int_{\boldsymbol{F_{1}}}P(\boldsymbol{F_{Y};F_{I}..;F_{1}|X})d{\boldsymbol{F_{I}}}\cdots
d{\boldsymbol{F_{1}}}\\\
&=P(\boldsymbol{F_{Y}|F_{I}})\int_{\boldsymbol{F_{I}}}P(\boldsymbol{F_{I}|F_{I-1}})d{\boldsymbol{F_{I}}}\cdots\int_{\boldsymbol{F_{i}}}P(\boldsymbol{F_{i}|F_{i-1}})d{\boldsymbol{F_{i}}}\cdots\int_{\boldsymbol{F_{1}}}P(\boldsymbol{F_{1}|X})d{\boldsymbol{F_{1}}}\end{split}$
(33)

where
$P(\boldsymbol{F_{i}|F_{i-1}})=\frac{1}{\boldsymbol{Z_{F_{i}}}}\text{exp}[-E_{\boldsymbol{F_{i}}}(\boldsymbol{f_{i-1}})]$
is a Boltzmann distribution. All the above integrations, e.g.,
$\int_{\boldsymbol{F_{i}}}P(\boldsymbol{F_{i}|F_{i-1}})d{\boldsymbol{F_{i}}}$,
are equal to 1, because the integral space, i.e., $\boldsymbol{F_{i}}$, are
entirely depends on $\boldsymbol{f_{i-1}}$. As a result,
$P(\boldsymbol{F_{Y}|X})=P(\boldsymbol{F_{Y}|F_{I}})$. Since
$\boldsymbol{f_{Y}}$ formulates a Boltzmann distribution, and its energy
function is
$E_{\boldsymbol{F_{Y}}}(\boldsymbol{f_{I}})=-f_{Y}(f_{I}(f_{I-1}\cdots
f_{1}(\boldsymbol{x})))$, we have

${\textstyle
P(\boldsymbol{F_{Y}|X})=\frac{1}{\boldsymbol{Z_{F_{Y}}}}\text{exp}[-f_{Y}(f_{I}(f_{I-1}\cdots
f_{1}(\boldsymbol{x})))]}$ (34)

Therefore $P(\boldsymbol{F_{Y}|X})$ is also a Boltzmann distribution, and the
energy function of which is also equivalent to the functionality of the entire
CNN.

###  6.3 The method for generating the synthetic dataset

Algorithm 1 The algorithm for generating the synthetic dataset

1:NIST dataset of handwritten digits by class

2:repeat

3: binarizing an image of NIST to obtain $\footnotesize\boldsymbol{z}$

4: extracting the central part of $\footnotesize\boldsymbol{z}$ to obtain
$\footnotesize\boldsymbol{z}_{c}$ with dimension $\footnotesize 64\times 64$

5: downsampling $\footnotesize\boldsymbol{z}_{c}$ to obtain
$\footnotesize\boldsymbol{z}_{cd}$ with dimension $\footnotesize 32\times 32$

6: extracting the edge of $\footnotesize\boldsymbol{z}_{cd}$ to obtain the
mask image $\footnotesize\boldsymbol{m}_{cd}$

7: decomposing $\footnotesize\boldsymbol{m}_{cd}$ into four parts, i.e.,
$\footnotesize\boldsymbol{m}_{\text{outside}}$,
$\footnotesize\boldsymbol{m}_{\text{outside-boundary}}$,
$\footnotesize\boldsymbol{m}_{\text{inside-boundary}}$, and
$\footnotesize\boldsymbol{m}_{\text{inside}}$.

8: sampling $\footnotesize\mathcal{N}(0,1024)$ to derive a random vector
$\footnotesize\boldsymbol{x}\in\boldsymbol{\mathbb{R}}^{1024\times 1}$

9: sorting $\footnotesize\boldsymbol{x}$ in the descending order to derive
$\footnotesize\boldsymbol{\hat{x}}$

10: decomposing $\footnotesize\boldsymbol{\hat{x}}$ into four parts, i.e.,
$\footnotesize\boldsymbol{\hat{x}}=\\{\boldsymbol{\hat{x}}_{\text{outside}},\boldsymbol{\hat{x}}_{\text{outside-
boundary}},\boldsymbol{\hat{x}}_{\text{inside-
boundary}},\boldsymbol{\hat{x}}_{\text{inside}}\\}$

11: placing each pixel of
$\footnotesize\boldsymbol{\hat{x}}=\\{\boldsymbol{\hat{x}}_{\text{outside}},\boldsymbol{\hat{x}}_{\text{outside-
boundary}},\boldsymbol{\hat{x}}_{\text{inside-
boundary}},\boldsymbol{\hat{x}}_{\text{inside}}\\}$ into the above masks.

12:until (20,000 synthetic images are generated)

13:The synthetic dataset

In this section, we present a novel approach to generate a synthetic dataset
obeying the Gaussian distribution based on the NIST
111https://www.nist.gov/srd/nist-special-database-19 dataset of handwritten
digits. The synthetic dataset has similar characteristics as the benchmark
MNIST dataset. It consists of 20,000 $32\times 32$ grayscale images in 10
classes (digits from 0 to 9), and each class has 1,000 training and 1,000
testing images. The approach to generate the synthetic dataset is summarized
in Algorithm 1, and Figure 8 visualize the spatial relation of
$\boldsymbol{\hat{x}}=\\{\boldsymbol{\hat{x}}_{\text{outside}},\boldsymbol{\hat{x}}_{\text{outside-
boundary}},\boldsymbol{\hat{x}}_{\text{inside-
boundary}},\boldsymbol{\hat{x}}_{\text{inside}}\\}$.

Figure 8:  The first row shows an original image, its edge, and the
corresponding synthetic image based on the original one. The second row uses
white pixels to show the spatial position of the four mask parts, i.e.,
$\boldsymbol{m}_{\text{outside}}$, $\boldsymbol{m}_{\text{outside-boundary}}$,
$\boldsymbol{m}_{\text{inside-boundary}}$, and
$\boldsymbol{m}_{\text{inside}}$. The third row shows the synthetic image
corresponding to each mask part, i.e.,
$\boldsymbol{\hat{x}}=\\{\boldsymbol{\hat{x}}_{\text{outside}},\boldsymbol{\hat{x}}_{\text{outside-
boundary}},\boldsymbol{\hat{x}}_{\text{inside-
boundary}},\boldsymbol{\hat{x}}_{\text{inside}}\\}$.

###  6.4 Proof: the equivalence between the gradient descent and the first
order approximation

In the context of deep learning, most learning algorithms belong to the
gradient descent algorithm. Given a
$\text{CNN}=\\{\boldsymbol{x;f_{1};...;f_{I};f_{Y}}\\}$ and a cost function
$H[\boldsymbol{f_{Y}},P(\boldsymbol{Y|X})]$, where $P(\boldsymbol{Y|X})$ is
the true distribution of the training label given the corresponding training
dataset. Let $\boldsymbol{\theta}$ be the parameters of the CNN, the gradient
descent aims to optimize $\boldsymbol{\theta}$ by minimizing
$H[\boldsymbol{f_{Y}},P(\boldsymbol{Y|X})]$ [28]. We typically update
$\boldsymbol{\theta}$ iteratively to derive $\boldsymbol{\theta}^{*}$, which
can be expressed as

$\boldsymbol{\theta}_{t+1}=\boldsymbol{\theta}_{t}-\alpha\nabla_{\theta_{t}}H[\boldsymbol{f_{Y}},P(\boldsymbol{Y|X})]$
(35)

where $\nabla_{\theta_{t}}H[\boldsymbol{f_{Y}},P(\boldsymbol{Y|X})]$ is the
Jacobian matrix of $H[\boldsymbol{f_{Y}},P(\boldsymbol{Y|X})]$ with respect to
$\boldsymbol{\theta}_{t}$ at the $t$th iteration, and $\alpha>0$ denotes a
constant indicating the learning rate. Since $P(\boldsymbol{Y|X})$ is
constant, we denote $H(\boldsymbol{f_{Y}})$ as
$H[\boldsymbol{f_{Y}},P(\boldsymbol{Y|X})]$ for simplifying the following
derivation.

Since the functions of all hidden layers are differentiable and the output of
a hidden layer is the input of its next layer, the Jacobian matrix of
$H(\boldsymbol{f_{Y}})$ with respect to the parameters of the $i$th hidden
layer, i.e., $\nabla_{\theta(i)}H(\boldsymbol{f_{Y}})$, can be expressed as
follows based on the chain rule.

$\nabla_{\boldsymbol{\theta}(i)}H(\boldsymbol{f_{Y}})=\nabla_{\boldsymbol{f_{Y}}}H(\boldsymbol{f_{Y}})\cdot\nabla_{\boldsymbol{f_{I}}}\boldsymbol{f_{Y}}\cdot\prod_{j=i+1}^{I}\nabla_{\boldsymbol{f_{j-1}}}\boldsymbol{f_{j}}\cdot\nabla_{\boldsymbol{\theta}(i)}\boldsymbol{f_{i}}$
(36)

where $\boldsymbol{\theta}(i)$ denote the parameters of the $i$th hidden
layer.

Table 2: One iteration of the backpropagation training procedure for the CNN

Layer |  Gradients update $\nabla_{\boldsymbol{\theta}_{t}(i)}H(\boldsymbol{f_{Y}})$ |  | Parameters and activations update |
---|---|---|---|---
$\boldsymbol{f_{Y}}$ | ${\scriptstyle\nabla_{\boldsymbol{f_{Y}}}H(\boldsymbol{f_{Y}})\nabla_{\boldsymbol{\theta}_{t}(Y)}\boldsymbol{f_{Y}}}$ | $\downarrow$ |  ${\scriptstyle\boldsymbol{\theta}_{t+1}(Y)=\boldsymbol{\theta}_{t+1}(Y)-\alpha[\nabla_{\boldsymbol{\theta}_{t}(Y)}H(\boldsymbol{f_{Y}})]}$, ${\scriptstyle\boldsymbol{f_{Y}}(\boldsymbol{f_{2}},\boldsymbol{\theta}_{t+1}(Y))}$ | $\uparrow$
$\boldsymbol{f_{3}}$ | ${\scriptstyle\nabla_{\boldsymbol{f_{Y}}}H(\boldsymbol{f_{Y}})\nabla_{\boldsymbol{f_{3}}}\boldsymbol{f_{Y}}\nabla_{\boldsymbol{\theta}_{t}(3)}\boldsymbol{f_{3}}}$ | $\downarrow$ |  ${\scriptstyle\boldsymbol{\theta}_{t+1}(3)=\boldsymbol{\theta}_{t+1}(3)-\alpha[\nabla_{\boldsymbol{\theta}_{t}(3)}H(\boldsymbol{f_{Y}})]}$, ${\scriptstyle\boldsymbol{f_{3}}(\boldsymbol{f_{2}},\boldsymbol{\theta}_{t+1}(3))}$ | $\uparrow$
$\boldsymbol{f_{2}}$ | ${\scriptstyle\nabla_{\boldsymbol{f_{Y}}}H(\boldsymbol{f_{Y}})\nabla_{\boldsymbol{f_{3}}}\boldsymbol{f_{Y}}\nabla_{\boldsymbol{f_{2}}}\boldsymbol{f_{3}}\nabla_{\boldsymbol{\theta}_{t}(2)}\boldsymbol{f_{2}}}$ | $\downarrow$ |  ${\scriptstyle\boldsymbol{\theta}_{t+1}(2)=\boldsymbol{\theta}_{t+1}(2)-\alpha[\nabla_{\boldsymbol{\theta}_{t}(2)}H(\boldsymbol{f_{Y}})]}$, ${\scriptstyle\boldsymbol{f_{2}}(\boldsymbol{f_{1}},\boldsymbol{\theta}_{t+1}(2))}$ | $\uparrow$
$\boldsymbol{f_{1}}$ | ${\scriptstyle\nabla_{\boldsymbol{f_{Y}}}H(\boldsymbol{f_{Y}})\nabla_{\boldsymbol{f_{3}}}\boldsymbol{f_{Y}}\nabla_{\boldsymbol{f_{2}}}\boldsymbol{f_{3}}\nabla_{\boldsymbol{f_{1}}}\boldsymbol{f_{2}}\nabla_{\boldsymbol{\theta}_{t}(1)}\boldsymbol{f_{1}}}$ | $\downarrow$ |  ${\scriptstyle\boldsymbol{\theta}_{t+1}(1)=\boldsymbol{\theta}_{t+1}(1)-\alpha[\nabla_{\boldsymbol{\theta}_{t}(1)}H(\boldsymbol{f_{Y}})]}$, ${\scriptstyle\boldsymbol{f_{1}}(\boldsymbol{x},\boldsymbol{\theta}_{t+1}(1))}$ | $\uparrow$
$\boldsymbol{x}$ | — |  | — |

  * •

The uparrow and downarrow indicate the order of gradients and
parameters(activations) update, respectively.

Based on Equation 35 and Equation 36, $\boldsymbol{\theta}(i)$ can be learned
by the gradient descent method as

$\boldsymbol{\theta}_{t+1}(i)=\boldsymbol{\theta}_{t}(i)-\alpha[\nabla_{\theta_{t}(i)}H(\boldsymbol{f_{Y}})]$
(37)

Table 2 summarizes the backpropagation training procedure for the CNN in
Figure 6.

If an arbitrary function $\boldsymbol{f}$ is differentiable at point
$\boldsymbol{p}^{*}$ in ${\mathcal{R}}^{N}$ and its differential is
represented by the Jacobian matrix $\nabla_{\boldsymbol{p}}\boldsymbol{f}$,
the first order approximation of $\boldsymbol{f}$ near the point
$\boldsymbol{p}$ can be formulated as

$\boldsymbol{f}(\boldsymbol{p})-\boldsymbol{f}(\boldsymbol{p}^{*})=(\nabla_{\boldsymbol{p}^{*}}\boldsymbol{f})\cdot(\boldsymbol{p}-\boldsymbol{p}^{*})+o(||\boldsymbol{p}-\boldsymbol{p}^{*}||)$
(38)

where $o(||\boldsymbol{p}-\boldsymbol{p}^{*}||)$ is a quantity that approaches
zero much faster than $||\boldsymbol{p}-\boldsymbol{p}^{*}||$ approaches zero.

Since the gradient descent can be viewed as the first order approximation [2],
updating the activations of the hidden layers of the CNN (Figure 6) during the
training procedure can be approximated as

$\begin{split}\boldsymbol{f_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t+1}(2)]&\approx\boldsymbol{f_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t}(2)]+(\nabla_{\boldsymbol{\theta}_{t}(2)}\boldsymbol{f_{2}})\cdot[\boldsymbol{\theta}_{t+1}(2)-\boldsymbol{\theta}_{t}(2)]\\\
\boldsymbol{f_{1}}[\boldsymbol{x},\boldsymbol{\theta}_{t+1}(1)]&\approx\boldsymbol{f_{1}}[\boldsymbol{x},\boldsymbol{\theta}_{t}(1)]+(\nabla_{\boldsymbol{\theta}_{t}(1)}\boldsymbol{f_{1}})\cdot[\boldsymbol{\theta}_{t+1}(1)-\boldsymbol{\theta}_{t}(1)]\\\
\end{split}$ (39)

where $\boldsymbol{f_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t}(2)]$
denote the activations of the second hidden layer based on the parameters
learned in the $t$th iteration, i.e., $\boldsymbol{\theta}_{t}(2)$, given the
activations of the first hidden layer, i.e., $\boldsymbol{f_{1}}$. The
definitions of
$\boldsymbol{f_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t+1}(2)]$ and
$\boldsymbol{f_{1}}(\boldsymbol{x},\boldsymbol{\theta}_{t}(1))$ are the same
as $\boldsymbol{f_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t}(2)]$.

Since $\boldsymbol{f_{2}}$ has $Q$ convolutional channels,
i.e.,$\boldsymbol{f_{2}}=\\{\boldsymbol{f^{q}_{2}}=\sigma_{2}(\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2}\cdot\boldsymbol{1})\\}_{q=1}^{Q}$,
the parameters of a single convolutional channel can be expressed as
$\boldsymbol{\theta^{q}}(2)=\\{\boldsymbol{S^{(q,1)}_{2}},\cdots,\boldsymbol{S^{(q,N)}_{2}},b^{q}_{2}\\}$.
As a result, $\nabla_{\boldsymbol{\theta^{q}}_{t}(2)}\boldsymbol{f^{q}_{2}}$
can be expressed as

$\begin{split}\nabla_{\boldsymbol{\theta^{q}}_{t}(2)}\boldsymbol{f^{q}_{2}}=(\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}})\cdot[\boldsymbol{f^{1}_{1}};\cdots,\boldsymbol{f^{N}_{1}},1]^{T}\end{split}$
(40)

where
$\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}}=\frac{\partial\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t}(2)]}{\partial\sigma_{2}}$.
If we only consider a single convolutional channel $\boldsymbol{f^{q}_{2}}$,
its activations updating can be approximated as

$\begin{split}\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta^{q}}_{t+1}(2)]\approx\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta^{q}}_{t}(2)]+(\nabla_{\boldsymbol{\theta^{q}}_{t}(2)}\boldsymbol{f^{q}_{2}})\cdot[\boldsymbol{\theta^{q}}_{t+1}(2)-\boldsymbol{\theta^{q}}_{t}(2)]\\\
\end{split}$ (41)

Substituting
$(\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}})\cdot[\boldsymbol{f^{1}_{1}};\cdots,\boldsymbol{f^{N}_{1}},1]^{T}$
for $\nabla_{\boldsymbol{\theta^{q}}_{t}(2)}\boldsymbol{f^{q}_{2}}$ in
Equation 41, we can derive

$\begin{split}\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t+1}(2)]\approx&\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t}(2)]+(\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}})\cdot[\boldsymbol{f^{1}_{1}};\cdots,\boldsymbol{f^{N}_{1}},1]^{T}\cdot\boldsymbol{\theta^{q}}_{t+1}(2)\\\
&-(\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}})\cdot[\boldsymbol{f^{1}_{1}};\cdots,\boldsymbol{f^{N}_{1}},1]^{T}\cdot\boldsymbol{\theta^{q}}_{t}(2)\end{split}$
(42)

Assuming
$\boldsymbol{\theta^{q}}_{t+1}(2)=\\{\boldsymbol{S^{(q,1)}_{2}},\cdots,\boldsymbol{S^{(q,N)}_{2}},b^{q}_{2}\\}$
and
$\boldsymbol{\theta^{q}}_{t}(2)=\\{\boldsymbol{S^{(q,1)}_{2*}},\cdots,\boldsymbol{S^{(q,N)}_{2*}},b^{q}_{2*}\\}$,
thus
$[\boldsymbol{f^{1}_{1}};\cdots,\boldsymbol{f^{N}_{1}},1]^{T}\cdot\boldsymbol{\theta^{q}}_{t+1}(2)=\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2}\cdot\boldsymbol{1}$
because a convolution operation is linear. As a result, Equation 42 can be
expressed as

$\begin{split}\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t+1}(2)]\approx&(\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}})\cdot[\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2}\cdot\boldsymbol{1}]\text{
(First order approximation)}\\\
&+\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t}(2)]-(\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}})\cdot[\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2*}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2*}\cdot\boldsymbol{1}]\text{
(Error)}\\\ \end{split}$ (43)

Equation 43 indicates that
$\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t+1}(2)]$ can
be expressed as its first order approximation with an error component based on
the activation in the previous iteration, i.e.,
$\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t}(2)]$. Since
$\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}}=\frac{\partial\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t}(2)]}{\partial\sigma_{2}}$
is only related to $\boldsymbol{f_{1}}$ and the parameters in the $t$th
training iteration, i.e, $\boldsymbol{\theta}_{t}(2)$, it can be regarded as a
constant. Also note that the error component do not contain any parameters in
the $(t+1)$th training iteration. In summary,
$\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t+1}(2)]$ can
be reformulated as

$\begin{split}\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t+1}(2)]&\approx\boldsymbol{C_{1}}\cdot[\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2}\cdot\boldsymbol{1}]+\boldsymbol{C_{2}}\\\
\end{split}$ (44)

where $\boldsymbol{C_{1}}=\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}}$ and
$\boldsymbol{C_{2}}=\boldsymbol{f^{q}_{2}}[\boldsymbol{f_{1}},\boldsymbol{\theta}_{t}(2)]-(\nabla_{\sigma_{2}}\boldsymbol{f^{q}_{2}})\cdot[\sum_{n=1}^{N}\boldsymbol{S^{(q,n)}_{2*}}\circ\boldsymbol{f^{n}_{1}}+b^{q}_{2*}\cdot\boldsymbol{1}]$.
Similarly, the activations in the first hidden layer $\boldsymbol{f_{1}}$ also
can be formulated as the first order approximation in the context of the
gradient descent learning algorithm.

###  6.5 Proof: the equivalence between a convolutional channel and a matrix
multiplication

Figure 9:  The equivalence between a convolution channel and a matrix
multiplication.

Assuming a $\text{CNN}=\\{\boldsymbol{x;f_{1};...;f_{I};f_{Y}}\\}$ has two
adjacent convolutional layers $\boldsymbol{f_{i}}$ and $\boldsymbol{f_{i+1}}$.
The first convolutional layer has $N$ channels, i.e.
$\boldsymbol{f_{i}}=\\{\boldsymbol{f^{n}_{i}}\\}_{n=1}^{N}$, and the second
convolutional layer has $K$ channels, i.e.
$\boldsymbol{f_{i+1}}=\\{\boldsymbol{f^{k}_{i+1}}\\}_{k=1}^{K}$. Therefore,
the dimension of the convolutional filters $\boldsymbol{S_{i+1}}$ in
$\boldsymbol{f_{i+1}}$ are $W\times H\times N\times K$, i.e., the dimension of
a single filter connecting the $\boldsymbol{f^{n}_{i}}$ channel and the
$\boldsymbol{f^{k}_{i+1}}$ channel are
$D(\boldsymbol{S^{(k,n)}_{i+1}})=W\times H$. As a result, the
$\boldsymbol{f^{k}_{i+1}}$ channel can be formulated as

$\boldsymbol{f^{k}_{i+1}}=\sigma_{i+1}(\sum_{n=1}^{N}\boldsymbol{S^{(k,n)}_{i+1}}\circ\boldsymbol{f^{n}_{i}}+b^{k}_{i+1}\cdot\boldsymbol{1})$
(45)

where $\sigma_{i+1}$ is the activation function and $b_{i+1}$ is the bias for
the $(i+1)$th channel.

We can regard the convolutional channel $\boldsymbol{f^{k}_{i+1}}$ as a matrix
multiplication. The weights of the matrix depend on the convolutional filter
and the spatial location of the convolutional output, and the output of the
matrix multiplication is the vectorized output of the convolutional channel.

Without considering the activation function, the equivalence between a
convolutional channel and a matrix multiplication can be explained in Figure
9, in which the first convolutional layer has two channels, i.e.,
$\boldsymbol{f_{i}}=\\{\boldsymbol{x};\boldsymbol{\hat{x}}\\}$, where
$\boldsymbol{x}=\\{x_{1},\cdots,x_{9}\\}$ and
$\boldsymbol{\hat{x}}=\\{\hat{x}_{1},\cdots,\hat{x}_{9}\\}$, and the $k$th
channel of the second convolutional layer is
$\boldsymbol{f^{k}_{i+1}}=\\{z_{1},z_{2},z_{3},z_{4}\\}$. In addition, the
dimension of the convolutional filter $\boldsymbol{S_{i+1}}$ is $2\times
2\times 2\times 1$, where $\boldsymbol{S^{(k,1)}_{i+1}}=[A,B;C,D]$ and
$\boldsymbol{S^{(k,2)}_{i+1}}=[E,F;G,H]$. Figure 9 shows that the output of
$\boldsymbol{f^{k}_{i+1}}$ can be expressed as a matrix multiplication.

If we consider the activation function and make a precise statement,
$\boldsymbol{f^{k}_{i+1}}$ can be formulated as

$\boldsymbol{f^{k}_{i+1}}=\sigma_{i+1}(\boldsymbol{W^{k}_{i+1}}\cdot\boldsymbol{f_{i}}+b^{k}_{i+1}\cdot\boldsymbol{1})$
(46)

where $\boldsymbol{W^{k}_{i+1}}$ is the matrix corresponding to all the
convolutional filters in the $k$th channel.