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):
 - Summary The paper is well written and proves how deep, wide, fully connected NNs are equivalent to GPs in the limit. This result, which was well known for single-layer NNs, is now extended to the multilayer case. Although there was already previous work suggesting GP this behavior, there was no formal proof under the specific conditions presented here. The convergence to a GP is also verified experimentally on some toy examples. - Relevance The result itself does not feel very novel because variants of it were already available. Unfortunately, although making other researchers aware of this is worthy, the application of this result seems limited, since in fact it describes and lets us know more about a regime that we would rather avoid, rather than one we want to exploit. Most of the applications of deep learning benefit from strong structured priors that cannot be represented as a GP. This is properly acknowledged in the paper. The lack of practical relevance combined with the not-groundbreaking novelty of the result makes this paper less appealing. - Other comments Page 6: "It does mean however that our empirical study does not extend to larger datasets where such inference is prohibitively expensive (...) prior dominated problems are generally regarded as an area of strength for Bayesian approaches and in this context our results are directly relevant." Although that argument can hold for datasets that are large in terms of amount of data points, it doesn*t for datasets that are large in terms of number of dimensions. The empirical study could have used very high-dimensional datasets with comparatively low amounts of training data. That would maintain a regime were the prior does matter but and better show the generality of the results. Page 6: "We use rectified linear units and correct the variances to avoid a loss of prior variance as depth is increased as discussed in Section 3" Are you sure this is discussed in Section 3? Page 4: "This is because for finite H the input activations do not have a multivariate normal distribution". Can you elaborate on this? Since we are interested in the infinite limit, why is this a problem?