# ON THE IMPLICIT BIASES OF ARCHITECTURE & GRADIENT DESCENT

**Anonymous authors**
Paper under double-blind review

ABSTRACT

Do neural networks generalise because of bias in the functions returned by gradient
descent, or bias already present in the network architecture?

_¿Por qué no los dos?_

This paper finds that while typical networks that fit the training data already generalise fairly well, gradient descent can further improve generalisation by selecting
networks with a large margin. This conclusion is based on a careful study of the
behaviour of infinite width networks trained by Bayesian inference and finite width
networks trained by gradient descent. To measure the implicit bias of architecture,
new technical tools are developed to both analytically bound and consistently
_estimate the average test error of the neural network–Gaussian process (NNGP)_
posterior. This error is found to be already better than chance, corroborating the
findings of Valle-Pérez et al. (2019) and underscoring the importance of architecture. Going beyond this result, this paper finds that test performance can be
substantially improved by selecting a function with much larger margin than is typical under the NNGP posterior. This highlights a curious fact: minimum a posteriori
functions can generalise best, and gradient descent can select for those functions. In
summary, new technical tools suggest a nuanced portrait of generalisation involving
both the implicit biases of architecture and gradient descent.

1 INTRODUCTION

Following an influential paper by Zhang et al. (2017), the basic question of why neural networks
_generalise is generally regarded to be open. The authors demonstrate a surprising fact: deep learning_
generalises even when the neural network is expressive enough to represent functions that do not
generalise. In turn, this implies that the theory of Vapnik & Chervonenkis (1971)—based on uniform
convergence of train error to population error—does not explain generalisation in neural networks.

Several hypotheses have since been proposed to fill this theoretical vacuum. One prominent hypothesis
states that, while most neural network solutions do not generalise well, there is an implicit bias in
the kinds of functions returned by gradient descent (Soudry et al., 2018). In sharp contrast, a second
hypothesis states that the solution space of a neural network is dominated by simple functions, while
the complex kinds of functions that overfit are relatively rare (Valle-Pérez et al., 2019).

This latter hypothesis dovetails with a particular “PAC-Bayesian” theorem of McAllester (1998),
which bounds the average population error of all classifiers consistent with a training sample. If this
bound is small, then the complex functions that overfit must indeed be rare. In more technical terms,
the PAC-Bayesian theorem can provide a meaningful certificate of generalisation even for machine
learning models with an infinite Vapnik-Chervonenkis dimension, or an arbitrarily large number of
parameters, by properly accounting for the measure of those very complex functions.

This paper takes a more nuanced position between these two hypotheses. While the average population error of all neural networks that fit a training sample is found to be good (and certifiably
_good by the PAC-Bayesian theorem), it is still possible for certain networks with special properties_
to perform substantially better than average (and likewise substantially better than the PAC-Bayes
bound). Moreover, gradient descent may be used to specifically target these special networks. This
subtly counters a position put forward by Mingard et al. (2021), which suggests that gradient descent
may be well-modelled as sampling randomly from a particular Bayesian posterior distribution.


-----

To support these claims, a careful study of the behaviour of both infinite width and finite width neural
networks is conducted. In particular, this paper makes the following technical contributions:

Section 3 Implicit Bias of Architecture. A purely analytical PAC-Bayes bound (Theorem 2) on the
population error of the neural network–Gaussian process (NNGP) posterior for binary
classification is derived. The bound furnishes an interpretable measure of model complexity
that depends on both the architecture and the training data. This exact analytical bound
improves upon an approximate computational approach due to Valle-Pérez et al. (2019).

Section 4 Testing the Bound. The new bound is found to be both non-vacuous and correlated with the
test error of finite width multilayer perceptrons trained by gradient descent. This provides
supporting evidence for the important role of architecture in generalisation, as put forward
by Valle-Pérez et al. (2019). Still, a gap exists with the performance of gradient descent.

Section 5 Implicit Bias of Gradient Descent. Going further beyond the work of Valle-Pérez et al.
(2019), a new theoretical tool (Theorem 3) is developed to enable consistent estimation
of the average error of the NNGP posterior on a given holdout set. The average is found
to be significantly worse than the holdout performance of Gaussian process draws with
large margin. The experiment is repeated for finite width neural networks trained by
gradient descent, and the same qualitative phenomenon persists. This finding demonstrates
the ability of gradient descent to select large margin functions with vanishing posterior
probability that nonetheless generalise significantly better than the posterior average.

2 RELATED WORK

**Margin-based generalisation theory** A rich body of work explores the connection between margin
and generalisation. For instance, Bartlett et al. (2017) propose a margin-based complexity measure
for neural networks derived via Rademacher complexity analysis, and Neyshabur et al. (2018) derive
a similar result via PAC-Bayes analysis. Neyshabur et al. (2019) test these bounds experimentally,
finding them to be vacuous and to scale poorly with network width. Biggs & Guedj (2021) further
develop this style of bound. Margin-based PAC-Bayes bounds go back at least to the work of Herbrich
(2001), who derived such a bound for linear classification. The standard idea is that a solution with
large margin implies the existence of nearby solutions with the same training error (but perhaps
smaller margin), facilitating a more spread out PAC-Bayes posterior. This style of margin-based
PAC-Bayes bound does not provide guidance on whether the original large margin classifier should
generalise better or worse than the additional nearby classifiers included in the posterior.

**Non-vacuous bounds for neural networks** Aside from PAC-Bayes, many styles of generalisation
bound for neural networks are vacuous (Dziugaite & Roy, 2017). Even many PAC-Bayes bounds
are vacuous (Achille & Soatto, 2018; Foret et al., 2021) due to their choice of PAC-Bayes prior.
Dziugaite & Roy (2017) construct a non-vacuous weight space PAC-Bayes bound by iteratively
optimising both the PAC-Bayes prior and posterior by gradient descent. Wu et al. (2020) extend
this approach to involve Hessians. Meanwhile, Valle-Pérez et al. (2019) instantiate a non-vacuous
PAC-Bayes bound in the function space of neural networks via the NNGP correspondence (Neal,
1994; Lee et al., 2018). The evaluation of this bound involves an iterative and statistically inconsistent
expectation-propagation approximation. Unlike this paper’s purely analytical bound in Theorem 2,
these non-vacuous bounds are all evaluated by iterative computational procedures.

**Implicit bias of gradient descent** Much research has gone into the role that gradient descent plays
in neural network generalisation. For instance, Zhang et al. (2017) suggest that gradient descent may
converge to solutions with special generalisation properties, while Wilson et al. (2017) suggest that
gradient descent may have a more favourable implicit bias than the Adam optimiser. Azizan et al.
(2021) make a related study in the context of mirror descent. Keskar et al. (2017) point to an implicit
bias present in mini-batch stochastic gradient descent, but Geiping et al. (2021) argue against the
importance of stochasticity. Meanwhile, Soudry et al. (2018) observe how gradient descent combined
with certain loss functions converges to max-margin solutions; Wei et al. (2018) make a similar
observation. On a different tack, Mingard et al. (2021), Valle-Pérez & Louis (2020) and Valle-Pérez
et al. (2019) argue that the implicit bias of neural architecture is more important than that of gradient
descent for understanding generalisation, and so long as gradient descent does not select solutions
pathologically, generalisation should occur.


-----

_f_ (x ; ·)

**solution space VS**

volume P[VS]

_f_ (x[(1)]; ·

**function space FS**

volume P[FS] = 1

**weight space Ω**

volume P[Ω] = 1

**solution space VS**

volume P[VS]

Figure 1: The solution space VS of a learning task denotes the subset of classifiers that attain zero
error on a dataset S. While the solution space can have a complicated geometry when described in
weight space, in function space the solution space of a binary classification task is just an orthant.

3 IMPLICIT BIAS OF ARCHITECTURE

This section derives an analytical bound on the population error of an infinitely wide neural network
averaged over all weight configurations that attain 0% train error—in other words, the average
population error of the NNGP posterior. Since this quantity assesses the performance of all solutions
rather than just those returned by gradient descent, it is intended to measure the implicit bias of
architecture. In Section 4, the bound is found to be non-vacuous for multilayer perceptrons. The
source of a substantial gap with the results of gradient descent training is investigated in Section 5,
and ultimately attributed to an additional implicit bias of gradient descent: namely margin.

3.1 PRELIMINARY NOTATION

Consider a training dataset S of n input–label pairs: S = {(x[(1)], y[(1)]), ..., (x[(][n][)], y[(][n][)])}. The inputs
_x[(][i][)]_ _∈_ R[d][0] are embedded in Euclidean d0-space, while to simplify the analysis the paper shall restrict
to binary labels y[(][i][)] _∈{±1}. A classifier f shall depend on a weight vector w ∈_ Ω, where Ω denotes
the weight space. In particular, for an input x the prediction is given by sign f (x; w).

One is interested in the relationship between the train error εS of a classifier on a dataset S to the
_population error εD on a data distribution D. For a weight vector w ∈_ Ω, these are defined as:

_n_

_εS(w) := [1]_ I sign f (x[(][i][)], w) = y[(][i][)][i] ; _εD(w) := P(x,y)_ _D [sign f_ (x; w) = y] . (1)

_n_ _̸_ _∼_ _̸_

_i=1_

X h

The solution space VS := {w ∈ Ω _| εS(w) = 0} denotes those classifiers that attain zero train error._


3.2 GEOMETRY OF SOLUTIONS IN FUNCTION SPACE

The function space FS is defined to be the set of outputs that the classifier f can realise on dataset S:

_FS :=_ _f_ (x[(1)], w), ..., f (x[(][n][)], w) _w ∈_ Ω _⊂_ R[n]. (2)
n  o

In weight space the geometry of solutions can be arbitrarily complicated, but in function space the
solutions occupy the single orthant picked out by the binary training labels.

To measure the volume of the solution space, one can define a measure P on weight space. It
is convenient to enforce P[Ω] = 1 so that P is a probability measure. In the absence of a better
alternative, P is usually set to a multivariate Gaussian or uniform distribution. The measure P on
weight space induces a measure P on function space via the relation:

P[F ⊂FS] := P _w ∈_ Ω _f_ (x[(1)], w), ..., f (x[(][n][)], w) _∈_ _F_ _._ (3)
n   o

The volume of solutions is denoted P[VS], and can be computed either in weight or function space.
In function space, P[VS] is an orthant probability. The situation is visualised in Figure 1.


-----

3.3 PAC-BAYES THEORY

PAC-Bayes relates the volume of solutions to their average population error. The following result
was derived by Valle-Pérez et al. (2019) as a corollary of a theorem due to Langford & Seeger (2001).
The result is similar in form to a theorem of McAllester (1998).
**Theorem 1 (A PAC-Bayesian theorem). First, fix a probability measure P over the weight space Ω**
_of a classifier. Let S denote a training set of n datapoints sampled iid from the data distribution D_
_and let VS denote the corresponding solution space. Consider the population error 0 ≤_ _εD(w) ≤_ 1
_of weight setting w ∈_ Ω, and its average over the solution space εD(VS) := Ew∼P[εD(w) | w ∈ _VS]._
_Then, for a proportion 1 −_ _δ of draws of the training set S,_

1 ln P[V1S ] [+ ln][ 2]δ[n]
_εD(VS)_ ln _._ (4)
_≤_ 1 _εD(VS)_ _n_ 1

_−_ _[≤]_ _−_

For large n, the ln 2n/δ term is negligible and the result says that the population error averaged over
solutions is less than the ratio of the negative log volume of solutions to the number of training points.

3.4 VOLUME OF SOLUTIONS VIA GAUSSIAN ORTHANT PROBABILITIES

Since infinitely wide neural networks induce a Gaussian measure on function space, computing the
volume of the solution space P[VS] amounts to computing a Gaussian orthant probability.

In more detail, let P denote a measure on the weight space Ω of an infinitely wide neural network
satisfying the conditions of the NNGP correspondence—see Theorem 4 in Appendix A.3 for an
example. Then for a weight vector w ∼ P, the network outputs on a training set S are distributed:

_f_ (x[(1)]; w), ..., f (x[(][n][)]; w) ∼N (0, Σ), (5)

with covariance Σij := Ew∼P[f (x[(][i][)]; w)f (x[(][j][)]; w)].

Therefore, under the NNGP correspondence, the volume of solutions computed in function space is
just the Gaussian probability of the orthant picked out by the binary training labels:

P[VS] = Pϕ (0,Σ)[sign ϕ1 = y[(1)], ..., sign ϕn = y[(][n][)]]. (6)
_∼N_
To facilitate estimating and bounding this probability, this paper has derived the following lemma.
**Lemma 1 (Gaussian orthant probability). For a covariance matrix Σ ∈** R[n][×][n], and a binary vector
_y ∈{±1}[n], let p denote the corresponding Gaussian orthant probability:_
_p := Pϕ_ (0,Σ)[sign(ϕ) = y]. (7)
_∼N_

_Letting I denote the n_ _×_ _n identity matrix, ⊙_ _the elementwise product and |·| the elementwise absolute_
_value, then p may be equivalently expressed as:_

_n_

_p = [1]_ e[−] 2[1] [(][y][⊙|][u][|][)][T] [(] _√det ΣΣ[−][1]−I)(y⊙|u|)[i]_ =: e[−C][0][(Σ][,y][)], (8)

2[n][ E][u][∼N][ (0][,][I][)]

_and p may be bounded as follows:h_


_p_ _n2_ _[−]_ _√[n]_
_≥_ 2[1][n][ e]


det Σ[( [1]2 _[−]_ _π[1]_ [) tr(Σ][−][1][)+][ 1]π _[y][T][ Σ][−][1][y][] ]=: [1]_ _n2_ _[−C][1][(Σ][,y][)]._ (9)

2[n][ e]


The proof is given in Appendix A.2. Equation 8 yields an unbiased Monte Carlo estimator of Gaussian
orthant probabilities, and Inequality 9 yields a lower bound. The complexity measures C0 and C1 are
defined for later use.

To gain intuition about the lemma, observe that 1/2[n] is the orthant probability for an isotropic
_n_
Gaussian. Depending on the degree of anisotropy _√det ΣΣ[−][1]_ _−_ I inherent in the covariance

matrix Σ, Equation 8 captures how the orthant probability may either be exponentially amplified or
suppressed compared to 1/2[n].

As an aside, using Inequality 9 to lower bound Equation 6 has a Bayesian interpretation. Since
the volume of solutions may be written P[VS] = Ω [d][P][(][w][)][ I][[][w][ ∈] _[V][S][]][, it may be interpreted as the]_

_Bayesian evidence for the network architecture under the likelihood function I[w_ _VS]. A Bayesian_

R _∈_

would then refer to Inequality 9 as an evidence lower bound.

The following generalisation bound is a basic consequence of Theorem 1, Equation 6 and Lemma 1:


-----

**Theorem 2 (Upper bound on the average population error of an infinitely wide neural network).**
_First, fix a probability measure P over the weight space Ω_ _of an infinitely wide neural network. Let S_
_denote a training set of n datapoints sampled iid from the data distribution D, let y ∈{±1}[n]_ _denote_
_the binary vector of training labels, and let VS denote the corresponding solution space. Consider the_
_population error 0_ _εD(w)_ 1 of weight setting w Ω, and its average over the solution space
_≤_ _≤_ _∈_
_εD(VS) := Ew∼P[εD(w) | w ∈_ _VS]. Let Σ denote the NNGP covariance matrix (Equation 5). Then,_
_for a proportion 1 −_ _δ of draws of the training set S,_

1 _δ_ _n5_ [+][ C][1][(Σ][, y][) + ln][ 2]δ[n]
_εD(VS)_ ln _,_ (10)
_≤_ 1 _εD(VS)_ _n_ 1 _≤_ _n_ 1

_−_ _[≤C][0][(Σ][, y] −[) + ln][ 2][n]_ _−_

_where the complexity measures C0 and C1 are defined in Lemma 1._

Since the complexity measure C1 is an analytical function of the NNGP covariance matrix Σ and the
binary vector of training labels y, Theorem 2 is an analytical generalisation bound for the NNGP
posterior. For large n, the result simplifies to: εD(VS) ⪅ _C0(Σ, y)/n ≤_ 1/5 + C1(Σ, y)/n. So the
bound depends on the ratio of complexity measures C0 and C1 to the number of datapoints n.

To gain further intuition about Theorem 2, consider two special cases. First, suppose that the neural
architecture induces no correlation between any pair of distinct data points, such that Σ = I. Then
_εD(VS) ⪅_ _C0(I, y)/n = ln 2 ≈_ 0.7 and the bound is worse than chance. This corresponds to
_pure memorisation of the training labels. Next, suppose that the neural architecture induces strong_
intra-class correlations and strong inter-class anti-correlations, such that Σij = yiyj. Although this Σ
is singular, it may be seen directly that P[VS] = 1/2. Then by Theorem 1, εD(VS) ⪅ ln 2/n which
is much better than chance for large n. This corresponds to pure generalisation from the training
labels. Interpolating between these two extremes would suggest that a good neural architecture would
impose a prior on functions with strong intra-class and weak inter-class correlations.

4 TESTING THE BOUND

This section compares the generalisation bound for infinite width networks (Theorem 2) to the
performance of finite width multilayer perceptrons (MLPs) trained by gradient descent. The bound is
found to be non-vacuous and correlated with the effects of varying depth and dataset complexity. Still,
there is a substantial gap between the bound and gradient descent, which is investigated in Section 5.

Three modified versions of the MNIST handwritten digit dataset (LeCun et al., 1998) of varying
“hardness” were used in the experiments, as detailed in Figure 2 (top left). MLPs were trained with L
layers and W hidden units per layer. Specifically, each 28px × 28px input image x was flattened
to lie in R[784], and normalised to satisfy ∥x∥2 = _√784. The networks consisted of an input layer_

in R[784][×][W], (L − 2) layers in R[W][ ×][W], and an output layer in R[W][ ×][1]. The nonlinearity ϕ was set
to ϕ(z) := _√2 · max(0, z). For this architecture, the width W →∞_ kernel is the compositional

_arccosine kernel described in Theorem 4 in Appendix A.3. For the finite width networks, the training_
loss was set to square loss using the binary training labels as regression targets, and the networks
were trained for 100 epochs using the Nero optimiser (Liu et al., 2021) with an initial learning rate
of 0.01 decayed by a factor of 0.9 every epoch, and a mini-batch size of min(50, training set size)
data points. The final train error was 0% in all reported experiments. All bounds were computed
with a failure probability of δ = 0.01, and C0 was estimated using 10[6] Monte-Carlo samples. All
experiments were run on one NVIDIA Titan RTX GPU.

The generalisation bound of Theorem 2 was first compared across three datasets of varying complexity.
The network architecture was set to a depth L = 7 MLP, and the bound was computed via MonteCarlo estimation of 0. The results are shown in Figure 2 (top right). The bound was found to reflect
_C_
the relative hardness of the datasets. For random labels, the bound was vacuous as desired. Next,
the generalisation bound was compared against the empirical performance of a depth L = 7 MLP
trained on the decimal digits dataset. The results are shown in Figure 2 (bottom left). While loose
compared to the holdout error of the finite width network, the bound is still non-vacuous. Finally,
the effect of varying network depth was investigated on the decimal digits dataset. Two depths
were compared: L = 2 and L = 7. The results are shown in Figure 2 (bottom right). The bounds
(computed via Monte-Carlo estimation of 0) appear to predict the relative holdout performance of
_C_
the two architectures at finite width. These results corroborate those of Valle-Pérez & Louis (2020)
but without the use of the statistically inconsistent expectation-propagation approximation.


-----

**MNIST Variant** **Inputs** **Labels**

random labels _{0, 1, ..., 9}_ coin flip
decimal digits _{0, 1, ..., 9}_ parity
binary digits _{0, 1}_ parity


Comparing to gradient descent Varying the number of layers

0.2 W=5000 0.2 L=7, W= bound

Varying the dataset complexity

0.6

0.4

Error

Random labels

0.2 Chance

Decimal digits
Binary digits

0.0

100 200 300 400 500

Number of training examples

Varying the number of layers

0.6

0.4

Chance

Error L=2, W= bound

L=2, W=5000

0.2 L=7, W= bound

L=7, W=5000

0.0

100 200 300 400 500

Number of training examples

Comparing to gradient descent

0.6

0.4

Chance

Error W= bound

W= estimator

0.2 W=5000

W=10000

0.0

100 200 300 400 500

Number of training examples


Figure 2: Testing Theorem 2—the bound on the average population error of the NNGP posterior.
The bound is found to be non-vacuous and correlated with the effects of varying network depth and
dataset complexity. For all curves, the mean and range are plotted over three global random seeds. A
substantial gap is visible between the bound and the results of gradient descent training, which is
investigated in Section 5.

Top left: The three datasets used in the experiments, listed in order of hardness. For random labels
there is no meaningful relationship between image and label, so generalisation is impossible. Binary
_digits is easier than decimal digits because each class is less diverse._

Top right: Comparing generalisation bounds from Theorem 2 for datasets of varying hardness. The
curves are computed by Monte-Carlo estimation of C0 for an infinite width depth 7 MLP. The ordering
of the bounds reflects the dataset hardness. For random labels the bound is rightfully vacuous.

Bottom left: Comparing generalisation bounds from Theorem 2 to the results of training finite
width (W = 5000 and W = 10000) networks by gradient descent. The bounds are computed by
both Monte-Carlo estimation of C0 (referred to as W = ∞ estimator) and exact computation of C1
(referred to as W = ∞ bound). The comparison is made for depth 7 MLPs on the decimal digits
dataset. The exact bound computed via C1 is looser than the C0 bound computed via Monte-Carlo
estimation but only slightly, and both are non-vacuous above 150 datapoints. While the bounds follow
the same trend as the results of training finite width networks, a substantial gap is visible.

Bottom right: Comparing generalisation bounds from Theorem 2 to the results of training finite
width networks by gradient descent on decimal digits, for networks of varying depth. The bounds are
computed by exact computation of 1. The ordering of the bounds matches the finite width results,
_C_
but again a substantial gap is visible between the bounds and the results of gradient descent training.


-----

5 IMPLICIT BIAS OF GRADIENT DESCENT

In Section 3, a bound was derived on the average population error of all infinitely wide neural
networks that fit a certain training set. Since the bounded quantity measures the performance of all
solutions rather than just those returned by gradient descent, it is intended to assess the implicit bias
of architecture. In Section 4, this bound was tested and found to be non-vacuous. Still a substantial
gap was found between the bound and the performance of finite width networks trained by gradient
descent. This gap could arise for several potential reasons:

i) slackness in the bounding technique;

ii) a difference between infinite width and finite width neural networks;

iii) an additional implicit bias of gradient descent.

This section tests these various possibilities, ultimately concluding that gradient descent does have an
important additional implicit bias: the ability to control the margin of the returned network.

5.1 AN EXACT FORMULA FOR THE AVERAGE HOLDOUT ERROR OF SOLUTIONS

To investigate whether slackness in the bounding technique is responsible for the gap, an additional
theoretical tool was developed: a formula for the holdout error of a binary classifier averaged over the
solution space. In contrast to Theorem 2 which gives an upper bound on population error, Theorem 3
is an equality and therefore does not suffer from slackness. The proof is given in Appendix A.1.

**Theorem 3 (Average holdout error of a binary classifier). First, fix a probability measure P over the**
_weight space Ω_ _of a binary classifier. Let S denote a training set and let VS denote the solution space._
_For a holdout set T of m datapoints, consider the holdout error 0 ≤_ _εT (w) ≤_ 1 of weight setting
_w ∈_ Ω, and its average over the solution space εT (VS) := Ew∼P[εT (w) | w ∈ _VS]. Then:_

P[VS (x, _y)]_

_εT (VS) = [1]_ _∪_ _−_ _._ (11)

_m_ P[VS]

(x,yX)∈T

In words, the average holdout error over solutions equals the reduction in the volume of solutions
when the training set is augmented with a negated holdout point, averaged over the holdout set. For an
infinitely wide neural network, Equation 11 involves computing a sum of ratios of Gaussian orthant
probabilities. These probabilities can be consistently estimated by Equation 8 in Lemma 1.

5.2 THREE-WAY COMPARISON OF HOLDOUT ERROR

Armed with Theorem 3, this subsection makes a three-way comparison between the holdout error of
infinite width networks averaged over solutions, the holdout error of finite width networks averaged
over solutions, and the holdout error of finite width networks trained by gradient descent.

To compute the holdout error of finite width networks averaged over solutions, weight vectors were
randomly sampled and all non-solutions were discarded. To make this process computationally
tractable, a very small training set was used consisting of only 5 samples from binary digits. A
holdout set of 50 datapoints was used, and the network was set to a 7-layer MLP. The results were:

average holdout error at infinite width: 0.337 ± 0.001;
average holdout error at width 10,000: 0.33 ± 0.01;
gradient descent holdout error at width 10,000: 0.178 ± 0.007.

Based on these results, three comments are in order. First, the close agreement between the average
holdout error at finite and infinite width suggests that the infinite width limit may not be responsible
for the significant gap observed in Section 4. Second, the significant gap between the gradient descent
holdout error and the holdout error averaged over solutions suggests slackness in the bounding
technique in Theorem 2 may also not be the main culprit. Third, the significant gap between the
gradient descent holdout error and the holdout error averaged over solutions suggests that gradient
descent does have an extra implicit bias. The next subsection attempts to diagnose this implicit bias.


-----

Holdout error along the label ray

10 2 10 1 10[0] 10[1] 10[2]

|Chan Aver Netw Netw NNGP|ce age holdout error over orthant ork holdout error at scale ork train error at scale holdout error at scale|
|---|---|


2 1 [0] [1] [2]

Label scale



0.6

0.4


0.2

0.0

|Col1|f(x(2); ·)|
|---|---|
|label ray|f(x(1);|
|||




Varying the label scale at depth 2

0.6

0.4

Error ChanceW= estimator

W=5000, label scale 0.1

0.2 W=5000, label scale 1.0

0.0

100 200 300 400 500

Number of training examples


Varying the label scale at depth 7

0.6

0.4

Error

0.2

0.0

100 200 300 400 500

Number of training examples


Figure 3: Exploring generalisation performance along the label ray defined in Section 5.3. Functions
further along the ray are found to generalise substantially better than average, for both NNGP draws
and finite width networks trained by gradient descent.

Top left: Schematic diagram illustrating the label ray in function space. While all functions in the grey
shaded orthant attain zero train error, the darker shaded solutions have small margin and thus—under
a Gaussian process model—holdout predictions are expected to be driven by random fluctuations.

Top right: Holdout error along the label ray of the small learning task described in Section 5.2.
Results are shown for both networks trained by gradient descent and NNGP posterior samples. The
label scale α is defined in Section 5.3 and measures the distance of the function along the label
ray. NNGP posterior samples were generated by sampling from a Gaussian distribution with mean
and covariance given by Equations 12 and 13. Networks trained by gradient descent were found by
training with the loss function given in Equation 14. Shading shows the standard deviation over 100
random intialisations or posterior samples. Deep into function space (large α) the holdout error is
significantly better than for functions close to the origin (small α). The average holdout error over
the orthant, as computed in Section 5.2, is also plotted. Large margin networks and NNGP posterior
samples both significantly outperform the orthant average. Finally, a gap is visible between large
margin NNGP posterior samples and large margin networks trained by gradient descent—this may be
due to an additional and undiagnosed implicit bias of gradient descent.

Bottom: Repeating experiments from Section 4 with a smaller label scale to sanity check the findings.
For a depth 2 and a depth 7 MLP on the decimal digits dataset, the experiment from Section 4 was
repeated with the loss function set to Lα=0.1 as defined in Equation 14. The infinite width W = ∞
curve shows the Theorem 2 bound estimated via 0. The other two curves show the results of
_C_
networks trained by gradient descent with Lα=0.1 and Lα=1. Despite all networks attaining 0% train
error, holdout error was significantly worse for networks trained using _α=0.1. Also, the networks_
_L_
trained using Lα=0.1 exhibit a substantially smaller gap with the Theorem 2 upper bound.


-----

5.3 DIAGNOSIS: MARGIN


This subsection finds that gradient descent has an important implicit bias in determining the margin
of the returned network. This conclusion is made by studying the generalisation error of functions
along the label ray in function space—depicted in Figure 3 (top left). For a training set S with a
vector of binary labels y 1, a function α-far along the label ray refers to the point αy _S._
_∈{±_ _}[n]_ _∈F_

For the case of a zero mean NNGP, the predictive distribution on a holdout set T conditioned on
training labels α-far along the label ray is given by Bishop (2006, Chapter 2.3.1):

posterior mean = α × ΣT SΣ[−]SS[1][y,] (12)

posterior covariance = ΣT T − ΣT SΣ[−]SS[1][Σ][ST][,] (13)

where ΣT S is the covariance between holdout and train inputs, and ΣST, ΣSS and ΣT T are defined
analogously. For large α, predictions are driven by the posterior mean, whereas for small α, they are
driven by random fluctuations. So one expects that letting α →∞ should improve holdout error.

For the case of gradient descent training, functions α-far along the label ray can be returned by
minimising the following loss function:

_n_

_α(W_ ) := [1] _f_ (x[(][i][)]; W ) _αy[(][i][)][][2]_ _._ (14)
_L_ _n_ _−_

_i=1_

X 

A subtle but important point is that for gradient descent training with Nero (Liu et al., 2021) the
norms of the weight matrices are constrained, thus α controls a properly normalised notion of margin.

An experiment was performed to measure the holdout error of networks α-far along the label ray—for
both NNGP draws and neural networks selected by Nero—for α ranging from 10[−][2] to 10[2]. As can
be seen in Figure 3 (top right), varying α appears to directly control the holdout error, despite all
solutions attaining 0% train error. Moreover, large α solutions significantly outperform the average
holdout error over the solution space. This suggests that gradient descent possesses a significant
implicit bias in its ability to control margin—going beyond the implicit bias of architecture.

To sanity check this result, experiments from Section 4 were repeated with the loss function _α=0.1_
_L_
replacing _α=1. As can be seen in Figure 3 (bottom), the holdout performance was significantly_
_L_
diminshed at α = 0.1, as was the gap with the PAC-Bayes bound from Theorem 2.

6 DISCUSSION AND CONCLUSION

This paper has explored the separate implicit biases of architecture and gradient descent. Section 3
derived an analytical generalisation bound on the NNGP posterior, Section 4 found this bound to be
non-vacuous, while Section 5 showed that large margin functions substantially outperform the bound.

The findings in this paper support the importance of the implicit bias of architecture: reproducing the
findings in Valle-Pérez & Louis (2020) but with improved technical tools, the average generalisation
performance of architecture was already found to be good. But the implicit bias of gradient descent
was also found to be important. In particular, in contrast to an assumption made by Valle-Pérez &
Louis (2020) that gradient descent “samples the zero-error region close to uniformly”, it was found
that gradient descent can be used to target zero-error functions with large margin. This also subtly
counters a proposal by Mingard et al. (2021) that gradient descent acts like a “Bayesian sampler”.
In particular, gradient descent can target functions with margin α →∞ for which the Bayesian
posterior probability → 0. And indeed, these minimum a posteriori functions seem to generalise best.

One direction of future work suggested by these results is an improvement to the function space
PAC-Bayes theory to account for margin. While margin-based PAC-Bayes bounds do already exist
(Herbrich, 2001; Langford & Shawe-Taylor, 2003; Neyshabur et al., 2018), these bounds operate in
weight space and further do not seem to imply an advantage to selecting the max-margin classifier
over the other classifiers included in the PAC-Bayes posterior.

Ultimately, a generalisation theory that properly accounts for the various implicit biases of deep
learning could provide a more principled basis for both neural architecture design as well as the design
of new regularisation schemes. It is hoped that the new analytical results as well as experimental
insights included in this paper contribute a step towards reaching that goal.


_α(W_ ) := [1]
_L_ _n_


_i=1_


-----

REFERENCES

Alessandro Achille and Stefano Soatto. Emergence of invariance and disentanglement in deep
representations. Journal of Machine Learning Research, 2018.

Navid Azizan, Sahin Lale, and Babak Hassibi. Stochastic mirror descent on overparameterized
nonlinear models. IEEE Transactions on Neural Networks and Learning Systems, 2021.

Peter L. Bartlett, Dylan J. Foster, and Matus Telgarsky. Spectrally-normalized margin bounds for
neural networks. In Neural Information Processing Systems, 2017.

Felix Biggs and Benjamin Guedj. On margins and derandomisation in PAC-Bayes. arXiv:2107.03955,
2021.

Christopher M. Bishop. Pattern Recognition and Machine Learning. Springer-Verlag, 2006.

Youngmin Cho and Lawrence Saul. Kernel methods for deep learning. In Neural Information
_Processing Systems, 2009._

Gintare Karolina Dziugaite and Daniel M. Roy. Computing nonvacuous generalization bounds for
deep (stochastic) neural networks with many more parameters than training data. In Uncertainty in
_Artificial Intelligence, 2017._

Pierre Foret, Ariel Kleiner, Hossein Mobahi, and Behnam Neyshabur. Sharpness-aware minimization
for efficiently improving generalization. In International Conference on Learning Representations,
2021.

Jonas Geiping, Micah Goldblum, Phillip E. Pope, Michael Moeller, and Tom Goldstein. Stochastic
training is not necessary for generalization. arXiv:2109.14119, 2021.

Ralf Herbrich. Learning Kernel Classifiers: Theory and Algorithms. MIT Press, 2001.

Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Peter
Tang. On large-batch training for deep learning: Generalization gap and sharp minima. In
_International Conference on Learning Representations, 2017._

John Langford and Matthias Seeger. Bounds for averaging classifiers. Technical report, Carnegie
Mellon University, 2001.

John Langford and John Shawe-Taylor. PAC-Bayes & margins. In Neural Information Processing
_Systems, 2003._

Yann LeCun, Corinna Cortes, and Christopher J.C. Burges. MNIST handwritten digit database, 1998.

Jaehoon Lee, Jascha Sohl-Dickstein, Jeffrey Pennington, Roman Novak, Sam Schoenholz, and
Yasaman Bahri. Deep neural networks as Gaussian processes. In International Conference on
_Learning Representations, 2018._

Yang Liu, Jeremy Bernstein, Markus Meister, and Yisong Yue. Learning by turning: Neural
architecture aware optimisation. In International Conference on Machine Learning, 2021.

David A. McAllester. Some PAC-Bayesian theorems. In Conference on Computational Learning
_Theory, 1998._

Chris Mingard, Guillermo Valle-Pérez, Joar Skalse, and Ard A. Louis. Is SGD a Bayesian sampler?
Well, almost. Journal of Machine Learning Research, 2021.

Radford M. Neal. Bayesian Learning for Neural Networks. Ph.D. thesis, Department of Computer
Science, University of Toronto, 1994.

Behnam Neyshabur, Srinadh Bhojanapalli, and Nathan Srebro. A PAC-Bayesian approach to
spectrally-normalized margin bounds for neural networks. In International Conference on Learning
_Representations, 2018._


-----

Behnam Neyshabur, Zhiyuan Li, Srinadh Bhojanapalli, Yann LeCun, and Nathan Srebro. The role
of over-parametrization in generalization of neural networks. In International Conference on
_Learning Representations, 2019._

Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro. The implicit
bias of gradient descent on separable data. Journal of Machine Learning Research, 2018.

Guillermo Valle-Pérez and Ard A. Louis. Generalization bounds for deep learning. arXiv:2012.04115,
2020.

Guillermo Valle-Pérez, Chico Q. Camargo, and Ard A. Louis. Deep learning generalizes because
the parameter–function map is biased towards simple functions. In International Conference on
_Learning Representations, 2019._

Aad W. van der Vaart. Asymptotic Statistics. Cambridge University Press, 1998.

Vladimir N. Vapnik and Alexey Ya. Chervonenkis. On the uniform convergence of relative frequencies
of events to their probabilities. Theory of Probability & Its Applications, 1971.

Colin Wei, J. Lee, Qiang Liu, and Tengyu Ma. On the margin theory of feedforward neural networks.
_arXiv:1810.05369, 2018._

Ashia C. Wilson, Rebecca Roelofs, Mitchell Stern, Nathan Srebro, and Benjamin Recht. The marginal
value of adaptive gradient methods in machine learning. In Neural Information Processing Systems,
2017.

Yikai Wu, Xingyu Zhu, Chenwei Wu, Annie Wang, and Rong Ge. Dissecting Hessian: Understanding
common structure of Hessian in neural networks. arXiv:2010.04261, 2020.

Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals. Understanding deep learning requires rethinking generalization. In International Conference on Learning
_Representations, 2017._


-----

APPENDIX A PROOFS

A.1 FUNCTION SPACE ORTHANT ARITHMETIC

The following PAC-Bayesian theorem is due to Valle-Pérez et al. (2019):
**Theorem 1 (A PAC-Bayesian theorem). First, fix a probability measure P over the weight space Ω**
_of a classifier. Let S denote a training set of n datapoints sampled iid from the data distribution D_
_and let VS denote the corresponding solution space. Consider the population error 0 ≤_ _εD(w) ≤_ 1
_of weight setting w ∈_ Ω, and its average over the solution space εD(VS) := Ew∼P[εD(w) | w ∈ _VS]._
_Then, for a proportion 1 −_ _δ of draws of the training set S,_

1 ln P[V1S ] [+ ln][ 2]δ[n]
_εD(VS)_ ln _._ (4)
_≤_ 1 _εD(VS)_ _n_ 1

_−_ _[≤]_ _−_

_Proof. The first inequality is a basic property of logarithms. The second inequality follows from_
Theorem 3 of Langford & Seeger (2001), by setting the prior measure to P (·) = P(·) and the posterior
measure to the conditional Q( ) = P( _VS). Under these settings, the average training error rate over_

_·_ _·|_ 1
the posterior is zero and KL(Q|P ) = ln P[VS ] [.]

**Theorem 3 (Average holdout error of a binary classifier). First, fix a probability measure P over the**
_weight space Ω_ _of a binary classifier. Let S denote a training set and let VS denote the solution space._
_For a holdout set T of m datapoints, consider the holdout error 0 ≤_ _εT (w) ≤_ 1 of weight setting
_w ∈_ Ω, and its average over the solution space εT (VS) := Ew∼P[εT (w) | w ∈ _VS]. Then:_

P[VS (x, _y)]_

_εT (VS) = [1]_ _∪_ _−_ _._ (11)

_m_ P[VS]

(x,yX)∈T

_Proof. The result follows by interchanging the order of expectations, rewriting the expectation of an_
indicator variable as a probability and finally applying the definition of conditional probability:


_εT (VS) := Ew∼P[εT (w) | w ∈_ _VS] = m[1]_

= [1]

_m_

= [1]

_m_

= [1]

_m_

A.2 GAUSSIAN ORTHANT PROBABILITIES


Ew∼P [I[sign f (x, w) = −y] | w ∈ _VS]_
(x,yX)∈T

P [sign f (x, w) = −y | w ∈ _VS]_
(x,yX)∈T

P [sign f (x, w) = −y and w ∈ _VS]_

P [w _VS]_

(x,yX)∈T _∈_

P[VS (x, _y)]_
_∪_ _−_ _._

P[VS]

(x,yX)∈T


**Lemma 1 (Gaussian orthant probability). For a covariance matrix Σ ∈** R[n][×][n], and a binary vector
_y ∈{±1}[n], let p denote the corresponding Gaussian orthant probability:_

_p := Pϕ_ (0,Σ)[sign(ϕ) = y]. (7)
_∼N_

_Letting I denote the n_ _×_ _n identity matrix, ⊙_ _the elementwise product and |·| the elementwise absolute_
_value, then p may be equivalently expressed as:_


det ΣΣ[−][1]−I)(y⊙|u|)[i] =: e[−C][0][(Σ][,y][)], (8)


_p = [1]_ e[−] 2[1] [(][y][⊙|][u][|][)][T] [(]

2[n][ E][u][∼N][ (0][,][I][)]
h

_and p may be bounded as follows:_


_p_ _n2_ _[−]_ _√[n]_
_≥_ 2[1][n][ e]


det Σ[( [1]2 _[−]_ _π[1]_ [) tr(Σ][−][1][)+][ 1]π _[y][T][ Σ][−][1][y][] ]=: [1]_ _n2_ _[−C][1][(Σ][,y][)]._ (9)

2[n][ e]


-----

_Proof. The orthant probability may first be expressed using the probability density function of the_
multivariate Normal distribution as follows:


e[−] 2[1] _[ϕ][T][ Σ][−][1][ϕ]_ dϕ.
_y⊙ϕ≥0_


_p =_


(2π)[n] det Σ


_y_ _ϕ_ 2n
By the change of variables u = 2√n _⊙det Σ_ [or equivalently][ ϕ][ =] _√det Σ(y ⊙_ _u), the orthant probability_

may be expressed as:


e[−] [1]2 [(][y][⊙][u][)][T n]√
_u≥0_


det ΣΣ[−][1](y⊙u) du


_p =_


(2π)[n]


= [1]

2[n]


= [1] 1 2 [(][y][⊙|][u][|][)][T n]√det ΣΣ[−][1](y⊙|u|) du

2[n] (2π)[n] R[n][ e][−] [1]

Z

_n_

= [1] p e[−] 2[1] [(][y][⊙|][u][|][)][T] [(] _√det ΣΣ[−][1]−I)(y⊙|u|)[i]_ _,_

2[n][ E][u][∼N][ (0][,][I][)]
h

where the second equality follows by symmetry, and the third equality follows by inserting a factor of
e[−][u][2][/][2]e[+][u][2][/][2] = 1 into the integrand.


(2π)[n]


Next, by Jensen’s inequality,

_n_

_p_ 2 [E][u][∼N][ (0][,][I][)][[][(][y][⊙|][u][|][)][T] [(] _√det ΣΣ[−][1]−I)(y⊙|u|)]_
_≥_ 2[1][n][ e][−] [1]

_n_

= [1] 2 _ij_ [E][u][∼N][ (0][,][I][)][[][y][i][y][j] _[|][u][i][||][u][j]_ _[|][(]_ _√det ΣΣ[−]ij[1][−][δ][ij][)]]_

2[n][ e][−] [1] P

_n_ _n_

= [1] 2 [[] _i[(]_ _√det ΣΣ[−]ii[1][−][1][)][+][ 2]π_ _i≠_ _j_ _[y][i][y][j]_ _√det ΣΣ[−]ij[1][]]_

2[n][ e][−] [1] P P

_n_

= [1] 2 [[] _i[(][(1][−]_ _π[2]_ [)][ n]√det ΣΣ[−]ii[1][−][1][)][+][ 2]π _ij_ _[y][i][y][j]_ _√det ΣΣ[−]ij[1][]]_

2[n][ e][−] [1] P P

_n_

= [1] 2 [[][(1][−] _π[2]_ [)][ n]√det Σ tr(Σ[−][1])−n+ _π[2]_ _√det Σy[T]_ Σ[−][1]y]

2[n][ e][−] [1]

= [1] _n2_ _[−]_ _√[n]_ det Σ[( [1]2 _[−]_ _π[1]_ [) tr(Σ][−][1][)+][ 1]π _[y][T][ Σ][−][1][y][]]._

2[n][ e]


The third equality follows by noting Eu∼N (0,I) |ui||ui| = Eu∼N (0,1) u[2] = 1, while for i ̸= j,
Eu∼N (0,I) |ui||uj| = [Eu∼N (0,1) |u|][2] = 2/π.

A.3 NEURAL NETWORKS AS GAUSSIAN PROCESSES

The essence of the following lemma is due to Neal (1994). The lemma will be used in the proof of
Theorem 4.
**Lemma 2 (NNGP correspondence). For the neural network layer given by Equation 18, consider**
_randomly sampling the weight matrix W_ [(][l][)]. If the following hold:

_(i) for every x ∈_ R[d][0] _, the activations ϕ_ _z1[(][l][−][1)](x)_ _, ..., ϕ_ _zd[(][l]l[−]−1[1)][(][x][)]_ _are iid with finite first_
_and second moment;_    

_(ii) the weights Wij[(][l][)]_ _[are drawn iid with zero mean and finite variance;]_

_(iii) for any random variable z with finite first and second moment, ϕ(z) also has finite first and_
_second moment;_

_then, in the limit that dl−1 →∞, the following also hold:_

_(1) for every x ∈_ R[d][0], the activations ϕ _z1[(][l][)][(][x][)]_ _, ..., ϕ_ _zd[(][l]l[)][(][x][)]_ _are iid with finite first and_
_second moment;_    

_(2) for any collection of k inputs x[(1)], ..., x[(][k][)], the distribution of the ith pre-activations_
_zi[(][l][)][(][x][(1)][)][, ..., z]i[(][l][)][(][x][(][k][)][)][ is jointly Normal.]_


-----

While condition (i) may seem non-trivial, notice that the lemma propagates this condition to the next
layer via entailment (1). This means that provided condition (i) holds for ϕ(z[(1)](x)) at the first layer,
then recursive application of the lemma implies that the network’s pre-activations are jointly Normal
_at all layers via entailment (2)._

_Proof of Lemma 2. To establish entailment (1), consider the dl-dimensional vector Z1_ :=
_z1[(][l][)][(][x][)][, ..., z]d[(][l]l[)][(][x][)]_ . Observe that Z1 satisfies:
h i


_dl−1_

_j=1_

X


_W1[(]j[l][)][ϕ]_ _zj[(][l][−][1)](x)_ _, ..., Wd[(]l[l]j[)]_ _[ϕ]_ _zj[(][l][−][1)](x)_ _._ (15)
   i


_Z1 =_


_dl−1_


By conditions (i) and (ii), the summands in Equation 15 are iid random vectors with zero mean,
and any two distinct components of the same vector summand have the same variance and zero
covariance. Then by the multivariate central limit theorem (van der Vaart, 1998, p. 16), in the limit
that dl 1, the components of Z1 are Gaussian with a covariance equal to a scaled identity
matrix. In particular, the components of− _→∞_ _Z1 are iid with finite first and second moment. Applying_
condition (iii) then implies that the same holds for ϕ(Z1). This establishes entailment (1).

To establish entailment (2), consider instead the _k-dimensional_ vector _Z2_ :=
_zi[(][l][)][(][x][(1)][)][, ..., z]i[(][l][)][(][x][(][k][)][)]_ . Observe that Z2 satisfies:
h i


_dl−1_

_j=1_

X


_Wij[(][l][)][ϕ]_ _zj[(][l][−][1)](x[(1)])_ _, ..., Wij[(][l][)][ϕ]_ _zj[(][l][−][1)](x[(][k][)])_ _._ (16)
   i


_Z2 =_


_dl−1_


Again by combining conditions (i) and (ii), the summands in Equation 16 are iid random vectors with
finite mean and finite covariance. Then as dl 1, the distribution of Z2 is jointly Normal—again
by the multivariate central limit theorem. This establishes entailment (2).− _→∞_

The essence of the following theorem appears in a paper by Lee et al. (2018), building on the work of
Cho & Saul (2009). The theorem and its proof are included for completeness.
**Theorem 4 (NNGP for relu networks). Consider an L-layer MLP defined recursively via:**

1 _d0_
_zi[(1)](x) =_ _√d0_ _j=1_ _Wij[(1)][x][j][,]_ (17)

X


_dl−1_

_Wij[(][l][)][ϕ]_ _zj[(][l][−][1)](x)_ _,_ (18)
_j=1_

X  


_zi[(][l][)][(][x][) =]_


_dl−1_


_where x ∈_ R[d][0] _denotes an input, z[(][l][)](x) ∈_ R[d][l] _denotes the pre-activations at the lth layer, W_ [(][l][)]
_denotes the weight matrix at the lth layer, and ϕ denotes the nonlinearity._

_Set the output dimension dL = 1 and set the nonlinearity to a scaled relu ϕ(z) :=_ _√2 · max(0, z)._

_Suppose that the weight matrices W_ [(1)], ..., W [(][L][)] _have entries drawn iid N_ (0, 1), and consider any
_collection of k inputs x[(1)], ..., x[(][k][)]_ _each with Euclidean norm_ _d0._

_[√]_

_If d1, ..., dL_ 1 _, the distribution of outputs z[(][L][)](x[(1)]), ..., z[(][L][)](x[(][k][)])_ R induced by random
_sampling of the weights is jointly Normal with mean zero and covariance:−_ _→∞_ _∈_


_xT x′_

E _z[(][L][)](x)z[(][L][)](x[′])_ = h _..._ _h_ ; (19)
_◦_ _◦_ _d0_
h i _L−1 times_  

_where h(t) :=_ _π[1]_ _√1_ _t[2]_ + t (π arccos t) _._ | {z }

_−_ _·_ _−_

 

_Proof. Condition (ii) of Lemma 2 holds at all layers for iid standard Normal weights, and condition_
(iii) holds trivially for the scaled relu nonlinearity. Provided one can establish condition (i) for
the first layer activations ϕ _z1[(1)][(][x][)]_ _, ..., ϕ_ _zd[(1)]1_ [(][x][)], then condition (i) will hold at all layers by
   


-----

recursive application of Lemma 2, thus establishing joint Normality of the pre-activations at all layers
(including the network outputs). But condition (i) holds at the first layer, since it is quick to check
by Equation 17 that for any x satisfying ∥x∥2 = _[√]d0, the pre-activations z1[(1)][(][x][)][, ..., z]d[(1)]1_ [(][x][)][ are iid]
_N (0, 1), and ϕ preserves both iid-ness and finite-ness of the first and second moment._

Since the pre-activations at any layer are jointly Normal, all that remains is to compute their first and
second moments. For the ith hidden unit in the lth layer, the first moment E[zi[(][l][)][(][x][)] = 0][. This can]
be seen by taking the expectation of Equation 18 and using the fact that the Wij[(][l][)] [are independent of]
the previous layer’s activations and have mean zero.

Since the pre-activations zi[(][l][)][(][x][)][ and][ z]i[(][l][)][(][x][′][)][ are jointly Normal with mean zero, their distribution is]
completely described by their covariance matrix Σl(x, x[′]), defined by:

_ρl(x, x[′]) := E_ _zi[(][l][)][(][x][)][z]i[(][l][)][(][x][′][)]_
h i

Σl(x, x[′]) := _ρρll((x, xx, x[′]))_ _ρρll((xx, x[′], x[′][′]))_ _,_
 

where the hidden unit index i is unimportant since hidden units in the same layer are identically
distributed.

The theorem statement will follow from an effort to express Σl(x, x[′]) in terms of Σl−1(x, x[′]),
and then recursing back through the network. By Equation 18 and independence of the Wij[(][l][)][, the]
covariance ρl(x, x[′]) may be expressed as:

_ρl(x, x[′]) = E_ _ϕ_ _zj[(][l][−][1)](x)_ _ϕ_ _zj[(][l][−][1)](x[′])_ _,_ (20)

where j indexes an arbitrary hidden unit in theh  (l − 1)th layer. To make progress, it helps to first i
evaluate:

2[]
_ρl(x, x) = E_ _ϕ_ _zj[(][l][−][1)](x)_ = [1]2
   _[·][ 2][ ·][ ρ][l][−][1][(][x, x][)][,]_

which follows by the definition of ϕ and symmetry of the Gaussian expectation around zero. Then,
by recursion:
_ρl(x, x) = ρl−1(x, x) = ... = ρ1(x, x) = 1,_
where the final equality holds because the first layer pre-activations are iid N (0, 1) by Equation 17.
Therefore, the covariance Σl 1 at layer l 1 is just:
_−_ _−_

Σl−1(x, x[′]) = ρl−1(1x, x[′]) _ρl−1(1x, x[′])_ _,_

Equation 20 may now be used to express ρl(x, x[′]) in terms of ρl−1(x, x[′]). Dropping the (x, x[′])
indexing for brevity:
_ρl = Eu,v_ (0,Σl 1) [ϕ (u) ϕ (v)]
_∼N_ _−_

1
= du dv exp _uv._

_π_ 1 − _ρ[2]l−1_ ZZu,v≥0 − _[u][2][ −]2(1[2][ρ] −[l][−]ρ[1][uv][2]l−1[ +][)]_ _[ v][2]_ 

By making the substitutionq _ρl−1 = cos θ, this integral becomes equivalent to_ _π[1]_ _[J][1][(][θ][)][ as expressed in]_

Equation 15 of Cho & Saul (2009). Substituting in the evaluation of this integral (Cho & Saul, 2009,
Equation 6), one obtains:
_ρl(x, x[′]) = h(ρl−1(x, x[′])),_ (21)

where h(t) := _π[1]_ _√1_ _t[2]_ + t (π arccos t) .

_−_ _·_ _−_

All that remains is to evaluate _ρ1(x, x[′]). Since E_ _Wij[(1)][W][ (1)]ik_ = δjk, this is given by:
h i

_ρ1(x, x[′]) := E_ _zi[(1)](x)zi[(1)](x[′])_
h i


_d0_

E _Wij[(1)][W][ (1)]ik_ _xjx[′]k_ [=][ x][T][ x][′]

_d0_

_j,k=1_

X h i


= [1]

_d0_


The proof is completed by combining this expression with the recurrence relation in Equation 21.


-----

