Title: NON-GAUSSIAN PROCESSES AND NEURAL NETWORKS AT FINITE WIDTHS

Abstract: Gaussian processes are ubiquitous in nature and engineering. A case in point is a class of neural networks in the infinite-width limit, whose priors correspond to Gaussian processes. Here we perturbatively extend this correspondence to finitewidth neural networks, yielding non-Gaussian processes as priors. The methodology developed herein allows us to track the flow of preactivation distributions by progressively integrating out random variables from lower to higher layers, reminiscent of renormalization-group flow. We further develop a perturbative prescription to perform Bayesian inference with weakly non-Gaussian priors.

1. INCEPTION:

Gaussian processes model many phenomena in the physical world. A prime example is Brownian motion (Brown, 1828) , modeled as the integral of Gaussian-distributed bumps exerted on a pointlike solute (Einstein, 1905) . The theory of elementary particles (Weinberg, 1995) also becomes a Gaussian process in the free limit where interactions between particles are turned off, and manybody systems as complex as glasses come to be Gaussian in the infinite-dimensional, mean-field, limit (Parisi & Zamponi, 2010) . In the context of machine learning, Neal (1996) pointed out that a class of neural networks give rise to Gaussian processes in the infinite-width limit, which can perform exact Bayesian inference from training to test data (Williams, 1997) . They occupy a corner of theoretical playground wherein the karakuri of neural networks is scrutinized (Lee et al., 2018; Matthews et al., 2018; Jacot et al., 2018; Chizat et al., 2018; Lee et al., 2019; Geiger et al., 2019) .

In reality, Gaussian processes are but mere idealizations. Brownian particles have finite-size structure, elementary particles interact, and many-body systems respond nonlinearly. In order to understand rich phenomena exhibited by these real systems, Gaussian processes rather serve as starting points to be perturbed around. Indeed many edifices in theoretical physics are build upon the successful treatment of non-Gaussianity, with a notable example being renormalization-group flow (Kadanoff, 1966; Wilson, 1971; Weinberg, 1996; Goldenfeld, 2018) . In the quest to elucidate behaviors of real neural networks away from the infinite-width limit, it is thus natural to wonder if the similar treatment of non-Gaussianity yields equally elegant and powerful machinery.

Here we set out on this program, perturbatively treating finite-width corrections to neural networks. Prior distributions of outputs are obtained through progressively integrating out preactivation of neurons layer by layer, yielding non-Gaussian priors. Intriguingly, intermediate recursion relations and their derivation resemble renormalization-group flow (Goldenfeld, 2018; Mehta & Schwab, 2014) . Such a recursive approach further enables us to treat finite-width corrections on Bayesian inference and their regularization effects, with arbitrary activation functions.

The rest of the paper is structured as follows. In Section 2 we review and set up basic concepts. Our master recursive formulae (R1,R2,R3) are derived in Section 3, which control the flow of preactivation distributions from lower to higher layers. After an interlude with concrete examples in Section 4, we extend the Gaussian-process Bayesian inference method to non-Gaussian priors in Section 5 and use the resulting scheme to study inference of neural networks at finite widths. We conclude in Section 6 with dreams.

2. TO INFINITY AND BEYOND:

In this paper we study real finite-width neural networks in the regime where the number of neurons in hidden layers is asymptotically large whereas input and output dimensions are kept constant.

2.1. GAUSSIAN PROCESSES AND NEURAL NETWORKS AT INFINITE WIDTHS:

Let us focus on a class of neural networks termed multilayer perceptrons, with model parameters,

θ = b ( ) i , W ( ) i,j

, and an activation function, σ. For each input, x ∈ R n0 , a neural network outputs a vector, z(x; θ) = z (L) ∈ R n L , recursively defined as sequences of preactivations through

EQUATION

z

EQUATION

Following Neal (1996), we assume priors for biases and weights given by independent and identically distributed Gaussian distributions with zero means, E b

( ) i = E W ( )

i,j = 0, and variances

EQUATION

EQUATION

Higher moments are then obtained by Wick's contractions (Wick, 1950; Zee, 2010) . For instance,

EQUATION

For those unfamiliar with Wick's contractions and connected correlation functions (a.k.a. cumulants), a pedagogical review is provided in Appendix A as our formalism heavily relies on them.

In the infinite-width limit where n 1 , n 2 , . . . , n L-1 → ∞ (but finite n 0 and n L ), it has been argued -with varying degrees of rigor (Neal, 1996; Lee et al., 2018; Matthews et al., 2018) -that the prior distribution of outputs is governed by the Gaussian process with a kernel

EQUATION

and all the higher moments given by Wick's contractions. In particular, there exists a recursive formula that lets us evaluate this kernel for any pair of inputs (Lee et al., 2018) [c.f. Equation (R1)]. Importantly, once the values of the kernel are evaluated for all the pairs of , 1997; Williams & Rasmussen, 2006) [c.f. Equation (GPM) ]. This should be contrasted with stochastic gradient descent (SGD) optimization (Robbins & Monro, 1951) , through which typically a single estimate for the optimal model parameters of the posterior, θ , is obtained and used to predict outputs for test inputs; Bayesian inference instead marginalizes over all model parameters, performing an ensemble average over the posterior distribution (MacKay, 1995).

N D = N R + N E input data, {x α } α=1,...,

2.2. BEYOND INFINITY:

We shall now study real finite-width neural networks in the regime n 1 , . . . , n L-1 ∼ n 1. 1 At finite widths, there are corrections to Gaussian-process priors. In other words, a whole tower of 1 Note that input and output dimensions, n0 and nL, are arbitrary. To be precise, defining n1, . . . , nL-1 ≡ µ1n, . . . , µL-1n, we send n

1 while keeping C ( ) b , C ( ) W =1,...,L

, µ1, . . . , µL-1, n0, and nL constants, and compute the leading 1/n corrections. In particular it is crucial to keep the number of outputs nL constant in order to consistently perform Bayesian inference within our approach.

nontrivial preactivation correlation functions beyond the kernel,

EQUATION

collectively dictate the distribution of preactivations. Our aim is to trace the flow of these distributions progressively and cumulatively all the way up to the last layer whereat Bayesian inference is executed. More specifically, we shall inductively and self-consistently show that two-point preactivation correlation functions take the form 2

G ( )α1,α2 i1,i2 = δ i1,i2 K α1,α2 ( ) + 1 n -1 S α1,α2 ( ) + O 1 n 2 (KS)

and connected four-point preactivation correlation functions

G ( )α1,α2,α3,α4 i1,i2,i3,i4 connected (V) ≡G ( )α1,α2,α3,α4 i1,i2,i3,i4 -G ( )α1,α2 i1,i2 G ( )α3,α4 i3,i4 -G ( )α1,α3 i1,i3 G ( )α2,α4 i2,i4 -G ( )α1,α4 i1,i4 G ( )α2,α3 i2,i3 = 1 n -1 δ i1,i2 δ i3,i4 V (α1,α2),(α3,α4) ( ) + δ i1,i3 δ i2,i4 V (α1,α3),(α2,α4) ( ) + δ i1,i4 δ i2,i3 V (α1,α4),(α2,α3) ( ) + O 1 n 2 ,

and higher cumulants are all suppressed by O 1 n 2 . 3 Here the Gaussian-process core kernel K α1,α2

( )

and the self-energy correction S α1,α2

( )

are symmetric under the exchange of sample indices α 1 ↔ α 2 and the four-point vertex

V (α1,α2),(α3,α4) ( )

is symmetric under α 1 ↔ α 2 , α 3 ↔ α 4 , and (α 1 , α 2 ) ↔ (α 3 , α 4 ). At the first layer the preactivation distribution is exactly Gaussian for any finite widths and hence Equations (KS) and (V) are trivially satisfied, with

EQUATION

= 0 , and

V (α1,α2),(α3,α4) (1) = 0 . (R0)

Obtained in Section 3 are the recursive formulae that link these core kernel, self-energy, and fourpoint vertex at the -th layer to those at the ( + 1)-th layer while in Section 5 these tensors at the last layer = L are used to yield the leading 1/n correction for Bayesian inference at finite widths.

2.3. RELATED WORK:

Our Schwinger operator approach is orthogonal to the replica approach by Cohen et al. (2019) and, unlike the planar diagrammatic approach by Dyer & Gur-Ari (2019) , applies to general activation functions, made possible by accumulating corrections layer by layer rather than dealing with them all at once. More substantially, in contrast to these previous approaches, we here study finite-width effects on Bayesian inference and find that the renormalization-group picture naturally emerges.

3. DISTRIBUTIONAL FLOW:

As auxiliary objects in recursive steps, let us introduce activation correlation functions

EQUATION

)

Our basic strategy is to establish relations ZIGZAG) zigzagging between sets of preactivation correlation functions and sets of activation correlation functions, keeping track of leading finite-width corrections. Below, relations G ( ) → H ( ) are obtained by integrating out preactivations while relations H ( ) → G ( +1) are obtained by integrating out biases and weights. At first glance the algebra in this paper may look horrifying but repeated applications of Wick's contractions are all there is to it. The results are summarized in Section 3.2.

EQUATION

2 In the main text we place tildes on objects that depend only on sample indices α's in order to distinguish them from those that depend both on sample indices α's and neuron indices i's.

3 Given that the means of biases and weights are zero, G ( )α 1 ,...,αm i 1 ,...,im = 0 for all odd m.

3.1. ZIGZAG RELATIONS FOR PREACTIVATION AND ACTIVATION CORRELATION FUNCTIONS:

Integrating over the Gaussian biases and weights at 's connections yield the relations that link activation correlations H ( ) to preactivation correlations G ( +1) at the next layer. Recalling Equations (KS) and (V), trivial Wick's contractions yield

K α1,α2 ( +1) + 1 n S α1,α2 ( +1) + O 1 n 2 = C ( +1) b + C ( +1) W   1 n n j=1 H ( )α1,α2 j,j   and (9) V (α1,α2),(α3,α4) ( +1) C ( +1) W 2 + O 1 n = 1 n n j,k=1 H ( )α1,α2,α3,α4 j,j,k,k -H ( )α1,α2 j,j H ( )α3,α4 k,k

.( 10)

The remaining task is to relate preactivation correlations G ( ) to activation correlations H ( ) within the same layer, which will complete the zigzag relation (ZIGZAG) for these correlation functions. 4Once the sorcery of Wick's contractions and connected correlation functions is mastered, it is simple to derive the following combinatorial hack (Appendix A.4): viewing prior preactivations

z ( ) ≡ z α i ≡ z ( ) i (x α ) i=1,...,n ;α=1,...,ND

as a random (n N D )-dimensional vector and defining the Gaussian integral with the kernel z

α1 i1 z α2 i2 K ( ) = K ( )α1,α2 i1,i2 ≡ δ i1,i2 K α1,α2 ( )

, the prior average

E F[z ( ) ] = F[z ( ) ] K ( ) + 1 n -1 F[z ( ) ]O S [z ( ) ] + F[z ( ) ]O V [z ( ) ] K ( ) + O 1 n 2 (HACK)

for any function F. Here the operators O S [z ( ) ] and O V [z ( ) ] capture 1/n corrections due to selfenergy and four-point vertex, respectively, and are defined as

O S [z ( ) ] ≡ 1 2 α1,α2 S ( ) α1,α2 n i=1 z α1 i z α2 i -n K α1,α2 ( )

and (OS):

EQUATION

,

where the sample indices are lowered by using the inverse core kernel as a metric, meaning

S ( ) α1,α2 ≡ α 1 ,α 2 K -1 ( ) α1,α 1 K -1 ( ) α2,α 2 S α 1 ,α 2 ( )

and ( 11)

EQUATION

)

Using the above hack, we can evaluate the activation correlations by straightforward algebra with -you guessed it -Wick's contractions. As the Gaussian integral is diagonal in the neuron index i, we just need to disentangle cases with repeated and unrepeated neuron indices. The solution for this exercise is in Appendix B: this is the most cumbersome algebra in the paper and the ability to perform it certifies the graduation from the magical school of Wick's crafts and wizardly cumulants.

3.2. MASTER RECURSIVE FLOW EQUATIONS:

Denoting the Gaussian integral with the core kernel zα1 zα2

K ( ) = K α1,α2 ( )

for a single-neuron random vector z( ) ≡ {z α } α=1,...,ND , and plugging in results of Appendix B into Equations ( 9) and ( 10), we arrive at our master recursion relations

K α1,α2 ( +1) = C ( +1) b + C ( +1) W σ(z α1 )σ(z α2 ) K ( ) , (R1) V (α1,α2),(α3,α4) ( +1) = C ( +1) W 2 σ(z α1 )σ(z α2 )σ(z α3 )σ(z α4 ) K ( ) (R2) -σ(z α1 )σ(z α2 ) K ( ) σ(z α3 )σ(z α4 ) K ( ) + 1 4 n n -1 α 1 ,α 2 ,α 3 ,α 4 V ( ) (α 1 ,α 2 ),(α 3 ,α 4 ) σ(z α1 )σ(z α2 )(z α 1 zα 2 -K α 1 ,α 2 ( ) ) K ( ) × σ(z α3 )σ(z α4 )(z α 3 zα 4 -K α 3 ,α 4 ( ) ) K ( )

, and

S α1,α2 ( +1) = n n -1 C ( +1) W 1 2 α 1 ,α 2 S ( ) α 1 ,α 2 σ(z α1 )σ(z α2 )(z α 1 zα 2 -K α 1 ,α 2 ) K ( ) (R3) + 1 8 α 1 ,α 2 ,α 3 ,α 4 V ( ) (α 1 ,α 2 ),(α 3 ,α 4 ) σ(z α1 )σ(z α2 ) × zα 1 zα 2 zα 3 zα 4 -2z α 1 zα 2 K α 3 ,α 4 -4z α 1 zα 3 K α 2 ,α 4 + K α 1 ,α 2 K α 3 ,α 4 + 2 K α 1 ,α 3 K α 2 ,α 4 K ( )

.

For = 1, a special note about the ratio n n -1 is in order: even though n 0 stays constant while n 1 1, the terms proportional to that ratio are identically zero due to the complete Gaussianity (R0).

The preactivation distribution in the first layer (R0) sets the initial condition for the flow from lower to higher layers dictated by these recursive equations. Once recursed up to the last layer = L, the resulting distribution of outputs z = z (L) can be succinctly encoded by the probability distribution

p[z] = e -H[z] dz e -H[z ] (D0) with the potential H[z] = H 0 [z] + H 1 [z] + O( 2 ) where ≡ 1 n L-1 1, H 0 [z] = 1 2 α1,α2 K -1 (L) α1,α2 n L i=1 z α1 i z α2 i

, and (D1)

EQUATION

)

By now, the reader should be deriving this through Wick's contractions without solicitation. It is important to note that n L is constant and thus H 1 [z] can consistently be treated perturbatively.5 

4. INTERLUDE: EXAMPLES:

The recursive relations obtained above can be evaluated numerically (Lee et al., 2018) [or sometimes analytically for ReLU (Cho & Saul, 2009) ], which is a perfectly adequate approach: at the leading order it involves four-dimensional Gaussian integrals at most. Here, continuing the theme of wearing out Wick's contractions, we develop an alternative analytic method that works for any polynomial activations (Liao & Poggio, 2017) , providing another perfectly cromulent approach.

For a general polynomial activation of degree p, σ(z) = p k=0 a k z k , the nontrivial term in Equation (R1) can be expanded as

EQUATION

)

Each term can then be evaluated by Wick's contractions and the same goes for all the terms in Equations (R2) and (R3). Below and in Appendix C, we illustrate this procedure with simple examples.

4.1. DEEP LINEAR NETWORKS:

When the activation function is linear, σ(z) = z, multilayer perceptrons go under the aweinspiring name of deep linear networks (Saxe et al., 2013) . Setting C

( ) b = 0 and C ( ) W = 1 for simplicity, our recursion relations reduce to K α1,α2 ( +1) = K α1,α2 ( ) , V (α1,α2),(α3,α4) ( +1) = K α1,α3 ( ) K α2,α4 ( ) + K α1,α4 ( ) K α2,α3

( )

+ n n -1 V (α1,α2),(α3,α4) ( )

, and S α1,α2

EQUATION

. Solving them yields the layer-independent core kernel and zero self-energy

EQUATION

=:

x α1 • x α2 n 0 and S α1,α2

EQUATION

and the linearly layer-dependent four-point vertex

EQUATION

K α2,α4

(1)

EQUATION

EQUATION

.

(

EQUATION

It succinctly reproduces the result that can be obtained through planar diagrams in this special setup (Dyer & Gur-Ari, 2019) . Quadratic activation (Li et al., 2018) is worked out in Appendix C.1.

4.2. RELU WITH SINGLE INPUT:

The recursion relations simplify drastically for the case of a single input, N D = 1, as worked out in detail in Appendix C.2. For instance, for rectified linear unit (ReLU) activation with C

( ) b = 0 and C ( ) W = 2, we obtain the layer-independent core kernel, zero self-energy, and the four-point vertex

EQUATION

EQUATION

)

Interestingly, as for deep linear networks, the factor (1/n ) appears again. This factor has also been found by Hanin & Rolnick (2018) , which provides guidance for network architectural design through its minimization. We generalize this factor for monomial activations in Appendix C.2.1

4.3. EXPERIMENTAL VERIFICATION: OUTPUT DISTRIBUTIONS FOR A SINGLE INPUT:

Here we put our theory to the test. For concreteness, take a single black-white image of handwritten digits with 28-by-28 pixels (i.e. n 0 = 784) from the MNIST dataset (LeCun et al., 1998) without preprocessing, set depth L = 3, bias variance C ( ) b = 0, weight variance C ( ) W = C W , and widths (n 0 , n 1 , n 2 , n 3 ) = (784, n, 2n, 1), and use activations σ(z) = z (linear) with C W = 1 and max(0, z) (ReLU) with C W = 2. In Figure 1 , for each width-parameter n of the hidden layers we record the prior distribution of outputs over 10 6 instances of Gaussian weights and compare it with the theoretical prediction -obtained by cranking the knob from the initial condition (R0) through the recursion relations (R1-R3) to the distribution (D0-D2). The prior distribution becomes increasingly non-Gaussian as networks narrow and the deviation from the Gaussian-process prior is correctly captured by our theory. Higher-order perturbative calculations are expected to systematically improve the quality -and extend the range -of the agreement. Additional experiments are performed in Appendix C.3, which further corroborates our theory.

Figure 1 : Comparison between theory and experiments for prior distributions of outputs for a single input. The agreement between our theoretical predictions (smooth thick lines) and experimental data (rugged thin lines) is superb, correctly capturing the initial deviations from Gaussian processes at n = ∞ (black), all the way down to n ∼ 10 for linear activation and to n ∼ 30 for ReLU activation.

5. BAYESIAN INFERENCE:

Let us take off from the terminal point of Section 3: we have obtained the recursive equations (R0-R3) for the Gaussian-process kernel and the leading finite-width corrections and codified them in the weakly non-Gaussian prior distributions p[z] (D0-D2) of outputs z ≡ β i } β=1,...,NR;i=1,...,n L . We shall develop a formalism to infer outputs for test inputs a lá Bayes, perturbatively extending the textbook by Williams & Rasmussen (2006) . For field theorists, our calculation is just a tree-level background field calculation (Weinberg, 1996) in disguise.

z α i ≡ z (L) i (x α ) i=1,...,

Taking the liberty of notations, we let the number of input-data arguments dictate the summation over sample indices α inside the potential H, and denote the joint probabilities

EQUATION

)

Given the training targets y R , the posterior distribution of test outputs are given by Bayes' rule:

p [z E |y R ] = p[y R , z E ] p[y R ] = dz R e -H[z R ] dz R dz E e -H[z R ,z E ] e -(H[yR,zE]-H[yR]) .

(Bayes)

The leading Gaussian-process contributions can be segregated out through the textbook manipulation (Williams & Rasmussen, 2006) [c.f. Appendix D.1]: denoting the full Gaussian-process kernel in the last layer as

EQUATION

and the Gaussian-process posterior mean prediction as

y GP E γ i ≡ β K ER K -1 RR γ β (y R ) β i , (GPM)

and defining a fluctuation (z E )

γ i ≡ y GP E γ i + (δz E ) γ i and a matrix K ∆ ≡ K EE -K ER K -1 RR K RE , H 0 [y R , z E ] -H 0 [z R ] = 1 2 i γ1,γ2 (δz E ) γ1 i K -1 ∆ γ1,γ2 (δz E ) γ2 i . (GP∆)

For any function F, its expectation over the Bayesian posterior (Bayes) then turns into

EQUATION

where the deviation kernel

(δz E ) γ1 i1 (δz E ) γ2 i2 K∆ ≡ δ i1,i2 K γ1,γ2

∆ and the normalization factor

EQUATION

In particular the mean posterior output is given by

(y E ) γ i ≡ dz E (z E ) γ i p [z E |y R ] = y GP E γ i + N (δz E ) γ i e -H1[yR,y GP E +δzE] K∆ (22) = y GP E γ i -(δz E ) γ i H 1 y R , y GP E + δz E K∆ + O( 2 ) .

Stringing together

φ α i ≡ [(y R ) β i , y GP E γ i ],

recalling Equation (D2) for H 1 , and using Wick's contractions for one last time, the mean prediction becomes

y GP E γ i (NGPM) + α1,γ1 K γ,γ1 ∆ φ α1 i S γ1,α1 - α2,α3 V (γ1,α2),(α1,α3) K α2,α3 + γ2,γ3 V (γ1,γ2),(α1,γ3) K γ2,γ3 ∆ + n L 2 γ2,γ3 V (α1,γ1),(γ2,γ3) K γ2,γ3 ∆ + α2,α3 V (γ1,α1),(α2,α3)   - n L 2 K α2,α3 + 1 2 j φ α2 j φ α3 j   .

With additional manipulations in Appendix D, this expression is simplified into the actionable form that is amenable to use in practice. For illustration, there, we also show a simple preliminary experiment, which indicates the 1/n regularization effect for sufficiently large width n and small amount of training data N R . This is in line with expectations that finite widths ameliorate overfitting and that non-Gaussian priors increase the expressivity of neural functions, but additional large-scale extensive experiments would be desirable in the future.

6. DREAMS:

In this paper, we have developed the perturbative formalism that captures the flow of preactivation distributions from lower to higher layers. The spiritual resemblance between our recursive equations and renormalization-group flow equations in high-energy and statistical physics is highly appealing. It would be exciting to investigate the structure of fixed points away from the Gaussian asymptopia (Schoenholz et al., 2016) and fully realize the dream articulated by Mehta & Schwab (2014) beyond their limited example of a mapping between two antiquated techniques -the audacious hypothesis that neural networks wash away microscopic irrelevancies and extract relevant features.

In addition we have developed the perturbative Bayesian inference scheme universally applicable whenever prior distributions are weakly non-Gaussian, and have applied it to the specific cases of neural networks at finite widths. In light of finite-width regularization effects, it would be prudent to revisit the empirical comparison between SGD optimization and Bayesian inference at finite widths (Lee et al., 2018; Novak et al., 2019) , especially for convolutional neural networks.

Finally, given surging interests in SGD dynamics within the large-width regime (Jacot et al., 2018; Chizat et al., 2018; Lee et al., 2019; Cohen et al., 2019; Dyer & Gur-Ari, 2019) , it would be natural to adapt our formalism for investigating corrections to neural tangent kernels, and even nonperturbatively aspire to capture a phase transition out of lazy-learning into feature-learning regimes.

A WICKED TRICKS:

Welcome to the magical school of Wick's crafts and wizardly cumulants. Here is all you need to know in order to follow the calculations in the paper. In the main text, Wick's contractions are used trivially for integrating out biases and weights as straightforward applications of Appendix A.1 while they are used more nontrivially for integrating out preactivations, with concepts of cumulants reviewed in Appendix A.2 and A.3, culminating in a hack derived in Appendix A.4. For most parts, we shall forget about the neuron index i for pedagogy and put them back in at the very end.

A.1 WICK'S CONTRACTIONS:

For Gaussian-distributed variables z = {z α } α=1,...,N with a kernel K α,α , moments

z α1 z α2 • • • z αm K ≡ dze -H0[z] z α1 z α2 • • • z αm dze -H0[z] with H 0 [z] ≡ 1 2 N α,α =1 z α K -1 α,α z α .

(S1) For any odd m such moments identically vanish. For even m, Isserlis-Wick's theorem states that

z α1 z α2 • • • z αm K = all pairing K α k 1 ,α k 2 • • • K α k m-1 ,α km (S2)

where the sum is over all the possible pairings of m variables,

(k 1 , k 2 ), . . . , (k m-1 , k m ). In general, there are (m -1)!! = (m -1) • (m -3) • • • 1 such pairings.

For a proof, see for example Zee (2010) .

In order to understand and use the theorem, it is instructive to look at a few examples:

EQUATION

and

z α1 z α2 z α3 z α4 z α5 z α6 K (S5) = K α1,α2 K α3,α4 K α5,α6 + K α1,α3 K α2,α4 K α5,α6 + K α1,α4 K α2,α3 K α5,α6 + K α1,α2 K α3,α5 K α4,α6 + K α1,α3 K α2,α5 K α4,α6 + K α1,α5 K α2,α3 K α4,α6 + K α1,α2 K α5,α4 K α3,α6 + K α1,α5 K α2,α4 K α3,α6 + K α1,α4 K α2,α5 K α3,α6 + K α1,α5 K α3,α4 K α2,α6 + K α1,α3 K α5,α4 K α2,α6 + K α1,α4 K α5,α3 K α2,α6 + K α5,α2 K α3,α4 K α1,α6 + K α5,α3 K α2,α4 K α1,α6 + K α5,α4 K α2,α3 K α1,α6 .

A.2 CONNECTED CORRELATIONS Given general (not necessarily Gaussian) random variables, connected correlation functions are defined inductively through

E [z α1 z α2 • • • z αm ] (S6) ≡ E [z α1 z α2 • • • z αm ] connected + all subdivisions E z α k [1] 1 • • • z α k [1] ν 1 connected • • • E z α k [s] 1 • • • z α k [s] νs connected

where the sum is over all the possible subdivisions of m variables into s > 1 clusters of sizes (ν 1 , . . . , ν s ) as (k

[1] 1 , . . . , k [1] ν1 ), . . . , (k [s] 1 , . . . , k [s]

νs ). In order to understand the definition, it is again instructive to look at a few examples. Assuming that all the odd moments vanish,

E [z α1 z α2 ] = E [z α1 z α2 ] connected and (S7) E [z α1 z α2 z α3 z α4 ] = E [z α1 z α2 z α3 z α4 ] connected (S8) +E [z α1 z α2 ] connected E [z α3 z α4 ] connected +E [z α1 z α3 ] connected E [z α2 z α4 ] connected +E [z α1 z α4 ] connected E [z α2 z α3 ] connected .

Rearranging them in particular yields

E [z α1 z α2 z α3 z α4 ] connected (S9) = E [z α1 z α2 z α3 z α4 ] -E [z α1 z α2 ] E [z α3 z α4 ] -E [z α1 z α3 ] E [z α2 z α4 ] -E [z α1 z α4 ] E [z α2 z α3 ] .

If these examples do not suffice, here is yet another example to chew on:

E [z α1 z α2 z α3 z α4 z α5 z α6 ] = E [z α1 z α2 z α3 z α4 z α5 z α6 ] connected (S10) +E [z α1 z α2 ] connected E [z α3 z α4 ] connected E [z α5 z α6 ] connected + [14 other (2, 2, 2) subdivisions] +E [z α1 z α2 z α3 z α4 ] connected E [z α5 z α6 ] connected + [14 other (4, 2) subdivisions]

and hence

E [z α1 z α2 z α3 z α4 z α5 z α6 ] connected (S11) = E [z α1 z α2 z α3 z α4 z α5 z α6 ] -{E [z α1 z α2 z α3 z α4 ] E [z α5 z α6 ] + [14 other (4, 2) subdivisions]} +2 {E [z α1 z α2 ] E [z α3 z α4 ] E [z α5 z α6 ] + [14 other (2, 2, 2) subdivisions]} .

We emphasize that these are just renderings of the definition (S6). The power of this definition will be illustrated in the next two subsections.

A.3 HIERARCHICAL CLUSTERING:

We often encounter situations with the hierarchy

E [z α1 • • • z αm ] connected = O( m-2 2 ) (S12)

where 1 is a small perturbative parameter and here again odd moments are assumed to vanish. Often comes with the hierarchical structure is the asymptotic limit → 0 where

E [z α1 z α2 ] = K α1,α2 + S α1,α2 + O( 2 ) (S13)

with the Gaussian kernel K α1,α2 at zero and the leading self-energy correction S α1,α2 . Let us also denote the leading four-point vertex

E [z α1 z α2 z α3 z α4 ] connected = V α1,α2,α3,α4 + O( 2 ) . (S14)

For instance this hierarchy holds for weakly-coupled field theories -from which we are importing names such as self-energy, vertex, and metric -and, in this paper, such hierarchical structure is inductively shown to hold for prior preactivations z ( ) with = 1 n -1 in the regime n 1 , . . . , n L-1 ∼ n 1. Note that, by definition, K α1,α2 and S α1,α2 are symmetric under α 1 ↔ α 2 and V α1,α2,α3,α4 is symmetric under permutations of (α 1 , α 2 , α 3 , α 4 ). 6With the review of connected correlation functions passed us, we can now readily see that

E [z α1 • • • z αm ] (CLUSTER) =E [z α1 z α2 ] • • • E [z αm-1 z αm ] + {[(m -1)!! -1] other pairings} + E [z α1 z α2 z α3 z α4 ] connected E [z α5 z α6 ] • • • E [z αm-1 z αm ] + m 4 × (m -5)!! -1 other (4, 2, 2, . . . , 2) clusterings + O( 2 ) . =K α1,α2 • • • K αm-1,αm + {[(m -1)!! -1] other pairings} + S α1,α2 K α3,α4 • • • K αm-1,αm + m 2 × (m -3)!! -1 other such clusterings + V α1,α2,α3,α4 K α5,α6 • • • K αm-1,αm + m 4 × (m -5)!! -1 other (4, 2, 2, . . . , 2) clusterings + O( 2 ) . = z α1 • • • z αm K + S α1,α2 z α3 • • • z αm K + m 2 -1 other self-energy contractions + V α1,α2,α3,α4 z α5 • • • z αm K + m 4 -1 other vertex contractions + O( 2 ) .

where in the last equality Wick's theorem was used backward.

A.4 COMBINATORIAL HACK So far we have reviewed the standard technology of Wick's contractions, connected correlation functions, and all that. Here is one sorcery, which lets us magically evaluate E [z α1 • • • z αm ] by mindless repetitions of Wick's contractions. Throughout, we shall assume the hierarchical structure (S12), 7 which is the inductive assumption in the main text, and start from its consequence (CLUSTER).

Below, let us use the inverse kernel K -1 α1,α2 as a metric to lower indices:

S α1,α2 ≡ α 1 ,α 2 K -1 α1,α 1 K -1 α2,α 2 S α 1 ,α 2 and (S15) V α1,α2,α3,α4 ≡ α 1 ,...,α 4 K -1 α1,α 1 • • • K -1 α4,α 4 V α 1 ,α 2 ,α 3 ,α 4 . (S16) First note that z α1 • • • z αm   α 1 ,α 2 S α 1 ,α 2 z α 1 z α 2   K = α 1 ,α 2 S α 1 ,α 2 z α 1 z α 2 K z α1 • • • z αm K +2 z α1 z α 1 K z α2 z α 2 K z α3 • • • z αm K + m 2 -1 other (α 1 , α 2 ) =   α 1 ,α 2 S α 1 ,α 2 K α 1 ,α 2   z α1 • • • z αm K +2S α1,α2 z α3 • • • z αm K + m 2 -1 other (α 1 , α 2 )

7 More precisely, we shall only use the weaker assumption that 2) for m ≥ 6 along with Equations (S13) and (S14).

E [z α 1 • • • z αm ] connected = O(

where the symmetry α 1 ↔ α 2 of S α1,α2 was used. Hence, defining

O S [z] ≡ 1 2 α 1 ,α 2 S α 1 ,α 2 z α 1 z α 2 -K α 1 ,α 2 , (OS')

we obtain, for one term in Equation (CLUSTER),

S α1,α2 z α3 • • • z αm K + m 2 -1 other (α 1 , α 2 ) (S17) = z α1 • • • z αm O S [z] K . (S18)

The similar algebraic exercise renders the other term in Equation (CLUSTER) to be

V α1,α2,α3,α4 z α5 • • • z αm K + m 4 -1 other (α 1 , α 2 , α 3 , α 4 ) (S19) = z α1 • • • z αm O V [z] K (S20) with O V [z] ≡ 1 24 α 1 ,...,α 4 V α 1 ,α 2 ,α 3 ,α 4 z α 1 z α 2 z α 3 z α 4 -6z α 1 z α 2 K α 3 ,α 4 + 3K α 1 ,α 2 K α 3 ,α 4 . (OV')

In summary, for any function

F[z] of random variables z α E {F[z]} = F[z] K + F[z]O S [z] K + F[z]O V [z] K + O( 2 ) . (HACK')

In order to get the expressions used in the main text at the -th layer, we need only to put back neuron indices i by replacing α → (α, i), identify = 1 n -1 , and use the inductive assumptions (KS)

E z α1 i1 z α2 i2 = δ i1,i2 K α1,α2 ( ) + 1 n -1 S α1,α2 ( ) + O 1 n 2 and (V) E z α1 i1 z α2 i2 z α3 i3 z α4 i4 connected = 1 n -1 δ i1,i2 δ i3,i4 V (α1,α2),(α3,α4) ( ) + δ i1,i3 δ i2,i4 V (α1,α3),(α2,α4) ( ) +δ i1,i4 δ i2,i3 V (α1,α4),(α2,α3) ( ) + O 1 n 2 .

The operators in Equations (OS') and (OV') then become

EQUATION

i.e., the operators in Equations (OS) and (OV) in the main text.

B FULL CONDENSED PROOF:

In this Appendix, we provide a full inductive proof for one of the main claims in the paper, streamlined in the main text. Namely, we assume at the -th layer that Equations (KS) and (V) hold and that all the higher-point connected preactivation correlation functions are of order O 1 n 2 -which are trivially true at = 1 -and prove the same for the ( + 1)-th layer. We assume the full mastery of Appendix A or, conversely, this section can be used to test the mastery of wicked tricks. 

EQUATION

Studiously disentangling cases with different numbers of repetitions in neuron indices (j 1 , . . . , j k ), we notice that at order O 1 n , terms without repetition or with only one repetition contribute, finding

1 n k n j1,...,j k =1 H ( )α1,α2,...,α 2k-1 ,α 2k j1,j1,...,j k ,j k (S22) = σ(z α1 )σ(z α2 ) K ( ) • • • σ(z α 2k-1 )σ(z α 2k ) K ( ) + 1 n σ(z α1 )σ(z α2 )σ(z α3 )σ(z α4 ) K ( ) -σ(z α1 )σ(z α2 ) K ( ) σ(z α3 )σ(z α4 ) K ( ) × σ(z α5 )σ(z α6 ) K ( ) • • • σ(z α 2k-1 )σ(z α 2k ) K ( ) + k 2 -1 others + 1 n -1 σ(z α1 1 )σ(z α2 1 ) • • • σ(z α 2k-1 k )σ(z α 2k k ) {O S [z] + O V [z]} K ( ) +O 1 n 2

where we used the inductive hierarchical assumption at the -th layer, i.e., its consequence (HACK) and denoted a single-neuron random vector z( ) = {z α } α=1,...,ND and the Gaussian integral with the core kernel zα1 zα2

K ( ) = K α1,α2 ( )

. Plugging in expressions (OS,OV) for operators O S,V [z],

σ(z α1 1 )σ(z α2 1 ) • • • σ(z This expression simplifies drastically through the identity

α 2k-1 k )σ(z α 2k k ) O S [z] K ( ) (S23) = 1 2 α 1 ,α 2 S α 1 ,α 2 σ(z α1 )σ(z α2 ) zα 1 zα 2 -K α 1 ,α 2 K ( ) × σ(z α3 )σ(z α4 ) K ( ) • • • σ(z α 2k-1 )σ(z α 2k ) K ( ) + [(k -1) others]

EQUATION

which can be checked explicitly, recalling K ∆ ≡ K EE -K ER K -1 RR K RE . Incidentally, this identity can also be used to prove Equation (GP∆). Now equipped with this identity, recalling φ 

:

The nontrivial parts of the inductive proof for Equations (KS) and (V) are to show (i) that the right-hand side of Equation (10) is finite as n → ∞, (ii) that the leading contribution of Equation (9) is the Gaussianprocess kernel, and (iii) that higher-point connected preactivation correlation functions are all suppressed by O 1 n 2 , all of which are verified in obtaining the recursive equations. See Appendix B for a full proof.

If nL were of order n 1, the potential H would become a large-n vector model, for which we would have to sum the infinite series of bubble diagrams(Moshe & Zinn-Justin, 2003).

In the main text the connected four-point preactivation correlation functions are symmetric under the permutations of four (sample, neuron) indices, {(α1, i1), (α2, i2), (α3, i3), (α4, i4)}.

Figures and Tables:



FigureS1: Comparison between theory and experiments for prior distributions of outputs for a single input. Our theoretical predictions (smooth thick lines) and experimental data (rugged thin lines) agree, correctly capturing the initial deviations from the Gaussian processes (black, n = ∞), at least down to n = n with n ∼ 10 for linear cases, n ∼ 30 for ReLU cases and depth L = 2 quadratic case, and n ∼ 100 for depth L = 3 quadratic case. This also illustrates that nonlinear activations quickly amplify non-Gaussianity.

we notice that φ i,β1 = β K -1 e. with one or both of training components (β 2 , β 3 ) replaced by test components γ] vanish. Equation the mean prediction. Equations (NGPM') and (NGPM") withφ i,β1 = β K -1 RR β1,β (y R ) β iare actionable, i.e., easy to program.

ND , consisting of N R training inputs with target outputs and N E test inputs with unknown targets, we can perform exact Bayesian inference to yield mean outputs as predictions for N E test data (Williams

Table content: <table/>

Let us now divide N D inputs into N R training and N E test inputs as {x α } α=1,...,ND = {x β

Table content: <table><tr><td>with ≡ 1 n L-1</td><td colspan="2">1. Examples in Section 4 illustrate that finite-width corrections stay perturba-</td></tr><tr><td colspan="2">tive typically when depth width</td><td>1.</td></tr></table>

Peer review (iclr2020):
The paper proposes a method based on Wick*s contractions that can integrate out intermediate layers of a Bayesian neural network with finite width. Similar to the connection between GP and infinite-width NN, the marginalization of intermediate layers could allow Bayesian inference for the network. This is a very technical paper with mainly the derivation of the recursion relations and the posterior mean for Bayesian inference. The math is tantalizing. But the main concern is how effective it*s gonna be when being used in a neural network inference problem. The simulated example shows that there could be a mismatch between theoretical predictions and experimental data when using ReLU. How much effect it*s gonna have in a prediction task? The recursion relations need to be evaluated numerically for most of cases. Then what*s the computational difference between the proposed method and general MC based integration? I would suggest the authors move some derivations to the appendix and add more experimental analyses and discussion about the model. Like what*s the trade-off between approximation accuracy, running time, predictive performance etc, like other papers did. The current format is not mature enough to get published I think.