Title: GAUSSIAN PROCESS BEHAVIOUR IN WIDE DEEP NEURAL NETWORKS

Abstract: Whilst deep neural networks have shown great empirical success, there is still much work to be done to understand their theoretical properties. In this paper, we study the relationship between Gaussian processes with a recursive kernel definition and random wide fully connected feedforward networks with more than one hidden layer. We exhibit limiting procedures under which finite deep networks will converge in distribution to the corresponding Gaussian process. To evaluate convergence rates empirically, we use maximum mean discrepancy. We then exhibit situations where existing Bayesian deep networks are close to Gaussian processes in terms of the key quantities of interest. Any Gaussian process has a flat representation. Since this behaviour may be undesirable in certain situations we discuss ways in which it might be prevented. 1

1. INTRODUCTION:

Deep feedforward neural networks have emerged as an essential component of modern machine learning. As such there has been significant research effort in trying to understand the theoretical properties of such models. One important branch of such research is the study of random networks. By assuming a probability distribution on the network parameters, a distribution is induced on the input to output function that such networks encode. This has proved important in the study of initialisation and learning dynamics (Schoenholz et al., 2017) and expressivity (Poole et al., 2016) . It is, of course, essential in the study of Bayesian priors on networks (Neal, 1996) . The Bayesian approach makes little sense if prior assumptions are not understood, and distributional knowledge can be essential in finding good posterior approximations.

Since we typically want our networks to have high modelling capacity, it is natural to consider limit distributions of networks as they become large. Whilst distributions on deep networks are generally challenging to work with exactly, the limiting behaviour can lead to more insight. Further, as we shall see, networks used in the literature may be very close to this behaviour.

The seminal work in this area is that of Neal (1996) , which showed that under certain conditions random neural networks with one hidden layer converge to a Gaussian process. The question of the type of convergence is non-trivial and part of our discussion. Historically this result was a significant one because it provided a connection between flexible Bayesian neural networks and Gaussian processes (Williams, 1998; Rasmussen & Williams, 2006) 1.1 OUR CONTRIBUTIONS We extend the theoretical understanding of random fully connected networks and their relationship to Gaussian processes. In particular, we prove a rigorous result (Theorem 1) on the convergence of certain finite networks with more than one hidden layer to Gaussian processes.

Further, we empirically study the distance between finite networks and their Gaussian process analogues by using maximum mean discrepancy (Gretton et al., 2012) as a distance measure. We find that Bayesian deep networks from the literature can exhibit predictions that are close to Gaussian processes. To demonstrate this, we systematically compare exact Gaussian process inference with 'gold standard' MCMC inference for Bayesian neural networks.

Our work is of relevance to the theoretical understanding of neural network initialisation and dynamics. It is also important in the area of Bayesian deep networks because it demonstrates that Gaussian process behaviour can arise in more situations of practical interest than previously thought. If this behaviour is desired then Gaussian process inference (exact and approximate) should also be considered. In some scenarios, the behaviour may not be desired because it implies a lack of a hierarchical representation. We therefore highlight promising ideas from the literature to prevent such behaviour.

1.2. RELATED WORK:

The case of random neural networks with one hidden layer was studied by Neal (1996) . Cho & Saul (2009) provided analytic expressions for single layer kernels including those corresponding to a rectified linear unit (ReLU). They also studied recursive kernels designed to 'mimic computation in large, multilayer neural nets'. As discussed in Section 3 they arrived at the correct kernel recursion through an erroneous argument. Such recursive kernels were later used with empirical success in the Gaussian process literature (Krauth et al., 2017) , with a similar justification to that of Cho and Saul. The first case we are aware of using a Gaussian process construction with more than one hidden layer is the work of Hazan & Jaakkola (2015) . Their contribution is similar in content to Lemma 1 discussed here, and the work has had increasing interest from the kernel community (Mitrovic et al., 2017) . Recent work from Daniely et al. (2016) uses the concept of 'computational skeletons' to give concentration bounds on the difference in the second order moments of large finite networks and their kernel analogue, with strong assumptions on the inputs. The Gaussian process view given here, without strong input assumptions, is related but concerns not just the first two moments of a random network but the full distribution. As such the theorems we obtain are distinct. A less obvious connection is to the recent series of papers studying deep networks using a mean field approximation (Poole et al., 2016; Schoenholz et al., 2017) . In those papers a second order approximation gives equivalent behaviour to the kernel recursion. By contrast, in this paper the claim is that the behaviour emerges as a consequence of increasing width and is therefore something that needs to be proved. Another surprising connection is to the analysis of self-normalizing neural networks (Klambauer et al., 2017) . In their analysis the authors assume that the hidden layers are wide in order to invoke the central limit theorem. The premise of the central limit theorem will only hold approximately in layers after the first one and this theoretical barrier is something we discuss here. An area that is less related than might be expected is that of 'Deep Gaussian Processes' (DGPs) (Damianou & Lawrence, 2013 ). As will be discussed in Section 6, narrow intermediate representations mean that the marginal behaviour is not close to that of a Gaussian process. Duvenaud et al. (2014) offer an analysis that largely applies to DGPs though they also study the Cho and Saul recursion with the motivating argument from the original paper.

2. THE DEEP WIDE LIMIT:

We consider a fully connected network as shown in Figure 1 . The inputs and outputs will be real valued vectors of dimension M and L respectively. The network is fully connected. The initial step and recursion are standard. The initial step is:

f (1) i (x) = M j=1 w (1) i,j x j + b (1) i .

(1) Figure 1 : In this paper we consider fully connected feedforward networks with more than one hidden layer. We call the pre-nonlinearity an activation and post-nonlinearity an activity. As the network becomes increasingly wide the distribution of the marginal distributions of the activations at each layer and of the output will become close to a Gaussian process in a sense described in the text.

We make the functional dependence on x explicit in our notation as it will help clarify what follows.

For a network with D hidden layers the recursion is, for each µ = 1, . . . , D, g

i (x) = φ(f

EQUATION

EQUATION

so that f (D+1) (x) is the output of the network given input x. φ denotes the non-linearity. In all cases the equations hold for each value of i; i ranges between 1 and H µ in Equation (2), and between 1 and H µ+1 in Equation (3) except in the case of the final activation where the top value is L. The network could of course be modified to be probability simplex-valued by adding a softmax at the end.

A distribution on the parameters of the network will be assumed. Conditional on the inputs, this induces a distribution on the activations and activities. In particular we will assume independent normal distributions on the weights and biases

EQUATION

)

We will be interested in the behaviour of this network as the widths H µ becomes large. The weight variances for µ ≥ 2 will be scaled according to the width of the network to avoid a divergence in the variance of the activities in this limit. As will become apparent, the appropriate scaling is

EQUATION

)

The assumption is that Ĉ(µ) w will remain fixed as we take the limit. Neal (1996) analysed this problem for D = 1, showing that as H 1 → ∞, the values of f

(2) i (x), the output of the network in this case, converge to a certain multi-output Gaussian process if the activities have bounded variance.

Since our approach relies on the multivariate central limit theorem we will arrange the relevant terms into (column) vectors to make the linear algebra clearer. Consider any two inputs x and x and all output functions ranging over the index i. We define the vector f (2) (x) of length L whose elements are the numbers f

(2) i (x). We define f (2) (x ) similarly. For the weight matrices defined by w (µ) i,j for fixed µ we use a 'placeholder' index • to return column and row vectors from the weight matrices. In particular w

(1) j,• denotes row j of the weight matrix at depth 1. Similarly, w

(2)

•,j denotes column j at depth 2. The biases are given as column vectors b (1) and b (2) . Finally we concatenate the two vectors f (2) (x) and f (2) (x ) into a single column vector F (2) of size 2L. The vector in question takes the form

EQUATION

•,j φ(w

(1)

j,• x + b (1) j ) w (2) •,j φ(w (1) j,• x + b (1) j ) (7)

The benefit of writing the relation in this form is that the applicability of the multivariate central limit theorem is immediately apparent. Each of the vector terms on this right hand side is independent and identically distributed conditional on the inputs x and x . By assumption, the activities have bounded variance. The scaling we have chosen on the variances is precisely that required to ensure the applicability of the theorem. Therefore as H becomes large F (2) converges in distribution to a multivariate normal distribution. The limiting normal distribution is fully specified by its first two moments. Defining γ ∼ N (0,

EQUATION

w I M ), the moments in question are:

EQUATION

Note that we could have taken a larger set of input points to give a larger vector F and again we would conclude that this vector converged in distribution to a multivariate normal distribution. More formally, we can consider the set of possible inputs as an index set. A set of consistent finite dimensional Gaussian distributions on an index set corresponds to a Gaussian process by the Kolmogorov extension theorem. The Gaussian process in question is a distribution over functions defined on the product σ-algebra, which has the relevant finite dimensional distributions as its marginals.

In the case of a multivariate normal distribution a set of variables having a covariance of zero implies that the variables are mutually independent. Looking at Equation ( 9), we see that the limiting distribution has independence between different components i, j of the output. Combining this with the recursion (2), we might intuitively suggest that the next layer also converges to a multivariate normal distribution in the limit of large H µ . Indeed we state the following lemma, which we attribute to Hazan & Jaakkola (2015) :

Lemma 1 (Normal recursion). If the activations of a previous layer are normally distributed with moments:

EQUATION

Then under the recursion (2) and as H → ∞ the activations of the next layer converge in distribution to a normal distribution with moments

EQUATION

)

where K is a 2 × 2 matrix containing the input covariances.

Unfortunately the lemma is not sufficient to show that the joint distribution of the activations of higher layers converge in distribution to a multivariate normals. This is because for finite H the input activations do not have a multivariate normal distribution -this is only attained (weakly or in distribution) in the limit. It could be the case that the rate at which the limit distribution is attained affects the distribution in subsequent layers. We are able to offer the following theorem rigorously:

Theorem 1. Consider a Bayesian deep neural network of the form in Equations (1) and (2) using ReLU activation functions. Then there exist strictly increasing width functions

h µ : N → N such that H 1 = h 1 (n), . . . , H D = h D (n)

, and for any countable input set (x[i]) ∞ i=1 , the distribution of the output of the network converges in distribution to a Gaussian process as n → ∞.

A proof is included in the appendix. We conjecture that a more general theorem will hold. In particular we expect that the width functions h µ can be taken to be the identity and that the nonlinearity can be extended to monotone functions with well behaved tails. Our conjecture is based on the intuition from Lemma 1 and from our experiments, in which we always take the width function to be the identity.

3 SPECIFIC KERNELS UNDER RECURSION Cho & Saul (2009) suggest a family of kernels based on a recurrence designed to 'mimic computation in large, multilayer neural nets'. It is therefore of interest to see how this relates to deep wide Gaussian processes. A kernel may be associated with a feature mapping Φ(x) such that K(x, x ) = Φ(x) • Φ(x ). Cho and Saul define a recursive kernel through a new feature mapping by compositions such as Φ(Φ(x)). However this cannot be a legitimate way to create a kernel because such a composition represents a type error. There is no reason to think the output dimension of the function Φ matches the input dimension and indeed the output dimension may well be infinite. Nevertheless, the paper provides an elegant solution to a different task: it derives closed form solution to the recursion from Lemma 1 (Hazan & Jaakkola, 2015) for the special case

EQUATION

)

where Θ is the Heaviside step function. Specifically, the recursive approach of Cho & Saul (2009) can be adapted by using the fact that u z for z ∼ N (0, LL ) is equivalent in distribution to (L u) ε with ε ∼ N (0, I), and by optionally augmenting u to incorporate the bias. Since r = 1 corresponds to rectified linear units we apply this analytic kernel recursion in all of our experiments.

4. MEASURING CONVERGENCE USING MAXIMUM MEAN DISCREPANCY:

In this section we use the kernel based two sample tests of Gretton et al. (2012) to empirically measure the similarity of finite random neural networks to their Gaussian process analogues. The maximum mean discrepancy (MMD) between two distributions P and Q is defined as:

EQUATION

where H denotes a reproducing kernel Hilbert space and || • || H denotes the corresponding norm. It gives the biggest possible difference between expectations of a function under the two distributions under the constraint that the function has Hilbert space norm less than or equal to one. We used the unbiased estimator of squared MMD given in Equation (3) of Gretton et al. (2012) .

In this experiment and all those that follow we take weight variance parameters Ĉ(µ) w = 0.8 and bias variance C b = 0.2. We took 10 standard normal input points in 4 dimensions and pass them through 2000 independent random neural networks drawn from the distribution discussed in this paper. This was then compared to 2000 samples drawn from the corresponding Gaussian process distribution. The experiment was performed with different numbers of hidden layers and numbers of units per hidden layer. We repeated each experiment 20 times which allows us to reduce variance in our results and give a simple estimate of measurement error. The experiments use an RBF kernel for the MMD estimate with lengthscale 1/2. In order to help give an intuitive sense of the distances involved we also include a comparison between two Gaussian processes with isotropic RBF kernels using the same MMD distance measure. The kernel length scales for this pair of 'calibration' Gaussian processes are taken to be l and 2l, where the characteristic length scale l = √ 8 is chosen to be sensible for the standard Normal input distribution on the four dimensional space.

The results of the experiment are shown in Figure 2 . We see that for each fixed depth the network converges towards the corresponding Gaussian process as the width increases. For the same number of hidden units per layer, the MMD distance between the networks and their Gaussian process analogue becomes higher as depth increases. The rate of convergence to the Gaussian process is slower as the number of hidden layers is increased.

5. COMPARING BAYESIAN DEEP NETWORKS TO GAUSSIAN PROCESSES:

In this section we compare the behaviour of finite Bayesian deep networks of the form considered in this paper with their Gaussian process analogues. If we make the networks wide enough the agreement will be very close. It is also of interest, however, to consider the behaviour of networks actually used in the literature, so we use 3 hidden layers and 50 hidden units which is typical of the networks used by Hernández-Lobato & Adams (2015) . Fully connected Bayesian deep networks with finite variance priors on the weights have also been considered in other works (Graves, 2011; Hernández-Lobato et al., 2016; Blundell et al., 2015) , though the specific details vary. We use rectified linear units and correct the variances to avoid a loss of prior variance as depth is increased. Our general strategy was to compare exact Gaussian process inference against expensive 'gold standard' Markov Chain Monte Carlo (MCMC) methods. We choose the latter because used correctly it works well enough to largely remove questions of posterior approximation quality from the calculus of comparison. It does mean however that our empirical study does not extend to datasets which are large in terms of number of data points or dimensionality, where such inference is challenging. We therefore sound a note of caution about extrapolating our empirical finite network conclusions too confidently to this domain. On the other hand, lower dimensional, prior dominated problems are generally regarded as an area of strength for Bayesian approaches and in this context our results are directly relevant.

We computed the posterior moments by the two different methods on some example datasets. For the MCMC we used Hamiltonian Monte Carlo (HMC) (Neal, 2010) updates interleaved with elliptical slice sampling (Murray et al., 2010) . We considered a simple one dimensional problem and a two dimensional real valued embedding of the four data point XOR problem. We see in Figures 3 and 4 (left) that the agreement in the posterior moments between the Gaussian process and the Bayesian deep network is very close.

A key quantity of interest in Bayesian machine learning is the marginal likelihood. It is the normalising constant of the posterior distribution and gives a measure of the model fit to the data. For a Bayesian neural network, it is generally very difficult to compute, but with care and computational time it can be approximated using Hamiltonian annealed importance sampling (Sohl-Dickstein & Culpepper, 2012) . The log-importance weights attained in this way constitute a stochastic lower bound on the marginal likelihood (Grosse et al., 2015) . Figure 4 (right) shows the result of such an experiment compared against the (extremely cheap) Gaussian process marginal likelihood computation on the XOR problem. The value of the log-marginal likelihood computed in the two different ways agree to within a single nat which is negligible from a model selection perspective (Grosse et al., 2015) .

Predictive log-likelihood is a measure of the quality of probabilistic predictions given by a Bayesian regression method on a test point. To compare the two models we sampled 10 standard normal train and test points in 4 dimensions and passed them through a random network of the type under study to get regression targets. We then discarded the true network parameters and compared the predictions of posterior inference between the two methods. We also compared the marginal predictive distributions of a latent function value. Figure 5 shows the results. We see that the correspondence in predictive log-likelihood is close but not exact. Similarly the marginal function values are close to those of a Gaussian process but are slightly more concentrated. 

6. AVOIDING GAUSSIAN PROCESS BEHAVIOUR:

When using deep Bayesian neural networks as priors, the emergence of Gaussian priors raises important questions in the cases where it is applicable, even if one sets aside questions of computational tractability. It has been argued in the literature that there are important cases where kernel machines with local kernels will perform badly (Bengio et al., 2005) . The analysis applies to the posterior mean of a Gaussian process. The emergent kernels in our case are hyperparameter free. Although they do not meet the strict definition of what could be considered 'local' the fact remains that any Gaussian process with a fixed kernel does not use a learnt hierarchical representation. Such representations are widely regarded to be essential to the success of deep learning. There is relevant literature here on learning the representation of a standard, usually structured, network composed with a Gaussian process (Wilson et al., 2016a; b; Al-Shedivat et al., 2017) . This differs from the assumed paradigm of this paper, where all model complexity is specified probabilistically and we do not assume convolutional, recurrent or other problem specific structure. Within this paradigm, the question therefore arises as to what can be done to avoid marginal Gaussian process behaviour if it is not desired. Speaking loosely, to stop the onset of the central limit theorem and the approximate analogues discussed in this paper one needs to make sure that one or more of its conditions is far from being met. Since the chief conditions on the summands are independence, bounded variance and many terms, violating these assumptions will remove Gaussian process behaviour. Deep Gaussian processes (Damianou & Lawrence, 2013) are not close to standard Gaussian processes marginally because they are typically used with narrow intermediate layers. It can be challenging to choose the precise nature of these narrow layers a priori. Neal (1996) suggests using networks with infinite variance in the activities. With a single hidden layer and correctly scaled, these networks become alpha stable processes in the wide limit. Neal also discusses variants that destroy independence by coupling weights. Our results about the emergence of Gaussian processes even with more than one hidden layer mean these ideas are of considerable interest going forward.

7. CONCLUSIONS:

Studying the limiting behaviour of distributions on feedforward networks has been a fruitful avenue for understanding these models historically. In this paper we have extended the state of knowledge about the wide limit, including for networks with more than one hidden layer. In particular, we have exhibited limit sequences of networks that converge in distribution to Gaussian processes with a certain recursively defined kernel. Our empirical study using MMD suggests that this behaviour is exhibited in a variety of models of size comparable to networks used in the literature. This led us to juxtapose finite Bayesian neural networks with their Gaussian process analogues, finding that the agreement in terms of key predictors is close empirically. If this Gaussian process behaviour is desired then exact and approximate inference using the analytic properties of Gaussian processes should be considered as an alternative to neural network inference. Since Gaussian processes have an equivalent flat representation then in the context of deep learning the behaviour may well not be desired and steps should be taken to avoid it.

We view these results as a new opportunity to further the understanding of neural networks in the work that follows. Initialisation and learning dynamics are crucial topics of study in modern deep learning which require that we understand random networks. Bayesian neural networks should offer a principled approach to generalisation but this relies on successfully approximating a clearly understood prior. In illustrating the continued importance of Gaussian processes as limit distributions, we hope that our results will further research in these broader areas.

A PROOF OF MAIN THEOREM A.1 STATEMENT OF THEOREM AND NOTATION:

In this section, we provide a proof of the main theorem of the paper, which we begin by recalling.

Theorem 1. Consider a Bayesian deep neural network of the form in Equations (1) and (2) using ReLU activation functions. Then there exist strictly increasing width functions h µ : N → N such that H 1 = h 1 (n), . . . , H D = h D (n), and for any countable input set (x[i]) ∞ i=1 the distribution of the output of the network converges in distribution to a Gaussian process as n → ∞.

The theorem is proven via use of the propositions that follow below. The broad structure of the proof is to use a particular variant of the Berry-Esseen inequality to upper bound how far each layer is from a multivariate normal distribution, and then to inductively propagate these inequalities through the network, leading to a bound on the distance between the output of the network for a collection of input points, and a multivariate Gaussian distribution. These notions will be made precise below. We begin in Section A.2 by stating the propositions that will be used in the proof of Theorem 1, but first establish notation that will be used in the remainder of the appendix.

Given a finite set of inputs x[1], . . . , x[n] ∈ R M , we will write:

• f (µ) (x) for the random variables (f (µ) (x[i])) n i=1 collectively taking values in R nHµ ; • f (µ) j (x) for the random variables (f (µ) j (x[i])) n i=1 collectively taking values in R n ; • g (µ) (x) for the random variables (g (µ) (x[i])) n i=1 collectively taking values in R nHµ ; • g (µ) j (x) for the random variables (g (µ) j (x[i])) n i=1 collectively taking values in R n ;

Throughout, if U 1 , U 2 are random variables taking in values in some Euclidean space R d , we will define d(U 1 , U 2 ) = sup

A⊆R d A convex |P(U 1 ∈ A) -P(U 2 ∈ A)| .

Note that convergence of a sequence of random variables in this metric implies convergence in distribution.

We will also consider multivariate normal distributions (Z (µ) j (x[i])|j = 1, . . . , H µ , i = 1, . . . , n) with covariance matrices of block diagonal form, such that

Cov(Z (µ) k (x[a]), Z (µ) l (x[b])) = 0 for distinct k, l ∈ {1, . . . , H µ } , for all x[a], x[b] .

To avoid writing this in full every time it is required, we will refer to this condition as blockwise independence with respect to the index j. We will avoid specification of all covariance values, deferring to the expression (13) given in the main paper. Finally, to simplify notation, we will assume that the network output is one-dimensional. Our proof trivially extends to arbitrary finite output dimension where the limiting distribution is a coordinate-wise independent multivariate GP. 

(f (µ-1) (x), Z (µ-1) (x)) ≤ 2 -((D+1)-(µ-1))-n D k=µ-1 H k ε , holds, where Z (µ-1) (x) = (Z (µ-1) j (x[i])|j = 1, . . . , H µ-1 , i = 1, . . . , n

) is mean-zero multivariate normal, with blockwise independence with respect to the index j. Then we have

d(f (µ) (x), Z (µ) (x)) ≤ 2 -((D+1)-µ)-n D k=µ H k ε , where Z (µ) (x) = (Z (µ) j (x[i])|j = 1, . . . , H µ , i = 1, . . . , n

) is mean-zero multivariate normal, with blockwise independence with respect to the index j.

We need a mildly modified version of this theorem to deal with iid random vectors X 1 , . . . , X n with non-identity covariance matrices. To this end, suppose that Σ is the (full-rank) covariance matrix of each X i , with decomposition Σ = RR , for some invertible matrix R; R can be obtained, for example, by using Cholesky or Schur decomposition. The random variables R -1 X 1 , . . . , R -1 X n are then iid, mean zero and with identity covariance matrices, so we may apply Theorem 2 to obtain:

sup A⊆R d A convex |P(R -1 S n ∈ A) -P(Y ∈ A)| ≤ 400d 1/4 β √ n , where β = E R -1 X i 3 . Now note that this is equivalent to sup A⊆R d A convex |P(S n ∈ RA) -P(RY ∈ RA)| ≤ 400d 1/4 β √ n , noting that RY ∼ N (0, Σ).

Since R is invertible, and recalling the definition of the distance d above, this is exactly equivalent to:

EQUATION

which is the variant of Bentkus' result we will require in the sequel.

We apply this bound to the sum

Hµ-1 j=1 g (µ-1) j (x) ⊗ w (µ) •,j

Noting that the summands indexed by j are iid by assumption, with the expected third moment norm featuring in the Berry-Esseen inequality upper-bounded by β ≤ C Hµ,n , for some constant C Hµ,n depending on H µ and n, but independent of H µ-1 (finiteness of C Hµ,n follows from Lemma 4).

As a consequence, we have the following bound:

d   Hµ-1 j=1 g (µ-1) j (x) ⊗ w (µ) •,j , Z (x)   ≤ 400C Hµ,n (nH µ ) 1/4 / H µ-1 ,

where Z (x) = (Z j (x[i])|j = 1, . . . , H µ , i = 1, . . . , n) is mean-zero multivariate normal, with blockwise independence with respect to the index j. We wish to demonstrate that this is less than or equal to 2 -(D-(µ-1))-n D k=µ H k ε when H D is sufficiently large. This is equivalent to showing that where Z(x) is mean-zero multivariate normal, with the same block-diagonal covariance structure as described for Z (x) above, and 1 n ∈ R n is a vector of 1's.

EQUATION

Proof of Lemma 3. Let A ⊆ R nHµ be an arbitrary convex set. First, observe that we have    is exactly that the vector φ( f (µ-1) (x)) lies in the preimage of the convex set A -1 n ⊗ b (µ) under the linear map w (µ) , which is again a convex set. Secondly, observe that for the specific ReLU nonlinearity φ, if C is an arbitrary convex set, then {(f (µ-1) (x)|φ(f (µ-1) (x)) ∈ C} may be written as the disjoint union of at most 2 nHµ-1 convex sets: Applying the assumed bound in the statement of the lemma to each of these sets, we obtain |P(φ(f (µ-1) (x)) ∈ C) -P(φ( f (µ-1) (x)) ∈ C)| ≤ 2 nHµ-1 ε .

{f (µ-1) (x)|φ(f (µ-1) (x)) ∈ C} = (φ -1 (C) ∩ {t ∈ R nHµ-1 |t i ≥ 0 ∀i})∪

Substituting this bound into the conditional probability (17) yields |P(f (µ) (x) ∈ A) -P( f (µ) (x) ∈ A)| ≤ 2 nHµ-1 ε .

Since A was an arbitrary convex set, the proof is complete.

Proof of Lemma 4. Note that by independence of g (µ-1) j (x) from w (µ)

•,j

we have that each X j has mean zero and covariance Σ ⊗ = Σ ⊗ Ĉ(µ) w I where Σ is the covariance matrix of g (µ-1) j (x). By standard properties of the Kronecker product, the Schur decomposition of Σ ⊗ is (QΛQ T )⊗( Ĉ(µ) w I) where QΛQ T is the Schur decomposition of Σ. Simple algebraic manipulation yields:

E Q ⊗ Λ -1/2 ⊗ Q T ⊗ ( g (µ-1) j (x) ⊗ w (µ)

•,j ) 3 = E (QΛ -1/2 Q T g (µ-1) j (x)) ⊗ (( Ĉ(µ) w ) -1/2 w (µ)

•,j ) 3 = E QΛ -1/2 Q T g (µ-1) j (x) 3 E ( Ĉ(µ) w ) -1/2 w (µ)

•,j

3 .

Notice that the random variable ( Ĉ(µ) w ) -1/2 w (µ)

•,j

follows the R Hµ -dimensional standard normal distribution, and thus its squared norm follows the chi-squared distribution with H µ degrees of freedom, which is also known as the Gamma(H µ /2, 1/2) distribution. Exponentiating to the power of 3/2 and taking the expectation, we obtain:

E ( Ĉ(µ)

w ) -1/2 w (µ)

•,j 3 = 2 3/2 Γ((H µ + 3)/2) Γ(H µ /2) .

:

Code for the experiments in the paper can be found at https://github.com/ widedeepnetworks/widedeepnetworks

Figures and Tables:

Figure2: A comparison of finite random neural networks to their corresponding Gaussian process analogue using an (RBF) kernel estimator of the squared maximum mean discrepancy (MMD). The results are consistent with the emergence of Gaussian process behaviour as the networks become wide. The red dashed line is for calibration and denotes the squared MMD between two Gaussian processes with isotropic RBF kernels and length scales l and 2l where l = √ 8 is the characteristic length scale of the input space (see text).

Figure 3: A comparison between Bayesian posterior inference in a Bayesian deep neural network and posterior inference in the analogous Gaussian process. The neural network has 3 hidden layers and 50 units per layer. The lines show the posterior mean and two σ credible intervals.

Figure 4: A comparison between posterior inference for a Gaussian process and a Bayesian deep network for a real value embedding of the XOR function. Left and centre: The two posterior means. The mean absolute different between the two posterior estimate grids is 0.027. Right: Kernel density estimate of the log weights from annealed importance sampling on a Bayesian deep network compared to the analogous Gaussian process marginal likelihood shown by the vertical line. The neural network has 3 hidden layers and 50 units per layer.

SUPPORTING RESULTS Proposition 1. Let ε > 0, and x[1], . . . , x[n] ∈ R M . Let µ ∈ {2, . . . , D + 1}, and let H k = 2 H 2 k+1 for k = 1, . . . , D -1. Then for H D sufficiently large, suppose the condition d

1, the left-hand side converges to 0 as H D increases (using the bound obtained for C Hµ,n in Lemma 4), so for all H D sufficiently large, we obtaind , Z(x)   ≤ 2 -((D+1)-(µ-1))-n D k=µ H k ε ,

-1 n ⊗ b (µ) w (µ) , b (µ) that for fixed w (µ) and b (µ) -1 n ⊗ b (µ) w (µ) , b(µ)

nHµ-1 |t I < 0, ∃y ∈ C s.t. y I c = t I c , y I = 0} .

Peer review (iclr2018):
 In part 1, the authors introduce motivation for studying wide neural networks and summarize related work. In part 2, they present a theorem (main theoretical result) stating that under conditions on the weight priors, the output function of a multi-layer neural network (conditionally to a given input) weakly converges to a gaussian process as the size of the hidden layers go to infinity. remark on theorem 1: This result generalizes a result proven in 2015 stating that the normality of a layer propagates to the next as the size of the first layer goes to infinity. The result stated in this paper is proven by bounding the gap between the output distribution and the corresponding gaussian process, and by propagating this bound across layers (appendix). In part 3, the authors discuss the choice of a nonlinearity function that enables easy computation of the kernels introduced in the covariance matrix of the limit normal distribution. Their choice lands on ReLU. In part 4, the focus is on the speed of the convergence presented in theorem 1. Experiments are conducted to show how the distance (maximum mean disrepancy) between the output distribution and its theoretical gaussian process limit vary when the sizes of the hidden layers increase. The results show that the convergence (in MMD) happens consistently, although it is slower when the number of hidden layers gets bigger. In part 5, the authors compare the distributions (finite Bayesian deep networks and their analogues Gaussian processes) in yet another way: by studying their agreement in terms of inference. For this purpose, the authors chose several crieteria: the first two moments of the posterior, the log marginal likelihood and the predictive log-likelihood. The authors judge that the distributions agree on those criteria, but do not provide further analysis. In part 6, now that It has been shown that the output distributions of Bayesian neural nets do not only weakly converge to Gaussian processes but also behave similarly in terms of inference, the authors discuss ways to avoid the gaussian process behaviour. Indeed, it seems that Gaussian processes with a fixed kernel cannot learn hierarchical representations, that are essential in deep learning. The idea to avoid the Gaussian process behaviour is to contradict one of the hypothesis of the CLT (so that it does not hold anymore), either by controlling the size of intermediate layers, by using networks with infinite variance in the activities, or by choosing non-independent weights. In part 7, it is concluded that the result that has been proven for size of layers going to infinity (Theorem 1) seems to empirically be verified on finite networks similar to those used in the literature. This can be used to simplify inference in cases were the gaussian process behaviour is desired, and opens questions on how to avoid this behaviour the rest of the time. Pros: The authors line of thought of the authors is overall quite easy to follow. The main theoretical convergence result is stated early on, and the remaining of the article is dedicated to observing this result empirically from different angles (MMD, inference, predictive capability..). The last part contains a discussion concerning the extent to which it is actually a desired or a undesired result in classical deep learning use-cases, and the authors provide intuitive conditions under which the convergence would not hold. The stated theorem is a clear improvement on the past literature and is promising in a context where multi-layers neural networks are more and more studied. Finally, the work is well documented. Cons: I have a some concerns with the main result (Theorem 1) and found that some of the notations / formulas were not very clear. Concerns with Theorem 1: * at the end of the proof of Lemma 2, H_.mu is to be chosen large enough in order to get the .epsilon bound of the statement. However, I think that H_.mu is constrained by the statement of Proposition 2, not to be larger than a constant times 2^(H_{.mu+1}). Isn*t that a problem? * In the proof of Lemma 4, it looks like matrix .Psi, from the schur decomposition of .tilde f, actually depends on H_{.mu-2}, thus making .psi_max depend on it too, as well as the final .beta bound, which would contradict the statement that it depends only on n and H_{.mu}. Could you please double check? Unclear statements/notations: * end of page 3, notations are not entirely consist with previous notations * I do not understand which distribution is assumed on epsilon and gamma when taking the expectancy in equation (9). * the notation x^(i) (in the theorem and the proof notably) could be changed, for the ^(i) index refers to the depth of the layer in the rest of the notations, and is here surprisingly referring to a set of observations. * the statement of Theorem 1: * I would change "for a countable input set" to "for any countable input set", if this holds true. * does not say that the width has to go to infinity for the convergence to happen, which goes a bit in contradiction with the adjective "wide". However, the authors say that in practice, they use the identity as width function. * I understood that the conclusion of part 3 was that the expectation of eq (9) was elegantly computable for certain non-linearity (including ReLU). However I don*t see the link with the "recursive kernel" idea (maybe it*s just the way to do the computation described in Cho&Saul(2009) ?) Some places where it appears that there are minor mistakes: * 7th line from the bottom of page 3, the vector f^{(2)}(x) contains f_i^{(1)}(x) but should contain f_i^{(2)}(x) * last display of page 3: change x and x*, and indicate upper limit of the sum * please double check variances C_w and/or .hat{C}_w appearing in equations in (9) and (13). * line 2 of second paragraph after equations (8) and (9). The authors refer to equation (8) concerning the independence of the components of the output. I think they rather wanted to refer to (9). Same for first sentence before eq (14). * middle of page 12: matrix LY should be RY.