1 Introduction
This work substantially extends the work of Matthews et al. (2018)
published at ICLR 2018. 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 this 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 the 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, finite 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 nontrivial and part of our discussion. Historically this result was significant because it provided a connection between flexible Bayesian neural networks and Gaussian processes (Williams, 1998; Rasmussen and 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 2.4
) on the convergence of certain sequences of finite fully connected networks with more than one hidden layer to Gaussian processes. The number of hidden layers can be any fixed number. The sizes of the hidden layers must strictly increase for each network in the sequence although the different hidden layers are allowed to grow at different rates. The weights are assumed to be independent normally distributed with their variances sensibly scaled as the network grows following the prescription of
Neal (1996). The nonlinearities are assumed to obey the ‘linear envelope’ Condition 2.3 which all commonly used nonlinearities do in fact obey. Since these are the only assumptions on the sequence of networks it will be seen that the result is a meaningfully general one.Further, we empirically study the distance between finite networks and their Gaussian process analogues by using maximum mean discrepancy (MMD, Gretton et al., 2012)
as a distance measure. We then systematically compare exact Gaussian process inference with ‘gold standard’ MCMC inference for finite Bayesian neural networks. Of the six datasets we consider, five show close agreement between the two models. Owing to the computational burden of the MCMC algorithms, the problems we can study by this method are constrained in terms of their network size, the data dimensionality and the number of data points. Nevertheless our results suggest that some experiments in the literature studied under the banner of Bayesian deep learning would have given very similar results to a Gaussian process with the appropriate kernel. A practical recommendation following from our study is that the Bayesian deep learning community should routinely compare their results to Gaussian processes with the kernels studied in this paper.
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 addition to standard techniques for inference in Bayesian deep learning. In some scenarios, the behaviour may not be desired because it implies a lack of a hierarchical representation and a Gaussian statistical assumption. 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 and 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 and Jaakkola (2015). Their contribution is similar in content to Lemma 2.3 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 selfnormalizing 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 and Lawrence, 2013). As will be discussed in Section 7, narrow intermediate representations mean that the marginal behaviour of DGPs 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.Simultaneously with the submission of the previous version of our paper to ICLR 2018, at the same conference venue, Lee et al. (2018)
released a paper that has overlap with our own. There are however some important differences. Empirically, whilst we compare finite Bayesian neural networks, using ‘gold standard’, asymptotically exact, sampling and MMD, to their Gaussian process analogues, Lee et al. compare finite neural networks trained with stochastic gradient descent (SGD) to Gaussian processes instead. The latter comparison to SGD is suggestive that this optimization method mimics Bayesian inference – an idea that has been receiving increasing attention
(Welling and Teh, 2011; Mandt et al., 2017; Smith and Le, 2018). This is of particular importance because typically SGD is still more scalable than traditional Markov Chain based methods, enabling Lee et al. to consider some relatively large datasets. The empirical comparison of Lee et al. is therefore particularly intriguing and we hope it will lead to follow up work. Therefore, whilst there is overlap, the two papers also have independent value going forward. Theoretically, there are also important differences. Lee et al. give an argument for the Gaussian process limit, although importantly this depends on sequentially taking the number of units in each successive layer to be infinite. The proof we give here concerns the case where the layers grow simultaneously, which is arguably more relevant in practice. Note that we show a precise type of convergence – namely ‘weak convergence’ also called ‘convergence in distribution’. Owing to the challenging nature of obtaining a full rigorous proof, the earlier version of this paper
(Matthews et al., 2018) did not achieve full generality either. We needed to assume specific growth rates for the sizes of the hidden layers, and the ReLU nonlinearity. What follows here removes these assumptions and thus resolves the conjecture made in the earlier version of this work in the affirmative. The new proof method, placing a particular emphasis on exchangeability, may well be of use more generally.2 The deep wide limit
2.1 The result for one hidden layer
We consider a fully connected network as shown in Figure 1
. The inputs and outputs will be realvalued vectors of dimension
and respectively. The network is fully connected. The initial step and recursion are standard. The initial step is:(1) 
We make the functional dependence on explicit in our notation as it will help clarify what follows. For a network with hidden layers the recursion is, for each ,
(2)  
(3) 
so that is the output of the network given input . denotes the nonlinearity. In all cases the equations hold for each value of ; ranges between and in Equation (2), and between and in Equation (3) except in the case of the final activation where the top value is . The network could of course be modified to be probability simplexvalued 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
(4)  
(5) 
We will be interested in the behaviour of this network as the widths becomes large. The weight variances for 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
(6) 
The assumption is that will remain fixed as we take the limit. Neal (1996) analysed this problem for , showing that as , the values of , the output of the network in this case, converge to a certain multioutput 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 and and all output functions ranging over the index . We define the vector of length whose elements are the numbers . We define similarly. For the weight matrices defined by for fixed we use a ‘placeholder’ index to return column and row vectors from the weight matrices. In particular denotes row of the weight matrix at depth . Similarly, denotes column at depth . The biases are given as column vectors and . Finally we concatenate the two vectors and into a single column vector of size . The vector in question takes the form
(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 and . 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, and is also in line with most commonly used initialisation strategies in practice. Therefore as becomes large converges in distribution to a multivariate normal distribution. The limiting normal distribution is fully specified by its first two moments. Defining , the moments in question are:
(8)  
(9) 
Note that we could have taken a larger set of input points to give a larger vector 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. What we have shown is that for any finite index set the distribution over functions converges to a multivariate normal. If we consider these limiting multivariate normals they obey a consistency property under marginalization. This means that the limiting distributions can be used to define a Gaussian process by the Kolmogorov extension theorem.
2.2 Definition of weak convergence of random functions
There are some important technical issues here that are not discussed in the original work of Neal (1996). In some sense, the convergence of the finitedimensional distributions is enough if we want to answer questions about finite events, just as many of the uses of Gaussian processes within machine learning (Rasmussen and Williams, 2006) can be expressed in terms of finitedimensional multivariate normal distributions. The reader who is content with restricting their attention to such a case may safely omit the rest of this subsection.
Given a consistent set of finitedimensional marginals, the Kolmogorov extension theorem ensures the existence of an underlying infinitedimensional object – a distribution over functions. If we want to make precise mathematical statements about convergence to this object some care is needed.
Firstly, the Kolmogorov theorem ensures the existence of a distribution which is uniquely defined on a specific algebra, namely the product algebra. The algebra defines which events we can assign probabilities to. If we try to consider events outside the algebra then the rules governing probability distributions (c.f. measures) can break down. Secondly, in abstract spaces, the definition of convergence in distribution is necessarily with respect to some topology. In everything that follows we will assume that this topology is generated by a metric. We also assume that the index set of the stochastic process is countably infinite. We use the metric :
(10) 
This metric metrises the product topology of the product of countably many copies of with the usual Euclidean topology (Dashti and Stuart, 2013). For such a countable index set, it is sufficient (Billingsley, 1999, p. 19)
to prove weak convergence of the finitedimensional marginals of the process to the corresponding multivariate Gaussian random variables. This is not generally the case if we remove the assumption of a countable index set
(Billingsley, 1999, p. 19).The restriction to countably infinite index sets means that phenomena that depend on uncountably many indices such as continuity, boundedness and differentiability are not covered by our theory. There is literature extending measures on the product algebra of an uncountable index set using, for instance, the Kolmogorov continuity theorem. One could then consider proving convergence with respect to the topology in question. We do not do this in this paper but it could certainly be of interest.
2.3 The recursion lemma and the linear envelope property
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 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 .
This will indeed be the case assuming that the nonlinearity does not induce heavy tail behaviour. We give an assumption on the nonlinearity that will be used throughout the sequel:
[Linear envelope property for nonlinearities] A nonlinearity is said to obey the the linear envelope property if there exist such that the following inequality holds
(11) 
The majority of commonly used nonlinearities, including the sigmoid, ReLU, ELU, and SeLU nonlinearities have the linear envelope property. Intuitively the linear bounds on the nonlinearity stop it from inducing heavy tail behaviour when a random variable is passed through it. An exponential nonlinearity would not have this property. We could indeed craft a nonlinearity that is designed to violate the linear envelope property and give heavy tail behaviour. Consider, for example, the composition of the Gaussian cumulative density function (CDF) followed by the Cauchy inverse CDF. Passing a standard normal variate through such a function would, by construction, give a Cauchy distributed variable, which has an undefined mean. Whilst it may not be the most general assumption possible for what will follow, the linear envelope assumption rules in most practically used nonlinearities and, as we shall see rules out all nonlinearities for which our theory does not hold.
Next we state the following lemma, which we attribute to Hazan and Jaakkola (2015):
[Normal recursion] If the activations of a previous layer are normally distributed with moments:
(12)  
(13) 
Then under the recursion (2) and as the activations of the next layer converge in distribution to a normal distribution with moments
(14)  
(15) 
where is a 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 normal. This is because for finite
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.Therefore the proof of our main result will require considerably more technical machinery then would be suggested by the recursion in Lemma 2.3. We discuss the more general result in the next section.
2.4 Convergence for more than one hidden layer
In order to state our theorem we will need one more definition, namely that of a width function:
[Width functions] For a given fixed input , a width function at depth specifies the number of hidden units at depth .
For a given fixed input , the set of width functions together fully specify a shape for a fully connected network. In this way, the countable sequence of natural numbers specifies a countable sequence of fully connected networks. We will be interested in the case where each of the width functions tends to infinity. Note that this includes the case of taking the width functions to be the identity, which gives the case where each hidden layer has the same number of hidden units and tends jointly to infinity rather than taking the limit in sequence. We are now ready to state the main theorem.
[] Consider a random deep neural network of the form in Equations (1) and (2) with a continuous nonlinearity obeying the linear envelope condition 2.3. Then for all sets of strictly increasing width functions and for any countable input set , the distribution of the output of the network converges in distribution to a Gaussian process as . The Gaussian process has mean function zero and covariance function is given by the recursion Lemma 2.3. The convergence in distribution in the statement of the theorem is to be understood in relation to the topology induced by the metric described in Expression (10). Note the generality allowed by the statement in terms of width functions. We could for instance have width functions growing at very different rates, such as . The special case in which all width functions are the identity is most common in other papers on fully connected networks and is used in the majority of our experiments. It is sufficiently important that we state it as a corollary.
Consider a random deep neural network of the form in Equations (1) and (2) with a continuous nonlinearity obeying the linear envelope condition 2.3 and with common number of hidden units for each hidden layer . Then for any countable input set , the distribution of the output of the network converges in distribution to a Gaussian process as . The Gaussian process has mean function zero and covariance function is as in the recursion Lemma 2.3.
We postpone the proof of the main theorem until Section 6. We next look at specific instances of the implied covariance function.
3 Specific kernels under recursion
Cho and 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 such that . Cho and Saul define a recursive kernel through a new feature mapping by compositions such as . 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 2.3 (Hazan and Jaakkola, 2015) for the special case
(16) 
where is the Heaviside step function. Specifically, the recursive approach of Cho and Saul (2009) can be adapted by using the fact that for is equivalent in distribution to with , and by optionally augmenting to incorporate the bias. Since 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 and is defined as:
(17) 
where denotes a reproducing kernel Hilbert space and
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 where required in what follows, we take weight variance parameters and bias variance . We took standard normal input points in dimensions and pass them through independent random neural networks drawn from the distribution discussed in this paper. This was then compared to samples drawn from the corresponding Gaussian process marginal distribution. The experiment was performed with different numbers of hidden layers, different choices of monotonic width functions (which will be described in the sequel), and network sequence index as described in Definition 2.4. We repeated each experiment 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 . 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 and , where the characteristic length scale is chosen to be sensible for the standard normal input distribution on the four dimensional space. Note that there are multiple different uses for kernels in this experiment. The first use is to estimate MMD, the second is for the covariance function of the calibration Gaussian processes and the third use is for the covariance function of the limit Gaussian process. The first and second cases both happen to use the RBF kernel with various length scales, but they should not be confused.
We investigated three choices of strictly increasing width functions, all of which meet the assumptions required by Theorem 2.4 for convergence in distribution to the corresponding Gaussian process. The identity width function corresponds to the case where all hidden layers are the same size and may be directly identified with the width of the network. To test a broader variety of the predictions made by the theory we introduced two other width function specifications. What we call the largest last width function is given by:
(18) 
For example, in a three hidden layer neural network, with , starting from the layer closest to the inputs, we would have hidden layer sizes . The largest first width function is given by:
(19) 
For example in a three hidden layer neural network, with , starting from the layers closest to the inputs, we would have have hidden layer sizes . For both the largest first and largest last width functions the sequence index may be identified with the width of the narrowest hidden layer.
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. Unsurprisingly, since the corresponding networks will have strictly more units, both the largest last and largest first width functions converge faster than the identity width function. The largest last width function seems to converge slightly faster than the largest last width function with respect to this metric. The comparison is more interesting in this case since these two width functions have similar numbers of units. All of the results are consistent with the predictions of Theorem 2.4.
5 Empirical Comparison of 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. For expectations of bounded continuous functions, 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. Fully connected Bayesian deep networks with finite variance priors on the weights have been considered in several recent works (Graves, 2011; HernándezLobato and Adams, 2015; Blundell et al., 2015; HernándezLobato et al., 2016), though the specific details vary. From a Bayesian perspective, the previous section could be interpreted as using MMD as a similarity metric between priors. By constrast, in this section we will compare data dependent quantities that are typically used in Bayesian modelling practice.
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, priordominated problems are generally regarded as an area of strength for Bayesian approaches and in this context our results are directly relevant.
We use hidden layers and hidden units which is typical of the smaller Bayesian neural networks used by HernándezLobato and Adams (2015). HernándezLobato and Adams (2015) also use the variance scaling of Neal (1996)
on their normally distributed weights and give a hierarchical treatment of the hyperparameters. Note that much larger networks have been used in the literature. For example,
Blundell et al. (2015) use as many as units per layer, though they use a two component scale mixture of Gaussians for the weight prior. This would require an extension of our theory to nonGaussian weight distributions for our results to be strictly applicable. Our modest choice of hidden units per layer is partly also motivated by necessity. For larger networks the MCMC would be prohibitively slow.The experiments are divided into those with fixed hyperparameters and those where the hyperparameters are learnt. The hyperparameters are specifically the noise variance, the raw weight variance and the bias variance . The latter two hyperparameters are shared across layers. The fixed hyperparameter experiments are the comparison most directly relevant to the theory presented here. However we found that as we moved to larger datasets both the neural network prior and the Gaussian process prior were often misspecified to an extent that made the results practically uninteresting. Since we were already computationally constrained by the neural network MCMC, we adopted the pragmatic solution of using the type II maximum likelihood parameter estimate of the Gaussian process model for both the neural network and Gaussian process priors. Although the number of hyperparameters is small, this technically adds dependency, so the fixedhyperparameter experiments are complementary.
5.1 Experiments with fixed hyperparameters
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 onedimensional regression problem and a two dimensional realvalued embedding of the four data point XOR problem. To distinguish this from a later larger embedding we term this the small XOR dataset. 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 (SohlDickstein and Culpepper, 2012). The logimportance 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 small XOR problem. The value of the logmarginal 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 loglikelihood 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 standard normal train and test points in 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 loglikelihood is close but not exact. Similarly the marginal function values are close to those of a Gaussian process but are slightly more concentrated.
. 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
hidden layers and units per layer.5.2 Experiments with learnt hyperparameters
As described above, in this section we compare neural networks and the corresponding Gaussian process on larger datasets using hyperparameters for both models that are taken from the learnt Gaussian process kernel, estimated using type II maximum likelihood.
We made a comparison for the data point Snelson dataset, a regression benchmark commonly used in the sparse Gaussian process literature (Snelson and Ghahramani, 2005). Figure 6 shows that the agreement is very close.
Next we made a comparison for a larger embedding of the real valued XOR function which we term the smooth XOR dataset, to distinguish it from the small XOR dataset above. In detail, we have:
(20) 
where and are chosen so that and . One hundred input points are sampled from a standard normal distribution and Gaussian noise of variance is added to the outputs. In order to allow better visualisation of the posterior we take test points along two linear cross sections as shown in Figure 7. This allows us to plot the two posteriors along the crosssections in a manner similar to a one dimensional regression problem. Figure 7 shows the results. We can see that there is again close agreement between the Bayesian neural network posterior and that of the Gaussian process.
Finally, we make a comparison on the Delft yacht hydrodynamics dataset. The task is to predict the residuary resistance per unit weight of displacement for a yacht hull based on six relevant attributes. We randomly partition the data into 100 training instances and 208 test instances. The data has very low noise. To make it a more challenging task for probabilistic modelling we add Gaussian noise of variance . We evaluate per test data point hold out log likelihood for both the Gaussian process and the neural network and the marginal posterior on a randomly selected test function value. The results are shown in Figure 8. The results indicate that on this dataset the Bayesian deep network and the Gaussian process do not make similar predictions. Of the two, the Bayesian neural network achieves significantly better log likelihoods on average, indicating that a finite network performs better than its infinite analogue in this case.
5.3 Summary and discussion of Bayesian posterior comparison
A summary of the datasets studied for comparing posteriors is given in Table 1. Of the datasets studied, the Bayesian neural network showed close agreement with the Gaussian process on five of the six datasets according to the various metrics used, the exception being the yacht dataset. It is notable that the yacht dataset has the highest dimensionality of those considered.
As already noted, our comparison method is computationally expensive as a result of the goldstandard MCMC algorithms used for Bayesian neural network inference. This means we are restricted to relatively small, low dimensional datasets. This caveat is particularly important in light of the yacht data results. On the other hand, we were also limited in the size of finite network we could consider for the same computational reason. As already discussed, the hidden unit networks we use are on the small end of the range of networks that have been studied in the literature (HernándezLobato and Adams, 2015), compared to values as high as used in other works (Blundell et al., 2015). We would, of course, expect that where the model matches the assumption of our theory the agreement would become closer as the number of hidden units increases. As a result of the empirical analysis in Section 4, we would predict more difference if the number of hidden layers was substantially increased, though this has been relatively rare in the existing Bayesian literature thus far.
Bringing these considerations together, it seems likely that some experiments in the literature studied under the banner of Bayesian deep learning would have given very similar results to a Gaussian process with the correct kernel. In the case where the two true posteriors are close, but the posterior approximation for the neural network is significantly worse than any approximation required for the Gaussian process, it would be expected that the Gaussian process would perform better. It should again be noted that the Bayesian neural network experiments were significantly slower than those conducted using the Gaussian process. The Snelson example took 44 hours on ten 3.2 GHz I7 CPU cores to obtain the two million samples required for the Bayesian neural network, where the Gaussian process took a matter of seconds.
Practically, we suggest that the Bayesian deep learning community routinely compare their results to Gaussian processes with the kernels studied here. This will be facilitated by the release of our covariance function code built on GPflow (Matthews et al., 2017). Such a convention would significantly increase our empirical knowledge of the phenomenon studied in this paper.
6 Proof of the main theorem
Let us first sketch the proof we will follow in this section. We first show that, for a countable set of inputs, the infinitedimensional convergence problem can be reduced to a set of onedimensional problems based on finite linear projections. When we examine these onedimensional projections, we find their structure involves a sum of terms we refer to as summands. For fixed width functions the summands are exchangeable, which leads us to consider central limit theorems for exchangeable arrays. A result of Blum et al. (1958) plays an essential role and requires certain moment conditions, that we show by induction through the layers of the network, starting nearest the input. There is a slight complication around the correct scaling of the summands to map onto the exchangeable central limit theorem, but this can be resolved with care.
We already pointed out in Section 2.2 that with a countable index set convergence with respect to the metric is equivalent to convergence of each finitedimensional marginal. The CramérWold device (Cramér and Wold, 1936) (Billingsley, 1986, p. 383) states that convergence of a sequence of finitedimensional vectors to some limit is equivalent to convergence on all possible linear projections to the corresponding realvalued random variable. Putting these two results together we obtain the following lemma.
[Convergence of finite linear projections] Consider a sequence of random functions taking values in each defined on a countable input set , with the sequence of functions indexed by . Let be a finite subset of the input set. Further, let . Then convergence in distribution of the sequence of random functions taking values in to a limiting random function with respect to the metric is equivalent to weak convergence of to the corresponding finite linear projection for every such and .
Therefore our task is reduced to that of proving convergence of a sequence of realvalued random variables to another realvalued random variable – a considerable simplification. In particular, we will leverage a theorem of Blum et al. (1958) on central limit theorems for exchangeable sequences.
It will be convenient to consider a sequence of ‘infinite width, finite fanout, networks’. By this we mean that the indices in the recursion (2) can be thought of as running over all natural numbers instead of just up to (hence infinite width). The limits of the sums in the recursion will retain the same finite values, which depend on the width functions evaluated at some (hence finite fanout). This makes only a superficial change because it adds extra copies of the same variables at each depth. For fixed , these extra variables will not effect the downstream distribution of the network. The change is however useful in the bookkeeping needed to prove convergence. We have defined a countable sequence of such networks because is a natural number.
It will also be useful to slightly rewrite the defining initialisation and recursion (2) from the more familiar form to one which is easier to manipulate:
(21) 
and:
(22)  
(23) 
where:
(24) 
This amounts to reparameterising the weights in terms of standard normals and making the previously mentioned infinite extension of the width variable . We reemphasize that neither step changes the distribution over final function values. With the aim of mapping onto Lemma 6 we make the following definitions:
[Projections and summands] The projections are defined in terms of a finite linear projection of the function values without biases:
(25) 
where is a finite set of tuples of data points and indices of prenonlinearities, with . is a vector parameterising the linear projection. The suffix indicates that the corresponding width functions are instantiated with input .
The summands are defined as:
(26) 
in order to ensure the summation relation
(27) 
The last relation follows from applying the definitions and rearranging the order of summation. Note the similarity between the definition of projections used here and in Lemma 6. We next show that the summands are exchangeable.
[Exchangeability of summands] For each fixed and , the countable sequence of summands are an exchangeable sequence with respect to the index .
To prove the lemma we use de Finetti’s theorem, which states that a sequence of random variables is exchangeable if and only if they are i.i.d. conditional on some set of random variables. It is therefore sufficient to exhibit such as set of random variables. To do this we apply the recursion. Removing some multiplicative constants we have:
(28)  
(29) 
with the convention that and for . Conditional on the finite set of random variables (where is the set of inputs points in ), the summands are independent and identically distributed.
Thus we are led to consider central limit theorems for sequences of exchangeable sequences. The work of Blum et al. (1958) will provide our starting point.
[CLT for exchangeable sequences (Blum et al., 1958)] For each positive integer let be an infinitely exchangeable process with mean zero, variance one, and finite absolute third moment. Define
(30) 
Then if the following conditions hold:
Then converges in distribution to a standard normal.
This is effectively a generalisation of the classical CLT from independent identically distributed variables to the more general class of exchangeable ones. We will need to address the fact that the theorem applies to unit variance variables and that we have nonidentity width functions. The next lemma adapts the work of Blum et al. to our specific requirements.
[Adapted CLT for sequences of exchangeable sequences] For each positive integer let be an infinitely exchangeable process with mean zero, finite variance , and finite absolute third moment. Suppose also that the variance has a limit . Define
(31) 
where is a strictly increasing function. Then if the following conditions hold:
Then converges in distribution to , where is interpreted as converging to .
We postpone the proof of Lemma 6 until Appendix A. Our next step will be to apply Lemma 6 to the projections and summands by showing they meet each condition. We first establish the existence of a limiting variance.
[Limiting variance] The limiting variance, defined as
(32) 
exists, where is the variance of the random variables , and has the value
(33) 
where is the Gram matrix implied by the recursion 2.3 without a bias correction on the final layer.
The proof of this Lemma can be found in Appendix B.1.
[Convergence in distribution of projections] As the projection converges in distribution to .
The full details of Lemma 6 are explained in Appendix B.1. Here we outline the key points of the approach. We apply Lemma 6 to the projections, using the fact that the summands are exchangeable for each and with the limiting variance derived in Lemma 6. Condition a) of Lemma 6 follows straightforwardly from the fact that the summands are uncorrelated. That Condition c) is fulfilled is intuitively reasonable given that we in fact expect this absolute third moment to tend to a constant. Condition c) will however still need to be shown carefully. This leaves Condition b). Convergence of the expectation of a sequence of random variables can be ensured if the sequence is uniformly integrable and the sequence converges in distribution (Billingsley, 1999). Thus the main work of Appendix B.1 is to prove these conditions in our case, by induction forwards through the network.
Lemma 6 shows consistency of convergence of the finite linear projections of the prebias function distribution with the stated Gaussian process. By Lemma 6, this is sufficient for convergence in distribution to the Gaussian process. As the biases are normally distributed it is straightforward to add them and get the final result. Therefore we are done.
7 Desirability of Gaussian process behaviour and methods to avoid it
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. The kernels considered in this paper have not been commonly used in the Gaussian process literature and warrant further analysis. It has been argued by previous authors 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 kernels considered in this paper do not meet the strict definition of what could be considered local, though the Euclidean inner product between two points is sufficient to compute the corresponding covariance. In any case, the fact remains that a 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. A complication to consider is when a hierarchical treatment of the model is taken, learning model hyperparameters. Typically only a few such hyperparameters are used and it seems unlikely this could offer the same benefits as full representation learning. Using significantly more hyperparameters would move the model beyond the scope of this paper. MacKay (2002, p. 547) famously reflected on what is lost when taking the Gaussian process limit of a single hidden layer network, remarking that Gaussian processes will not learn hidden features. Neal (1996, p. 43) makes similar comments and also expresses the hope that Bayesian neural networks could expand the range of probabilistic models beyond Gaussian processes. In light of the results in this paper for networks with more than one hidden layer these considerations are of considerable importance going forward.
There is literature on learning the representation of a standard, usually structured, network composed with a Gaussian process (Wilson et al., 2016a, b; AlShedivat 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 the paradigm considered here, 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 and 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. These alternatives each arguably have a mechanism to discover hierarchies of features. Again, given the convergence results for multiple hidden layer networks from this paper, there is now further motivation to study the nonGaussian alternatives as well.
8 Conclusions
Studying the limiting behaviour of distributions on feedforward neural networks has been a fruitful avenue for understanding these models historically. In this paper we have formalised and extended prior results by Neal (1996) to deep networks. In particular, we have shown that, under broad conditions, as we make the architecture increasingly wide, the implied random function converges in distribution to a Gaussian process. 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. In several cases there was close agreement, leading us to conclude that it is likely some results from the existing Bayesian deep learning literature would be very similar to those obtained with the corresponding Gaussian process model. We recommend that empirical investigation of Bayesian neural networks should routinely include comparison to their Gaussian process analogue. If 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 there may well be cases where the behaviour is not 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.
9 Acknowledgements
We wish to thank Neil Lawrence and Jascha SohlDickstein for helpful conversations. We also thank anonymous reviewers of previous versions for their insights. Alexander Matthews and Zoubin Ghahramani acknowledge the support of EPSRC Grant EP/N014162/1 and EPSRC Grant EP/N510129/1 (The Alan Turing Institute). Jiri Hron holds a Nokia CASE Studentship. Mark Rowland acknowledges support by EPSRC grant EP/L016516/1 for the Cambridge Centre for Analysis. Richard E. Turner is supported by Google as well as EPSRC grants EP/M0269571 and EP/L000776/1.
References
 AlShedivat et al. (2017) Maruan AlShedivat, Andrew G. Wilson, Yunus Saatchi, Zhiting Hu, and Eric P. Xing. Learning Scalable Deep Kernels with Recurrent Structure. Journal of Machine Learning Research (JMLR), 2017.

Bengio et al. (2005)
Yoshua Bengio, Olivier Delalleau, and Nicolas Le Roux.
The Curse of Dimensionality for Local Kernel Machines.
Technical Report 1258, Département d’informatique et recherche opérationnelle, Université de Montréal, 2005.  Billingsley (1986) Patrick Billingsley. Probability and Measure. John Wiley and Sons, second edition, 1986.
 Billingsley (1999) Patrick Billingsley. Convergence of Probability Measures. John Wiley & Sons Inc., Second edition, 1999.
 Blum et al. (1958) J. R. Blum, H. Chernoff, M. Rosenblatt, and H. Teicher. Central limit theorems for interchangeable processes. Canadian Journal of Mathematics, 10:222–229, 1958.
 Blundell et al. (2015) Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra. Weight Uncertainty in Neural Networks. International Conference on Machine Learning (ICML), 2015.
 Cho and Saul (2009) Youngmin Cho and Lawrence K. Saul. Kernel Methods for Deep Learning. Advances in Neural Information Processing Systems (NIPS), 2009.
 Cramér and Wold (1936) H. Cramér and H. Wold. Some theorems on distribution functions. Journal of the London Mathematical Society, s111(4):290–294, 1936.

Damianou and Lawrence (2013)
Andreas C. Damianou and Neil D. Lawrence.
Deep Gaussian Processes.
International Conference on Artificial Intelligence and Statistics (AISTATS)
, 2013.  Daniely et al. (2016) Amit Daniely, Roy Frostig, and Yoram Singer. Toward Deeper Understanding of Neural Networks: The Power of Initialization and a Dual View on Expressivity. Advances in Neural Information Processing Systems (NIPS), 2016.
 Dashti and Stuart (2013) M. Dashti and A. M. Stuart. The Bayesian Approach To Inverse Problems. ArXiv eprints, February 2013.
 Duvenaud et al. (2014) David Duvenaud, Oren Rippel, Ryan P. Adams, and Zoubin Ghahramani. Avoiding Pathologies in very Deep Networks. International Conference on Artificial Intelligence and Statistics (AISTATS), 2014.
 Graves (2011) Alex Graves. Practical Variational Inference for Neural Networks. Advances in Neural Information Processing Systems (NIPS), 2011.
 Gretton et al. (2012) Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola. A Kernel Twosample test. Journal of Machine Learning Research (JMLR), 2012.
 Grosse et al. (2015) Roger B. Grosse, Zoubin Ghahramani, and Ryan P. Adams. Sandwiching the marginal likelihood using bidirectional Monte Carlo. ArXiv eprints, November 2015.
 Hazan and Jaakkola (2015) Tamir Hazan and Tommi Jaakkola. Steps Toward Deep Kernel Methods from Infinite Neural Networks. ArXiv eprints, August 2015.

HernándezLobato and Adams (2015)
José M. HernándezLobato and Ryan P. Adams.
Probabilistic Backpropagation for Scalable Learning of Bayesian Neural Networks.
International Conference on Machine Learning (ICML), 2015.  HernándezLobato et al. (2016) José M. HernándezLobato, Yingzhen Li, Mark Rowland, Thang Bui, Daniel HernándezLobato, and Richard E. Turner. Blackbox alpha divergence minimization. International Conference on Machine Learning (ICML), 2016.
 Klambauer et al. (2017) Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter. Selfnormalizing neural networks. In Advances in Neural Information Processing Systems (NIPS). 2017.
 Krauth et al. (2017) Karl Krauth, Edwin V. Bonilla, Kurt Cutajar, and Maurizio Filippone. AutoGP: Exploring the capabilities and limitations of Gaussian Process models. Conference on Uncertainty in Artificial Intelligence (UAI), 2017.
 Lee et al. (2018) Jaehoon Lee, Yasaman Bahri, Roman Novak, Samuel S. Schoenholz, Jeffrey Pennington, and Jascha SohlDickstein. Deep Neural Networks as Gaussian Processes. International Conference on Learning Representations (ICLR), 2018.
 MacKay (2002) David J. C. MacKay. Information Theory, Inference & Learning Algorithms. Cambridge University Press, 2002.
 Mandt et al. (2017) S. Mandt, M. D. Hoffman, and D. M. Blei. Stochastic Gradient Descent as Approximate Bayesian Inference. ArXiv eprints, April 2017.

Matthews et al. (2017)
Alexander G. de G. Matthews, Mark van der Wilk, Tom Nickson, Keisuke Fujii,
Alexis Boukouvalas, Pablo LeónVillagrá, Zoubin Ghahramani, and
James Hensman.
GPflow: A Gaussian Process Library using TensorFlow.
Journal of Machine Learning Research, 18(40):1–6, 2017.  Matthews et al. (2018) Alexander G. de G. Matthews, Jiri Hron, Mark Rowland, Richard E. Turner, and Zoubin Ghahramani. Gaussian Process Behaviour in Wide Deep Neural Networks. In International Conference on Learning Representations (ICLR), 2018.
 Mitrovic et al. (2017) Jovana Mitrovic, Dino Sejdinovic, and Yee Whye Teh. Deep Kernel Machines via the Kernel Reparametrization Trick. In International Conference on Learning Representations (ICLR) Workshop Track, 2017.
 Murray et al. (2010) Iain Murray, Ryan P. Adams, and David J. C. MacKay. Elliptical Slice Sampling. International Conference on Artificial Intelligence and Statistics (AISTATS), 2010.
 Neal (1996) Radford M. Neal. Bayesian Learning for Neural Networks. Springer, 1996.
 Neal (2010) Radford M. Neal. MCMC using Hamiltonian Dynamics. Handbook of Markov Chain Monte Carlo, 2010.
 Poole et al. (2016) Ben Poole, Subhaneil Lahiri, Maithreyi Raghu, Jascha SohlDickstein, and Surya Ganguli. Exponential expressivity in Deep Neural Networks through Transient Chaos. Advances in Neural Information Processing Systems (NIPS), 2016.
 Rasmussen and Williams (2006) Carl E. Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2006.
 Schoenholz et al. (2017) Samuel S. Schoenholz, Justin Gilmer, Surya Ganguli, and Jascha SohlDickstein. Deep Information Propagation. International Conference on Learning Representations (ICLR), 2017.
 Smith and Le (2018) Samuel L. Smith and Quoc V. Le. A Bayesian perspective on generalization and stochastic gradient descent. International Conference on Learning Representations (ICLR), 2018.
 Snelson and Ghahramani (2005) Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudoinputs. Advances in Neural Information Processing Systems (NIPS), 2005.
 SohlDickstein and Culpepper (2012) Jascha SohlDickstein and Benjamin J. Culpepper. Hamiltonian Annealed Importance Sampling for partition function estimation. CoRR, abs/1205.1925, 2012.
 Vaart (1998) A. W. van der Vaart. Asymptotic Statistics. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 1998.
 Welling and Teh (2011) Max Welling and Yee Whye Teh. Bayesian learning via stochastic gradient langevin dynamics. International Conference on Machine Learning (ICML), 2011.
 Williams (1998) Christopher K. I. Williams. Computing with Infinite Networks. Advances in Neural Information Processing Systems (NIPS), 1998.
 Wilson et al. (2016a) Andrew G. Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P. Xing. Deep Kernel Learning. International Conference on Artificial Intelligence and Statistics (AISTATS), 2016a.
 Wilson et al. (2016b) Andrew G. Wilson, Zhiting Hu, Ruslan R. Salakhutdinov, and Eric P. Xing. Stochastic Variational Deep Kernel Learning. Advances in Neural Information Processing Systems (NIPS), 2016b.
Appendix A Adapting the exchangeable CLT of Blum et al. 1958.
This section gives further detail on our adaption of Theorem 6 to our specific needs. It states and proves an intermediate Lemma A and then using that lemma gives the postponed proof of Lemma 6.
[Variance adapted CLT for sequences of exchangeable sequences] For each positive integer let be an infinitely exchangeable process with mean zero, finite variance , and finite absolute third moment. Suppose also that the variance has a limit . Define
(34) 
Then if the following conditions hold:
Then converges in distribution to , where is interpreted as converging to .
[Proof of Lemma 6] Either or it does not. We deal with each case separately.
In the case where , we have:
(35)  
(36)  
(37) 
where we have used property (i) that the distinct elements in a row are uncorrelated. Now the proof can take a similar route to that used in proving the weak law of large numbers with finite variance. Chebyshev’s inequality we have:
(38) 
for all , So that:
(39) 
That is to say converges in probability to . In such a case of a constant target convergence in probability is equivalent to convergence in distribution.
In the case where there is always some such that for by the definition of a limit. Let us assume we are in this range of . Then the standardised values will obey the conditions of Theorem 6, which we now show term by term. We clearly have mean zero, unit variance and finite third moment, and conditions 1) and i) are identical. This leaves us to validate conditions 2) and 3). Starting from ii) we have:
(40)  
(41)  
(42) 
which clearly implies condition 2). Starting from condition iii) we have:
(43)  
(44) 
which implies condition 3). Since multiplication by a constant is a continuous function we have therefore showed that converges in distribution to .
Note that the sequence converges surely to . This certainly implies convergence in probability of the same sequence to zero. We can therefore invoke a general result on convergence of sequences that says if a sequence of random variables converges to and converges in probability to zero, then converges in distribution to (Vaart, 1998).
[Proof of Lemma 6] Lemma A applies to what are known as triangular arrays in the literature. This lemma is the generalisation to arrays that are not strictly triangular. To do this we embed the nontriangular array in a large triangular one. We fill the extra spaces with standard normal random variables. This gives an interleaved sequence. The terms we actually care about will obey the necessary conditions if conditions 1) 2) and 3) if a) b) and c) hold. The conditions 1) 2) and 3) will hold trivially for the standard normal rows. Thus the whole sequence converges in distribution. But since any subsequence also converges in distribution we get our required result.
Appendix B Details of the proof of Theorem 2.4
Here, we summarise the highlevel structure of the proof of Theorem 2.4. The argument is inductive, showing sequentially that the hidden units in each layer of the network converge in distribution; to avoid repetition, all mentions of convergence in distribution of infinitedimensional random variables in what follows are specifically with respect to the topology generated by the metric introduced in Section 2.2. The main part of the inductive argument is summarised in the following proposition.
For any , suppose that the collection of random variables converges in distribution as to a centred Gaussian with covariance function of the form given in Lemma 2.3. Then any finite linear combination (with finite and ) of prenonlinearities at the next layer also converges in distribution to a centred Gaussian of the form described in Lemma 6.
Note that the conclusion of Proposition B leads to the statement of Lemma 6. By the CramérWold device discussed in Section 6, the convergence of the finite linear projections established in Proposition B guarantees convergence of all finitedimensional marginal distributions. Adding in the independent bias terms yields convergence of finitedimensional marginals of the preactivations at layer
; this may be demonstrated via a standard argument using characteristic functions. Due to the remarks on weak convergence in Section
2.2, convergence in distribution of all finitedimensional marginals guarantees convergence in distribution of the full collection of random variables in the next layer, completing the inductive step.The proof of Theorem 2.4 is then concluded by observing that the prenonlinearities in the first hidden layer,
, have a fixed Gaussian distribution that does not depend on
.We thus turn our attention to proving Proposition B. The main idea is to use Lemma 6, taking each of the random variables (for ) appearing in the statement of the Lemma to be the summands appearing in the finite linear projections :
(45) 
Addressing the conditions of Lemma 6, we note that the exchangeability condition is provided by Lemma 6, the meanzero condition is immediate, the limiting variance condition is dealt with by Lemma 6. Condition a) of Lemma 6 holds trivially as the random variables and are meanzero and uncorrelated. The remaining conditions of Lemma 6 are dealt with through the following results; Lemma B deals with Condition b), whilst Lemma B deals with the growth of third absolute moments as required by Condition c).
[Convergence of ] Consider arbitrary and the corresponding set of random variables . Assume that the countably infinite vector of random variables converges in distribution to a centred Gaussian process with covariance specified by the recursion in Lemma 2.3 as . Then
[Bound on ] For arbitrary given , , and , with independent of . Thus .
Thus, all that remains to establish Theorem 2.4 is the proof of these intermediate lemmas; the proofs are given in the sections that follow.
b.1 Proofs of main lemmas and corollaries
Throughout this section, we simplify the notation by defining
where is the ^{th} member of the set . Without loss of generality, in what follows we will take to lighten notation.
To prove Lemma B, we need to know the value of where as defined in Lemma 6. Lemma B.1 combined with the inductive propagation of convergence in distribution verifies Lemma 6 and thus yields .
Consider arbitrary . Assume that the countably infinite vector of random variables converges in distribution to a centred Gaussian process with covariance specified by the recursion in Lemma 2.3 as . Then
Lemma 6 introduces which is the marginal covariance of the limiting Gaussian process without the bias term (c.f. the recursion in Lemma 2.3).
We use exchangeability of over the index to obtain
Hence the limit of is fully determined by the behaviour of as .
We can thus focus on individual entries of the expectation on the RHS of the above equation. For entry with , we have
Since and the collection converges in distribution as by assumption, we can use the continuity of and the continuous mapping theorem to deduce that the postnonlinearities are converging in distribution. Because the function is continuous, we can apply the continuous mapping theorem again to deduce that the twoway products of postnonlinearities are converging in distribution to the limit specified by the pushforward of the limiting multivariate normal distribution.
Theorem 3.5 in (Billingsley, 1999) tells us that the expectation
if the family of random variables indexed by is uniformly integrable. Uniform integrability is a corollary of Lemma B.2. Inspection of the recursion in Lemma 6 finishes the proof.
[Proof of Lemma B] Substituting for and , we have
(46) 
The expectation on the RHS can be rewritten as
Hence the limit of the LHS of Equation (46) is fully determined by the behaviour of , as . We can thus focus on individual entries of the expectation on the RHS of the above equation. For entry with , we have
(47) 