Title: PROBABILISTIC MODELING THE HIDDEN LAYERS OF DEEP NEURAL NETWORKS Anonymous

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.

1. INTRODUCTION:

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.

2. PRELIMINARIES:

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 between X and Y , and θ indicate the parameters of P θ (X, Y ).

A dataset D = {(x j , y j )|x j ∈ R M , y j ∈ R L } J j=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; f 1 ; ...; f I ; f Y } and trained by D, where (x, y) ∈ D are the input of the DNN and the corresponding training label, and f Y is an estimation of the true distribution P (Y |X). As a result, x ∼ P (X) and f Y ≈ P (Y |X). The meaning of f i 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, F i and F Y are the random variables corresponding to f i and f Y , respectively. In this work, we use the MLP (Figure 1 ) to derive most theoretical conclusions unless otherwise specified. 

3. BACKGROUND:

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

EQUATION

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 f k (x; θ k ) (Wang et al., 2013; Goodfellow et al., 2016) . Moreover, assuming a single potential function defines an expert, i.e., F k = exp[-f k (x; θ k )], we can regard MRFs as a PoE model, i.e., P (x; θ) = 1 (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.

Z(θ) k F k , where Z(θ) = k Z(θ k )

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.

4. MAIN RESULTS:



4.1. THE PARAMETERS OF DNNS DO NOT SATISFY THE I.I.D. ASSUMPTION:

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 } M m=1 ∼ P (α n ) in the MLP (Figure 1 ). Since all the parameters are i.i.d., different neurons should be uncorrelated, i.e., ∀n = 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 = 1 M M m=1 α mn is the sample mean of {α mn } M m=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. Third, we visualize all the parameters of the NN = {x; f ; f Y } 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..

4.2. ACTIVATIONS BEING I.I.D. IS NOT VALID FOR ALL THE HIDDEN LAYERS OF DNNS:

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 (A mn ), β nk ∼ P (W nk ), b 1n ∼ P (B 1n ), and b 2k ∼ P (B 2k ), thus we regard

G 1n = M m=1 A mn X m + B 1n as the random variable of g 1n (x) = M m=1 α mn • x m + b 1n , and G 2k = N n=1 W nk F 1n + B 2k as the random variable of g 2k (f 1 ) = N n=1 β nk • f 1n + b 2k .

Hence, the random variables of two arbitrary activations in f 1 and f 2 , e.g., f 1n and f 2k , can be expressed as

F 1n = σ 1 (G 1n ) and F 2k = σ 2 (G 2k ), respectively.

Given two random variables A and B, 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 {F 1n } N n=1 are i.i.d. and {f 1n } N n=1 ∼ P (F 1 ), we derive the following necessary conditions for the activations of the second hidden layer, i.e. ,{F 2k } K k=1 , being i.i.d..

EQUATION

σ 2 (•) is strictly increasing, differentiable, and invertible.

(

EQUATION

where W k = N n=1 W nk . Equation 3 and 4 are the necessary conditions for {G 2k } K k=1 being i.i.d.. Equation 5 constrains the activation function σ 2 (•) to guarantee that if {G 2k } K k=1 are i.i.d. then {F 2k } K k=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. 

E(X) Cov(An, A n ) E(An) f2 E(F1) Cov(W k , W k ) E(W k ) f Y E(F2) Cov(C l , C l ) E(C l ) {xm} M m=1 ∼ P (X), {f1n} N n=1 ∼ P (F1), {f 2k } K k=1 ∼ P (F2) An = M

m=1 Amn, where αmn ∼ P (Amn) Based on the empirical approach, we estimate the statistical measures (Table 1 ) to examine whether the activations of the three layers (i.e., f 1 , f 2 , and f Y ) 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 f 2 can be simplified as

W k = N n=1 W nk , where β nk ∼ P (W nk ) C l = K k=1 C kl ,

∀(k, k ) ∈ S 1 = {(k, k ) ∈ Z 2 |k = k , 1 ≤ k ≤ K, 1 ≤ k ≤ K}, we must have Cov(G 2k , G 2k ) = Cov(W k , W k )E 2 (F 1 ) = 0 (6) E(F 1 )E(W k ) = E(F 1 )E(W k ) (7)

First, Figure 3 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, f 1 = 0.436 and f 2 = 0.498 imply E(F 1 ) = 0 and E(F 2 ) = 0 even considering the estimation error. Therefore, the activations of the last two hidden layer, i.e., f 2 and f Y , 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.

A N × N f 1 W K × K f 2 C L × L

4.3. THE DISTRIBUTION OF A HIDDEN LAYER IS A GIBBS DISTRIBUTION:

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

EQUATION

)

where

Z Y = L l=1 exp(f yl ) is the partition function. Since f yl = K k=1 γ kl • f 2k + b yl , we have P (F Y ) = {P (f yl ) = 1 Z Y exp( K k=1 γ kl • f 2k + b yl )} L l=1 (9)

Based on the properties of the exponential function, i.e., exp(a

EQUATION

where

Z Y = Z Y /exp(b yl ). Since {exp(f 2k )} K k=1 are scalar, we can introduce a new partition function Z F2 = K k=1 exp(f 2k ) such that { 1 Z F 2 exp(f 2k )} K

k=1 become a probability measure. As a result, we can further reformulate P (F Y ) = {P (f yl )} L l=1 as a Product of Expert (PoE) model as

EQUATION

)

where

Z Y = Z Y /[exp(b yl ) • K k=1 [Z F2 ] γ kl ]

and each expert is defined as 1

Z F 2 exp(f 2k ).

It is noteworthy that all the experts { 1 Z F 2 exp(f 2k )} K k=1 form a probability measure and establish an exact one-to-one correspondence to all the neurons in the second hidden layer f 2 , i.e., {f 2k } K k=1 . Therefore, the distribution of f 2 can be expressed as

P (F 2 ) = {P (f 2k ) = 1 Z F2 exp(f 2k )} K k=1 (12)

Based on the definition of Gibbs distribution, Equation 8 and 12 show that f Y and f 2 formulate two multivariate Gibbs distributions and their energy functions are equivalent to the negative of the output nodes {f yl } L l=1 and the neurons {f 2k } K k=1 , respectively. In addition, Equation 11 implies that the connection between P (f yl ) and P (f 2k ) can be expressed as the PoE model, namely that P (f yl ) is equal to the product of all the experts {P (f 2k )} K k=1 with their respective weights {γ kl } K k=1 . Moreover, we prove that the probabilistic representation is valid for other hidden layers in the MLP. Since

{f 2k = σ 2 ( N n=1 β nk • f 1n + b 2k )} K

k=1 , P (F 2 ) can be expressed as

P (F 2 ) = {P (f 2k ) = 1 Z F2 exp[σ 2 ( N n=1 β nk • f 1n + b 2k )]} K k=1 (13)

Based on the equivalence between the gradient descent algorithm and the first order approximation (Battiti, 1992) , we prove that σ 2 (

EQUATION

)

where C 21 and C 22 only depend on the activations {f 1n } N n=1 in the previous training iteration, thus they can be regarded as constants and absorbed by β nk and b 2k . The proof for the approximation is included in Appendix G. Therefore, P (F 2 ) also can be modeled as a PoE model.

EQUATION

)

where

Z F2 = Z F2 /[exp(b 2k ) • N n=1 [Z F1 ] β nk ] and Z F1 = N n=1 exp(f 1n ).

Finally, the distribution of the first hidden layer f 1 is formulated as

EQUATION

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 (f 2k ) as a PoE model, in which all the experts { 1

Z F 1 exp(f 1n )} N n=1

are defined by all the neurons {f 1n } N n=1 . As a comparison, the GP explanation derives P (f 2k ) as a Gaussian distribution in the limit of infinite neurons, i.e., {f 1n } N →∞ n=1 , under the assumption that {f 1n } N n=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 f i with K convolutional channels, i.e.,

f i = {f k i = σ i ( N n=1 S (k,n) i • f n i-1 + b k i • 1)} K k=1

, formulate a MRF model as

EQUATION

where

f i-1 = {f n i-1 } N n=1 is the input, ϕ k i (f i-1 ) = N n=1 S (k,n) • f n i-1 + b k i • 1)

is the kth linear convolutional channel, S (k,n) is the kth convolutional filter applying to the nth input channel, b k i is the bias, σ i (•) is an non-linear operator, and the partition function

Z F i = exp[ K k=1 ϕ k (f i-1 )]df i-1 . The corresponding energy function E F i =-K k=1 ϕ k (f i-1 ).

In summary, the distribution of a hidden layer f i can be formulated as a Gibbs distribution P (F i ), especially the energy function E F i establish a connection between the deterministic function of the hidden layer and the Gibbs distribution P (F i ). The correspondence between Gibbs distributions and hidden layers can be summarized by Table 2 . 

f k i f i = {f k i } K k=1

Explanation energy function discrete Gibbs energy function MRF:

E F 1n =-f 1n {P (f 1n )} N n=1 E F i =-K k=1 ϕ k (f i-1 ) 1 Z F i exp[ K k=1 ϕ k i (f i-1 )]

We use Equation 16and Equation 17as examples to illustrate the correspondence.

Figure 4 : The top-left figure shows a 32×32 synthetic image as the input, which is sampled from the Gaussian distribution N (0, 1) and sorted in the primary diagonal direction by the descending order. The remaining figures visualize all the weights {αmn} 1024 m=1 of the seven neurons. In addition, the subscript in each figure, i.e., (gn(x), fn(x), P (fn)), denotes the linear output gn(x), activation fn(x), and Gibbs probability P (fn).

5. SIMULATIONS:

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.

5.1. THE DISTRIBUTION OF A FULLY CONNECTED LAYER IS A MULTIVARIATE DISCRETE GIBBS DISTRIBUTION:

A fully connected layer with N neurons, i.e., f = {f n } N n=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., E fn = -f n (x).

EQUATION

)

where f n (x) = σ[g n (x)], gn(x)= M m=1 αmn•xm+bn is a linear filter with α mn being the weights and b n being the bias, σ(•) is the activation function, and Z F = N n=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., g 1 (x) = 122.09, g 2 (x) = 113.70, and g 4 (x) = 140.20). However, Neuron 6 and Neuron 7 do not describe the input feature correctly, thus they derive large negative outputs (i.e., g 6 (x) = -134.08, g 7 (x) = -133.36). After applying the sigmoid function, we have

f 1 (x) = f 2 (x) = f 4 (x) = 1 and f 6 (x) = f 7 (x) = 0.

Finally, we derive the probability of different neurons, i.e., P (f 1 ) = P (f 2 ) = P (f 4 ) = 0.20 and P (f 6 ) = P (f 7 ) = 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. 

5.2. THE DISTRIBUTION OF A CONVOLUTIONAL LAYER IS A MRF MODEL:

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; f 1 ; f 2 ; f 3 ; f 4 ; f Y } for classifying the synthetic dataset. The CNN has two convolutional layers with ReLU (i.e., f 1 and f 3 ), two max-pooling layers (i.e., f 2 and f 4 ) and a softmax output layer f Y . The detailed information about the CNN is included in Appendix I.

Based on the MRF explanation (Equation 17), the distribution of f 1 can be formulated as

EQUATION

where

ϕ k 1 (x) = S k •x +b k 1 •1 is the kth linear convolutional channel, S k

is the linear convolutional filter and b k 1 is the bias. The energy function of P (F 1 ) is equivalent to

E F1 = - 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. 

{P (F n 1 ) = 1 Z F n 1 exp(f n 1 )} 5 n=1 .

The second row shows their respective distributions, and the KLDs denote their distances to the true distribution, i.e., KLD[P (X)||P (F n 1 )].

After the CNN is well trained, we randomly choose a testing image x as the input to derive P (F 1 ). Since x ∼ P (X), we can use the histogram of the testing image, i.e., P (x), to estimate P (X).

In addition, E F1 is a sufficient statistics for P (F 1 ) such that we can use the histogram of E F1 to estimate P (F 1 ). As a result, we can verify the probabilistic representation by calculating the distance between P (F 1 ) and P (X), i..e, KL[P (X)||P (F 1 )]. We visualize their distributions in Figure 6 , which shows that P (F 1 ) is very close to P (X), i.e., KL(P (X)||q(F 1 )) = 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 f 1 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.

6. DISCUSSION:

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.

6.1. SGD CAN BE EXPLAINED AS A VARIATIONAL INFERENCE:

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; f 1 ; ...; f I ; f Y }, we prove that the conditional distribution of f Y given x, i. 

EQUATION

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.

6.2. THE ENTIRE ARCHITECTURE OF DNNS CAN BE EXPLAINED AS A BAYESIAN NETWORK:

Since the output of the hidden layer f i in the DNN = {x; f 1 ; ...; f I ; f Y } is commonly the input of the next layer f i+1 , the DNN formulates a Markov chain

F 1 → • • • → F I → F Y .

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.

EQUATION

Richard & Lippmann (1991) prove that P (F Y |X) is an estimation of the true posterior distribution P (Y |X) ∝ P (X|Y )P (Y ) and we demonstrate that P (F 1 |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) .

6.3. A NOVEL REGULARIZATION ALGORITHM FOR DNNS:

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.

7. CONCLUSION:

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. 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 22and derive the correlation matrix for all the weights in every hidden layer.

EQUATION

where α n = 1 M M m=1 α mn is the sample mean of {α mn } M m=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 f 1 are uncorrelated because their correlation coefficients are very close to zero, i.e., |r(α n , α n )| ≈ 0, but many weights in f 2 and f Y are correlated to others, i.e., |r(β k , β k )| = 0 and |r(γ l , γ l )| = 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. 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.

A N × N f 1 (0.08) B K × K f 2 (0.26) C L × L f Y (0.

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 f 1 has 64 channels, f 1 has totally 3 × 64 = 192 filters. Hence, the dimension of the correlation matrix for f 1 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 f 3 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 f 1 and f 2 , 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 f 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 ; f Y } 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. 

, f2 = {σ2[g 2k (f1)]} K k=1 where g 2k (f1) = N n=1 β nk • f1n + b 2k . In addition, since f Y is softmax, f Y = { 1 Z Y exp[g yl (f2)]} L l=1 where g yl (f2) = K k=1 γ kl • f 2k + b yl , and the partition function is Z Y = L l=1 exp[g yl (f2)].

In the context of Bayesian probability, all the weights and the biases in the MLP (Figure 15 ) have prior distributions, i.e., α mn ∼ P (A mn ), 

β

G 2k = N n=1 W nk F 1n + B 2k as the random variable of g 2k (f 1 ) = N n=1 β nk • f 1n + b 2k .

Therefore, the random variables of arbitrary activations in f 1 and f 2 , e.g., f 1n and f 2k , can be expressed as

F 1n = σ 1 (G 1n ) and F 2k = σ 2 (G 2k ), respectively.

Based on the theorem (Appendix D) that functions of independent random variables are independent, we can derive that if {F 2k } K k=1 are independent and the activation function σ 1 (•) is invertible, then {G 2k } K k=1 are independent. In other words, in order to prove {F 2k } K k=1 are independent, we only need to prove {G 2k } K k=1 being independent as long as σ 2 (•) is invertible. The necessary condition for

EQUATION

Based on the definition of covariance, i.e., Cov(A, B) = E(AB) -E(A)E(B), we have

EQUATION

Based on the law of total expectation, i.e, E(AB) = E A (AE B|A (B)), we have

EQUATION

Since the statistical properties of the input {F 1n } N n=1 are unrelated to the weights β nk and biases b 2k , we have

E F1nF 1n |W nk W n k (F 1n F 1n ) = E F1nF 1n (F 1n F 1n ), E F1n|W nk B 2k (F 1n ) = E F1n (F 1n )

, and E F1n|W nk (F 1n ) = E F1n (F 1n ). As a result, we can derive that

E(W nk W n k • F 1n F 1n ) = E W nk W n k (W nk W n k )E F1nF 1n (F 1n F 1n ) E(W nk B 2k F 1n ) = E W nk B 2k (W nk B 2k )E F1n (F 1n ) E(W nk F 1n ) = E W nk (W nk )E F1n (F 1n ) (27)

Substituting Equation ( 27) for the corresponding expectations in Equation ( 25), we have

EQUATION

Assuming {F 1n } N n=1 are i.i.d. and F 1n ∼ P (F 1 ), we have

EQUATION

Substituting Equation ( 29) for the corresponding covariances in Equation ( 24), we have

EQUATION

Therefore, given σ 2 (•) is invertible and {F 1n } N n=1 are i.i.d., the necessary condition for {F 2k } K k=1 being independent can be expressed as

EQUATION

where

W k = N n=1 W nk and W k = N n=1 W nk .

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.

C PROOF: THE NECESSARY CONDITION FOR THE ACTIVATIONS OF HIDDEN LAYERS BEING IDENTICALLY DISTRIBUTED:

Given the same DNN drawn in Figure 15 , we assume σ 2 (•) being strictly increasing and differentiable, so σ 2 (•) is invertible and its inverse σ -1 2 (•) is also strictly increasing. The cumulative distribution function of F 2k can be expressed as

EQUATION

where φ(F 2k ≤ f ) is the probability of F 2k takes on a value less than or equal to f . Subsequently, we can obtain that

EQUATION

Equation ( 34) means that if the activation function σ 2 (•) is strictly increasing and differentiable, {F 2k } K k=1 being identically distributed implies {G 2k } K k=1 being identically distributed.

In the context of Bayesian probability,

EQUATION

Based on the law of total expectation and Equation 27, we have

EQUATION

Assuming {F 1n } N n=1 are i.i.d. and F 1n ∼ F 1 , we have E(F 1n ) = E(F 1 ). Hence, we can derive that

EQUATION

Based on the property of expectation, i.e., E(A) + E(B) = E(A + B), we have

EQUATION

Therefore, given σ 2 (•) is strictly increasing and differentiable and {F 1n } N n=1 are i.i.d., the necessary condition for {F 2k } K k=1 being identically distributed is that ∀(n, n ) ∈ S 1 , we have

EQUATION

where

W k = N n=1 W nk and W k = N n=1 W nk .

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.

D THEOREM: FUNCTIONS OF INDEPENDENT RANDOM VARIABLES ARE:



INDEPENDENT:

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 of A and B, respectively, we have

EQUATION

Based on the definition of independence, we can derive that

EQUATION

Based on the definition of preimage, we can derive that

EQUATION

Therefore, g(X) and h(Y ) are independent random variables.

E THE NECESSARY CONDITIONS FOR ACTIVATIONS BEING I.I.D. ARE VALID FOR CONVOLUTIONAL LAYERS:

Assuming a DNN = {x; f 1 ; ...; f I ; f Y } has two adjacent convolutional layers f i and f i+1 . The first convolutional layer has N channels, i.e. f i = {f n i } N n=1 , and the second convolutional layer has K channels, i.e. f i+1 = {f k i+1 } K k=1 . Therefore, the dimension of the convolutional filters S i+1 in f i+1 are W × H × N × K, i.e., the dimension of a single filter connecting the f n i channel and the

f k i+1 channel are D(S (k,n) i+1 ) = W × H.

As a result, the f k i+1 channel can be formulated as

EQUATION

where σ i+1 is the activation function and b i+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 f k i+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., f i = {x; x}, where x = {x 1 , • • • , x 9 } and x = {x 1 , • • • , x9 }, and the kth channel of the second convolutional layer is

f k i+1 = {z 1 , z 2 , z 3 , z 4 }. In addition, the dimension of the convolutional filter S i+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 f k i+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 matrix W , and the bias is b,

thus f k i+1 = W • x + b • 1.

If we consider the activation function and make a precise statement, f k i+1 can be formulated as

f k i+1 = σ i+1 (W k i+1 • f i + b k i+1 • 1) (44)

where W k i+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 f i+1 are i.i.d., the necessary conditions are

Cov(f k i+1 (m), f k i+1 (m )) = 0, E(f k i+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. 

) E(An) f2 E(F1) Cov(W k , W k ) E(W k ) f Y E(F2) Cov(C l , C l ) E(C l ) {xm} M m=1 ∼ P (X), {f1n} N n=1 ∼ P (F1), {f 2k } K k=1 ∼ P (F2) An = M

m=1 Amn, where αmn ∼ P (Amn) 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..

W k = N n=1 W nk , where β nk ∼ P (W nk ) C l = K k=1 C kl ,

F.1 THE MLP ON THE MNIST DATASET:

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 ) ∈ S 1 = {(k, k ) ∈ Z 2 |k = k , 1 ≤ k ≤ K, 1 ≤ k ≤ K}, we must have Cov(G 2k , G 2k ) = Cov(W k , W k )E 2 (F 1 ) = 0 (45) E(F 1 )E(W k ) = E(F 1 )E(W k ) (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(F 1 ), and E(F 2 ), we use the sample means x, f 1 , and f 2 , respectively. The definition of the sample mean is

x = 1 M M m=1 x m , f 1 = 1 N N n=1 f 2 (n), f 2 = 1 K K k=1 f 2 (k) (47)

After the MLP is well-trained, we derive the sample means as x = 0.131, f 1 = 0.436, and f 2 = 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(G 1n , G 1n ) = 0 and E(G 1n ) = E(G 1n ) = 0, thus f 1 satisfy the necessary conditions for activations being i.i.d..

However, f 1 = 0.436 and f 2 = 0.498 imply that E(F 1 ) = 0 and E(F 2 ) = 0 even considering the estimation error. Therefore, we need to estimate Cov(W k , W k ) and E(W k ) to examine whether the activations {f 2k } K k=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)} T t=1 to estimate the random variable W nk , thereby estimating Cov(W k , W k ) and E(W k ). More specifically, we use the sample correlation r(β k , β k ) to check whether Cov(W k , W k ) are close to zero and use the sample mean

EQUATION

)

where

β k = 1 T T t=1 β k (t) is the sample mean of {β k (t)} T t=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)} 20 t=1 . For an activation f 2k , it has 128 weights because its previous layer has 128 neurons. Hence, a single sample β k (t) = We use the same method to estimate Cov(A n , A n ), E(A n ), Cov(C l , C l ), and E(C l ).

Finally, we shows the absolute correlation matrices i.e., A N ×N , W K×K , and C L×L , in Figure 3 . Since many elements in the three matrices are far from zero, we conclude that their corresponding covariances, i.e, Cov(A n , A n ), Cov(W k , W k ), and Cov(C k , C k ), are not close to zero based on connection between correlation and covariance

Corr(W k , W k ) = Cov(W k , W k ) σ W k σ W k (49)

where σ W k > 0 is the standard derivation of W k .

Figure 17 visualizes the expectations, i.e., E(A n ), E(W k ), and E(C l ), for each layer in the MLP. We find that some expectations show great variations, especially {E(A n )} N n=1 and {E(C l )} L l=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.

F.2 THE CNN ON THE CIFAR-10 DATASET:

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

3×3×3×64 3×3×64×128 8192×1024 1024×256 256×10

Weight Dim.

M ×N N ×K K×L L×Q Q×R (27×64) (576×128) (8192×1024) (1024×256) (256×10) Weight PDF amn∼P (Amn) b nk ∼P (B nk ) c kl ∼P (C kl ) d lq ∼P (D lq ) fqr ∼P (Fqr )

Sum. PDF

An= M m=1 Amn B k = N n=1 B nk C l = K k=1 C kl Dq = L l=1 D lq Fr = Q q=1 Fqr Covariance Cov(An,A n ) Cov(B k ,B k ) Cov(C l ,C l ) Cov(Dq ,D q ) Cov(Fr ,F r )

Sum. Exp.

E(An) E(B k ) E(C l ) E(Dq ) E(Fr ) Activation Exp. E(X) E(F 1 ) E(F 2 ) E(F 3 ) E(F 4 )

Dim. is short for dimension. Sum. is short for Weights Summation. Exp. is short for Expectation.

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

EQUATION

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 f 1 are 3 × 3 × 3, thus a single convolutional channel in f 1 has 27 weights, so the weights dimension of f 1 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, f 1 = -707.074, f 2 = -6347.836, f 3 = 36.870, f 4 = 15.131, thus all the expectations cannot be zero except the input layer, i.e., E(X) = 0,

E(F 1 ) = 0, E(F 2 ) = 0, E(F 3 ) = 0, E(F 4 ) = 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(C l )} L l=1 and {E(D q )} Q q=1 , which means that most expectations of weights summations of different neurons in f 3 and f 4 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. et al., 1986) . We typically update θ iteratively to derive θ * , which can be expressed as

G THE EQUIVALENCE:

EQUATION

where

∇ θt H[f Y , P (Y |X)] is the Jacobian matrix of H[f Y , 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(f Y ) as H[f Y , 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 of H(f Y ) with respect to the parameters of the ith hidden layer, i.e., ∇ θ(i) H(f Y ), can be expressed as follows based on the chain rule.

EQUATION

where θ(i) denote the parameters of the ith hidden layer.

Based on Equation 52and Equation 53, θ(i) can be learned by the gradient descent method as

EQUATION

Table 8 summarizes the backpropagation training procedure for the MLP in Figure 15 . 

f Y ∇ f Y H(f Y )∇ θ t (Y ) f Y ↓ θ t+1 (Y )=θ t+1 (Y )-α[∇ θ t (Y ) H(f Y )], f Y (f 2 ,θ t+1 (Y )) ↑ f2 ∇ f Y H(f Y )∇ f 2 f Y ∇ θ t (2) f 2 ↓ θ t+1 (2)=θ t+1 (2)-α[∇ θ t (2) H(f Y )], f 2 (f 1 ,θ t+1 (2)) ↑ f1 ∇ f Y H(f Y )∇ f 2 f Y ∇ f 1 f 2 ∇ θ t (1) f 1 ↓ θ t+1 (1)=θ t+1 (1)-α[∇ θ t (1) H(f Y )], f 1 (x,θ t+1 (1)) ↑ x - -

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

If an arbitrary function f is differentiable at point p * in R N and its differential is represented by the Jacobian matrix ∇ p f , the first order approximation of f near the point p can be formulated as 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

EQUATION

EQUATION

where

f 2 [f 1 , θ 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., f 1 . The

EQUATION

where 56, we can derive

∇ σ2 f 2 = ∂f2[f1,θt(2)] ∂σ2 . Substituting (∇ σ2 f 2 ) • [f 1 ; 1] T for ∇ θt(2) f 2 in Equation

EQUATION

If we only consider a single neuron, such as

f 2k , in f 2 , θ t+1 (2k) = [β 1k ; • • • ; β N k ; b 2k ] and θ t (2k) = [β * 1k ; • • • ; β * N k ; b * 2k ], thus [f 1 ; 1] T • θ t+1 (2k) = N n=1 β nk • f 1n + b 2k .

As a result, Equation 58 with respect to a single neuron can be expressed as

EQUATION

Equation 59 indicates that f 2k [f 1 , θ 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.,

f 2k [f 1 , θ t (2k)]. Since ∇ σ2 f 2k = ∂f 2k [f1,θt(2)]

∂σ2 is only related to f 1 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, f 2k (f 1 , θ t+1 (2k)) can be reformulated as

EQUATION

where

C 1 = ∇ σ2 f 2k and C2 = f 2k (f 1 , θ t (2k)) -(∇ σ2 f 2k ) • [ N n=1 β * nk • f 1n + b * 2k ].

Similarly, the activations in the first hidden layer, i.e, f 1 , also can be formulated as the first order approximation in the context of the gradient descent learning algorithm. We flatten the output of f2 as the input of f Y .

H THE PROPOSED PROBABILISTIC REPRESENTATION HOLDS FOR CONVOLUTIONAL LAYERS:

This section proves that the proposed probabilistic representation is valid for convolutional layers.

In the above CNN, the input

x has Q channels x = {x q } Q q=1 (e.g., Q = 3 if x are color images), f 1 has N convolutional channels f 1 = {f n 1 } N n=1 , f 2 has K convolutional channels f 2 = {f k 2 } K k=1 .

The output layer f Y is the softmax with L nodes, thus its distribution can be formulated as 

EQUATION

EQUATION

)

where {f 2j } J j=1 is the flattened output of f 2 , γ jl is the weight of the edge between f 2j and f yl , and b yl denotes the bias. Since f 2 is a convolutional layer with K convolutional channels, it can be reformulated as {f 2j } J j=1 = {f k 2 } K k=1 . As a result, P (F Y ) can be reformulated as

P (F Y ) = {P (f yl ) = 1 Z F Y exp( K k=1 γ T kl • f k 2 + b yl )} L l=1 (63)

where γ kl is a subset of all the parameters {γ jl } J j=1 such that  and γ T kl is the transpose of γ kl . Therefore, P (F Y ) can be reformulated as

K k=1 γ T kl • f k 2 = J j=1 γ jl • f 2j ,

EQUATION

)

where Based on the element-wise matrix power, we can reformulate P (F Y ) as 

Z F Y = Z F Y /[

EQUATION

P (F Y |X) = • • • P (F Y ; F I ..; F 1 |X)dF I • • • dF 1 (79) Since P (F Y ; F I ..; F 1 |X) = P (F Y |F I ) • ... • P (F i+1 |F i ) • ... • P (F 1 |X)

and only F 2 depends on F 1 , we can derive

P (F Y |X) = • • • P (F Y |F I ) • ... • P (F 3 |F 2 )dF I • • • dF 2 • P (F 2 |F 1 )P (F 1 |X)dF 1 (80) Since P (F 2 |F 1 )P (F 1 |X)dF 1 = P (F 2 ; F 1 |X)dF 1 = P [F 2 |f 1 (X)],

where f 1 (•) is the deterministic function defined by the corresponding hidden layer, P (F Y |X) can be simplified as

P (F Y |X) = • • • P (F Y |F I ) • ... • P (F 3 |F 2 )P [F 2 |f 1 (X)]dF I • • • dF 2 (81)

Similarly, since only F 3 depends on F 2 , we can obtain

P (F Y |X) = • • • P (F Y |F I ) • ... • P (F 3 |f 2 (f 1 (X)))dF I • • • dF 3 (82)

If we iteratively apply the integral with respect to F 3 until F I , we can obtain

EQUATION

Since the output layer f Y is commonly a fully connected layer. we can formulate P (F Y |X) as a multivariate discrete Gibbs distribution based on Equation 8.

EQUATION

It is important to note that the energy function 

E DNN (f Y |x) = -f Y (f I (• • •f 2 (f 1 (x))• • •)), i.

(x) = exp[f Y (f I (• • •f 2 (f 1 (x))• • •))]df Y (f I (• • •f 2 (f 1 (x))• • •)).

For clarity, here we only consider integral no matter a Gibbs distribution is discrete or continuous, because f1

•••f I exp[-E DNN (f Y |x)

] is a special case of integral. All functions involving f i should include parameters θ fi in the conditioning set, but we omit for the same reason.

J.2 SGD CAN BE EXPLAINED AS A VARIATIONAL INFERENCE ALGORITHM:

As a dominant paradigm for Bayesian posterior inference P (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) .

EQUATION

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. 

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 (F Y |X) (Equation 84) and SGD aims to find P * (F Y |X) such that it can minimize the distance to the true posterior distribution P (Y |X).

J.3 SGD IS EQUIVALENT TO THE ENERGY MINIMIZATION OPTIMIZATION:

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 (F Y |X)] = -log[P (f Y (l * )|x)].

Based on the definition of P (F Y |X) (Equation 84), we can obtain

EQUATION

where f yl * is the l * th output node of the output layer f Y . In particular, since log[Z DNN (x)] is not related to the output f Y , 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

EQUATION

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. Based on the architecture of CNN in Table 9 , E CN N (f yl * |x) can be formulated as

EQUATION

In Figure 23 , the first row visualizes the variation of E CN N (f yl * |x) over 30 training epochs. Though there are certain fluctuations, E CN N (f yl * |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 (F Y |X) corresponding to the MLP is defined as

EQUATION

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.

K.2 EXPERIMENT: THE PARAMETERS LEARNED BY SGD ARE NOT I.I.D.

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 2to 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 f 1 and f 2 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 f 2 is huge, i.e., the weights of f 2 cannot be i.i.d. anymore. However, the increase for f 1 is very small, i.e., all the weights of f 1 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 f 1 and f 2 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 f 1 and f 2 are not i.i.d. anymore. In addition, the average correlation for f 3 and f 4 also have an increasing trend, but their increase are very small, i.e., the parameters of f 3 and f 4 still could be i.i.d.. Overall, we conclude that not all the learned parameters in the CNN are i.i.d.. Based on the MRF explanation for convolutional layers, the first convolutional layer f 1 formulates a MRF model as the prior distribution P (F 1 ), which can be expressed as

EQUATION

where ϕ k 1 (x) = S k •x +b k 1 •1 is the kth linear convolutional channel, S k is the linear convolutional filter and b k 1 is the bias. 

× 1 × 20 1 × 1 × 20 1 × 1 × 20 F Y f Y Output(softmax) 1 × 1 × 10 1 × 1 × 10 1 × 1 × 10 R.V.

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

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 (F 1 )] 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 f 1 , 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 (F 1 )] and the testing error. Specifically, the testing error decreases as the learned prior distribution P (F 1 ) approaching to the truly prior distribution P (X), and the testing error becomes stable when P (F 1 ) 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. 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 {x j } 20000 j=1 , we design a FoE model with 12 Gaussian Scale Mixture (GSM) experts, thus it can be expressed as

EQUATION

where 11 . It is noteworthy that the number of convolutional filters in f 1 and the filter dimension are the same as the above FoE model, thus we can use the well-trained FoE to initialize f 1 as a regularization for the CNN classifying the CIFAR-10 dataset based on Algorithm 2.

f GSM k [f k (x)] = 11 i=1 π ki • N (f k (x); 0,

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. Bayes prior indicates the proposed regularization and its testing error is 29.90%.

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. 

:

https://www.nist.gov/srd/nist-special-database-19

Figures and Tables:

Figure 1: A MLP = {x; f1; f2; f Y }, in which f1 has N neurons, i.e., f1 = {f1n = σ1[g1n(x)]} N n=1 , where g1n(x) = M m=1 αmn • xm + b1n is the nth linear filter, αmn is the weight of the edge between xm and f1n, b1n is the bias, and σ1 is a non-linear activation function. Similarly, f2 = {σ2[g 2k (f1)]} K k=1 where g 2k (f1) = N n=1 β nk • f1n + b 2k . In addition, since f Y is the softmax, f Y = { 1 Z Y exp[g yl (f2)]} L l=1 where g yl (f2) = K k=1 γ kl • f 2k + b yl , and the partition function is Z Y = L l=1 exp[g yl (f2)].

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.

Figure 2: Three absolute correlation matrices for the weights in f1, f2, and f Y . In addition, 0.08, 0.26, and 0.31 are the average of the absolute correlation coefficients for all the weights in f1, f2, and f Y , respectively.

shows the correlation matrices, i.e., A N ×N , W K×K , and C L×L , to check if their corresponding covariances, i.e, Cov(A n , A n ), Cov(W k , W k ), and Cov(C k , C k ), are close to zero.

Figure 3: Three absolute correlation matrices for showing whether Cov(An, A n ), Cov(W k , W k ), and Cov(C k , C k ) are close to zero.

Figure 5: The first row shows five synthetic images of handwritten digits, the second row visualizes their respective histograms, and the red curve indicates the Gaussian distribution N (0, 1024).

Figure 6: The red curve indicates the truly prior distribution P (X) = N (0, 1024). The blue curves are different histograms. (A) the synthetic image x is the input. (B) the histogram of x, i.e., P (x), and KL[P (X)||P (x)] = 0.62. (C) the histogram of E F 1 for estimating P (F1) and KL[P (X)||P (F1)] = 0.83.

e., P (F Y |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 * (F Y |X) such that it can minimize the distance between P (F Y |X) and the true posterior distribution P (Y |X) is determined by the training dataset D.

Figure 8: Three absolute correlation matrices for the weights in f1, f2, and f Y . In addition, 0.08, 0.26, 0.31 are the average of the absolute correlation coefficients for all the weights in f1, f2, and f Y , respectively.

Figure 9: A neural network with a single hidden layer NN = {x; f ; f Y }.

Figure11shows the absolute correlation matrices for every layer. We can find that most weights are correlated in f 1 and f 2 , 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 f 3 and f 4 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 Figure11, 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.

Figure12: Four synthetic images with different features. All images are sampled from N (0, 1). The only difference is that they are sorted in different diagonal directions by the ascending or descending orders. Image0 and Image1 are sorted in the primary diagonal direction by the ascending order and the descending order, respectively. Image2 and Image3 are sorted in the secondary diagonal direction by the ascending order and the descending order, respectively.

Figure 13: Eight neurons in the fully connected hidden layer. Their dimension is 32 × 32 Since the synthetic dataset is very simple, we set the hidden layer f have only eights neurons and the output layer f Y 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: The left figure visualizes the four linear filters in the output layer. The right figure visualizes the correlation matrix of the four filters in the output layer.

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

nk ∼ P (W nk ), b 1n ∼ P (B 1n ), and b 2k ∼ P (B 2k ), thus we regard G 1n = M m=1 A mn X m +B 1n as the random variable of g 1n (x) = M m=1 α mn •x m +b 1n , and

Figure16: The equivalence between a convolution channel and a fully connected layer.

where γ kl ∼ P (C kl ) F ACTIVATIONS BEING I.I.D. CANNOT BE VALID FOR ALL THE HIDDEN LAYERS OF DNNS

Figure 17: The expectation of the weights summation, i.e., E(An), E(W k ), and E(C l ), in each layer. The number at the top of the figure indicates the variance of the points drawn in the figure.

n=1 β nk (t). Based on the 20 samples, i.e., {β k (t)} 20 t=1 , we can calculate r(β k , β k ) to examine whether Cov(W k , W k ) are close to zero and derive the sample mean to estimate E(W k ).

Figure 18: Five absolute correlation matrices for the weight summations in each layer. The number at the top is the average of the absolute correlation coefficients for all the weights summations of the layer

BETWEEN THE GRADIENT DESCENT ALGORITHM AND THE FIRST ORDER APPROXIMATION In the context of deep learning, most learning algorithms belong to the gradient descent algorithm. Given a DNN = {x; f 1 ; ...; f I ; f Y } and a cost function H[f Y , 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[f Y , P (Y |X)] (Rumelhart

) where o(||p -p * ||) is a quantity that approaches zero much faster than ||p -p * || approaches zero.

Figure 20: The above CNN has two convolutional layers, i.e., f1, f2, and one softmax output layer, i.e., f Y .

where Z F Y = L l=1 exp(f yl ) denotes the partition function. Since {f yl = J j=1 γ jl • f 2j + b yl } L l=1 , we can derive

exp(b yl )]. Recall the element-wise matrix power, e.g., exp(a) b , where a = [a 1 ; a 2 ; a 3 ] and b = [b 1 ; b 2 ; b 3 ], we can derive that exp(a) = [exp(a 1 ); exp(a 2 ); exp(a 3 )] and exp(a) b = [exp(a 1 ) b1 ; exp(a 2 ) b2 ; exp(a 3 ) b3 ] = [exp(a 1 b 1 ); exp(a 2 b 2 ); exp(a 3 b 3 )]. As a result, exp(ab) = exp(a 1 b 1 + a 2 b 2 + a 3 b 3 ) = |a| exp(a) b , where |a| is the element numbers of a.

Figure 22: 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., moutside, moutside-boundary, minside-boundary, and minside. The third row shows the synthetic image corresponding to each mask part, i.e., x = {xoutside, xoutside-boundary, xinside-boundary, xinside}.

Given a DNN = {x; f 1 ; ...; f I ; f Y }, we choose the cross entropy H[P (Y |X), P (F Y |X)] as the loss function, where P (Y |X) is the true posterior distribution determined by the training dataset D and P (F Y |X) is the conditional distribution of f Y given x derived from the DNN.H[P (Y |X), P (F Y |X)] = H[P (Y |X)] + KL[P (Y |X) P (F Y |X)] [y(l)|x]log[P (f Y (l)|x)](86)whereH[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 (f Y |x) is the learned conditional distribution of f Y given the same x, e.g., P (f Y |x) could be [0.2, 0.7, 0.1, 0]. In the context of deep learning, SGD aims to find the P * (F Y |X) such that it can minimize the loss function H[P (Y |X), P (F Y |X)] given the training dataset D. P * (F Y |X) = arg min P (F Y |X)∈Q H[P (Y |X), P (F Y |X)] (87) Given the training dataset D, P (Y |X) is a constant, i.e., H[P (Y |X)] = 0, thus we can obtain P * (F Y |X) = arg min P (F Y |X)∈Q KL[P (Y |X) P (F Y |X)]

Figure 23: The first row shows the testing error of the CNN and the variation of ECNN (f yl * |x) given three testing datasets with different labels over 30 training epochs. The second row shows testing error of the MLP and the variation of EMLP (f yl * |x) given three testing datasets with different labels over 30 training epochs.

Figure 24: (A) shows the testing error over 15 training epochs. (B), (C), and (D) show the average correlation between different weights of three layers in the MLP over 15 training epochs.

Figure 25: shows the testing error over 30 training epochs. (B), (C), (D), (E) and (F) show the average correlation between different weights of five layers in the CNN over 30 training epochs.

5) + ReLU 11 × 11 × 20 11 × 11 × 20 11 × 11 ×

Figure 26: (A) shows the average KL[P (X)||P (F1)] overall all testing images. (B) shows 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 convolutional channels of f1.

) ), σ 2 b is a fixed base variance, δ = {δ(1), • • • , δ(11)} is a range of scales, f k (•) 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 = [f 1 , ..., f 12 ] and π = {π ki } k=12,i=11 k=1,i=1 . The architecture of the CNN for image recognition is summarized in Table

Figure 28: The average testing error of the CNN on the last 5 training epochs based on various regularizations.

Figure 29: The average testing error of the CNN on the last 5 training epochs based on various combined regularizations. Bayes prior + both indicates that we combine all three regularizations together.

The statistical measures for examining the necessary conditions

Table content: <table><tr><td>Layer Activation expectation</td><td>Covariance</td><td>Weight expectation</td></tr><tr><td>f1</td><td/><td/></tr></table>

where γ kl ∼ P (C kl ) 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.. 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 W nk and B 2k , thus we can use the sample covariance to estimate Cov(W k , W k ), Cov(W k , B 2k ), and etc.

Table content: <table/>

The correspondence between hidden layers and Gibbs distributions

Table content: <table><tr><td>Architecture</td><td>neuron f1n</td><td>fully connected layer f1 = {f1n} N n=1</td><td>conv. channel</td><td>conv. layer</td></tr></table>

The architecture of the MLP for MNIST classification

Table content: <table><tr><td>Layer</td><td>Description</td><td colspan="2">Dimension Filter dimension (M × N )</td></tr><tr><td>x</td><td>Input</td><td>784</td><td>-</td></tr><tr><td>f1</td><td>FC(sigmoid)</td><td>128</td><td>128 × 784</td></tr><tr><td>f2</td><td>FC(sigmoid)</td><td>64</td><td>64 × 128</td></tr><tr><td>f Y</td><td>Output(softmax)</td><td>10</td><td>10 × 64</td></tr><tr><td colspan="3">FC denotes the fully connected layer</td><td/></tr><tr><td colspan="4">Dimension: the number of nodes (i.e., neurons) in every layer.</td></tr><tr><td colspan="4">Filter dimension (M × N ): a layer has M filters with N × 1 dimension</td></tr><tr><td colspan="4">A THE PARAMETERS OF DNNS CANNOT SATISFY THE I.I.D. PRIOR</td></tr><tr><td colspan="2">A.1 A MLP ON THE MNIST DATASET</td><td/><td/></tr></table>

The architecture of the NN = {x; f ; f Y } for MNIST classification NEURAL NETWORK WITH A SINGLE HIDDEN LAYER ON THE MNIST DATASET Based on the seminal work(Neal, 1994), we design a neural network with a single fully connected hidden layer (Figure9) 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. Table4summarizes the architecture of the network. In addition, {α mn } M m=1 denote the weights of the nth linear filter in the fully conntected hidden layer and {β nl } N n=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 } 784 m=1 . The output layer has totally 10 linear filters and their dimension is 1 × 128, i.e., {β nl } 128 n=1 .

Table content: <table><tr><td>Layer</td><td colspan="2">Description</td><td colspan="4">Dimension Filter dimension (M × N )</td></tr><tr><td>x</td><td>Input</td><td/><td>784</td><td/><td/><td>-</td></tr><tr><td>f</td><td colspan="2">FC(sigmoid)</td><td>128</td><td/><td/><td>128 × 784</td></tr><tr><td>f Y</td><td colspan="2">Output(softmax)</td><td>10</td><td/><td/><td>10 × 128</td></tr><tr><td colspan="4">FC denotes the fully connected layer</td><td/><td/></tr><tr><td colspan="7">Dimension: the number of nodes (i.e., neurons) in every layer.</td></tr><tr><td colspan="7">Filter dimension (M × N ): a layer has M filters with 1 × N dimension</td></tr><tr><td colspan="4">A.2 A 𝑓 "" 𝑥 "</td><td>𝛽 ""</td><td/></tr><tr><td/><td>𝑥 %</td><td>⋮</td><td>𝑓 "' 𝛼 &amp;' 𝛼 %' 𝛼 "'</td><td>𝛽 '" 𝛽 /" 𝛽 ("</td><td>⋮</td><td>𝑓 ) (1)</td></tr><tr><td/><td colspan="2">⋮</td><td/><td/><td/><td>𝑓 ) (𝐿)</td></tr><tr><td/><td/><td/><td>𝑓 "(</td><td/><td/></tr><tr><td/><td>𝑥 &amp;</td><td/><td>𝑓 "/ ⋮</td><td/><td/></tr><tr><td/><td colspan="2">𝒙</td><td>𝒇</td><td/><td colspan="2">𝒇 𝒀</td></tr></table>

The architectures of the CNN for CIFAR-10 classification

Table content: <table><tr><td>Layer</td><td>Description</td><td>Output dimension</td><td>Filters dimension</td><td>Correlation matrix</td></tr><tr><td>x</td><td>Input</td><td>32 × 32 × 3</td><td>-</td><td>-</td></tr><tr><td>f1</td><td>Conv (3 × 3) + Maxpool</td><td>16 × 16 × 64</td><td>3 × 3 × 3 × 64</td><td>A192×192</td></tr><tr><td>f2</td><td>Conv (3 × 3) + Maxpool</td><td>8 × 8 × 128</td><td>3 × 3 × 64 × 128</td><td>B8192×8192</td></tr><tr><td>f3</td><td>FC + Sigmoid</td><td>1024</td><td>8192 × 1024</td><td>C1024×1024</td></tr><tr><td>f4</td><td>FC + Sigmoid</td><td>256</td><td>1024 × 256</td><td>D256×256</td></tr><tr><td>f Y</td><td>Output(softmax)</td><td>10</td><td>256 × 10</td><td>E10×10</td></tr><tr><td colspan="3">A.3 A CNN ON THE CIFAR-10 DATASET</td><td/><td/></tr></table>

The statistical measures for examining the necessary conditions

Table content: <table><tr><td colspan="2">Layer Activation expectation</td><td>Covariance</td><td>Weight expectation</td></tr><tr><td>f1</td><td>E(X)</td><td>Cov(An, A n</td></tr></table>

The statistical measures for examining the necessary conditions of each layer in the CNN

Table content: <table><tr><td/><td>f1</td><td>f2</td><td>f3</td><td>f4</td><td>f Y</td></tr><tr><td>Description</td><td>Conv Maxpool</td><td>Conv Maxpool</td><td>FC Sigmoid</td><td>FC Sigmoid</td><td>Output Softmax</td></tr><tr><td>Activation Dim.</td><td>16×16×64</td><td>8×8×128</td><td>1024×1</td><td>256×1</td><td>10×1</td></tr></table>

One iteration of the backpropagation training procedure for the MLP Layer Gradients update ∇ θ t (i) H(f Y ) Parameters and activations update

Table content: <table/>

The architectures of the CNN for synthetic image classificationThe architecture of the CNN for classifying the synthetic dataset is summarized in Table9.J THE EQUIVALENCE BETWEEN SGD AND VARIATIONAL INFERENCE J.1 THE DISTRIBUTION OF THE OUTPUT LAYER GIVEN THE INPUT IS A GIBBS DISTRIBUTIONIn Section 6.2, we demonstrate the entire architecture of a DNN = {x; f 1 ; ...; f I ; f Y } corresponds to a joint distribution P (F Y ; F I ..; F 1 |X). As a result, the conditional distribution of f Y given x, i.e., P (F Y |X), can be formulated

Table content: <table><tr><td colspan="2">R.V. Layer</td><td>Description</td><td>CNN</td></tr><tr><td>X</td><td>x</td><td>Input</td><td>32 × 32 × 1</td></tr><tr><td>F1</td><td>f1</td><td colspan="2">Conv (3 × 3) + ReLU 30 × 30 × 20</td></tr><tr><td/><td>f2</td><td>Maxpool</td><td>15 × 15 × 20</td></tr><tr><td>F2</td><td>f3</td><td colspan="2">Conv (5 × 5) + ReLU 11 × 11 × 60</td></tr><tr><td/><td>f4</td><td>Maxpool</td><td>5 × 5 × 60</td></tr><tr><td>F Y</td><td>f Y</td><td>Output(softmax)</td><td>1 × 1 × 10</td></tr><tr><td colspan="4">R.V. is the random variable of the hidden layer(s).</td></tr></table>

e., E DNN (•) is determined by the entire architecture of DNNs. In addition, since P (F Y |X) is a conditional distribution of F Y given X, the partition function Z DNN (x) should integrate all possible activations of the output layer f Y to guarantee P (F Y |X) being a valid conditional probability, i.e., Z DNN

Table content: <table/>

The architectures of CNNs for classifying the synthetic dataset

Table content: <table><tr><td colspan="2">R.V. Layer</td><td>Description</td><td>CNN1</td><td>CNN2</td><td>CNN3</td></tr><tr><td>X</td><td>x</td><td>Input</td><td>32 × 32 × 1</td><td>32 × 32 × 1</td><td>32 × 32 × 1</td></tr><tr><td>F1</td><td>f1 f2</td><td>Conv (3 × 3) +ReLU Maxpool</td><td>30 × 30 × 4 15 × 15 × 4</td><td colspan="2">30 × 30 × 12 30 × 30 × 36 15 × 15 × 12 15 × 15 × 36</td></tr></table>

The architectures of CNN for CIFAR-10 classification Testing error of CNN without any regularization (33.95%) Testing error of CNN with droupout (32.26%) Testing error of CNN with L2 norm (29.58%) Testing error of CNN with Bayes prior (29.90%)

Table content: <table><tr><td>Layer</td><td>Description</td><td>Output</td></tr><tr><td>x</td><td>Input</td><td>32 × 32 × 1</td></tr><tr><td>f1</td><td>Conv (3 × 3)</td><td>30 × 30 × 12</td></tr><tr><td>f2</td><td>ReLU</td><td>30 × 30 × 12</td></tr><tr><td>f3</td><td>Conv (5 × 5)</td><td>26 × 26 × 32</td></tr><tr><td>f4</td><td colspan="2">ReLU + Maxpool 13 × 13 × 32</td></tr><tr><td>f5</td><td>Conv (5 × 5)</td><td>9 × 9 × 128</td></tr><tr><td>f6</td><td>ReLU + Maxpool</td><td>4 × 4 × 128</td></tr><tr><td>f7</td><td>Fully connected</td><td>1 × 1 × 512</td></tr><tr><td>f8</td><td>ReLU</td><td>1 × 1 × 512</td></tr><tr><td>f9</td><td>Fully connected</td><td>1 × 1 × 32</td></tr><tr><td>f10</td><td>ReLU</td><td>1 × 1 × 32</td></tr><tr><td>f Y</td><td>Output(softmax)</td><td>1 × 1 × 10</td></tr></table>

Peer review (iclr2020):
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. The paper is well written and well postulated. > Fig 4: What is the information presented by each neuron? How would this have looked with the i.i.d. prior in place. > There are places in the paper where one must refer to the supplementary, for example sections H and J with the simulations. Do consider moving these crucial sections to the main paper. > One recurring thought I had when the authors bring up Bayesian Hierarchical model, is that most of the BHMs rely on i.i.d assumptions both in the prior space and with the observations. How would you stand by your claim of explaining a DNN*s layers to be modelled as a BHM?