Year: 2019

Abstract:
In this paper, we demonstrate that the parameters of Deep Neural Networks (DNNs) do not satisfy the i.i.d. prior assumption and activations being i.i.d. is not valid for all the hidden layers of DNNs, thus the Gaussian process cannot correctly model the distribution of the hidden layers. Alternatively, we propose a novel probabilistic representation for the hidden layers of DNNs in two aspects: (i) a hidden layer formulates a Gibbs distribution, of which the energy function is determined by the activations, and (ii) the statistical connection between two adjacent layers can be formulated as a product of experts model. Based on the proposed probabilistic representation, we further confirm the equivalence between the variational inference and the stochastic gradient descent, and prove the later is also equivalent to the energy minimization optimization. As a result, the entire architecture of DNNs can be viewed as a Bayesian network, in which the hidden layers close to the input formulate the prior distributions of the training dataset. In addition, we propose a new regularization approach to improve the generalization performance of DNNs by pre-training the hidden layers corresponding to the prior distributions. Simulation results validate the proposed theories.

heading: 1 INTRODUCTION
text: Recently, interpreting the hidden layers of Deep Neural Networks (DNNs) as the Gaussian Process (GP) has attracted a great deal of attention because it provides a novel probabilistic perspective to clarify the working mechanism of deep learning. Neal (1994) initially demonstrates the equivalence between GP and neural networks with a single fully connected layer in the limit of infinite neurons. Following this seminal work, Lee et al. (2018); Matthews et al. (2018) further extend the equivalence to neural networks with multiple hidden layers and Garriga-Alonso et al. (2018); Novak et al. (2018) establish the correspondence between the convolutional layers and GP.
It is important to note that the GP explanations rely on a fundamental probabilistic premise that the activations of a hidden layer are independent and identically distributed (i.i.d.). More specifically, based on the classical Central Limit Theorem (CLT), the prerequisite of an activation of a hidden layer approaching a Gaussian distribution is that all the activations of the previous layer are i.i.d.. Though the Lyapunov CLT could relax the probabilistic premise to be independence only, all the previous works do not thoroughly discuss the probabilistic premise except assuming an i.i.d. prior for all the parameters of DNNs. Due to its extreme importance but unclear status, it is necessary to thoroughly investigate the probabilistic premise for improving the interpretability of DNNs.
In this work, we demonstrate that the parameters of DNNs are correlated and depend on the training dataset, i.e., they do not satisfy the i.i.d. prior. Moreover, we demonstrate that activations being i.i.d. cannot hold for all the hidden layers of DNNs. In the context of Bayesian probability, we theoretically derive the necessary conditions for the activations of a hidden layer being i.i.d. given the assumption that the activations of the previous layer are i.i.d.. Subsequently, we experimentally show that typical DNNs, such as the Multilayer Perceptron (MLP) and the Convolutional Neural Networks (CNNs), cannot satisfy the necessary conditions. As a result, activations being i.i.d. is not valid for all hidden layers. In other words, GP with the i.i.d. assumption for the parameters of DNNs cannot correctly explain all the hidden layers of DNNs. In addition, some previous works show that GP is sensitive to the curse of dimensionality (Bengio et al., 2005; Hinton et al., 2012) and cannot clarify the hierarchical representation, an essential of deep learning (Matthews et al., 2018).
Alternatively, we propose a novel probabilistic representation for the hidden layers based on the Gibbs distribution (LeCun et al., 2006) and the Product of Experts (PoE) model (Hinton, 2002). The probabilistic representation explains a hidden layer in two aspects: (i) the distribution of a hidden layer can be formulated as a Gibbs distribution, of which the energy function is determined by the activations, and (ii) the connection between two adjacent layers can be modeled by a PoE model. Specifically, the distribution of a single neuron or hidden unit (e.g., convolutional channel) can be expressed as a PoE, in which all the experts are defined by the neurons or hidden units in the previous hidden layer. Compared to the GP explanations, the proposed probabilistic representation provides a more specific explanation for the hidden layers of DNNs without any assumption.
Since the output of a hidden layer is commonly the input of the next layer, DNNs form a Markov chain, thereby corresponding to a joint distribution. Moreover, we confirm the equivalence between the Stochastic Gradient Descent (SGD) and the variational inference (Mandt et al., 2017; Chaudhari & Soatto, 2018) and prove the former is also equivalent to the energy minimization optimization based on the proposed probabilistic representation. As a result, the entire architecture of DNNs can be explained as a Bayesian network (Nielsen & Jensen, 2009; Koski & Noble, 2011). In particular, we demonstrate that the hidden layers close to the input formulate prior distributions of the training dataset and propose a novel regularization approach to improve the generalization performance of DNNs through pre-training the hidden layers corresponding to prior distributions.

heading: 2 PRELIMINARIES
text: We assume X and Y are two random variables and Pθ(X,Y ) = P (Y |X)P (X) is an unknown joint distribution between X and Y , where P (X) describes the prior knowledge of X , P (Y |X) describes the statistical connection betweenX and Y , and θ indicate the parameters of Pθ(X,Y ).
A dataset D = {(xj ,yj)|xj ∈ RM ,yj ∈ RL}Jj=1 is composed of i.i.d. samples generated from Pθ(X,Y ). A neural network with I hidden layers is denoted as DNN = {x;f1; ...;fI ;fY } and trained by D, where (x,y) ∈ D are the input of the DNN and the corresponding training label, and fY is an estimation of the true distribution P (Y |X). As a result, x ∼ P (X) and fY ≈ P (Y |X). The meaning of fi is two-fold: (i) denoting all the neurons in the ith hidden layer, and (ii) representing the activations of the ith hidden layer. In addition, Fi and FY are the random variables corresponding to fi and fY , respectively. In this work, we use the MLP (Figure 1) to derive most theoretical conclusions unless otherwise specified.

heading: 3 BACKGROUND
text: As a fundamental probabilistic model, Gibbs distribution (a.k.a., energy based model, Boltzmann distribution, or renormalization group) formulates the dependence of random variables X through associating an energy to each dependence structure (Geman & Geman, 1984).
P (X;θ) = 1Z(θ)exp[−E(x;θ)] (1) where E(x;θ) is the energy function, θ denotes the parameters, and Z(θ) = ∫ x
exp[−E(x;θ)]dx is the partition function. Especially, a Gibbs distribution can be easily reformulated as various probabilistic models by redefining E(x;θ), e.g., it would become Markov Random Fields (MRFs) model if the energy function is defined as the summation of multiple potential functions, i.e., E(x;θ) = − ∑ k fk(x;θk) (Wang et al., 2013; Goodfellow et al., 2016). Moreover, assuming a single potential function defines an expert, i.e., Fk = exp[−fk(x;θk)], we can regard MRFs as a PoE model, i.e., P (x;θ) = 1Z(θ) ∏ k Fk, where Z(θ) = ∏ k Z(θk) (Hinton, 2002).
A Gibbs distribution has two appealing properties. First, E(x;θ) is the sufficient statistics of P (x;θ) because Z(θ) only relies on θ, i.e., P (x;θ) entirely depends on E(x;θ). The property establishes a connection between P (X;θ) and the deterministic function E(x;θ), which allow us to explain a deterministic hidden layer in a probabilistic way. Second, the energy minimization is the commonly used optimization for learning θ, i.e., θ∗ = arg minθ E(x;θ) (LeCun et al., 2006), which is helpful to explain the training procedure of DNNs, because the energy minimization can be implemented by the gradient descent algorithm as long as E(x;θ) is differentiable.
Numerous efforts have been devoted to explain DNNs from the viewpoint of Gibbs distribution. The classical work is the Restricted Boltzmann Machine (RBM) (Salakhutdinov & Hinton, 2009), which is a simplified Gibbs distribution. Specifically, RBM only considers the connections between visible units x and hidden units h but overlooks the internal connections within x and h, thus it cannot explain complex DNNs, e.g., CNNs. To the best of our knowledge, Mehta & Schwab (2014) first explain the distribution of hidden layers as a Gibbs distribution, and Yaida (2019) indirectly proves the distribution of fully connected layers is a Gibbs distribution. Lin et al. (2017) clarify some properties of DNNs based on the Gibbs distribution. However, all previous works still cannot derive an explicitly Gibbs explanation for complex hidden layers, e.g., convolutional layers.

heading: 4 MAIN RESULTS

heading: 4.1 THE PARAMETERS OF DNNS DO NOT SATISFY THE I.I.D. ASSUMPTION
text: We examine the statistical properties of the parameters of DNNs and show they are correlated. First, we derive a necessary condition for the parameters of DNNs satisfying the i.i.d. prior, Second, we experimentally show that typical DNNs cannot satisfy the necessary condition on the MNIST and CIFAR-10 datasets. Hence, the parameters of DNNs do not satisfy the i.i.d. assumption.
If all the parameters of DNNs are i.i.d., the parameters of each neuron are identically distributed, e.g., {αmn}Mm=1 ∼ P (αn) in the MLP (Figure 1). Since all the parameters are i.i.d., different neurons should be uncorrelated, i.e., ∀n 6= n′, Corr(αn,αn′) = 0. In particular, the i.i.d. assumption enables us to use the sample correlation r(αn,αn′) to estimate Corr(αn,αn′).
r(αn,αn′) = ∑M
m=1(αmn−αn)(αmn′−αn′ )√∑M m=1(αmn−αn) 2∑M m=1(αmn′−αn′ ) 2 ≈ Corr(αn,αn′) (2)
where αn = 1M ∑M m=1 αmn is the sample mean of {αmn}Mm=1. Therefore, taking into account the estimation error, |r(αn,αn′)| close to zero is the necessary condition for all parameters being i.i.d.. We examine the necessary condition for all parameters being i.i.d. in various DNNs. First, we use the MLP to classify the benchmark MNIST dataset. The information about the network architecture and training methods is included in Appendix A.1. After training the MLP well, we calculate the absolute sample correlations, namely, |r(αn,αn′)|, |r(βk,βk′)|, and |r(γl,γl′)|, and derive three absolute correlation matrices, i.e., A′N×N , B ′ K×K , and C ′ L×L, shown in Figure 2. We find that most |r(αn,αn′)| are close to zero, but most |r(βk,βk′)| and |r(γl,γl′)| are non-zero. That indicates the necessary condition cannot hold for all the parameters of the MLP.
Second, we construct a shallow neural network with a single fully connected hidden layer, namely, NN = {x;f ;fY }, for classifying the MNIST dataset based on the seminal work (Neal, 1994) in Appendix A.2, and construct a CNN for classifying the CIFAR-10 dataset in Appendix A.3. Using the same method in the MLP, we obtain the absolute sample correlations between the weights of different nodes in the both neural networks. The results in Appendix A.2 and Appendix A.3 show that the weights of DNNs do not satisfy the necessary condition.
Third, we visualize all the parameters of the NN = {x;f ;fY } for classifying a synthetic dataset in Appendix A.4. Since the synthetic dataset is much simpler than MNIST, we can directly show the correlation of different parameters and their dependences on the dataset by visualizing them. Overall, these simulations demonstrate that the parameters of DNNs are not i.i.d..

heading: 4.2 ACTIVATIONS BEING I.I.D. IS NOT VALID FOR ALL THE HIDDEN LAYERS OF DNNS
text: Although we demonstrate that the parameters of DNNs do not satisfy the i.i.d. assumption, we are still unclear whether the activations of a hidden layer are i.i.d.. In this section, we demonstrate that activations being i.i.d. is not valid for all the hidden layers of DNNs.
In the context of Bayesian probability, all the weights and biases of the MLP (Figure 1) have prior distributions, i.e., αmn ∼ P (Amn), βnk ∼ P (Wnk), b1n ∼ P (B1n), and b2k ∼ P (B2k), thus we regard G1n = ∑M m=1AmnXm +B1n as the random variable of g1n(x) = ∑M m=1 αmn · xm + b1n,
and G2k = ∑N n=1WnkF1n + B2k as the random variable of g2k(f1) = ∑N n=1 βnk · f1n + b2k. Hence, the random variables of two arbitrary activations in f1 and f2, e.g., f1n and f2k, can be expressed as F1n = σ1(G1n) and F2k = σ2(G2k), respectively.
Given two random variablesA andB, the necessary conditions for them being i.i.d. are that they are uncorrelated (i.e., Cov(A,B) = 0) and have the same expectation (i.e., E(A) = E(B)). Similarly, assuming {F1n}Nn=1 are i.i.d. and {f1n}Nn=1 ∼ P (F1), we derive the following necessary conditions for the activations of the second hidden layer, i.e. ,{F2k}Kk=1, being i.i.d..
∀(k, k′) ∈ S1 = {(k, k′) ∈Z2|k 6= k′, 1 ≤ k ≤ K, 1 ≤ k′ ≤ K}, we must have Cov(G2k, G2k′) = Cov(Wk,Wk′)E 2(F1) + Cov(Wk, B2k′)E(F1)
+ Cov(Wk′ , B2k)E(F1) + Cov(B2k, B2k′) = 0
(3)
E(F1)E(Wk) + E(B2k) = E(F1)E(Wk′) + E(B2k′) (4)
σ2(·) is strictly increasing, differentiable, and invertible. (5) where Wk = ∑N n=1Wnk. Equation 3 and 4 are the necessary conditions for {G2k}Kk=1 being i.i.d.. Equation 5 constrains the activation function σ2(·) to guarantee that if {G2k}Kk=1 are i.i.d. then {F2k}Kk=1 are i.i.d. as well. The necessary conditions are valid for arbitrary fully connected layers as long as properly changing the subscripts. The detailed derivations are included in Appendix B and Appendix C. Moreover, we prove that the necessary conditions also hold for convolutional layers based on the connection between fully connected layers and convolutional layers in Appendix E.
It is noteworthy that the i.i.d. assumption for the parameters of DNNs indeed satisfy the necessary conditions for activations being i.i.d.. Nevertheless, since the i.i.d. prior is not an appropriate prior for the parameters of DNNs, it cannot guarantee activations of a hidden layer being i.i.d..
The theoretical necessary conditions indicate an applicable approach to examine whether activations of a hidden layer being i.i.d. given the assumption that the activations of the previous layer are i.i.d.. More specifically, since {F1n}Nn=1 are i.i.d., we can use the sample mean f1 to estimate E(F1). Because the samples of the training dataset D are also assumed to be i.i.d., we can train the MLP multiple times to draw multiple independent observations of Wnk and B2k, thus we can use the sample covariance to estimate Cov(Wk,Wk′), Cov(Wk, B2k′), and etc.
Based on the empirical approach, we estimate the statistical measures (Table 1) to examine whether the activations of the three layers (i.e., f1, f2, and fY ) being i.i.d. in the MLP for classifying the MNIST dataset. Without loss of generality, we restrict all the layers of the MLP from using bias to decrease computation complexity, thus the necessary conditions for f2 can be simplified as
∀(k, k′) ∈ S1 = {(k, k′) ∈ Z2|k 6= k′, 1 ≤ k ≤ K, 1 ≤ k′ ≤ K}, we must have Cov(G2k, G2k′) = Cov(Wk,Wk′)E 2(F1) = 0 (6)
E(F1)E(Wk) = E(F1)E(Wk′) (7)
First, Figure 3 shows the correlation matrices, i.e., AN×N , WK×K , and CL×L, to check if their corresponding covariances, i.e, Cov(An, An′), Cov(Wk,Wk′), and Cov(Ck, Ck′), are close to zero. It shows that many correlations are far from zero, so their corresponding covariances are not close to zero either. Second, we use the sample mean, e.g., x, to estimate the expectation, e.g., E(X). Taking into account the estimation error, we can regard E(X) = 0 based on x = 0.131. However, f1 = 0.436 and f2 = 0.498 imply E(F1) 6= 0 and E(F2) 6= 0 even considering the estimation error. Therefore, the activations of the last two hidden layer, i.e., f2 and fY , cannot be independent because they cannot satisfy the first necessary condition (Equation 6). In addition, we show the necessary condition for activations being identically distributed also cannot be satisfied by the MLP in Appendix F.1. Furthermore, we demonstrate that the necessary conditions for activations being i.i.d. also cannot be satisfied in the CNN for classifying the CIFAR-10 dataset in Appendix F.2
Overall, we conclude that activations being i.i.d. is not valid for all the hidden layers of DNNs based on the theoretical necessary conditions and the experimental simulations on typical DNNs. In other words, we cannot use the CLT to establish the equivalence between GP and all the hidden layers of DNNs, thus GP can not correctly clarify all the hidden layers of DNNs. That necessitates a more general probabilistic representation for explaining the hidden layers of DNNs.

heading: 4.3 THE DISTRIBUTION OF A HIDDEN LAYER IS A GIBBS DISTRIBUTION
text: This section proves that the distribution of a hidden layer can be formulated as a Gibbs distribution and the statistical connection between two adjacent layers can be modeled as a PoE model.
Assuming the output layer in the MLP (Figure 1) is the softmax, its distribution can be expressed as
P (FY ) = {P (fyl) = 1
ZY exp(fyl)}Ll=1 (8)
where ZY = ∑L l=1 exp(fyl) is the partition function. Since fyl = ∑K k=1 γkl · f2k + byl, we have
P (FY ) = {P (fyl) = 1
ZY exp( K∑ k=1 γkl · f2k + byl)}Ll=1 (9)
Based on the properties of the exponential function, i.e., exp(a+ b) = exp(a) · exp(b) and exp(a · b) = [exp(b)]a, we can reformulate P (FY ) as
P (FY ) = {P (fyl) = 1
Z ′Y K∏ k=1 [exp(f2k)]γkl}Ll=1 (10)
where Z ′Y = ZY /exp(byl). Since {exp(f2k)}Kk=1 are scalar, we can introduce a new partition function ZF2 = ∑K k=1 exp(f2k) such that {
1 ZF2 exp(f2k)}Kk=1 become a probability measure. As a result, we can further reformulate P (FY ) = {P (fyl)}Ll=1 as a Product of Expert (PoE) model as
P (FY ) = {P (fyl) = 1
Z ′′Y K∏ k=1 [ 1 ZF2 exp(f2k)]γkl}Ll=1 (11)
where Z ′′Y = ZY /[exp(byl) · ∏K k=1[ZF2 ]
γkl ] and each expert is defined as 1ZF2 exp(f2k).
It is noteworthy that all the experts { 1ZF2 exp(f2k)} K k=1 form a probability measure and establish an exact one-to-one correspondence to all the neurons in the second hidden layer f2, i.e., {f2k}Kk=1. Therefore, the distribution of f2 can be expressed as
P (F2) = {P (f2k) = 1
ZF2 exp(f2k)}Kk=1 (12)
Based on the definition of Gibbs distribution, Equation 8 and 12 show that fY and f2 formulate two multivariate Gibbs distributions and their energy functions are equivalent to the negative of the output nodes {fyl}Ll=1 and the neurons {f2k}Kk=1, respectively. In addition, Equation 11 implies that the connection between P (fyl) and P (f2k) can be expressed as the PoE model, namely that P (fyl) is equal to the product of all the experts {P (f2k)}Kk=1 with their respective weights {γkl}Kk=1. Moreover, we prove that the probabilistic representation is valid for other hidden layers in the MLP. Since {f2k = σ2( ∑N n=1 βnk · f1n + b2k)}Kk=1, P (F2) can be expressed as
P (F2) = {P (f2k) = 1
ZF2 exp[σ2( N∑ n=1 βnk · f1n + b2k)]}Kk=1 (13)
Based on the equivalence between the gradient descent algorithm and the first order approximation (Battiti, 1992), we prove that σ2( ∑N n=1 βnk · f1n + b2k)] can be approximated as
σ2( N∑ n=1 βnk · f1n + b2k)] ≈ C21 · [ N∑ n=1 βnk · f1n + b2k] + C22 (14)
where C21 and C22 only depend on the activations {f1n}Nn=1 in the previous training iteration, thus they can be regarded as constants and absorbed by βnk and b2k. The proof for the approximation is included in Appendix G. Therefore, P (F2) also can be modeled as a PoE model.
P (F2) = {P (f2k) ≈ 1
Z ′′F2 N∏ n=1 [ 1 ZF1 exp(f1n)]βnk}Kk=1 (15)
where Z ′′F2 = ZF2/[exp(b2k) · ∏N n=1[ZF1 ] βnk ] and ZF1 = ∑N n=1 exp(f1n).
Finally, the distribution of the first hidden layer f1 is formulated as
P (F1) = {P (f1n) = 1
ZF1 exp(f1n)}Nn=1 (16)
Overall, we introduce a novel probabilistic representation for fully connected layers in two aspects. First, a fully connected layer formulates a multivariate Gibbs distribution (Equation 8, 12, and 16) and the energy functions is equivalent to the negative activation. Second, the connection between two adjacent fully connected layers can be formulated as a PoE model (Equation 11 and 15).
The proposed probabilistic representation has two advantages over the GP explanations. First, the FoE explanation does not require any assumption, e.g., we formulate P (f2k) as a PoE model, in which all the experts { 1ZF1 exp(f1n)} N n=1 are defined by all the neurons {f1n}Nn=1. As a comparison, the GP explanation derives P (f2k) as a Gaussian distribution in the limit of infinite neurons, i.e., {f1n}N→∞n=1 , under the assumption that {f1n}Nn=1 are i.i.d.. Second, the Gibbs distribution provides a more explicit explanation for the functionality of a fully connected layer. Specifically, a Gibbs distribution derives the probability of the multiple features defined by the neurons occurring in the input, but the covariance function of a GP cannot clearly explain the functionality of each neuron.
Moreover, we prove that a convolutional layer formulates a specific Gibbs distribution, namely the MRF model, in Appendix H. More specifically, a convolutional layer fi with K convolutional channels, i.e., fi = {fki = σi( ∑N n=1 S (k,n) i ◦ fni−1 + bki · 1)}Kk=1, formulate a MRF model as
P (Fi) = 1
ZFi K∏ k=1 exp[ϕki (fi−1)] = 1 ZFi exp[ K∑ k=1 ϕki (fi−1)] (17)
where fi−1 = {fni−1}Nn=1 is the input, ϕki (fi−1) = ∑N n=1 S
(k,n) ◦fni−1 + bki ·1) is the kth linear convolutional channel, S(k,n) is the kth convolutional filter applying to the nth input channel, bki is the bias, σi(·) is an non-linear operator, and the partition function ZFi= ∫ exp[ ∑K k=1 ϕ
k(fi−1)]dfi−1. The corresponding energy function EFi=− ∑K k=1 ϕ k(fi−1).
In summary, the distribution of a hidden layer fi can be formulated as a Gibbs distribution P (Fi), especially the energy function EFi establish a connection between the deterministic function of the hidden layer and the Gibbs distribution P (Fi). The correspondence between Gibbs distributions and hidden layers can be summarized by Table 2.

heading: 5 SIMULATIONS
text: This section uses two synthetic datasets to demonstrate the proposed probabilistic representation: (i) the distribution of a fully connected layer is a multivariate discrete Gibbs distribution and (ii) the distribution a convolutional layer is a MRF model.

heading: 5.1 THE DISTRIBUTION OF A FULLY CONNECTED LAYER IS A MULTIVARIATE DISCRETE GIBBS DISTRIBUTION
text: A fully connected layer with N neurons, i.e., f = {fn}Nn=1, formulates a multivariate discrete Gibbs distribution to describe the probability of the N features defined by their respective neurons occurring in the input. The energy function is equal to the negative activation, i.e., Efn = −fn(x).
P (F ) = {P (fn) = 1
ZF exp[fn(x)]}Nn=1 (18)
where fn(x) = σ[gn(x)], gn(x)= ∑M
m=1 αmn·xm+bn is a linear filter with αmn being the weights and bn being the bias, σ(·) is the activation function, and ZF=∑Nn=1 exp[fn(x)] is the partition function. Since benchmark datasets consist of very complex features and typical DNNs have too complicated architecture, it is very hard to directly demonstrate the Gibbs distribution explanation. Alternatively, we choose the neural network for classifying the synthetic dataset in Appendix A.4, because the neural network only has a single fully connected hidden layer with seven neurons, i.e., N = 7, and the synthetic dataset only has four features. After the network is trained well, we choose a synthetic image as the input and show all the information of the learned hidden layer in Figure 4.
Since the input image is sorted in the primary diagonal direction, the input feature can be described as high values above the secondary diagonal (the red line) and low values below the red line. Figure 4 shows that Neuron 1, Neuron 2 and Neuron 4 correctly describe the input feature, thus they derive large positive outputs (i.e., g1(x) = 122.09, g2(x) = 113.70, and g4(x) = 140.20). However, Neuron 6 and Neuron 7 do not describe the input feature correctly, thus they derive large negative outputs (i.e., g6(x) = −134.08, g7(x) = −133.36). After applying the sigmoid function, we have f1(x) = f2(x) = f4(x) = 1 and f6(x) = f7(x) = 0. Finally, we derive the probability of different neurons, i.e., P (f1) = P (f2) = P (f4) = 0.20 and P (f6) = P (f7) = 0.07, which indicates that the features defined by Neuron 1, Neuron 2 and Neuron 4 have high probability occurring in the input but the features defined by Neuron 6 and Neuron 7 have low probability occurring in the input.
Under review as a conference paper at ICLR 2020 Binary 2
Binary 1
Binary 8
Binary 0
Binary 5

heading: 5.2 THE DISTRIBUTION OF A CONVOLUTIONAL LAYER IS A MRF MODEL
text: Since the distributions of benchmark datasets are unknown, it is impossible to use them to verify the probabilistic explanation for convolutional layers. Alternatively, we generate a synthetic dataset obeying the Gaussian distribution N (0, 1024) based on the NIST dataset of handwritten digits. It consists of 20,000 32 × 32 grayscale images in 10 classes (digits from 0 to 9), and each class has 1,000 training images and 1,000 testing images. Figure 5 shows five synthetic images and their perspective histograms. The method for generating the synthetic dataset is included in Appendix I. We design a CNN = {x;f1;f2;f3;f4;fY } for classifying the synthetic dataset. The CNN has two convolutional layers with ReLU (i.e., f1 and f3), two max-pooling layers (i.e., f2 and f4) and a softmax output layer fY . The detailed information about the CNN is included in Appendix I.
Based on the MRF explanation (Equation 17), the distribution of f1 can be formulated as
P (F1) = 1
ZF1 20∏ k=1 exp[ϕk1(x)] = 1 ZF1 exp[ 20∑ k=1 ϕk1(x)] (19)
where ϕk1(x) = S k ◦x+bk1 ·1 is the kth linear convolutional channel, Sk is the linear convolutional filter and bk1 is the bias. The energy function of P (F1) is equivalent to EF1 = − ∑20 k=1 f k 1 .
In order to model a high dimensional dataset, e.g., the 32 × 32 synthetic image, we typically need a multivariate distribution with the same dimension, e.g., the dimension of the covariance function of a Gaussian process should be 1024 × 1024. However, most structured probabilistic graphical models, especially the MRF model, have two fundamental assumptions: stationary and Markov (Lyu & Simoncelli, 2007). 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 the histogram of samples to simulate the true distribution.
After the CNN is well trained, we randomly choose a testing image x as the input to derive P (F1). Since x ∼ P (X), we can use the histogram of the testing image, i.e., P (x), to estimate P (X). In addition, EF1 is a sufficient statistics for P (F1) such that we can use the histogram of EF1 to estimate P (F1). As a result, we can verify the probabilistic representation by calculating the distance between P (F1) and P (X), i..e, KL[P (X)||P (F1)]. We visualize their distributions in Figure 6, which shows that P (F1) is very close to P (X), i.e., KL(P (X)||q(F1)) = 0.83, which demonstrates the distribution of a convolutional layer can be modeled as a MRF model.
In addition, the experiment provides a new viewpoint to demonstrate the activations of different convolutional channels are not i.i.d.. Given another synthetic image, we visualize the output of five convolutional channels in f1 in Figure 7. It shows that the activations of different channels at the same position are correlated, i.e., they are not independent. We also estimate the distribution of each convolutional channel based on the histograms of their respective energy functions. It shows that the distributions of different channels are different, i.e., they are not identically distributed.

heading: 6 DISCUSSION
text: Based on the proposed probabilistic explanations for hidden layers, we discuss certain theoretical and practical topics of deep learning. First, we further confirm the equivalence between SGD and the variational inference. Second, we demonstrate that the entire architecture of DNNs can be explained as a Bayesian network, in which the hidden layers close to the input formulate prior distributions of the training dataset. Third, we propose a novel regularization method to improve the generalization performance of DNNs by pre-training the hidden layers corresponds to the prior distributions.

heading: 6.1 SGD CAN BE EXPLAINED AS A VARIATIONAL INFERENCE
text: Recent works prove that SGD works as a specific Bayesian posterior inference, i.e., the variational inference (Mandt et al., 2017; Chaudhari & Soatto, 2018). Here we further confirm the Bayesian explanation for SGD based on the proposed probabilistic representation. Briefly speaking, given a DNN = {x;f1; ...;fI ;fY }, we prove that the conditional distribution of fY given x, i.e., P (FY |X), is also a Gibbs distribution and its energy function is equivalent to the deterministic function formulated by the entire architecture of the DNN, i.e., the DNN describes a family of Gibbs distribution Q. Based on the definition of variational inference, SGD aims to find a Gibbs distribution P ∗(FY |X) such that it can minimize the distance between P (FY |X) and the true posterior distribution P (Y |X) is determined by the training datasetD.
P ∗(FY |X) = arg min P∈Q KL[P (Y |X)||P (FY |X)] (20)
Moreover, we prove that SGD is also equivalent to the energy minimization optimization, which in turn validates the Gibbs explanation for the hidden layer of DNNs. The proof is in Appendix J.
We justify the theoretical explanation for SGD based on the CNN classifying the synthetic dataset and the MLP classifying the MNIST dataset in Appendix K.1. In addition, we demonstrate that the learned parameters via SGD are not i.i.d. based on in Appendix K.2, which indirectly justifies that GP with the i.i.d. assumption for the parameters of DNNs cannot correctly clarify DNNs.

heading: 6.2 THE ENTIRE ARCHITECTURE OF DNNS CAN BE EXPLAINED AS A BAYESIAN NETWORK
text: Since the output of the hidden layer fi in the DNN = {x;f1; ...;fI ;fY } is commonly the input of the next layer fi+1, the DNN formulates a Markov chain F1→ · · · → FI → FY . As a result, the distribution of a hidden layer in the DNN should be regarded as a conditional distribution and the entire architecture of the DNN can be formulated as a joint distribution.
P (FY ;FI ..;F1|X) = P (FY |FI) · ... · P (Fi+1|Fi) · ... · P (F1|X) (21)
Richard & Lippmann (1991) prove that P (FY |X) is an estimation of the true posterior distribution P (Y |X) ∝ P (X|Y )P (Y ) and we demonstrate that P (F1|X) is a prior distribution of the training dataset and the equivalence between SGD and the variational inference, thus the entire architecture of DNNs can be explained as a Bayesian network. That provides a probabilistic explanation for the hierarchy property of DNNs, which cannot be clarified by GP (Matthews et al., 2018).

heading: 6.3 A NOVEL REGULARIZATION ALGORITHM FOR DNNS
text: Bayesian theory indicates that a better prior distribution corresponds to a better posterior distribution given the same likelihood distribution, thus pre-training the hidden layers corresponding to prior distributions should be an effective approach to improve the generalization performance of DNNs. Erhan et al. (2010) first demonstrate that the unsupervised pre-training is an effective regularization for DNNs through capturing the prior knowledge of the input. However, since they do not know which hidden layer corresponds to a prior distribution, they have to pre-training all the hidden layers, which restricts its wide application due to the computation complexity. In addition, Le et al. (2018) prove that incorporating the prior knowledge of the input results in uniform stability and provides a bound on generalization error for a specific DNN, i.e., the supervised auto-encoders.
Based on the proposed probabilistic representation, we first justify that a better prior distribution indeed has a positive correlation to the generalization performance in Appendix K.3. Subsequently, we propose a novel a novel regularization algorithm to improve the generalization performance of DNNs through pre-training the hidden layers only corresponding to prior distributions and show that it outperforms classical regularizations, e.g., dropout (Srivastava et al., 2014), for the image recognition task on the benchmark CIFAR-10 dataset in Appendix L.

heading: 7 CONCLUSION
text: In this paper, we demonstrate that the parameters of DNNs do not satisfy the i.i.d assumption and the activations being i.i.d. cannot be valid for all the hidden layers of DNNs. As a result, GP cannot correctly explain all the hidden layer of DNNs. Alternatively, we propose a novel probabilistic representation for DNNs in four aspects: (i) a hidden layer formulates a Gibbs distribution, of which the energy function is determined by the activations, (ii) the connection between two adjacent layers can be modeled by a PoE model, (iii) the entire architecture of DNNs can be explained as a Bayesian network, and (iv) the SGD can be explained as the Bayesian variational inference algorithm.
There are two general directions for future research. The theoretical direction is to examine the probabilistic representation in other popular DNNs, e.g., the residual networks, which can further verify the proposed probabilistic representation. The practical direction is to validate the proposed regularization in more complex DNNs and datasets.

heading: A THE PARAMETERS OF DNNS CANNOT SATISFY THE I.I.D. PRIOR

heading: A.1 A MLP ON THE MNIST DATASET
text: We design a MLP (Figure 1) to classify the benchmark MNIST dataset. The MLP has two hidden layers, each hidden layer is a fully connected layer with the sigmoid activation function, and the output layer is the softmax. Table 3 summarizes the architecture of the MLP.
After the MLP is well trained (its training and testing accuracy is 97.7% and 97.1%, respectively), we quantify the correlation of different weights in the same layer based on Equation 22 and derive the correlation matrix for all the weights in every hidden layer.
r(αn,αn′) = ∑M m=1(αmn − αn)(αmn′ − αn′)√∑M
m=1(αmn − αn) 2∑M m=1(αmn′ − αn′) 2 ≈ Corr(αn,αn′) (22)
where αn = 1M ∑M m=1 αmn is the sample mean of {αmn}Mm=1.
To examine if different weights are correlated, we only need to check whether their correlation coefficients are close to zero, so we take the absolute value of the correlation matrix element-wisely. Figure 8 visualizes the correlation matrices of each hidden layer, i.e., A′N×N , B ′ K×K , and C ′ L×L.
We find that most weights in f1 are uncorrelated because their correlation coefficients are very close to zero, i.e., |r(αn,αn′)| ≈ 0, but many weights in f2 and fY are correlated to others, i.e., |r(βk,βk′)| 6= 0 and |r(γl,γl′)| 6= 0. In addition, we derive the average of the absolute correlation coefficients of all the weights in each layer and show them at the top of Figure 8, which further validates that the weights of some hidden layer are correlated.
Therefore, the parameters of the MLP cannot satisfy the i.i.d. prior.

heading: A.2 A NEURAL NETWORK WITH A SINGLE HIDDEN LAYER ON THE MNIST DATASET
text: Based on the seminal work (Neal, 1994), we design a neural network with a single fully connected hidden layer (Figure 9) for classifying the benchmark MNIST dataset (LeCun et al., 1998). The input layer dimension is 1 × 784, the hidden layer has 128 neurons with the sigmoid activation function, and the output layer is the softmax with 1×10 dimension. Table 4 summarizes the architecture of the network. In addition, {αmn}Mm=1 denote the weights of the nth linear filter in the fully conntected hidden layer and {βnl}Nn=1 are the weights of the lth linear filter in the output layer. Therefore, the hidden layer has 128 linear filters with 1 × 784 dimension, i.e., {αmn}784m=1. The output layer has totally 10 linear filters and their dimension is 1× 128, i.e., {βnl}128n=1.
After the network is well trained (its training and testing accuracy is 97.4% and 96.9%, respectively), we calculate the absolute sample correlations, i.e., |r(αn,αn′)| and |r(βk,βk′)|, and derive two absolute correlation matrices, i.e., A′N×N and B ′ L×L, based on Equation 2. They are shown in Figure 10. We find that most |r(αn,αn′)| are close to zero, but |r(βk,βk′)| not. That indicates the parameters of the NN do not satisfy the necessary condition, thus they cannot satisfy the i.i.d. prior.

heading: A.3 A CNN ON THE CIFAR-10 DATASET
text: We design a CNN to classify the benchmark CIFAR-10 dataset (Krizhevsky, 2009). The CNN has two convolutional layers without activation function, two max-pooling layers for dimension reduction, and two fully connected layers with the sigmoid activation function. In addition, the output layer is the softmax. Table 5 summarizes the architecture of the CNN.
After we train the CNN well (its training and testing accuracy are 98.0% and 69.4%, respectively). We quantify the correlation of different weights in the same layer based on Equation 2, thereby deriving the absolute correlation matrix for all the weights in every hidden layer.
Since x has 3 channels and f1 has 64 channels, f1 has totally 3 × 64 = 192 filters. Hence, the dimension of the correlation matrix for f1 is 192 × 192. In particular, we flatten the output of the second convolutional layer, i.e., reshaping 8×8×128 to 1×8192, thus the filter dimension of f3 is 8192× 1024. Table 5 shows the dimension of the correlation matrices for all the layers in the CNN. Figure 11 shows the absolute correlation matrices for every layer. We can find that most weights are correlated in f1 and f2, i.e., the weights in the convolutional layers are correlated. In addition, the weighs in the output layer are also shown to be correlated. However, the weights in f3 and f4 are not correlated, because their correlation coefficients are close to zero. In other words, the weights in the fully connected layers are uncorrelated. Moreover, we derive the average of the absolute correlation coefficients of all the weights in each layer and show them at the top of Figure 11, which further validates the correlation of the weights in each layer. Overall, the parameters of the CNN cannot satisfy the necessary condition. Therefore, they cannot satisfy the i.i.d. prior.

heading: A.4 NEURAL NETWORKS WITH A SINGLE HIDDEN LAYER ON A SYNTHETIC DATASET
text: In order to further demonstrate that the parameters of DNNs cannot satisfy the i.i.d. assumption, we construct a neural network, i.e., NN = {x;f ;fY } shown in Figure 9, to classify a synthetic dataset. The synthetic dataset consists of 4,800 32× 32 grayscale images with 4 classes, and every synthetic image is sampled from the Gaussian distribution N (0, 1). To make the synthetic image have certain features that can be learned by neural networks, we sort the synthetic image in the primary or secondary diagonal directions by the ascending or descending orders. Thus, the synthetic dataset only has four features, i.e., four classes, and their corresponding images are drawn above.
There are two important reasons why we choose the synthetic dataset. Above all, compared to benchmark datasets containing too complex features, the synthetic dataset only has four simple features, thus we can easily show the correlation of different parameters by simply visualizing them. In addition, since a single synthetic image only has a single feature, we can easily examine the connection between the parameters of networks and the synthetic image.
Since the synthetic dataset is very simple, we set the hidden layer f have only eights neurons and the output layer fY have four output nodes. After training the network well, we show all the eight neurons in Figure 13. In contrast to the i.i.d. assumption, their weights show strong correlation in three aspects. First, weights show strong internal correlation within each neuron, i.e., if a weight has large magnitude, its neighbors have high probability to be large, and vice versa. Second, weights show strong external correlation between different neurons. For example, most weights of Neuron1 and Neuron5 at the same position show strong positive correlation, and most weights of Neuron3 and Neuron7 at the same position show strong negative correlation. Third, the weights show strong dependence on the training dataset. More specifically, since all the features of the synthetic dataset concentrate in the primary and secondary diagonal directions, most linear filters also demonstrate related features in the same directions.
Figure 14 visualizes the four linear filters in the output layer. Though they not show too strong internal correlation because their dimension is too small, their weights at the same position still show strong external correlation, e.g., the weights of Filter0 and Filter1 at the same position show strong negative correlation. Moreover, we calculate the correlation coefficients of all the four filters based on Equation 2 and derive the correlation matrix show in Figure 14. It shows that the correlation coefficient of Filter0 and Filter1 is close to−1, which further validates that they have strong negative correlation, hence the parameters of the output layer do not satisfy the i.i.d. prior.
Compared to the simulations using benchmark datasets, the synthetic dataset enable us to expose the correlation of different weights and their dependence on the training dataset in the more clearly way. Overall, we demonstrate that the parameters of neural networks are not i.i.d..
B PROOF: THE NECESSARY CONDITION FOR THE ACTIVATIONS OF HIDDEN LAYERS BEING INDEPENDENT
In the context of Bayesian probability, all the weights and the biases in the MLP (Figure 15) have prior distributions, i.e., αmn ∼ P (Amn), βnk ∼ P (Wnk), b1n ∼ P (B1n), and b2k ∼ P (B2k), thus we regardG1n = ∑M m=1AmnXm+B1n as the random variable of g1n(x) = ∑M m=1 αmn·xm+b1n,
and G2k = ∑N n=1WnkF1n + B2k as the random variable of g2k(f1) = ∑N n=1 βnk · f1n + b2k. Therefore, the random variables of arbitrary activations in f1 and f2, e.g., f1n and f2k, can be expressed as F1n = σ1(G1n) and F2k = σ2(G2k), respectively.
Based on the theorem (Appendix D) that functions of independent random variables are independent, we can derive that if {F2k}Kk=1 are independent and the activation function σ1(·) is invertible, then {G2k}Kk=1 are independent. In other words, in order to prove {F2k}Kk=1 are independent, we only need to prove {G2k}Kk=1 being independent as long as σ2(·) is invertible.
The necessary condition for {G2k}Kk=1 being independent is ∀(k, k′) ∈ S1 = {(k, k′) ∈ Z2|k 6= k′, 1 ≤ k ≤ K, 1 ≤ k′ ≤ K}, the covariance Cov(G2k, G2k′) = 0, which can be formulated as
Cov(G2k, G2k′) = Cov( N∑ n=1 WnkF1n +B2k, N∑ n′=1 Wn′k′F1n′ +B2k′) = 0 (23)
Based on Cov(X + Y,W + Z) = Cov(X,W ) + Cov(X,Z) + Cov(Y,W ) + Cov(Y,Z), we can derive that
Cov(G2k, G2k′) = N∑ n=1 N∑ n′=1 Cov(WnkF1n,Wn′k′F1n′) + N∑ n=1 Cov(WnkF1n, B2k′)
+ N∑ n′=1 Cov(Wn′k′F1n′ , B2k) + Cov(B2k, B2k′)
(24)
Based on the definition of covariance, i.e., Cov(A,B) = E(AB)− E(A)E(B), we have Cov(WnkF1n,Wn′k′F1n′) = E(WnkWn′k′ · F1nF1n′)− E(WnkF1n)E(Wn′k′F1n′)
Cov(WnkF1n, B2k′) = E(WnkB2k′F1n)− E(WnkF1n)E(B2k′) (25)
Based on the law of total expectation, i.e, E(AB) = EA(AEB|A(B)), we have
E(WnkWn′k′ · F1nF1n′) = EWnkWn′k′ (WnkWn′k′ · EF1nF1n′ |WnkWn′k′ (F1nF1n′)) E(WnkB2k′F1n) = EWnkB2k′ (WnkB2k′ · EF1n|WnkB2k′ (F1n))
E(WnkF1n) = EWnk(Wnk · EF1n|Wnk(F1n)) (26)
Since the statistical properties of the input {F1n}Nn=1 are unrelated to the weights βnk and biases b2k, we have EF1nF1n′ |WnkWn′k′ (F1nF1n′) = EF1nF1n′ (F1nF1n′), EF1n|WnkB2k′ (F1n) = EF1n(F1n), and EF1n|Wnk(F1n) = EF1n(F1n). As a result, we can derive that
E(WnkWn′k′ · F1nF1n′) = EWnkWn′k′ (WnkWn′k′)EF1nF1n′ (F1nF1n′) E(WnkB2k′F1n) = EWnkB2k′ (WnkB2k′)EF1n(F1n)
E(WnkF1n) = EWnk(Wnk)EF1n(F1n)
(27)
Substituting Equation (27) for the corresponding expectations in Equation (25), we have
Cov(WnkF1n,Wn′k′F1n′) = E(WnkWn′k′)E(F1nF1n′)− E(Wnk)E(Wn′k′)E(F1n)E(F1n′) Cov(WnkF1n, B2k′) = E(WnkB2k′)E(F1n)− E(Wnk)E(B2k′)E(F1n)
(28)
Assuming {F1n}Nn=1 are i.i.d. and F1n ∼ P (F1), we have E(F1nF1n′) = E(F1n)E(F1n′) = E2(F1). Hence, we can derive that
Cov(WnkF1n,Wn′k′F1n′) = Cov(Wnk,Wn′k′)E 2(F1)
Cov(WnkF1n, B2k′) = Cov(Wnk, B2k′)E(F1) (29)
Substituting Equation (29) for the corresponding covariances in Equation (24), we have
Cov(G2k, G2k′) = N∑ n=1 N∑ n′=1 Cov(Wnk,Wn′k′)E 2(F1) + N∑ n=1 Cov(Wnk, B2k′)E(F1)
+ N∑ n′=1 Cov(Wn′k′ , B2k)E(F1) + Cov(B2k, B2k′)
(30)
Based onCov(X + Y,W + Z) = Cov(X +W ) + Cov(X + Z) + Cov(Y +W ) + Cov(Y + Z), we can derive that
Cov(G2k, G2k′) = Cov( N∑ n=1 Wnk, N∑ n′=1 Wn′k′)E 2(F1) + Cov( N∑ n=1 Wnk, B2k′)E(F1)
+ Cov( N∑ n′=1 Wn′k′ , B2k)E(F1) + Cov(B2k, B2k′)
(31)
Therefore, given σ2(·) is invertible and {F1n}Nn=1 are i.i.d., the necessary condition for {F2k}Kk=1 being independent can be expressed as
Cov(G2k, G2k′) = Cov(Wk,Wk′)E 2(F1) + Cov(Wk, B2k′)E(F1)
+ Cov(Wk′ , B2k)E(F1) + Cov(B2k, B2k′) = 0 (32)
where Wk = ∑N n=1Wnk and Wk′ = ∑N n=1Wnk′ .
Overall, in order to guarantee the activations of a hidden layer are independent with each other in the context of Bayesian probability, the prior distributions of the weights and biases of the hidden layer should satisfy the above necessary condition given the inputs of the hidden layer are i.i.d. and the activation functions are invertible.

heading: C PROOF: THE NECESSARY CONDITION FOR THE ACTIVATIONS OF HIDDEN LAYERS BEING IDENTICALLY DISTRIBUTED
text: Given the same DNN drawn in Figure 15, we assume σ2(·) being strictly increasing and differentiable, so σ2(·) is invertible and its inverse σ−12 (·) is also strictly increasing. The cumulative distribution function of F2k can be expressed as
ΦF2k(f) = φ(F2k ≤ f) = φ(σ1(G2k) ≤ f) = φ(G2k ≤ σ−11 (f)) = ΦG2k(σ −1 1 (f))
(33)
where φ(F2k ≤ f) is the probability of F2k takes on a value less than or equal to f . Subsequently, we can obtain that
PF2k(f) = ∂ΦF2k(f) ∂f = ∂ΦG2k(σ
−1 2 (f))
∂f
= PG2k(σ −1 2 (f))
∂σ−12 (f)
∂f
(34)
Equation (34) means that if the activation function σ2(·) is strictly increasing and differentiable, {F2k}Kk=1 being identically distributed implies {G2k}Kk=1 being identically distributed.
In the context of Bayesian probability, G2k = ∑N n=1WnkF1n +B2k, we can derive that
E(G2k) = N∑ n=1 E(WnkF1n) + E(B2k) (35)
Based on the law of total expectation and Equation 27, we have
E(WnkF1n) = EWnk(Wnk)EF1n(F1n) (36)
Assuming {F1n}Nn=1 are i.i.d. and F1n ∼ F1, we have E(F1n) = E(F1). Hence, we can derive that
E(G2k) = E(F1) N∑ n=1 E(Wnk) + E(B2k) (37)
Based on the property of expectation, i.e., E(A) + E(B) = E(A+B), we have
E(G2k) = E(F1)E( N∑ n=1 Wnk) + E(B2k) (38)
Therefore, given σ2(·) is strictly increasing and differentiable and {F1n}Nn=1 are i.i.d., the necessary condition for {F2k}Kk=1 being identically distributed is that ∀(n, n′) ∈ S1, we have
E(F1)E(Wk) + E(B2k) = E(F1)E(Wk′) + E(B2k′) (39)
where Wk = ∑N n=1Wnk and Wk′ = ∑N n=1Wnk′ .
Overall, in order to guarantee the activations of a hidden layer are identically distributed, the prior distributions of the weights and the biases of the hidden layer should satisfy the above necessary condition in the context of Bayesian probability given the inputs of the hidden layer are i.i.d. and the activation functions are strictly increasing and differentiable.

heading: D THEOREM: FUNCTIONS OF INDEPENDENT RANDOM VARIABLES ARE INDEPENDENT
text: Let X and Y be independent random variables on a probability space (Ω,F ,P). Let g and h be real-valued functions defined on the codomains of X and Y , respectively. Then g(X) and h(Y ) are independent random variables.
Proof: Let A ⊆ R and B ⊆ R be the range of g and h, the joint distribution between g(X) and h(Y ) can be formulate as P(g(X) ∈ A, h(Y ) ∈ B). Let g−1(A) and h−1(B) denote the preimages ofA andB, respectively, we have
P(g(X) ∈ A, h(Y ) ∈ B) = P(X ∈ g−1(A),Y ∈ h−1(B)) (40)
Based on the definition of independence, we can derive that
P(g(X) ∈ A, h(Y ) ∈ B) = P(X ∈ g−1(A))P(Y ∈ h−1(B)) = P(g(X) ∈ A)P(h(Y ) ∈ B)
(41)
Based on the definition of preimage, we can derive that
P(g(X) ∈ A, h(Y ) ∈ B) = P(g(X) ∈ A)P(h(Y ) ∈ B) (42)
Therefore, g(X) and h(Y ) are independent random variables.

heading: E THE NECESSARY CONDITIONS FOR ACTIVATIONS BEING I.I.D. ARE VALID FOR CONVOLUTIONAL LAYERS
text: Assuming a DNN = {x;f1; ...;fI ;fY } has two adjacent convolutional layers fi and fi+1. The first convolutional layer hasN channels, i.e. fi = {fni }Nn=1, and the second convolutional layer has K channels, i.e. fi+1 = {fki+1}Kk=1. Therefore, the dimension of the convolutional filters Si+1 in fi+1 are W ×H ×N ×K, i.e., the dimension of a single filter connecting the fni channel and the fki+1 channel are D(S (k,n) i+1 ) = W ×H . As a result, the fki+1 channel can be formulated as
fki+1 = σi+1( N∑ n=1 S (k,n) i+1 ◦ f n i + b k i+1 · 1) (43)
where σi+1 is the activation function and bi+1 is the bias for the (i+ 1)th channel.
Since Garriga-Alonso et al. (2018) demonstrate the equivalence between a convolution operation and the matrix multiplication, we can regard the convolutional channel fki+1 as a fully connected layer. Specifically, the input of the virtual fully connected layer is the same as the input of the convolutional channel, the output of the virtual fully connected layer is the vectorized output of the convolutional channel, and the weights of the virtual fully connected layer depend on the convolutional filter and the spatial location of the output of the convolutional channel.
Without considering the activation function, the equivalence between a convolutional channel and a fully connected layer can be explained in Figure 16, in which the first convolutional layer has two channels, i.e., fi = {x; x̂}, where x = {x1, · · · , x9} and x̂ = {x̂1, · · · , x̂9}, and the kth channel of the second convolutional layer is fki+1 = {z1, z2, z3, z4}. In addition, the dimension of the convolutional filter Si+1 is 2×2×2×1, where S(k,1)i+1 = [A,B;C,D] and S (k,2) i+1 = [E,F ;G,H]. Figure 16 shows that the output of fki+1 can be expressed as a fully connected layer. The input of the fully connected layer is the vector x, the weights of each neuron is represented by each row of the matrixW , and the bias is b, thus fki+1 = W · x+ b · 1.
If we consider the activation function and make a precise statement, fki+1 can be formulated as
fki+1 = σi+1(W k i+1 · fi + bki+1 · 1) (44)
whereW ki+1 is the matrix corresponding to all the convolutional filters in the kth channel.
To derive the necessary conditions for two arbitrary different channels are i.i.d., we need to derive the necessary conditions for two arbitrary activations in two arbitrary different channels are i.i.d. For example, if the mth activation of the k channel and the m′th activation of the k′ channel in fi+1 are i.i.d., the necessary conditions are Cov(fki+1(m),f k′ i+1(m ′)) = 0, E(fki+1(m)) = E(f k′ i+1(m ′)).
Since we demonstrate the equivalence between a convolutional channel and a fully connected layer, the necessary conditions for activations being i.i.d. in a fully connected layer are also valid for a convolutional layer. In other words, we can also use the Equation 3, Equation 4, and Equation 5 to examine whether activations being i.i.d. is valid in a convolutional layer under the assumption that all activations being i.i.d. in the previous layer.

heading: F ACTIVATIONS BEING I.I.D. CANNOT BE VALID FOR ALL THE HIDDEN LAYERS OF DNNS
text: In this section, we demonstrate that activations being i.i.d. cannot be valid for all the hidden layers of DNNs through showing that the hidden layers of two typical DNNs, i.e, the MLP and the CNN, cannot satisfy the necessary conditions for activations being i.i.d..

heading: F.1 THE MLP ON THE MNIST DATASET
text: The architecture of the MLP for classifying the MNIST dataset is visualized in Figure 15, and its specific parameters are summarized in Table 3. Since the activation functions of all the hidden layers in the MLP are the sigmoid function, which is strictly increasing, differentiable, and invertible, we only need to examine the first two necessary conditions in the MLP.
Since we restrict all the hidden layers from using the bias, the first two necessary conditions for activations being i.i.d. can be simplified as
∀(k, k′) ∈ S1 = {(k, k′) ∈ Z2|k 6= k′, 1 ≤ k ≤ K, 1 ≤ k′ ≤ K}, we must have Cov(G2k, G2k′) = Cov(Wk,Wk′)E 2(F1) = 0 (45)
E(F1)E(Wk) = E(F1)E(Wk′) (46)
Based on the architecture of the MLP, Table 6 summarizes all the statistical measures for each layer of the MLP. Therefore, the key to examine if the layers of the MLP satisfy the necessary conditions is to estimate the statistical measures.
To estimate the expectations E(X), E(F1), and E(F2), we use the sample means x, f1, and f2, respectively. The definition of the sample mean is
x = 1
M M∑ m=1 xm,f1 = 1 N N∑ n=1 f2(n),f2 = 1 K K∑ k=1 f2(k) (47)
After the MLP is well-trained, we derive the sample means as x = 0.131,f1 = 0.436, and f2 = 0.498 given the benchmark MNIST testing dataset.
Taking into account the estimation error, we can regard E(X) = 0. As a result, we can derive that ∀(n, n′), Cov(G1n, G1n′) = 0 and E(G1n) = E(G1n′) = 0, thus f1 satisfy the necessary conditions for activations being i.i.d..
However, f1 = 0.436 and f2 = 0.498 imply that E(F1) 6= 0 and E(F2) 6= 0 even considering the estimation error. Therefore, we need to estimate Cov(Wk,Wk′) and E(Wk) to examine whether the activations {f2k}Kk=1 being i.i.d.. Since the training dataset is commonly viewed as i.i.d., we can use the same training dataset to train the MLP several times to derive multiple independent observation samples of {βnk(t)}Tt=1 to estimate the random variable Wnk, thereby estimating Cov(Wk,Wk′) and E(Wk).
More specifically, we use the sample correlation r(βk,βk′) to check whether Cov(Wk,Wk′) are close to zero and use the sample mean βk to estimate E(Wk).
r(βk,βk′) = ∑T t=1(βk(t)− βk)(wk′(t)− βk′)√∑T
t=1(βk(t)− βk) 2∑T t=1(βk′(t)− βk′) 2 ≈ Corr(Wk,Wk′) (48)
where βk = 1T ∑T t=1 βk(t) is the sample mean of {βk(t)}Tt=1 and βk(t) = ∑N n=1 βnk(t).
In order to make a tradeoff between the estimation precision and computation complexity, we use the same training dataset to train the MLP 20 times to derive the samples, i.e., {βnk(t)}20t=1. For an activation f2k, it has 128 weights because its previous layer has 128 neurons. Hence, a single sample βk(t) = ∑128 n=1 βnk(t). Based on the 20 samples, i.e., {βk(t)}20t=1, we can calculate r(βk,βk′) to examine whether Cov(Wk,Wk′) are close to zero and derive the sample mean to estimate E(Wk). We use the same method to estimate Cov(An, An′), E(An), Cov(Cl, Cl′), and E(Cl).
Finally, we shows the absolute correlation matrices i.e., AN×N , WK×K , and CL×L, in Figure 3. Since many elements in the three matrices are far from zero, we conclude that their corresponding covariances, i.e, Cov(An, An′), Cov(Wk,Wk′), and Cov(Ck, Ck′), are not close to zero based on connection between correlation and covariance
Corr(Wk,Wk′) = Cov(Wk,Wk′)
σWkσWk′ (49)
where σWk > 0 is the standard derivation of Wk.
Figure 17 visualizes the expectations, i.e., E(An), E(Wk), and E(Cl), for each layer in the MLP. We find that some expectations show great variations, especially {E(An)}Nn=1 and {E(Cl)}Ll=1. In other words, the MLP cannot guarantee that all the expectations of the weights summation are equivalent even though we take into account the estimation error. Therefore, the MLP cannot satisfy the necessary condition for activations being identically distributed.
Overall, we conclude that not all the hidden layers of the MLP can satisfy the necessary conditions for activations being i.i.d. under the assumption that the activations of the previous layer are i.i.d.

heading: F.2 THE CNN ON THE CIFAR-10 DATASET
text: The architecture of the CNN for classifying the CIFAR-10 dataset is in Table 5. In particular, the two convolutional layers do not have activations function, i.e., their activation functions can be viewed as y = x, which is strictly increasing, differentiable, and invertible. The sigmoid activation function is also strictly increasing, differentiable, and invertible, thus we only need to examine the first two necessary conditions in the CNN. To simplify derivation and decrease computation complexity, we assume the max-pooling layers are special activation functions and satisfy the third necessary condition for activations being i.i.d.
Similar as the MLP, we also restrict all the layers of the CNN from using the bias, thus the first two necessary conditions for activations being i.i.d. can be simplified as
∀(k, k′) ∈ S1 = {(k, k′) ∈ Z2|k 6= k′, 1 ≤ k ≤ K, 1 ≤ k′ ≤ K}, we must have Cov(G2k, G2k′) = Cov(Wk,Wk′)E 2(F1) = 0 (50)
E(F1)E(Wk) = E(F1)E(Wk′) (51)
Table 7 summarizes all the statistical measures for every layer of the CNN based on the architecture of the CNN in Table 5. The definitions of most random variables are similar as the MLP except the dimension of the weights. To check if two arbitrary convolutional channels are i.i.d., we need to know all the weights of two convolutional channels. Based on the equivalence between a convolutional channel and a fully connected layer (Appendix E), all the weights of a single convolutional channel consist of all the convolutional filters related to the convolutional channel. For example, the dimension of all the convolutional filters related to a convolutional channel in f1 are 3× 3× 3, thus a single convolutional channel in f1 has 27 weights, so the weights dimension of f1 is 27× 64. After training the CNN well (its training accuracy is 96.2% and testing accuracy is 67.9%), we use the same methods as the MLP to estimate the above statistical measures in Table 7.
For the activation expectations of each layer in the CNN, we have x = −0.023, f1 = −707.074, f2 = −6347.836, f3 = 36.870, f4 = 15.131, thus all the expectations cannot be zero except the input layer, i.e., E(X) = 0, E(F1) 6= 0, E(F2) 6= 0, E(F3) 6= 0, E(F4) 6= 0. For the covariances, we show all the correlation matrices of each layer of the CNN in Figure 18. We find that many correlation coefficients are far from zero in all the correlation matrices, thus their corresponding covariances are not close zero either. Based on the estimations for the activation expectation and the covariance between different weights summations in each layer, we can derive that activations being independent cannot be valid for all the layers of the CNN.
For the weights summation expectation, we visualize all the expectations for each layer in the CNN in Figure 19. We find that some expectations have great variations, especially {E(Cl)}Ll=1 and {E(Dq)}Qq=1, which means that most expectations of weights summations of different neurons in f3 and f4 cannot be equivalent. As a result, the activations of the two fully connected layers in the CNN cannot be identically distributed.
Overall, we conclude that not all the hidden layers of the CNN can satisfy the necessary conditions for activations being i.i.d. under the assumption that the activations of the previous layer are i.i.d.

heading: G THE EQUIVALENCE BETWEEN THE GRADIENT DESCENT ALGORITHM AND THE FIRST ORDER APPROXIMATION
text: In the context of deep learning, most learning algorithms belong to the gradient descent algorithm. Given a DNN = {x;f1; ...;fI ;fY } and a cost function H[fY , P (Y |X)], where P (Y |X) is the true distribution of the training label given the corresponding training dataset. Let θ be the parameters of the DNN, the gradient descent aims to optimize θ by minimizing H[fY , P (Y |X)] (Rumelhart et al., 1986). We typically update θ iteratively to derive θ∗, which can be expressed as
θt+1 = θt − α∇θtH[fY , P (Y |X)] (52)
where∇θtH[fY , P (Y |X)] is the Jacobian matrix of H[fY , P (Y |X)] with respect to θt at the tth iteration, and α > 0 denotes a constant indicating the learning rate. Since P (Y |X) is constant, we denote H(fY ) as H[fY , P (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 ofH(fY ) with respect to the parameters of the ith hidden layer, i.e., ∇θ(i)H(fY ), can be expressed as follows based on the chain rule.
∇θ(i)H(fY ) = ∇fY H(fY ) · ∇fIfY · I∏
j=i+1
∇fj−1fj · ∇θ(i)fi (53)
where θ(i) denote the parameters of the ith hidden layer.
Based on Equation 52 and Equation 53, θ(i) can be learned by the gradient descent method as
θt+1(i) = θt(i)− α[∇θt(i)H(fY )] (54)
Table 8 summarizes the backpropagation training procedure for the MLP in Figure 15.
If an arbitrary function f is differentiable at point p∗ in RN and its differential is represented by the Jacobian matrix∇pf , the first order approximation of f near the point p can be formulated as f(p)− f(p∗) = (∇p∗f) · (p− p∗) + o(||p− p∗||) (55) where o(||p− p∗||) is a quantity that approaches zero much faster than ||p− p∗|| approaches zero. Based on the equivalence between the gradient descent method and the first order approximation (Battiti, 1992), updating the activations of the hidden layers of the MLP in Figure 15 during the backpropagation training procedure can be approximated as
f2[f1,θt+1(2)] ≈ f2[f1,θt(2)] + (∇θt(2)f2) · [θt+1(2)− θt(2)] f1[x,θt+1(1)] ≈ f1[x,θt(1)] + (∇θt(1)f1) · [θt+1(1)− θt(1)]
(56)
where f2[f1,θt(2)] denote the activations of the second hidden layer based on the parameters learned in the tth iteration, i.e., θt(2), given the activations of the first hidden layer, i.e., f1. The definitions of f2[f1,θt+1(2)] and f1(x,θt(1)) are the same as f2[f1,θt(2)].
Because f2 has K neurons, i.e., f2 = {f2k = σ2( ∑N n=1 βnk · f1n + b2k)}Kk=1, and each neuron has N + 1 parameters, i.e., θ(2) = {[β1k; · · · ;βNk; b2k]}Kk=1, the dimension of ∇θt(2)f2 is equal to K × (N + 1) and ∇θt(2)f2 can be expressed as
∇θt(2)f2 = (∇σ2f2) · [f1; 1] T (57)
where∇σ2f2 = ∂f2[f1,θt(2)] ∂σ2 .
Substituting (∇σ2f2) · [f1; 1]T for ∇θt(2)f2 in Equation 56, we can derive f2[f1,θt+1(2)] ≈ f2[f1,θt(2)] + (∇σ2f2) · [f1; 1]T · θt+1(2)− (∇σ2f2) · [f1; 1]T · θt(2)
(58)
If we only consider a single neuron, such as f2k, in f2, θt+1(2k) = [β1k; · · · ;βNk; b2k] and θt(2k) = [β ∗ 1k; · · · ;β∗Nk; b∗2k], thus [f1; 1]T · θt+1(2k) = ∑N n=1 βnk · f1n + b2k. As a result, Equation 58 with respect to a single neuron can be expressed as
f2k[f1,θt+1(2k)] ≈ (∇σ2f2k) · [ N∑ n=1 βnk · f1n + b2k] (First order approximation)
+ f2k[f1,θt(2k)]− (∇σ2f2k) · [ N∑ n=1 β∗nk · f1n + b∗2k] (Error)
(59)
Equation 59 indicates that f2k[f1,θt+1(2k)] can be expressed as its first order approximation with an error component based on the activation in the previous iteration, i.e., f2k[f1,θt(2k)]. Since ∇σ2f2k = ∂f2k[f1,θt(2)] ∂σ2
is only related to f1 and the parameters in the tth training iteration, i.e, θt(2), it can be regarded as a constant. Also note the error component not contains any parameters in the (t+ 1)th training iteration. In summary, f2k(f1,θt+1(2k)) can be reformulated as
f2k(f1,θt+1(2k)) ≈ C1 · [ N∑ n=1 βnk · f1n + b2k] + C2 (60)
where C1 = ∇σ2f2k and C2 = f2k(f1,θt(2k))− (∇σ2f2k) · [ ∑N n=1 β ∗ nk ·f1n+b∗2k]. Similarly, the activations in the first hidden layer, i.e, f1, also can be formulated as the first order approximation in the context of the gradient descent learning algorithm.

heading: H THE PROPOSED PROBABILISTIC REPRESENTATION HOLDS FOR CONVOLUTIONAL LAYERS
text: This section proves that the proposed probabilistic representation is valid for convolutional layers. In the above CNN, the input x has Q channels x = {xq}Qq=1 (e.g., Q = 3 if x are color images), f1 has N convolutional channels f1 = {fn1 }Nn=1, f2 has K convolutional channels f2 = {fk2 }Kk=1. The output layer fY is the softmax with L nodes, thus its distribution can be formulated as
P (FY ) = {P (fyl) = 1
ZFY exp(fyl)}Ll=1 (61)
where ZFY = ∑L l=1 exp(fyl) denotes the partition function. Since {fyl = ∑J j=1 γjl ·f2j +byl}Ll=1, we can derive
P (FY ) = {P (fyl) = 1
ZFY exp( J∑ j=1 γjl · f2j + byl)}Ll=1 (62)
where {f2j}Jj=1 is the flattened output of f2, γjl is the weight of the edge between f2j and fyl, and byl denotes the bias. Since f2 is a convolutional layer with K convolutional channels, it can be reformulated as {f2j}Jj=1 = {fk2 }Kk=1. As a result, P (FY ) can be reformulated as
P (FY ) = {P (fyl) = 1
ZFY exp( K∑ k=1 γTkl · fk2 + byl)}Ll=1 (63)
where γkl is a subset of all the parameters {γjl}Jj=1 such that ∑K k=1 γ T kl · fk2 = ∑J j=1 γjl · f2j , and γTkl is the transpose of γkl. Therefore, P (FY ) can be reformulated as
P (FY ) = {P (fyl) = 1
Z ′FY K∏ k=1 exp(γTkl · fk2 )}Ll=1 (64)
where Z ′FY = ZFY /[exp(byl)]. Recall the element-wise matrix power, e.g., exp(a) b, where a = [a1; a2; a3] and b = [b1; b2; b3], we can derive that exp(a) = [exp(a1); exp(a2); exp(a3)] and exp(a)b = [exp(a1)b1 ; exp(a2)b2 ; exp(a3)b3 ] = [exp(a1b1); exp(a2b2); exp(a3b3)]. As a result, exp(ab) = exp(a1b1 + a2b2 + a3b3) = ∏ |a| exp(a) b, where |a| is the element numbers of a.
Based on the element-wise matrix power, we can reformulate P (FY ) as
P (FY ) = {P (fyl) ≈ 1
Z ′FY K∏ k=1 exp(γTkl · fk2 ) = 1 Z ′FY K∏ k=1 ∏ |fk2 | [exp(fk2 )] γkl}Ll=1 (65)
Moreover, we introduce a new partition function ZFk2 = ∑ fk2
exp(fk2 ) to guarantee that 1
Z Fk2
exp(fk2 ) is a probability measure. As a result, P (FY ) can be reformulated as
P (FY ) = {P (fyl) = 1
Z ′′FY K∏ k=1 ∏ |fk2 | [ 1 ZFk2 exp(fk2 )] γkl}Ll=1 (66)
where Z ′′FY = Z ′ FY / ∏K k=1[ZFk2 ]
γkl·|fk2 |. Overall, P (FY ) can be reformulated as a PoE model, in which each expert is defined by every convolutional channel in f2
P (fk2 ) = 1
ZFk2 exp(fk2 ) (67)
Since f2 has K convolutional channels, i.e., f2 = {fk2 }Kk=1, and each channel can be formulated as the summation of all the convolutional channels in f1, i.e.,
f2 = {fk2 = σ2( N∑ n=1 S (k,n) 2 ◦ fn1 + bk2 · 1)}Kk=1 (68)
where S(k,n)2 ◦ fn1 is the output of the kth convolutional filter applying into the nth channel of f1, bk2 is the bias, and σ2(·) is the activation function. Therefore, we have
P (fk2 ) = 1
ZFk2 exp[σ2( N∑ n=1 S (k,n) 2 ◦ fn1 + bk2 · 1)] (69)
Based on the equivalence between the gradient descent learning and the first order approximation (Appendix G), Equation 69 can be approximated as
P (fk2 ) ≈ 1
ZFk2 exp( N∑ n=1 S (k,n) 2 ◦ fn1 + bk2 · 1) (70)
Let us define a linear function as ϕk2(f1) = ∑N n=1 S (k,n) 2 ◦ fn1 + bk2 · 1, thus the distribution of the kth convolutional channel in f2, i.e., P (fk2 ), can be expressed as
P (fk2 ) ≈ 1
ZFk2 exp[ϕk2(f1)] (71)
If we regard all the linear filters {ϕk2(f1)}Kk=1 as potential functions for modeling signal structures of f1, we can formulate P (F2) as a specific Gibbs distribution, i.e., the Markov Random Fields (MRFs), which can be expressed as
P (F2) = 1
ZF2 K∏ k=1 P (fk2 ) = 1 ZF2 exp[ K∑ k=1 ϕk2(f1)] (72)
where ZF2 = ∑ f1 exp[ ∑K k=1 ϕ k 2(f1)] is the partition function for P (F2). In summary, we prove that the convolutional layer f2 formulates a MRF model to describe the statistical property of f1.
Here we prove that the convolutional layer f1 also formulates a MRF model to describe the statistical property of the input x. Based on Equation 70, we have
P (fk2 ) ≈ 1
Z ′ Fk2 exp( N∑ n=1 S (k,n) 2 ◦ fn1 ) (73)
where Z ′ Fk2 = ZFk2 /exp(b k 2).
Based on the equivalence between a convolutional channel and a fully connected layer proven in Appendix E, we can regard a convolutional operation as a matrix multiplication, which is visualized in Figure 21. As a result, we have
P (fk2 ) ≈ 1
Z ′ Fk2 exp( N∑ n=1 S (k,n) 2 ◦ fn1 ) = 1 Z ′ Fk2 exp( N∑ n=1 W (k,n) 2 · fn1 ) (74)
whereW (k,n)2 is a matrix corresponding to the convolution filter S (k,n) 2 .
Based on the element-wise matrix power, P (fk2 ) can be reformulated as
P (fk2 ) ≈ 1
Z ′ Fk2 N∏ n=1 exp(W (k,n)2 · fn1 ) = 1 Z ′ Fk2 N∏ n=1 exp(fn1 ) W (k,n) 2 (75)
Similarly, we introduce a new partition function ZFn1 = ∑ fn1
exp(fn1 ) to guarantee that 1
ZFn1 exp(fn1 ) is a probability measure. As a result, P (f k 2 ) can be reformulated as
P (fk2 ) = P (fyl) ≈ 1
Z ′′ Fk2 N∏ n=1 [ 1 ZFn1 exp(fn1 )] W (k,n) 2 (76)
where Z ′′ Fk2 = Z ′ Fk2 / ∏N n=1[ZFn1 ] W (k,n) 2 . Therefore, P (fk2 ) can be reformulated as a PoE model, in which all the experts are defined as P (fn1 ) = 1
ZFn1 exp(fn1 )
Since f1 has N convolutional channels
f1 = {fn1 = σ1( Q∑ q=1 S (n,q) 1 ◦ xq + bn1 · 1)}Nn=1 (77)
Based on the same derivation as Equations 70, 71, 72, P (F1) also can be expressed as a MRF model
P (F1) = 1
ZF1 N∏ n=1 P (fn1 ) = 1 ZF1 exp[ N∑ n=1 ϕn1 (x)] (78)
where ϕn1 (x) = ∑Q q=1 S (n,q) 1 ◦xq+bn1 ·1) and the partition functionZF1 = ∑ x exp[ ∑N n=1 ϕ n 1 (x)]. Overall, the distribution of a convolutional layer can be formulated as a MRF model.
Algorithm 1 The algorithm for generating the synthetic dataset input: NIST dataset of handwritten digits by class
1: repeat 2: binarizing an image of NIST to obtain z 3: extracting the central part of z to obtain zc with dimension 64× 64 4: downsampling zc to obtain zcd with dimension 32× 32 5: extracting the edge of zcd to obtain the mask image mcd 6: decomposing mcd into four parts, i.e., moutside, moutside-boundary, minside-boundary, and minside. 7: samplingN (0, 1024) to derive a random vector x ∈ R1024×1 8: sorting x in the descending order to derive x̂ 9: decomposing x̂ into four parts, i.e., x̂ = {x̂outside, x̂outside-boundary, x̂inside-boundary, x̂inside}
10: placing each pixel of x̂ = {x̂outside, x̂outside-boundary, x̂inside-boundary, x̂inside} into the above masks. 11: until (20,000 synthetic images are generated) output: The synthetic dataset

heading: I THE APPROACH TO GENERATE THE SYNTHETIC DATASET BASED ON NIST AND THE ARCHITECTURE OF THE CNN
text: In this section, we present a novel approach to generate a synthetic dataset obeying the Gaussian distribution based on the NIST 1 dataset of handwritten digits by class. The synthetic dataset has similar characteristics as the benchmark MNIST dataset. It consists of 20,000 32 × 32 grayscale images in 10 classes (digits from 0 to 9), and each class has 1,000 training and 1,000 testing images. More specifically, the approach to generate the synthetic dataset is summarized in Algorithm 1, and Figure 22 visualize the spatial relation of x̂ = {x̂outside, x̂outside-boundary, x̂inside-boundary, x̂inside}.
1https://www.nist.gov/srd/nist-special-database-19

heading: X x Input 32× 32× 1
text: The architecture of the CNN for classifying the synthetic dataset is summarized in Table 9.

heading: J THE EQUIVALENCE BETWEEN SGD AND VARIATIONAL INFERENCE

heading: J.1 THE DISTRIBUTION OF THE OUTPUT LAYER GIVEN THE INPUT IS A GIBBS DISTRIBUTION
text: In Section 6.2, we demonstrate the entire architecture of a DNN = {x;f1; ...;fI ;fY } corresponds to a joint distribution P (FY ;FI ..;F1|X). As a result, the conditional distribution of fY given x, i.e., P (FY |X), can be formulated as
P (FY |X) = ∫ · · · ∫ P (FY ;FI ..;F1|X)dFI · · · dF1 (79)
Since P (FY ;FI ..;F1|X) = P (FY |FI) · ... · P (Fi+1|Fi) · ... · P (F1|X) and only F2 depends on F1, we can derive
P (FY |X) = ∫ · · · ∫ P (FY |FI) · ... ·P (F3|F2)dFI · · · dF2 · ∫ P (F2|F1)P (F1|X)dF1 (80)
Since ∫ P (F2|F1)P (F1|X)dF1 = ∫ P (F2;F1|X)dF1 = P [F2|f1(X)], where f1(·) is the deterministic function defined by the corresponding hidden layer, P (FY |X) can be simplified as
P (FY |X) = ∫ · · · ∫ P (FY |FI) · ... · P (F3|F2)P [F2|f1(X)]dFI · · · dF2 (81)
Similarly, since only F3 depends on F2, we can obtain P (FY |X) = ∫ · · · ∫ P (FY |FI) · ... · P (F3|f2(f1(X)))dFI · · · dF3 (82)
If we iteratively apply the integral with respect to F3 until FI , we can obtain
P (FY |X) = P [FY |fI(· · ·f2(f1(X))· · ·)] (83)
Since the output layer fY is commonly a fully connected layer. we can formulate P (FY |X) as a multivariate discrete Gibbs distribution based on Equation 8.
P (FY |X) = 1
ZDNN(x) exp[fY (fI(· · ·f2(f1(x))· · ·))] (84)
It is important to note that the energy function EDNN(fY |x) = −fY (fI(· · ·f2(f1(x))· · ·)), i.e., EDNN(·) is determined by the entire architecture of DNNs. In addition, since P (FY |X) is a conditional distribution of FY givenX , the partition function ZDNN(x) should integrate all possible activations of the output layer fY to guarantee P (FY |X) being a valid conditional probability, i.e., ZDNN(x) = ∫ exp[fY (fI(· · ·f2(f1(x))· · ·))]dfY (fI(· · ·f2(f1(x))· · ·)).
For clarity, here we only consider integral no matter a Gibbs distribution is discrete or continuous, because ∑ f1···fI exp[−EDNN(fY |x)] is a special case of integral. All functions involving fi should include parameters θfi in the conditioning set, but we omit for the same reason.

heading: J.2 SGD CAN BE EXPLAINED AS A VARIATIONAL INFERENCE ALGORITHM
text: As a dominant paradigm for Bayesian posterior inferenceP (H|E) ∝ P (E|H)P (H), the variational inference converts the inference problem into an optimization problem, where the prior distribution P (H) is the probability of arbitrary hypothesis H with respect to the observation E, and P (E|H) is the likelihood distribution. Specifically, the variational inference posits a family of approximate distributions Q and aims to find a distribution Q∗(H) that minimizes the Kullback-Leibler (KL) divergence between P (H|E) and Q(H) (Blei et al., 2017).
Q∗(H) = arg min Q∈Q KL(P (H|E)||Q(H)) (85)
Notably, recent works demonstrate that SGD performs variational inference during training DNNs (Mandt et al., 2017; Chaudhari & Soatto, 2018). Following these works, we further confirm the Bayesian explanation for SGD based on the proposed probabilistic for DNNs.
Given a DNN = {x;f1; ...;fI ;fY }, we choose the cross entropy H[P (Y |X), P (FY |X)] as the loss function, where P (Y |X) is the true posterior distribution determined by the training datasetD and P (FY |X) is the conditional distribution of fY given x derived from the DNN.
H[P (Y |X), P (FY |X)] = H[P (Y |X)] + KL[P (Y |X) ‖ P (FY |X)]
= − L∑ l=1 P [y(l)|x]log[P (fY (l)|x)] (86)
where H[P (Y |X)] is the entropy of P (Y |X) and P (y(l)|x) is the true posterior distribution given the training data x. For example, if L = 4 and P (y|x) is [0, 1, 0, 0], so P [y(2)|x] = 1 and P [y(3)|x] = 0. P (fY |x) is the learned conditional distribution of fY given the same x, e.g., P (fY |x) could be [0.2, 0.7, 0.1, 0]. In the context of deep learning, SGD aims to find the P ∗(FY |X) such that it can minimize the loss function H[P (Y |X), P (FY |X)] given the training datasetD.
P ∗(FY |X) = arg min P (FY |X)∈Q H[P (Y |X), P (FY |X)] (87)
Given the training datasetD, P (Y |X) is a constant, i.e., H[P (Y |X)] = 0, thus we can obtain P ∗(FY |X) = arg min
P (FY |X)∈Q KL[P (Y |X) ‖ P (FY |X)] (88)
Comparing the Equation 88 and 85, we can confirm the equivalence between SGD and variational inference based on the proposed probabilistic explanations for DNNs. More specifically, the DNN defines a family of conditional Gibbs distribution P (FY |X) (Equation 84) and SGD aims to find P ∗(FY |X) such that it can minimize the distance to the true posterior distribution P (Y |X).

heading: J.3 SGD IS EQUIVALENT TO THE ENERGY MINIMIZATION OPTIMIZATION
text: In most deep learning applications, P (y|x) is commonly an one-hot vector, i.e., P (y|x) only has a single non-zero element. If we assume that l∗ is the index of the non-zero element of P (y|x), i.e., P [y(l∗)|x] = 1, we can simplify KL[P (Y |X) ‖ P (FY |X)] = −log[P (fY (l∗)|x)]. Based on the definition of P (FY |X) (Equation 84), we can obtain
KL[P (Y |X) ‖ P (FY |X)] = −log[P (fY (l∗)|x)] = −fyl∗(fI(· · ·f2(f1(x)) · · · )) + log[ZDNN(x)] = EDNN (fyl∗ |x) + log[ZDNN(x)]
(89)
where fyl∗ is the l∗th output node of the output layer fY . In particular, since log[ZDNN(x)] is not related to the output fY , it can be viewed as a constant with respect to the gradient calculation indicated by Table 8. Finally, combining Equation 88 and 89, we can obtain
P ∗(FY |X) = arg min P (FY |X)∈Q EDNN (fyl∗ |x) (90)
The above equation indicates that SGD is equivalent to the energy minimization optimization, which in turn validates the Gibbs explanation for the hidden layer of DNNs.

heading: K EXPERIMENTS VALIDATING THE THEORETICAL EXPLANATION FOR SGD

heading: K.1 EXPERIMENT: SGD IS EQUIVALENT TO THE ENERGY MINIMIZATION OPTIMIZATION
text: To validate the proof that SGD learns a P (FY |X) with the minimal energy from the family of Gibbs distribution determined by a DNN, we use SGD to train the CNN defined in Table 9 for classifying the synthetic dataset shown in Figure 5. Specifically, we use the entire training dataset to train the CNN over 30 epochs. During each training epoch, we calculate the energy function of P (FY |X) given a testing dataset with the same label. Since the testing dataset has the same label, we have the same index of the single non-zero element of P [y|x]. As a result, we can easily examine if the corresponding energy function ECNN (fyl∗ |x) is decreasing during the training procedure. Based on the architecture of CNN in Table 9, ECNN (fyl∗ |x) can be formulated as
ECNN (fyl∗ |x) = −fyl∗(f4(f3(f2(f1(x))))) (91) In Figure 23, the first row visualizes the variation ofECNN (fyl∗ |x) over 30 training epochs. Though there are certain fluctuations,ECNN (fyl∗ |x) indeed shows a decreasing trend. Due to limited space, we only visualize the energy minimization procedure given three labels, but the procedure is valid for any label. In addition, we demonstrate the equivalence based on the MLP (Table 3) for classifying the MNIST dataset. The energy function of P (FY |X) corresponding to the MLP is defined as
EMLP (fyl∗ |x) = −fyl∗(f2(f1(x))) (92) In Figure 23, the second row shows the same energy decreasing trend during training the MLP. Overall, we demonstrate the equivalence between SGD and the energy minimization optimization, which in turn validates the Gibbs explanation for DNNs.

heading: K.2 EXPERIMENT: THE PARAMETERS LEARNED BY SGD ARE NOT I.I.D.
text: To validate the parameters learned by SGD not i.i.d., we use SGD to train the MLP defined in Table 3 and the CNN defined in Table 5. Following the commonly used initialization approach, we initialize all the weights of the two DNNs by Gaussian random numbers, i.e., all the parameters can be viewed as i.i.d. before training. We use the absolute value of the sample correlation defined by Equation 2 to quantify the dependence between different weights during the training procedure.
Figure 24 visualizes the the variation of the average correlation between different weights of the three layers in the MLP. Initially, the average correlation for f1 and f2 are very small, because parameters are Gaussian random values. As the training procedure continues, we find that all the correlation values are increasing. Notably, the increase for f2 is huge, i.e., the weights of f2 cannot be i.i.d. anymore. However, the increase for f1 is very small, i.e., all the weights of f1 still could be i.i.d.. Therefore, not all the learned parameters in the MLP are i.i.d..
Figure 25 visualizes the variation of the average correlation between different weights of the five layers in the CNN. Similar as the weights in the MLP, the average correlation for f1 and f2 initially are very small, because they are initialized as Gaussian random values. As the training procedure continues, they have a sharp increasing and quickly achieve stationary, thus the parameters of f1 and f2 are not i.i.d. anymore. In addition, the average correlation for f3 and f4 also have an increasing trend, but their increase are very small, i.e., the parameters of f3 and f4 still could be i.i.d.. Overall, we conclude that not all the learned parameters in the CNN are i.i.d..

heading: K.3 EXPERIMENT: A BETTER PRIOR DISTRIBUTION INDICATES A BETTER GENERALIZATION PERFORMANCE
text: We design three CNNs = {x;f1;f2;f3;f4;f5;fY } in Table 10 to classify the synthetic dataset shown in Figure 5. We use the entire training dataset to train the CNNs over 30 epochs and use the entire testing dataset to examine the generalization performance of the CNNs.
Based on the MRF explanation for convolutional layers, the first convolutional layer f1 formulates a MRF model as the prior distribution P (F1), which can be expressed as
P (F1) = 1
ZF1
∏K k=1 exp[ϕ k 1(x)] =
1 ZF1
exp[ ∑K k=1 ϕ k 1(x)] (93)
where ϕk1(x) = S k ◦x+bk1 ·1 is the kth linear convolutional channel, Sk is the linear convolutional filter and bk1 is the bias.

heading: X x Input 32× 32× 1 32× 32× 1 32× 32× 1
text: Since the synthetic dataset obeys the Gaussian distribution P (X), we can measure the quality of the prior distribution by averaging KLD[P (X)||P (F1)] over all the testing images for every training epoch. In addition, we use the testing error to measure the generalization performance. Also note that the only difference between the tree CNNs is they have different K numbers of convolutional channels in f1, i.e., different CNNs have different prior distributions. Therefore, we can examine the effect of the prior distribution by checking the generalization performance of the three CNNs.
Figure 26 shows the learned prior distributions and the corresponding testing error over 30 training epochs. We can find that there is a positive correlation between KLD[P (X)||P (F1)] and the testing error. Specifically, the testing error decreases as the learned prior distribution P (F1) approaching to the truly prior distribution P (X), and the testing error becomes stable when P (F1) achieves stationary for all the three CNNs.
Moreover, we find that the CNN with more convolutional channels achieves the better generalization performance. In terms of Bayesian theory, that is because a convolutional layer with more convolutional channels can formulate a more powerful prior distribution, such that it can convey more prior knowledge of the dataset, thereby rendering better generalization performance. Overall, we justify that a better prior distribution indicates a better generalization performance.

heading: L A REGULARIZATION ALGORITHM BASED ON THE PROBABILISTIC REPRESENTATION

heading: L.1 THE PROPOSED REGULARIZATION ALGORITHM
text: Based on the probabilistic representation, the first hidden layer formulates a MRF model as a prior distribution of the training dataset and other hidden layers recursively serve as prior distributions for their respective next layer in a CNN. However, the backpropagation calculates the gradient in the backward direction, i.e., updating the gradient of a hidden layer corresponding to a prior distribution must go over all hidden layers behind it. In addition, the degradation problem could makes the learned prior distribution not precise enough to regularize DNNs and results in overfitting, especially in very deep neural networks (He et al., 2016).
In this section, we propose to a new regularization learning algorithm for DNNs, which includes two steps: (i) we learn the prior distributions of DNNs directly from the training dataset {xj}Jj=1 as the blue arrow shown in Figure 27 and (ii) we use the backpropagation to train the initialized DNNs by the learned prior distributions as the red arrow shown in Figure 27.
Since directly learning prior distribution from {xj}Jj=1 is less complicated than the backpropagation learning algorithm, we can relieve the effects of degradation problem and derive prior distributions more precisely, thereby achieving the better generalization performance of DNNs.
An important question need to answer is how many hidden layers of DNNs need to be pre-trained in the first step. Since most prior information is stored in the training dataset {xj}Jj=1, the hidden layer closest to the input of DNNs, i.e., the first hidden layer f1, should learn the most prior information from {xj}Jj=1 as long as f1 has powerful representation ability through including sufficient hidden units, e.g., convolutional filters. Therefore, the proposed regularization method only pre-train the first hidden layer f1 in order to decrease the computation complexity. Pre-training more hidden layers could obtain better generalization performance, which can be a direction for future research.
We choose the Field of Experts (FoE) model to directly learn a prior distribution of {xj}Jj=1 because of two reasons. Above all, FoE achieves great performance on learning prior knowledge of given datasets on various applications, e.g., image denoising (Roth & Black, 2005), image restoration Roth & Black (2008), and image inpainting (Schmidt et al., 2011). More importantly, FoE can be regarded as a Convolutional Neural Networks (CNN) with a single convolutional layer, hence we can easily use the learned FoE model to initialize the first hidden layer of CNNs as long as they have the same architecture. A FoE model can be formulated as
pFoE(x) = 1
ZFoE K∏ k=1 {f expertk [fk(x)]}, (94)
where fk(·) is a linear filter, f expertk is an expert function, and ZFoE = ∫ x ∏K k=1{f expert k [fk(x)]}dx is the partition function. We use the contrastive divergence learning method to train pFoE(x) under the Kullback-Leibler divergence (KLD) criterion (Carreira-Perpinan & Hinton, 2005). More detail of the training procedure is in (Roth & Black, 2005; 2008). After deriving a well-trained FoE model, we use the learned linear filters {fk}Kk=1 to initialize the first convolutional layer of CNNs. Finally, we use the backpropagation algorithm to train the initialized CNN for a specific application, such as image recognition. We summarize the proposed regularization learning algorithm in Algorithm 2.
Algorithm 2 Regularization learning algorithm input: D = {(xn,yn)}Nn=1, DNN = {x;f1; ...;fI ;fY }, and pFoE(x)
1: setup 2: initialize iteration k = 0, and training epochs K 3: specify the architecture of the DNN and pFoE(x) 4: regularization learning 5: train pFoE(x) given {xn}Nn=1 6: initialize f1 by the well-trained pFoE(x) 7: repeat 8: train the DNN given D 9: until (k = K)
output: the well-trained DNN for image recognition

heading: L.2 SIMULATIONS
text: In this section, we validate the proposed regularization learning algorithm for image recognition application on the CIFAR-10 benchmark dataset, which includes 70,000 32 × 32 color images in 10 classes, and each class has 6,000 training images and 1,000 testing images. Specifically, we randomly choose 2,000 training images of each class as a new training dataset and use the same testing dataset. Smaller training dataset means overfitting being more likely to occur, which can help us validate the proposed regularization learning method more convincingly. In addition, we covert color images to grayscale to simplify computation.
In order to directly learn regularization from {xj}20000j=1 , we design a FoE model with 12 Gaussian Scale Mixture (GSM) experts, thus it can be expressed as
pFoE(x) = 1 ZFoE ∏12 k=1{fGSMk [fk(x)]}, (95)
where fGSMk [fk(x)] = ∑11 i=1 πki · N (fk(x); 0, σ2b δ(i) ), σ 2 b is a fixed base variance, δ = {δ(1), · · · , δ(11)} is a range of scales, fk(·) is a 3 × 3 convolutional filter, and πki denotes the weight of each Gaussian distribution. The parameters of the FoE model is θ = {f,π}, where f = [f1, ..., f12] and π = {πki}k=12,i=11k=1,i=1 . The architecture of the CNN for image recognition is summarized in Table 11. It is noteworthy that the number of convolutional filters in f1 and the filter dimension are the same as the above FoE model, thus we can use the well-trained FoE to initialize f1 as a regularization for the CNN classifying the CIFAR-10 dataset based on Algorithm 2.
In order to make a comprehensive comparison, we choose two commonly used regularizations, i.e., the dropout (drop rate is 0.3) and the L2 norm, as references. We use the testing error to exam the generalization performances of different regularizations. Figure 28 visualizes that the proposed regularization method (abbr. Bayes Prior (BP)) successfully decreases the testing error. In particular, it outperforms dropout (32.26%) and achieve similar performance as L2 norm (29.58%). Therefore, this experiment validates the effectiveness of the proposed regularization learning method.
Moreover, the proposed regularization method can be easily combined with other regularizations methods to improve the generalization performance further. We combine two arbitrary aforementioned regularization methods together and examine their generalization performances under the same experimental conditions as above. Figure 29 shows that the combination of two regularizations together achieves better generalization performance than any single one. In particular, the combination of BP and L2 norm achieves the lowest testing error (28.60%). Also note that if we combine three regularizations together, we achieve the best generalization performance (27.69%). This experiment further validate the proposed regularization learning algorithm.

Review summaries:
The authors show that parameters of a DNN do not satisfy the i.i.d. prior assumption and that neural layer activations considered as i.i.d. are not valid assumptions for all hidden layers of the network. One can therefore not rightfully use GPs to describe the network’s hidden layers. The authors suggest formulating the neurons per layer as energy functions thereby rendering a hidden layer as a Gibbs distribution and the connection between adjacent hidden layers as a PoE model.