A theory of data variability in Neural Network Bayesian inference
Bayesian inference and kernel methods are well established in machine
learning. The neural network Gaussian process in particular provides a concept
to investigate neural networks in the limit of infinitely wide hidden layers by
using kernel and inference methods. Here we build upon this limit and provide a
field-theoretic formalism which covers the generalization properties of
infinitely wide networks. We systematically compute generalization properties
of linear, non-linear, and deep non-linear networks for kernel matrices with
heterogeneous entries. In contrast to currently employed spectral methods we
derive the generalization properties from the statistical properties of the
input, elucidating the interplay of input dimensionality, size of the training
data set, and variability of the data. We show that data variability leads to a
non-Gaussian action reminiscent of a ($\varphi^3+\varphi^4$)-theory. Using our
formalism on a synthetic task and on MNIST we obtain a homogeneous kernel
matrix approximation for the learning curve as well as corrections due to data
variability which allow the estimation of the generalization properties and
exact results for the bounds of the learning curves in the case of infinitely
many training data points.
Javed Lindner, David Dahmen, Michael Krämer, Moritz Helias
# A theory of data variability in Neural Network Bayesian inference

###### Abstract

Bayesian inference and kernel methods are well established in machine learning. The neural network Gaussian process in particular provides a concept to investigate neural networks in the limit of infinitely wide hidden layers by using kernel and inference methods. Here we build upon this limit and provide a field-theoretic formalism which covers the generalization properties of infinitely wide networks. We systematically compute generalization properties of linear, non-linear, and deep non-linear networks for kernel matrices with heterogeneous entries. In contrast to currently employed spectral methods we derive the generalization properties from the statistical properties of the input, elucidating the interplay of input dimensionality, size of the training data set, and variability of the data. We show that data variability leads to a non-Gaussian action reminiscent of a \(\varphi^{3}+\varphi^{4}\)-theory. Using our formalism on a synthetic task and on MNIST we obtain a homogeneous kernel matrix approximation for the learning curve as well as corrections due to data variability which allow the estimation of the generalization properties and exact results for the bounds of the learning curves in the case of infinitely many training data points.

## I Introduction

Machine learning and in particular deep learning continues to influence all areas of science. Employed as a scientific method, explainability, a defining feature of any scientific method, however, is still largely missing. This is also important to provide guarantees and to guide educated design choices to reach a desired level of accuracy. The reason is that the underlying principles by which artificial neural networks reach their unprecedented performance are largely unknown. There is, up to date, no complete theoretical framework which fully describes the behavior of artificial neural networks so that it would explain the mechanisms by which neural networks operate. Such a framework would also be useful to support architecture search and network training.

Investigating the theoretical foundations of artificial neural networks on the basis of statistical physics dates back to the 1980s. Early approaches to investigate neural information processing were mainly rooted in the spin-glass literature and included the computation of the memory capacity of the perceptron, path integral formulations of the network dynamics [36], and investigations of the energy landscape of attractor network [1; 10; 11].

As in the thermodynamic limit in solid state physics, modern approaches deal with artificial neural networks (ANN) with an infinite number of hidden neurons to simplify calculations. This leads to a relation between ANNs and Bayesian inference on Gaussian processes [29; 39], known as the Neural Network Gaussian Process (NNGP) limit: The prior distribution of network outputs across realizations of network parameters here becomes a Gaussian process that is uniquely described by its covariance function or kernel. This approach has been used to obtain insights into the relation of network architecture and trainability [30; 34]. Other works have investigated training by gradient descent as a means to shape the corresponding kernel [15]. A series of recent studies also captures networks at finite width, including adaptation of the kernel due to feature learning effects [4; 20; 44; 27; 33; 45]. Even though training networks with gradient descent is the most abundant setup, different schemes such as Bayesian Deep Learning [26] provide an alternative perspective on training neural networks. Rather than finding the single-best parameter realization to solve a given task, the Bayesian approach aims to find the optimal parameter distribution.

In this work we adopt the Bayesian approach and investigate the effect of variability in the training data on the generalization properties of wide neural networks. We do so in the limit of infinitely wide linear and non-linear networks. To obtain analytical insights, we apply tools from statistical field theory to derive approximate expressions for the predictive distribution in the NNGP limit.The remainder of this work is structured in the following way: In Section II we describe the setup of supervised learning in shallow and deep networks in the framework of Bayesian inference and we introduce a synthetic data set that allows us to control the degree of pattern separability, dimensionality, and variability of the resulting overlap matrix. In Section III we develop the field theoretical approach to learning curves and its application to the synthetic data set as well as to MNIST [17]: Section III.1 presents the general formalism and shows that data variability in general leads to a non-Gaussian process. Here we also derive perturbative expressions to characterize the posterior distribution of the network output. We first illustrate these ideas on the simplest but non-trivial example of linear Bayesian regression and then generalize them first to linear and then to non-linear deep networks. We show results for the synthetic data set to obtain interpretable expressions that allow us to identify how data
variability affects generalization; we then illustrate the identified mechanisms on MNIST. In Section IV we summarize our findings, discuss them in the light of the literature, and provide an outlook.

## II Setup

In this background section we outline the relation between neural networks, Gaussian processes, and Bayesian inference. We further present an artificial binary classification task which allows us to control the degree of pattern separation and variability and test the predictive power of the theoretical results for the network generalization properties.

### Neural networks, Gaussian processes and Bayesian inference

The advent of parametric methods such as neural networks is preceded by non-parametric approaches such as Gaussian processes. There are, however, clear connections between the two concepts which allow us to borrow from the theory of Gaussian processes and Bayesian inference to describe the seemingly different neural networks. We will here give a short recap on neural networks, Bayesian inference, Gaussian processes, and their mutual relations.

Figure 1: **Field theory of generalization in Bayesian inference.****a)** A binary classification task, such as distinguishing pairs of digits in MNIST, can be described with help of an overlap matrix \(K^{x}\) that represents similarity across the \(c=c_{1}+c_{2}\) images of the training set of two classes, \(1\) and \(2\) with \(D_{1}\) and \(D_{2}\) samples respectively. Entries of the overlap matrix are heterogeneous. Different drawings of \(c\) example patterns each lead to different realizations of the overlap matrix; the matrix is stochastic. We here describe the matrix elements by a correlated multivariate Gaussian. **b)** The data is fed through a feed-forward neural network to produce an output \(y\). In the case of infinitely wide hidden layers and under Gaussian priors on the network weights, the output of the network is a Gaussian process with the kernel \(K^{y}\), which depends on the network architecture and the input kernel \(K^{x}\). **c)** To obtain statistical properties of the posterior distribution, we compute its disorder-averaged moment generating function \(\overline{Z}(J,l_{*})\) diagrammatically. **d)** The leading-order contribution from the homogeneous kernel \(\langle y^{*}\rangle_{0}\) is corrected by \(\langle y^{*}\rangle_{1}\) due to the variability of the overlaps; both follow as derivatives of \(\overline{Z}(J,l_{*})\). **e)** Comparing the mean network output on a test point \(\langle y^{*}\rangle\), the zeroth order theory \(\langle y_{*}\rangle_{0}\) (blue dashed), the first-order approximation in the data-variability \(\langle y_{*}\rangle_{0+1}\) (blue-red dashed) and empirical results (black crosses) as a function of the amount of training data (learning curve) shows how variability in the data set limits the network performance and validates the theory.
Background: Neural Networks

In general a feed forward neural network maps inputs \(x_{\alpha}\in\mathbb{R}^{N_{\mathrm{dim}}}\) to outputs \(y_{\alpha}\in\mathbb{R}^{N_{\mathrm{out}}}\) via the transformations

\[h_{\alpha}^{(l)} =\mathbf{W}^{(l)}\phi^{(l)}\left(h_{\alpha}^{(l-1)}\right)\quad \text{with}\quad h_{\alpha}^{0}=\mathbf{V}x_{\alpha}\,,\] \[y_{\alpha} =\mathbf{U}\phi^{(L+1)}\left(h_{\alpha}^{(L)}\right)\,, \tag{1}\]

where \(\phi^{(l)}(x)\) are activation functions, \(\mathbf{V}\in\mathbb{R}^{N_{h}\times N_{\mathrm{dim}}}\) are the read-in weights, \(N_{\mathrm{dim}}\) is the dimension of the input, \(\mathbf{W}^{(l)}\in\mathbb{R}^{N_{h}\times N_{h}}\) are the hidden weights, \(N_{h}\) denotes the number of hidden neurons, and \(\mathbf{U}\in\mathbb{R}^{N_{\mathrm{out}}\times N_{h}}\) are the read-out weights. Here \(l\) is the layer index \(1\leq l\leq L\) and \(L\) the number of layers of the network; we here assume layer- independent activation functions \(\phi^{(l)}=\phi\). The collection of all weights are the model parameters \(\Theta=\left\{\mathbf{V},\mathbf{W}^{(1)},\ldots,\mathbf{W}^{(L)},\mathbf{U}\right\}\). The goal of training a neural network in a supervised manner is to find a set of parameters \(\hat{\Theta}\) which reproduces the input-output relation \(\left(x_{\mathrm{tr},\alpha},y_{\mathrm{tr},\alpha}\right)_{1\leq\alpha\leq D}\) for a set of \(D\) pairs of inputs and outputs as accurately as possible, while also maintaining the ability to generalize. Hence one partitions the data into a training set \(\mathcal{D}_{\mathrm{tr}}\), \(\left|\mathcal{D}_{\mathrm{tr}}\right|=D\), and a test-set \(\mathcal{D}_{\mathrm{test}}\), \(\left|\mathcal{D}_{\mathrm{test}}\right|=D_{\mathrm{test}}\). The training data is given in the form of the matrix \(\mathbf{x}_{\mathrm{tr}}\in\mathbb{R}^{N_{\mathrm{dim}}\times N_{\mathrm{tr}}}\) and \(\mathbf{y}_{\mathrm{tr}}\in\mathbb{R}^{N_{\mathrm{out}}\times N_{\mathrm{tr}}}\). The quality of how well a neural network is able to model the relation between inputs and outputs is quantified by a task-dependent loss function \(\mathcal{L}\left(\Theta,x_{\alpha},y_{\alpha}\right)\). Starting with a random initialization of the parameters \(\Theta\), one tries to find an optimal set of parameters \(\hat{\Theta}\) that minimizes the loss \(\sum_{\alpha=1}^{D}\mathcal{L}\left(\Theta,x_{\mathrm{tr},\alpha},y_{\mathrm{tr },\alpha}\right)\) on the training set \(\mathcal{D}_{\mathrm{tr}}\). The parameters \(\hat{\Theta}\) are usually obtained through methods such as stochastic gradient descent. The generalization properties of the network are quantified after the training by computing the loss \(\mathcal{L}\left(\hat{\Theta},x_{\alpha},y_{\alpha}\right)\) on the test set \(\left(x_{\mathrm{test},\alpha},y_{\mathrm{test},\alpha}\right)\in\mathcal{D}_{ \mathrm{test}}\), which are data samples that have not been used during the training process. Neural networks hence provide, by definition, a parametric modeling approach, as the goal is to a find an optimal set of parameters \(\hat{\Theta}\).

#### ii.2.2 Background: Bayesian inference and Gaussian processes

The parametric viewpoint in Section II.1.1 which yields a point estimate \(\hat{\Theta}\) for the optimal set of parameters can be complemented by considering a Bayesian perspective [26, 21, 25]: For each network input \(x_{\alpha}\), the network equations (1) yield a single output \(y\left(x_{\alpha}|\Theta\right)\). One typically considers a stochastic output \(y\left(x_{\alpha}|\Theta\right)+\xi_{\alpha}\) where the \(\xi_{\alpha}\) are Gaussian independently and identically distributed (i.i.d.) with variance \(\sigma_{\mathrm{reg}}^{2}\)[38]. This regularization allows us to define the probability distribution \(p\left(y|x_{\alpha},\Theta\right)=\left\langle\delta\left[y_{\alpha}-y(x_{ \mathrm{tr},\alpha}|\Theta)-\xi_{\alpha}\right]\right\rangle_{\xi_{\alpha}}= \mathcal{N}\left(y_{\alpha};\,y(x_{\alpha}|\Theta),\sigma_{\mathrm{reg}}^{2}\right)\). An alternative interpretation of \(\xi_{\alpha}\) is a Gaussian noise on the labels. Given a particular set of the network parameters \(\Theta\) this implies a joint distribution \(p\left(\mathbf{y}|\mathbf{x}_{\mathrm{tr}},\Theta\right):=\prod_{\alpha=1}^{D} \left\langle\delta\left[y_{\alpha}-y(x_{\mathrm{tr},\alpha}|\Theta)-\xi_{ \alpha}\right]\right\rangle_{\xi_{\alpha}}=\prod_{\alpha=1}^{D}p(y_{\alpha}|x_{ \alpha},\Theta)\) of network outputs \(\left\{y_{\alpha}\right\}_{1\leq\alpha\leq D}\), each corresponding to one network input \(\left\{x_{\mathrm{tr},\alpha}\right\}_{1\leq\alpha\leq D}\). One aims to use the training data \(\mathcal{D}_{\mathrm{tr}}\) to compute the posterior distribution for the weights \(\mathbf{V},\mathbf{W}^{(1)},\ldots,\mathbf{W}^{(L)},\mathbf{U}\) by conditioning on the network outputs to agree to the desired training values. Concretely, we here assume as a prior for the model parameters that the parameter elements \(V_{ij},W_{ij}^{(l)},U_{ij}\) are i.i.d. according to centered Gaussian distributions \(V_{ij}\sim\mathcal{N}\left(0,\sigma_{v}^{2}/N_{\mathrm{dim}}\right)\), \(W_{ij}^{(l)}\sim\mathcal{N}\left(0,\sigma_{w}^{2}/N_{h}\right)\), and \(U_{ij}\sim\mathcal{N}\left(\sigma_{w}^{2}/N_{\mathrm{out}}\right)\).

The posterior distribution of the parameters \(p\left(\Theta|\mathbf{x}_{\mathrm{tr}},\mathbf{y}_{\mathrm{tr}}\right)\) then follows from Bayes' theorem as

\[p\left(\Theta|\mathbf{x}_{\mathrm{tr}},\mathbf{y}_{\mathrm{tr}}\right)=\frac{p \left(\mathbf{y}_{\mathrm{tr}}|\mathbf{x}_{\mathrm{tr}},\Theta\right)\,p\left( \Theta\right)}{p\left(\mathbf{y}_{\mathrm{tr}}|\mathbf{x}_{\mathrm{tr}}\right)}\,, \tag{2}\]

withthe likelihood \(p\left(\mathbf{y}_{\mathrm{tr}}|\mathbf{x}_{\mathrm{tr}},\Theta\right)\), the weight prior \(p\left(\Theta\right)\) and the model evidence \(p\left(\mathbf{y}_{\mathrm{tr}}|\mathbf{x}_{\mathrm{tr}}\right)=\int d\Theta\,p \left(\mathbf{y}_{\mathrm{tr}}|\mathbf{x}_{\mathrm{tr}},\Theta\right)p\left(\Theta\right)\), which provides the proper normalization. The posterior parameter distribution \(p\left(\Theta|\mathbf{x}_{\mathrm{tr}},\mathbf{y}_{\mathrm{tr}}\right)\) also determines the distribution of the network output \(y_{\star}\) corresponding to a test-point \(x_{\ast}\) by marginalizing over the parameters \(\Theta\)

\[p\left(y_{\ast}|x_{\ast},\mathbf{x}_{\mathrm{tr}},\mathbf{y}_{\mathrm{tr}}\right) =\int d\Theta\,p\left(y_{\ast}|x_{\ast},\Theta\right)p\left(\Theta| \mathbf{x}_{\mathrm{tr}},\mathbf{y}_{\mathrm{tr}}\right)\,, \tag{3}\] \[=\frac{p\left(y_{\ast},\mathbf{y}_{\mathrm{tr}}|x_{\ast},\mathbf{x }_{\mathrm{tr}}\right)}{p\left(\mathbf{y}_{\mathrm{tr}}|\mathbf{x}_{\mathrm{tr}} \right)}\,. \tag{4}\]

One can understand this intuitively: The distribution in (2) provides a set of viable parameters \(\Theta\) based on the training data. An initial guess for the correct choice of parameters via the prior \(p\left(\Theta\right)\) is refined, based on whether the choice of parameters accurately models the relation of the training-data, which is encapsulated in the likelihood \(p\left(\mathbf{y}_{\mathrm{tr}}|\mathbf{x}_{\mathrm{tr}},\Theta\right)\). This viewpoint of Bayesian parameter selection is also equivalent to what is known as Bayesian deep learning [26]. The distribution \(p\left(y_{\ast},\mathbf{y}_{\mathrm{tr}}|x_{\ast},\mathbf{x}_{\mathrm{tr}}\right)\) describes the joint network outputs for all training points and the test point. In the case of wide networks, where \(N_{h}\rightarrow\infty\), [29, 39] showed that the distribution of network outputs \(p\left(y_{\ast},\mathbf{y}_{\mathrm{tr}}|x_{\ast},\mathbf{x}_{\mathrm{tr}}\right)\) approaches a Gaussian process \(y\sim\mathcal{N}\left(0,K^{y}\right)\), where the covariance \(\left\langle y_{\alpha}y_{\beta}\right\rangle=K_{\alpha\beta}^{y}\) is also denoted as the kernel. This is beneficial, as the inference for the network output \(y_{\ast}\) for a test point \(x_{\ast}\) then also follows a Gaussian distribution with mean and covariance given by [32]

\[\left\langle y_{\ast}\right\rangle =K_{\ast\alpha}^{y}\left(K^{y}\right)_{\alpha\beta}^{-1}\,y_{ \mathrm{tr},\beta} \tag{5}\] \[\left\langle\left(y_{\ast}-\left\langle y_{\ast}\right\rangle \right)^{2}\right\rangle =K_{\ast\ast}^{y}-K_{\ast\alpha}^{y}\left(K^{y}\right)_{\alpha\beta}^{-1}\,K_{ \beta\ast}^{y}\,, \tag{6}\]
where summation over repeated indices is implied. There has been extensive research in relating the outputs of wide neural networks to Gaussian processes [3, 19, 29] including recent work on corrections due to finite-width effects \(N_{h}\gg 1\)[27, 20, 2, 44, 35, 45, 44].

### Our contribution

A fundamental assumption of supervised learning is the existence of a joint distribution \(p(x_{\mathrm{tr}},y_{\mathrm{tr}})\) from which the set of training data as well as the set of test data are drawn. In this work we follow the Bayesian approach and investigate the effect of variability in the training data on the generalization properties of wide neural networks. We do so in the kernel limit of infinitely wide linear and non-linear networks. Variability here has two meanings: First, for each drawing of \(D\) pairs of training samples \((x_{\mathrm{tr},\alpha},y_{\mathrm{tr},\alpha})_{1\leq\alpha\leq D}\) one obtains a \(D\times D\) kernel matrix \(K^{y}\) with heterogeneous entries; so in a single instance of Bayesian inference, the entries of the kernel matrix vary from one entry to the next. Second, each such drawing of \(D\) training data points and one test data point \((x_{*},y_{*})\) leads to a different kernel \(\{K^{y}_{\alpha\beta}\}_{1\leq\alpha,\beta\leq D+1}\), which follows some probabilistic law \(K^{y}\sim p(K^{y})\).

Our work builds upon previous results for the NNGP limit to formalize the influence of such stochastic kernels. We here develop a field theoretic approach to systematically investigate the influence of the underlying kernel stochasticity on the generalization properties of the network, namely the learning curve, the dependence of \(\langle y_{*}\rangle\) on the number of training samples \(D=|\mathcal{D}_{\mathrm{tr}}|\). As we assume Gaussian i.i.d. priors on the network parameters, the output kernel \(K^{y}_{\alpha\beta}\) solely depends on the network architecture and the input overlap matrix

\[K^{x}_{\alpha\beta}=\sum_{i=1}^{N_{\mathrm{dim}}}x_{\alpha i}x_{\beta i}\quad x _{\alpha},x_{\beta}\in\mathcal{D}_{\mathrm{tr}}\cup\mathcal{D}_{\mathrm{test}}\,, \tag{7}\]

with \(\alpha,\beta=1...D+1\). We next define a data model which allows us to approximate the probability measure for the data variability.

### Definition of a synthetic data set

To investigate the generalization properties in a binary classification task, we introduce a synthetic stochastic binary classification task. This task allows us to control the statistical properties of the data with regard to the dimensionality of the patterns, the degree of separation between patterns belonging to different classes, and the variability in the kernel. Moreover, it allows us to construct training-data sets \(\mathcal{D}_{\mathrm{tr}}\) of arbitrary sizes and we will show that the statistics of the resulting kernels is indeed representative for more realistic data sets such as MNIST.

The data set consists of pattern realizations \(x_{\alpha}\in\{-1,1\}^{N_{\mathrm{dim}}}\) with dimension \(N_{\mathrm{dim}}\) even. We denote the entries \(x_{\alpha,i}\) of this \(N_{\mathrm{dim}}\)-dimensional vector for data point \(\alpha\) as pixels that randomly take either of two values \(x_{\alpha,i}\in\{-1,1\}\) with respective probabilities \(p(x_{\alpha,i}=1)\) and \(p(x_{\alpha,i}=-1)\) that depend on the class \(c(\alpha)\in\{1,2\}\) of the pattern realization and whether the index \(i\) is in the left half ( \(i\leq N_{\mathrm{dim}}/2\)) or the right half ( \(i>N_{\mathrm{dim}}/2\)) of the pattern: For class \(c(\alpha)=1\) each pixel \(x_{\alpha,1\leq i\leq N_{\mathrm{dim}}}\) is realized independently as a binary variable as

\[x_{\alpha,i} =\begin{cases}1&\text{with }p\\ -1&\text{with }(1-p)\end{cases}\quad\text{for }i\leq\frac{N_{\mathrm{dim}}}{2}\,, \tag{8}\] \[x_{\alpha,i} =\begin{cases}1&\text{with }(1-p)\\ -1&\text{with }p\end{cases}\quad\text{for }i>\frac{N_{\mathrm{dim}}}{2}\,. \tag{9}\]

For a pattern \(x_{\alpha}\) in the second class \(c(\alpha)=2\) the pixel values are distributed independently of those in the first class with a statistics that equals the negative pixel values of the first class, which is \(P\left(x_{\alpha i}\right)=P\left(-x_{\beta i}\right)\) with \(c(\beta)=1\). The number of training samples \(N_{\mathrm{dim}}\) is \(N_{\mathrm{dim}}\). The number of training samples \(N_{\mathrm{dim}}\) is \(N_{\mathrm{dim}}\).
\(1\) and \(c(\alpha)=2\). There are two limiting cases for \(p\) which illustrate the construction of the patterns: In the limit \(p=1\), each pattern \(x_{\alpha}\) in \(c=1\) consists of a vector, where the first \(N_{\text{dim}}/2\) pixels have the value \(x_{\alpha i}=1\), whereas the second half consists of pixels with the value \(x_{\alpha,i}=-1\). The opposite holds for patterns in the second class \(c=2\). This limiting case is shown in Figure 2 (right column). In the limit case \(p=0.5\) each pixel assumes the value \(x_{\alpha,i}=\pm 1\) with equal probability, regardless of the pattern class-membership or the pixel position. Hence one cannot distinguish the class membership of any of the training instances. This limiting case is shown in Figure 2 (left column). If \(c(\alpha)=1\) we set \(y_{\text{tr},\alpha}=-1\) and for \(c(\alpha)=2\) we set \(y_{\text{tr},\alpha}=1\).

We now investigate the description of this task in the framework of Bayesian inference. The hidden variables \(h_{\alpha}^{0}\) (1) in the input layer under a Gaussian prior on \(V_{ij}\overset{\text{i.i.d.}}{\sim}\mathcal{N}\left(0,\sigma_{v}^{2}/N_{\text {dim}}\right)\) follow a Gaussian process with kernel \(K^{(0)}\) given by

\[K_{\alpha\beta}^{0} =\langle h_{\alpha}^{0}h_{\beta}^{0}\rangle_{V\sim\mathcal{N} \left(0,\frac{\sigma_{v}^{2}}{N_{\text{dim}}}\right)}\,, \tag{10}\] \[=\frac{\sigma_{v}^{2}}{N_{\text{dim}}}\,\sum_{i=1}^{N_{\text{dim }}}x_{\alpha i}x_{\beta i}\,. \tag{11}\]

Separability of the two classes is reflected in the structure of this input kernel \(K^{0}\) as shown in Figure 2: In the cases with \(p=0.8\) and \(p=1\) one can clearly distinguish blocks; the diagonal blocks represent intra-class overlaps, the off-diagonal blocks inter-class overlaps. This is not the case for \(p=0.5\), where no clear block-structure is visible. In the case of \(p=0.8\) one can further observe that the blocks are not as clear-cut as in the case \(p=1\), but rather noisy, similar to \(p=0.5\). This is due to the probabilistic realization of patterns, which induces stochasticity in the blocks of the input kernel \(K^{0}\) (10). To quantify this effect, based on the distribution of the pixel values (9) we compute the distribution of the entries of \(K^{0}\) for the binary classification task. The mean of the overlap elements \(\mu_{\alpha\beta}\) and their covariances \(\Sigma_{(\alpha\beta)(\gamma\delta)}\) are defined via

\[\mu_{\alpha\beta} =\left\langle K_{\alpha\beta}^{0}\right\rangle\,, \tag{12}\] \[\Sigma_{(\alpha\beta)(\gamma\delta)} =\left\langle\delta K_{\alpha\beta}^{0}\,\delta K_{\gamma\delta }^{0}\right\rangle\,,\] (13) \[\delta K_{\alpha\beta}^{0} =K_{\alpha\beta}^{0}-\mu_{\alpha\beta}\,, \tag{14}\]

where the expectation value \(\left\langle\cdot\right\rangle\) is taken over drawings of \(D\) training samples each. By construction we have \(\mu_{\alpha\beta}=\mu_{\beta\alpha}\). The covariance is further invariant under the exchange of \((\alpha,\beta)\leftrightarrow(\gamma,\delta)\) and, due to the symmetry of \(K_{\alpha\beta}^{0}=K_{\beta\alpha}^{0}\), also under swapping \(\alpha\leftrightarrow\beta\) and \(\gamma\leftrightarrow\delta\) separately. In the artificial task-setting, the parameter \(p\), the pattern dimensionality \(N_{\text{dim}}\), and the variance \(\sigma_{v}^{2}/N_{\text{dim}}\) of each read-in weight \(V_{ij}\) define the elements of \(\mu_{\alpha\beta}\) and \(\Sigma_{(\alpha\beta)(\gamma\delta)}\), which read

\[\mu_{\alpha\beta} =\sigma_{v}^{2}\begin{cases}1&\alpha=\beta\\ u&c_{\alpha}=c_{\beta}\\ -u&c_{\alpha}\neq c_{\beta}\end{cases}\,,\] \[\Sigma_{(\alpha\beta)(\alpha\beta)} =\frac{\sigma_{v}^{4}}{N_{\text{dim}}}\kappa,\] \[\Sigma_{(\alpha\beta)(\alpha\delta)} =\frac{\sigma_{v}^{4}}{N_{\text{dim}}}\begin{cases}\nu&\text{ for }\begin{cases}c_{\alpha}=c_{\beta}=c_{\delta}\\ c_{\alpha}\neq c_{\beta}=c_{\delta}\\ -\nu&\text{ for }\begin{cases}c_{\alpha}=c_{\beta}\neq c_{\delta}\\ c_{\alpha}=c_{\delta}\neq c_{\beta}\end{cases}\end{cases}\,,\] \[\text{with}\quad\kappa :=1-u^{2}\,,\] \[\nu :=u\left(1-u\right),\] \[u :=4p(p-1)+1\,. \tag{15}\]

In addition to this, the tensor elements of \(\Sigma_{(\alpha\beta)(\gamma\delta)}\) are zero for the following index combinations because we fixed the value of \(K_{\alpha\alpha}^{0}\) by construction:

\[\Sigma_{(\alpha\beta)(\gamma\delta)} =0\quad\text{with}\quad\alpha\neq\beta\neq\gamma\neq\delta\,,\] \[\Sigma_{(\alpha\alpha)(\beta\gamma)} =0\quad\text{with}\quad\alpha\neq\beta\neq\gamma\,,\] \[\Sigma_{(\alpha\alpha)(\beta\beta)} =0\quad\text{with}\quad\alpha\neq\beta\,,\] \[\Sigma_{(\alpha\alpha)(\alpha\beta)} =0\quad\text{with}\quad\alpha\neq\beta\,,\] \[\Sigma_{(\alpha\alpha)(\alpha\alpha)} =0\quad\text{with}\quad\alpha\neq\beta\,. \tag{16}\]

The expressions for \(\Sigma_{(\alpha\beta)(\alpha\beta)}\) and \(\Sigma_{(\alpha\beta)(\alpha\delta)}\) in (15) show that the magnitude of the fluctuations are controlled through the parameter \(p\) and the pattern dimensionality \(N_{\text{dim}}\): The covariance \(\Sigma\) is suppressed by a factor of \(1/N_{\text{dim}}\) compared to the mean values \(\mu\). Hence we can use the pattern dimensionality \(N_{\text{dim}}\) to investigate the influence of the strength of fluctuations. As illustrated in Figure 1a, the elements \(\Sigma_{(\alpha\beta)(\alpha\beta)}\) denote the variance of individual entries of the kernel, while \(\Sigma_{(\alpha\beta)(\alpha\gamma)}\) are covariances of entries across elements of a given row \(\alpha\), visible as horizontal or vertical stripes in the color plot of the kernel.

Equation (15) implies, by construction, a Gaussian distribution of the elements \(K_{\alpha\beta}^{0}\) as it only provides the first two cumulants. One can show that the higher-order cumulants of \(K_{\alpha\beta}^{0}\) scale sub-leading in the pattern dimension and are hence suppressed by a factor \(\mathcal{O}\left(1/N_{\text{dim}}\right)\) compared to \(\Sigma_{(\alpha\beta)(\gamma\delta)}\).

## III Results

In this section we derive the field theoretic formalism which allows us to compute the statistical properties of the inferred network output in Bayesian inference with a stochastic kernel. We show that the resulting process is non-Gaussian and reminiscent of a \(\varphi^{3}+\varphi^{4}\)-theory. Specifically, we compute the mean of the predictive distribution
of this process conditioned on the training data. This is achieved by employing systematic approximations with the help of Feynman diagrams.

Subsequently we show that our results provide an accurate bound on the generalization capabilities of the network. We further discuss the implications of our analytic results for neural architecture search.

### Field theoretic description of Bayesian inference

#### iii.1.1 Bayesian inference with stochastic kernels

In general, a network implements a map from the inputs \(x_{\alpha}\) to corresponding outputs \(y_{\alpha}\). In particular a model of the form (1) implements a non-linear map \(\psi:\mathbb{R}^{N_{\text{dim}}}\to\mathbb{R}^{N_{h}}\) of the input \(x_{\alpha}\in\mathbb{R}^{N_{\text{dim}}}\) to a hidden state \(h_{\alpha}\in\mathbb{R}^{N_{h}}\). This map may also involve multiple hidden-layers, biases and non-linear transformations. The read-out weight \(\mathbf{U}\in\mathbb{R}^{1\times N_{h}}\) links the scalar network output \(y_{\alpha}\in\mathbb{R}\) and the transformed inputs \(\psi\left(x_{\alpha}\right)\) with \(1\leq\alpha\leq D_{\text{tot}}=D+D_{\text{test}}\) which yields

\[y_{\alpha}=\mathbf{U}\,\psi\left(x_{\alpha}\right)+\xi_{\alpha}\,, \tag{17}\]

where \(\xi_{\alpha}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,\sigma_{\text{reg}}^{2})\) is a regularization noise in the same spirit as in [38]. We assume that the prior on the read-out vector elements is a Gaussian \(\mathbf{U}_{i}\overset{\text{i.i.d.}}{\sim}\mathcal{N}\left(0,\sigma_{u}^{2}/ N_{h}\right)\). The distribution of the set of network outputs \(y_{1\leq\alpha\leq D_{\text{tot}}}\) is then in the limit \(N_{h}\to\infty\) a multivariate Gaussian [29]. The kernel matrix of this Gaussian is obtained by taking the expectation value with respect to the read-out vector, which yields

\[\left\langle y_{\alpha}\,y_{\beta}\right\rangle_{\mathbf{U}} =:K_{\alpha\beta}^{y}=\sigma_{u}^{2}\,K_{\alpha\beta}^{\psi}+ \delta_{\alpha\beta}\,\sigma_{\text{reg}}^{2}\,, \tag{18}\] \[K_{\alpha\beta}^{\psi} =\frac{1}{N_{h}}\,\sum_{i=1}^{N_{h}}\psi_{i}\left(x_{\alpha} \right)\,\psi_{i}\left(x_{\beta}\right)\,. \tag{19}\]

The kernel matrix \(K_{\alpha\beta}^{y}\) describes the covariance of the network's output and hence depends on the kernel matrix \(K_{\alpha\beta}^{\psi}\). The additional term \(\delta_{\alpha\beta}\,\sigma_{\text{reg}}^{2}\) acts as a regularization term, which is also known as a ridge regression [14] or Tikhonov regularization [41]. In the context of neural networks one can motivate the regularizer \(\sigma_{\text{reg}}^{2}\) by using the \(L^{2}\)-regularization in the readout layer. This is also known as weight decay [12]. Introducing the regularizer \(\sigma_{\text{reg}}^{2}\) is necessary to ensure that one can properly invert the matrix \(K_{\alpha\beta}^{y}\).

Different drawings of sets of training data \(\mathcal{D}_{\text{tr}}\) lead to different realizations of kernel matrices \(K^{\psi}\) and \(K^{y}\). The network output \(y_{\alpha}\) hence follows a multivariate Gaussian with a stochastic kernel matrix \(K^{y}\). A more formal derivation of the Gaussian statistics, including an argument for its validity in deep neural networks, can be found in [19]. A consistent derivation using field theoretical methods and corrections in terms for the width of the hidden layer \(N_{h}\) for deep and recurrent networks has been presented in [35].

In general, the input kernel matrix \(K^{0}\) (10) and the output kernel matrix \(K^{y}\) are related in a non-trivial fashion, which depends on the specific network architecture at hand. From now on we make an assumption on the stochasticity of \(K^{0}\) and assume that the input kernel matrix \(K^{0}\) is distributed according to a multivariate Gaussian

\[K^{0}\sim\mathcal{N}\left(\mu,\,\Sigma\right)\,, \tag{20}\]

where \(\mu\) and \(\Sigma\) are given by (12) and (13), respectively.

In the limit of large pattern dimensions \(N_{\text{dim}}\gg\)1 this assumption is warranted for the kernel matrix \(K^{0}\). This structure further assumes, that the overlap statistics are unimodal, which is indeed mostly the case for data such as MNIST. Furthermore we assume that this property holds for the output kernel matrix \(K^{y}\) as well and that we can find a mapping from the mean \(\mu\) and covariance \(\Sigma\) of the input kernel to the mean \(m\) and covariance \(C\) of the output kernel \(\left(\mu_{\alpha\beta},\,\Sigma_{\left(\alpha\beta\right)\left(\gamma\delta \right)}\right)\to\left(m_{\alpha\beta},\,C_{\left(\alpha\beta\right)\left( \gamma\delta\right)}\right)\) so that \(K^{y}\) is also distributed according to a multivariate Gaussian

\[K^{y}\sim\mathcal{N}(m,\,C). \tag{21}\]

For each realization \(K_{\alpha\beta}^{y}\), the joint distribution of the network outputs \(y_{1\leq\alpha\leq D_{\text{tot}}}\) corresponding to the training and test data points \(\mathbf{x}\) follow a multivariate Gaussian

\[p\left(\mathbf{y}|\mathbf{x}\right)\sim\mathcal{N}\!\left(0,K^{y}\right). \tag{22}\]

The kernel allows us to compute the conditional probability \(p\left(y_{*}|\mathbf{x}_{\text{tr}},\mathbf{y}_{\text{tr}},x_{*}\right)\) (3) for a test point \(\left(x_{*},y_{*}\right)\in\mathcal{D}_{\text{test}}\) conditioned on the data from the training set \(\left(\mathbf{x}_{\text{tr}},\mathbf{y}_{\text{tr}}\right)\in\mathcal{D}_{ \text{tr}}\). This distribution is Gaussian with mean and variance given by (5) and (6), respectively. It is our goal to take into account that \(K^{0}\) is a stochastic quantity, which depends on the particular draw of the training and test data set \(\left(\mathbf{x}_{\text{tr}},\mathbf{y}_{\text{tr}}\right)\in\mathcal{D}_{ \text{tr}},\left(x_{*},y_{*}\right)\in\mathcal{D}_{\text{test}}\). The labels \(\mathbf{y}_{\text{tr}},y_{*}\) are, by construction, deterministic and take either one of the values \(\pm 1\). In the following we investigate the mean of the predictive distribution on the number of training samples, which we call the learning curve. A common assumption is that this learning curve is rather insensitive to the very realization of the chosen training points. Thus we assume that the learning curve is self-averaging.. The mean computed for a single draw of the training data is hence expected to agree well to the average over many such drawings. Under this assumption it is sufficient to compute the data-averaged mean inferred network output, which reduces to computing the disorder-average of the following quantity

\[\left\langle y_{*}\right\rangle_{K^{y}}=\left\langle K_{*\alpha}^{y}\left[K^{ y}\right]_{\alpha\beta}^{-1}\right\rangle_{K^{y}}y_{\beta}\,. \tag{23}\]
To perform the disorder average and to compute perturbative corrections, we will follow these steps

* construct a suitable dynamic moment-generating function \(Z_{K^{y}}(l_{*})\) with the source term \(l_{*}\),
* propagate the input stochasticity to the network output \(K^{0}_{\alpha\beta}\to K^{y}_{\alpha\beta}\),
* disorder-average the functional using the model \(K^{y}_{\alpha\beta}\sim\mathcal{N}(m_{\alpha\beta},\,C_{(\alpha\beta)(\gamma \delta)})\),
* and finally perform the computation of perturbative corrections using diagrammatic techniques.

#### ii.2.2 Constructing the dynamic moment generating function \(Z_{K^{y}}(l_{*})\)

Our ultimate goal is to compute learning curves. Therefore we want to evaluate the disorder averaged mean inferred network output (23). Both the presence of two correlated random matrices and the fact that one of the matrices appears as an inverse complicate this process. One alternative route is to define the moment-generating function

\[Z(l_{*}) =\int dy_{*}\exp(l_{*}y_{*})p\left(y_{*}|x_{*},\mathbf{x}_{\text{ tr}},\mathbf{y}_{\text{tr}}\right)\,, \tag{24}\] \[=\frac{\int dy_{*}\exp(l_{*}y_{*})p\left(y_{*},\mathbf{y}_{\text{ tr}}|x_{*},\mathbf{x}_{\text{tr}}\right)}{p\left(\mathbf{y}_{\text{tr}}| \mathbf{x}_{\text{tr}}\right)}\,,\] (25) \[=:\frac{\mathcal{Z}(l_{*})}{\mathcal{Z}(0)}\,, \tag{26}\]

with joint Gaussian distributions \(p\left(y_{*},\mathbf{y}_{\text{tr}}|x_{*},\mathbf{x}_{\text{tr}}\right)\) and \(p\left(\mathbf{y}_{\text{tr}}|\mathbf{x}_{\text{tr}}\right)\) that each can be readily averaged over \(K^{y}\). Equation (23) is then obtained as

\[\left\langle y_{*}\right\rangle_{K^{y}}=\frac{\partial}{\partial l_{*}}\left\langle \frac{\mathcal{Z}(l_{*})}{\mathcal{Z}(0)}\right\rangle_{K^{y}}\Bigg{|}_{l_{*}= 0}\,. \tag{27}\]

A complication of this approach is that the numerator and denominator co-fluctuate. The common route around this problem is to consider the cumulant-generating function \(W(l_{*})=\ln\mathcal{Z}(l_{*})\) and to obtain \(\left\langle y_{*}\right\rangle_{K^{y}}=\frac{\partial}{\partial l_{*}}\left\langle W (l_{*})\right\rangle_{K^{y}}\), which, however, requires averaging the logarithm. This is commonly done with the replica trick [8; 23].

We here follow a different route to ensure that the disorder-dependent normalization \(\mathcal{Z}(0)\) drops out and construct a dynamic moment generating function [5]. Our goal is hence to design a dynamic process where a time dependent observable is related to our mean-inferred network output \(y_{*}\). We hence define the linear process in the auxiliary variables \(q_{\alpha}\)

\[\frac{\partial q_{\alpha}(t)}{\partial t}=-K^{y}_{\alpha\beta}\,q_{\beta}(t)+ y_{\alpha}\,, \tag{28}\]

for \((x_{\alpha},y_{\alpha})\in\mathcal{D}_{\text{tr}}\). From this we see directly that \(q_{\alpha}(t\to\infty)=\left[K^{y}\right]^{-1}_{\alpha\beta}y_{\beta}\) is a fixpoint. The fact that \(K^{y}_{\alpha\beta}\) is a covariance matrix ensures that it is positive semi-definite and hence implies the convergence to a fixpoint. We can obtain (5) \(\left\langle y_{*}\right\rangle=K^{y}_{\alpha\alpha}\left[K^{y}\right]^{-1}_{ \alpha\beta}y_{\beta}\) from (28) as a linear readout of \(q_{\alpha}(t\to\infty)\) with the matrix \(K^{y}_{\alpha\alpha}\). Using the Martin-Siggia-Rose-deDominicis-Janssen formalism [13; 16; 22; 37] one can express this as the first derivative of the moment generating function \(Z_{K^{y}}(l_{*})\) in frequency space

\[Z_{K^{y}}(l_{*}) =\int\mathcal{DQ}\tilde{Q}\exp\left(S(Q,\tilde{Q},l_{*})\right)\,, \tag{29}\] \[S(Q,\tilde{Q},l_{*}) =\tilde{Q}^{\top}_{\alpha}\big{(}-i\omega\mathbb{I}+K^{y}\big{)}_ {\alpha\beta}\,Q_{\beta}\] (30) \[-\tilde{Q}(\omega=0)_{\alpha}y_{\alpha}\] \[+l_{*}K^{y}_{\alpha,\alpha}Q_{\alpha}(\omega=0)\,, \tag{31}\]

where \(\tilde{Q}^{\top}_{\alpha}\left(\cdots\right)Q_{\beta}=\frac{1}{2\pi}\int d \omega\tilde{Q}_{\alpha}(\omega)\left(\cdots\right)Q_{\beta}(-\omega)\). As \(Z_{K^{y}}(l_{*})\) is normalized such that \(Z_{K^{y}}(0)=1\quad\forall\,K^{y}\), we can compute (23) by evaluating the derivative of the disorder-averaged moment-generating function \(\overline{Z}(l_{*})\)

\[\overline{Z}(l_{*}) =\left\langle\int\mathcal{D}\{Q,\tilde{Q}\}\exp\left(S(Q,\tilde{Q },l_{*})\right)\right\rangle_{K^{y}}\,, \tag{32}\] \[\left\langle y_{*}\right\rangle_{K^{y}} =\frac{\partial\overline{Z}(l_{*})}{\partial l_{*}}\bigg{|}_{l_{*}= 0}\,. \tag{33}\]

By construction the distribution of the kernel matrix entries \(K^{y}_{\alpha\beta}\) is a multivariate Gaussian (20). In the following we will treat the stochasticity of \(K^{y}_{\alpha\beta}\) perturbatively to gain insights into the influence of input stochasticity.

ii.2.3 Perturbative treatment of the disorder averaged moment generating function \(\overline{Z}(l_{*})\)

To compute the disorder averaged mean-inferred network output (23) we need to compute the disorder average of the dynamic moment generating function \(\overline{Z}(l_{*})\) and its derivative at \(l_{*}=0\). Due to the linear appearance of \(K^{y}\) in the action (30) and the Gaussian distribution for \(K^{y}\) (21) we can do this directly and obtain the action

\[\overline{Z}(l_{*}) =\int\mathcal{DQ}\mathcal{D}\tilde{Q}\,\left\langle\exp\left(S \right)\right\rangle_{K^{y}}\,,\] \[=\int\mathcal{DQ}\mathcal{D}\tilde{Q}\,\exp\left(\overline{S} \right)\,, \tag{34}\] \[\overline{S}(Q,\tilde{Q},l_{*}) =\tilde{Q}^{\top}\left(-i\omega\mathbb{I}+m\right)Q\] \[-\tilde{Q}^{0}_{\eta}\,y_{\eta}\] \[+l_{*}m_{*\epsilon}Q^{0}_{\epsilon}\] \[+\frac{1}{2}\tilde{Q}^{\top}_{\alpha}Q_{\beta}\,C_{(\alpha\beta)( \gamma\delta)}\,\tilde{Q}^{\top}_{\gamma}Q_{\delta}\] \[+l_{*}C_{(\epsilon\alpha)(\beta)}(\beta)^{0}_{\alpha}\tilde{Q}^{ \top}_{\beta}Q_{\gamma}\,, \tag{35}\]
with \(Q^{0}:=Q(\omega=0)\) and \(\tilde{Q}^{0}:=\tilde{Q}(\omega=0)\).

As we ultimately aim to obtain corrections for the mean inferred network output \(\langle y_{*}\rangle\), we utilize the action in (35) and established results from field theory to derive the leading order terms as well as perturbative corrections diagrammatically. The presence of the variance and covariance terms in (35) introduces corrective factors, which cannot appear in the 0th-order approximation, which corresponds to the homogeneous kernel that neglects fluctuations in \(K^{y}\) by setting \(C_{(\alpha\beta)(\gamma\delta)}=0\). This will provide us with the tools to derive an asymptotic bound for the mean inferred network output \(\langle y^{*}\rangle\) in the case of an infinitely large training data set. This bound is directly controlled by the variability in the data. We provide empirical evidence for our theoretical results for linear, non-linear, and deep-kernel-settings and show how the results could serve as indications to aid neural architecture search based on the statistical properties of the underlying data set.

iii.2.4 Field theoretic elements to compute the mean inferred network output \(\langle y_{*}\rangle\)

The field theoretic description of the inference problem in form of an action (35) allows us to derive perturbative expressions for the statistics of the inferred network output \(\langle y_{*}\rangle_{K^{y}}\) in a diagrammatic manner. This diagrammatic treatment for perturbative calculations is a powerful tool and is standard practice in statistical physics [46], data analysis and signal reconstruction [7], and more recently in the investigation of artificial neural networks [6].

Comparing the action (35) to prototypical expressions from classical statistical field theory such as the \(\varphi^{3}+\varphi^{4}\) theory[13; 46] one can similarly associate the elements of a field theory:

* \(-\tilde{Q}^{0}_{\alpha}y_{\alpha}\doteq\) is a monopole term
* \(l_{*}m_{*\epsilon}Q^{0}_{\epsilon}\doteq\) is a source term
* \(\Delta_{\alpha\beta}:=(i\omega\mathbb{I}-m)^{-1}_{\alpha\beta}\doteq\) is a propagator that connect the fields \(Q_{\alpha}(\omega),\tilde{Q}_{\beta}\,(-\omega)\)
* \(l_{*}C_{(\ast\alpha)(\beta\gamma)}Q^{0}_{\alpha}\tilde{Q}^{\top}_{\beta}Q_{ \gamma}\doteq\) is a three-point vertex
* \(\frac{1}{2}\tilde{Q}^{\top}_{\alpha}Q_{\beta}\,C_{(\alpha\beta)(\gamma\delta) }\,\tilde{Q}^{\top}_{\gamma}Q_{\delta}\doteq\) is a four-point vertex.

The following rules for Feynman diagrams simplify calculations:

1. To obtain corrections to first order in \(C\sim\mathcal{O}\left(1/N_{\text{dim}}\right)\), one has to compute all diagrams with a single vertex (three-point or four-point) [13]. This approach assumes that the interaction terms \(C_{(\alpha\beta)(\gamma\delta)}\) that stem from the variability of the data are small compared to the mean \(m_{\alpha,\beta}\). In the case of strong heterogeneity one cannot use a conventional expansion in the number of vertices \(C_{(\alpha\beta)(\gamma\delta)}\) and would have to resort to other methods.
2. Vertices, source terms, and monopoles have to be connected with one another using the propagator \(\Delta_{\alpha\beta}=(i\omega\mathbb{I}-m)^{-1}_{\alpha\beta}\) which couple \(Q_{\alpha}(\omega)\) and \(\tilde{Q}_{\beta}(-\omega)\) which each other.
3. We only need diagrams with a single external source term \(l_{*}\) because we seek corrections to the mean-inferred network output. Because the source \(l_{*}\) couples to the \(\omega=0\) component \(Q^{0}\) of the field \(Q\), propagators to these external legs are evaluated at \(\omega=0\), thus replacing \((i\omega\mathbb{I}-m)^{-1}_{\alpha\beta}\rightarrow-(m^{-1})_{\alpha\beta}\).
4. The structure of the integrals appearing in the four-point and three-point vertices containing \(C_{(\alpha\beta)(\gamma\delta)}\) with contractions by \(\Delta_{\alpha,\beta}\) or \(\Delta_{\gamma,\delta}\) within a pair of indices \((\alpha\beta)\) or \((\gamma\delta)\) yield vanishing contributions; such diagrams are known as closed response loops [13]. This is because the propagator \(\Delta_{\alpha,\beta}(t-s)\) in time domain vanishes for \(t=s\), which corresponds to the integral \(\int d\omega\,\Delta_{\alpha,\beta}(\omega)\) over all frequencies \(\omega\).
5. As we have frequency conservation at the vertices in the form \(\frac{1}{2}\tilde{Q}^{\top}_{\alpha}Q_{\beta}\,C_{(\alpha,\beta)(\gamma,\delta) }\,\tilde{Q}^{\top}_{\gamma}Q_{\delta}\) and since by point 4. above we only need to consider contractions by \(\Delta_{\beta\gamma}\) or \(\Delta_{\delta\alpha}\) by attaching the external legs all frequencies are constrained to \(\omega=0\), so also propagators within a loop are replaced by \(\Delta_{\alpha\beta}=(i\omega\mathbb{I}-m)^{-1}_{\alpha\beta}\rightarrow-(m^{-1} )_{\alpha\beta}\).

These rules directly yield that the corrections for the disorder averaged mean-inferred network to first order in \(C_{(\alpha\beta)(\gamma\delta)}\) can only include the diagrams

\[\langle y_{*}\rangle \doteq\] \[+\ \mathcal{O}\left(C^{2}\right) \tag{36}\]

which translate to our main result

\[\langle y_{*}\rangle_{0+1} =m_{*\alpha}m^{-1}_{\alpha\beta}\,y_{\beta}\] \[+m_{*\epsilon}m^{-1}_{c\alpha}C_{(\alpha\beta)(\gamma\delta)}m^{- 1}_{\beta\gamma}\,m^{-1}_{\delta\rho}\,y_{\rho}\] \[-C_{(\ast\alpha)(\beta\gamma)}m^{-1}_{\alpha\beta}m^{-1}_{\gamma \delta}\,y_{\delta}+\mathcal{O}\left(C^{2}\right)\,. \tag{37}\]

We here define the first line as the zeroth-order approximation \(\langle y_{*}\rangle_{0}:=m_{*\alpha}m^{-1}_{\alpha\beta}\,y_{\beta}\), which has the same form as (5), and the latter two lines as perturbative corrections \(\langle y_{*}\rangle_{1}=\mathcal{O}\left(C\right)\) which are of linear order in \(C\).
Evaluation of expressions for block-structured overlap matrices

To evaluate the first order correction \(\left\langle y_{*}\right\rangle_{1}\) in (37) we make use of the fact that Bayesian inference is insensitive to the order in which the training data are presented. We are hence free to assume that all training samples of one class are presented en bloc. Moreover, supervised learning assumes that all training samples are drawn from the same distribution. As a result, the statistics is homogeneous across blocks of indices that belong to the same class. The propagators \(-m_{\alpha\beta}^{-1}\) and interaction vertices \(C_{(\alpha\beta)(\gamma\delta)}\) and \(C_{(*\alpha)(\beta\gamma)}\), correspondingly, have a block structure. To obtain an understanding how variability of the data and hence heterogeneous kernels affect the ability to make predictions, we consider the simplest yet non-trivial case of binary classification where we have two such blocks.

In this section we focus on the overlap statistics given by the artificial data set described in Section II.3. This data set entails certain symmetries. Generalizing the expressions to a less symmetric task is straightforward. For the classification task, with two classes \(c_{\alpha}\in\{1,2\}\), the structure for the mean overlaps \(\mu_{\alpha\beta}\) and their covariance \(\Sigma_{(\alpha\beta)(\gamma\delta)}\) at the read-in layer of the network given by (15) are inherited by the mean \(m_{\alpha\beta}\) and the covariance \(C_{(\alpha\beta)(\gamma\delta)}\) of the overlap matrix at the output of the network. In particular, all quantities can be expressed in terms of only four parameters \(a\), \(b\), \(K\), \(v\) whose values, however, depend on the network architecture and will be given for linear and non-linear networks below. For four indices \(\alpha,\beta,\gamma,\delta\) that are all different

\[m_{\alpha\alpha} =a\,,\] \[m_{\alpha\beta} =\begin{cases}b&c_{\alpha}=c_{\beta}\\ -b&c_{\alpha}\neq c_{\beta}\end{cases}\,,\] \[C_{(\alpha\alpha),(\gamma\delta)} =0\,,\] \[C_{(\alpha\beta)(\alpha\beta)} =K\,,\] \[C_{(\alpha\beta)(\alpha\delta)} =\begin{cases}v&c_{\alpha}=c_{\beta}=c_{\delta};\quad c_{\alpha} \neq c_{\beta}=c_{\delta}\\ -v&c_{\alpha}=c_{\beta}\neq c_{\delta};\quad c_{\alpha}=c_{\delta}\neq c_{ \beta}\end{cases}\,. \tag{38}\]

This symmetry further assumes that the network does not have biases and utilizes point-symmetric activation functions \(\phi(x)\) such as \(\phi(x)=\text{erf}(\text{x})\). In general, all tensors are symmetric with regard to swapping \(\alpha\leftrightarrow\beta\) as well as \(\gamma\leftrightarrow\delta\) and the tensor \(C_{(\alpha\beta)(\gamma\delta)}\) is invariant under swaps of the index-pairs \((\alpha\beta)\leftrightarrow(\gamma\delta)\). We further assume that the class label for class 1 is \(y\) and that the class label for class 2 is \(-y\). In subsequent calculations and experiments we consider the prediction for the class \(y=-1\).

This setting is quite natural, as it captures the presence of differing mean intra- and inter-class overlaps. Further \(K\) and \(v\) capture two different sources of variability. Whereas \(K\) is associated with the presence of i.i.d. distributed variability on each entry of the overlap matrix separately, \(v\) corresponds to variability stemming from correlations between different patterns. Using the properties in (38) one can evaluate (37) for the inference of test-points \(*\) within class \(c_{1}\) on a balanced training set with \(D\) samples explicitly to

\[\left\langle y_{*}\right\rangle_{0} =Dgy\,, \tag{39}\] \[\left\langle y_{*}\right\rangle_{1} =vg\hat{y}\left(q_{1}+3q_{2}\right)\left(D^{3}-3D^{2}+2D\right)\] \[+Kg\hat{y}\left(q_{1}+q_{2}\right)\left(D^{2}-D\right)\] \[-v\hat{y}\left(q_{1}+q_{2}\right)\left(D^{2}-D\right)\] \[+\mathcal{O}\left(C_{(\alpha\beta)(\gamma\delta)}^{2}\right)\quad \text{for}\quad*\in c_{1} \tag{40}\]

with the additional variables

\[g =\frac{b}{(a-b)+bD}\,,\] \[q_{2} =-\frac{1}{(a-b)+bD}\,,\] \[q_{1} =\frac{1}{a-b}+q_{2}\,,\] \[\hat{y} =\frac{y}{(a-b)+bD}\,,\] \[g =\frac{b}{(a-b)+bD}\,, \tag{41}\]

which stem from the analytic inversion of block-matrices. Carefully treating the dependencies of the parameters in (41) and (40), one can compute the limit \(D\gg 1\) and show that the \(\mathcal{O}(1)\)-behavior of (40) for test points \(*\in c_{1}\) for the zeroth-order approximation, \(\lim_{D\to\infty}\left\langle y_{*}\right\rangle_{0}:=\left\langle y_{*} \right\rangle_{0}^{(\infty)}\), and the first-order correction, \(\lim_{D\to\infty}\left\langle y_{*}\right\rangle_{1}:=\left\langle y_{*} \right\rangle_{1}^{(\infty)}\), is given by

\[\left\langle y_{*}\right\rangle_{0}^{(\infty)} =y\;, \tag{42}\] \[\left\langle y_{*}\right\rangle_{1}^{(\infty)} =\frac{y}{(a-b)b}\left((K-4v)-v\frac{a-b}{b}\right)\;. \tag{43}\]

This result implies that regardless of the amount of training data \(D\), the lowest value of the limiting behavior is controlled by the data variability represented by \(v\) and \(K\). Due to the symmetric nature of the task setting, neither the limiting behavior (43) nor the original expression (40) explicitly show the dependence on the relative number of training samples in the two respective classes \(c_{1,2}\). This is due to the fact that the task setup in (38) is symmetric. In the case of asymmetric statistics this behavior changes. Moreover, the difference between variance \(a\) and covariance \(b\) enters the expression in a non-trivial manner

Using those results, we will investigate the implications for linear, non-linear, and deep kernels using the artificial data set, Section II.3, as well as real-world data.
### Applications to linear, non-linear and deep non-linear NNGP kernels

#### iii.2.1 Linear Kernel

Before going to the non-linear case, let us investigate the implications of (40) and (43) for a simple one-layer linear network. We assume that our network consists of a read-in weight \(\mathbf{V}\in\mathbb{R}^{1\times N_{\mathrm{dim}}};\mathbf{V}_{i}\sim\mathcal{N }\left(0,\sigma_{v}^{2}/N_{\mathrm{dim}}\right)\), which maps the \(N_{\mathrm{dim}}\) dimensional input vector to a one-dimensional output space. Including a regularization noise, the output hence reads

\[y_{\alpha}=\mathbf{V}x_{\alpha}+\xi_{\alpha}\,. \tag{44}\]

In this particular case the read-in, read-out, and hidden weights in the general setup (1) coincide with each other. Computing the average with respect to the weights \(\mathbf{V}\) yields the kernel

\[K_{\alpha\beta}^{y}=\left\langle y_{\alpha}y_{\beta}\right\rangle_{\mathbf{V}} =K_{\alpha\beta}^{0}+\delta_{\alpha\beta}\,\sigma_{\mathrm{reg}}^{2}\,, \tag{45}\]

where \(K_{\alpha\beta}^{0}\) is given by (10); it is hence a rescaled version of the overlap of the input vectors and the variance of the regularization noise.

We now assume that the matrix elements of the input-data overlap (45) are distributed according to a multivariate Gaussian (20).

As the mean and the covariance of the entries \(K_{\alpha\beta}^{y}\) are given by the statistics (15) we evaluate (40) and (43) with

\[a^{(\mathrm{Lin})} =\sigma_{v}^{2}+\sigma_{\mathrm{reg}}^{2}\,,\] \[b^{(\mathrm{Lin})} =\sigma_{v}^{2}\,u\,,\] \[K^{(\mathrm{Lin})} =\frac{\sigma_{v}^{4}}{N_{\mathrm{dim}}}\left(1-u^{2}\right),\] \[v^{(\mathrm{Lin})} =\frac{\sigma_{v}^{4}}{N_{\mathrm{dim}}}\,u\left(1-u\right),\] \[u :=4p(p-1)+1\,. \tag{46}\]

The asymptotic result for the first order correction, assuming that \(\sigma_{v}^{2}\neq 0\), can hence be evaluated, assuming \(p\neq 0.5\), as

\[\left\langle y_{*}\right\rangle_{1}^{(\infty)}= \frac{y_{1}\sigma_{v}^{2}\frac{\left(1-u\right)}{N_{\mathrm{dim} }}}{\left(\sigma_{v}^{2}\left(1-u\right)+\sigma_{\mathrm{reg}}^{2}\right)u} \left(-2u-\frac{\sigma_{\mathrm{reg}}^{2}}{\sigma_{v}^{2}}\right)\,. \tag{47}\]

Using this explicit form of \(\left\langle y_{*}\right\rangle_{1}^{(\infty)}\) one can see

* as \(u\in[0,1]\) the corrections are always negative and hence provide a less optimistic estimate for the generalization compared to the zeroth-order approximation;
* in the limit \(\sigma_{v}^{2}\rightarrow\infty\) the regularizer in (47) becomes irrelevant and the matrix inversion becomes unstable.
* taking \(\sigma_{v}^{2}\to 0\) yields a setting where constructing the limiting formula (47) is not useful, as all relevant quantities (40) like \(g,v,K\to 0\) vanish; hence the inference yields zero which is consistent with our intuition: \(\sigma_{v}^{2}\to 0\) implies that only the regularizer decides, which is unbiased with regards to the class membership of the data. Hence the kernel cannot make any prediction which is substantially informed by the data.

Figure (3) shows that the zeroth-order approximation\(\left\langle y^{*}\right\rangle_{0}\), even though it is able to capture some dependence on the amount of training data, is indeed too optimistic and predicts a mean-inferred network output closer to its negative target value \(y=-1\) than numerically obtained. The first-order correction on the other hand is able to reliably predict the results. Furthermore the limiting results \(D\rightarrow\infty\) match the numerical results for different task settings \(p\). These limiting results are consistently higher than the zeroth-order approximation \(\left\langle y_{*}\right\rangle_{0}\) and depend on the level of data variability. Deviations of the empirical results from the theory in the case \(p=0.6\) could stem from the fact that for \(p=0.5\) the fluctuations are maximal and our theory assumes small fluctuations.

#### iii.2.2 Non-Linear Kernel

We will now investigate how the non-linearities \(\phi\) present in typical network architectures (1) influence our results for the learning curve (40) and (43).

As the ansatz in Section III.1 does not make any assumption, apart from Gaussianity, on the overlap-matrix \(K^{y}\), the results presented in Section III.1.5 are general. One can use the knowledge of the statistics of the overlap matrix in the read-in layer \(K^{0}\) in (15) to extend the result (40) to both non-linear and deep feed-forward neural networks.

As in Section III.2.1 we start with the assumption that the input kernel matrix is distributed according to a multivariate Gaussian: \(K_{\alpha\beta}^{0}\sim\mathcal{N}(\mu_{\alpha\beta},\Sigma_{(\alpha\beta)( \gamma\delta)})\). In the non-linear case, we consider a read-in layer \(\mathbf{V}\in\mathbb{R}^{N_{h}\times N_{\mathrm{dim}}};\mathbf{V}_{i,j}\sim \mathcal{N}\left(0,\sigma_{v}^{2}/N_{\mathrm{dim}}\right)\), which maps the inputs to the hidden-state space and a separate read-out layer \(\mathbf{W}\in\mathbb{R}^{1\times N_{h}};\mathbf{W}_{i}\sim\mathcal{N}\left(0, \sigma_{w}^{2}/N_{h}\right)\), obtaining a neural network with a single hidden layer

\[h_{\alpha}^{(0)} =\mathbf{V}x_{\alpha}\,,\] \[y_{\alpha} =\mathbf{W}\phi\left(h_{\alpha}^{(0)}\right)+\xi_{\alpha}\,, \tag{48}\]

and network kernel
\[\left\langle y_{\alpha}y_{\beta}\right\rangle_{\mathbf{V},\mathbf{W}} =\frac{\sigma_{w}^{2}}{N_{h}}\sum_{i=1}^{N_{h}}\left\langle\phi \left(h_{\alpha i}^{(0)}\right)\phi\left(h_{\beta i}^{(0)}\right)\right\rangle_{ \mathbf{V}}\] \[+\delta_{\alpha\beta}\,\sigma_{\text{reg}}^{2}\,. \tag{49}\]

As we consider the limit \(N_{h}\to\infty\), one can replace the empirical average \(\frac{1}{N_{h}}\sum_{i=1}^{N_{h}}...\) with a distributional average \(\frac{1}{N_{h}}\sum_{i=1}^{N_{h}}...\to\left\langle...\right\rangle_{\mathbf{h} ^{(0)}}\)[18, 30]. This yields the following result for the kernel matrix \(K_{\alpha\beta}^{y}\) of the multivariate Gaussian

\[K_{\alpha\beta}^{y}\underset{N_{h}\to\infty}{\rightarrow}\sigma_ {w}^{2}\Big{\langle}\phi\left(h_{\alpha}^{(0)}\right)\phi\left(h_{\beta}^{(0)} \right)\Big{\rangle}_{\mathbf{h}^{(0)},\mathbf{V}}+\delta_{\alpha\beta}\, \sigma_{\text{reg}}^{2}\,. \tag{50}\]

The expectation over the hidden states \(h_{\alpha}^{(0)},h_{\beta}^{(0)}\) is with regard to the Gaussian

\[\left(\begin{array}{c}h_{\alpha}^{(0)}\\ h_{\beta}^{(0)}\end{array}\right)\sim\mathcal{N}\left(\left(\begin{array}{c }0\\ 0\end{array}\right),\left(\begin{array}{cc}K_{\alpha\alpha}^{0}&K_{\alpha \beta}^{0}\\ K_{\beta\alpha}^{0}&K_{\beta\beta}^{0}\end{array}\right)\right)\,, \tag{51}\]

with the variance \(K_{\alpha\alpha}^{0}\) and the covariance \(K_{\alpha\beta}^{0}\) given by (10). Evaluating the Gaussian integrals in (50) is analytically possible in certain limiting cases [3, 40]. For an erf-activation function, as a prototype of a saturating activation function, this average yields

\[\left\langle\phi^{2}\left(h_{\alpha}^{(0)}\right)\right\rangle_{ \mathbf{h}^{(0)}}= \frac{4}{\pi}\arctan\left(\sqrt{1+4K_{\alpha\alpha}^{0}}\right)-1\,, \tag{52}\] \[\left\langle\phi\left(h_{\alpha}^{(0)}\right)\phi\left(h_{\beta}^ {(0)}\right)\right\rangle_{\mathbf{h}^{(0)}}= \frac{2}{\pi}\arcsin\left(\frac{2K_{\alpha\beta}^{0}}{1+2K_{ \alpha\alpha}^{0}}\right)\,. \tag{53}\]

We use that the input kernel matrix \(K^{0}\) is distributed as \(K_{\alpha\beta}^{0}\sim\mathcal{N}(\mu_{\alpha\beta},\Sigma_{(\alpha\beta)( \gamma\delta)})\). Equation (50) hence provides information on how the mean overlap \(m_{\alpha\beta}\) changes due to the application of the non-linearity \(\phi(\cdot)\), fixing the parameters \(a\), \(b\), \(K\), \(v\) of the general form (38) as

\[a^{\text{(Non-lin)}} =K_{\alpha\alpha}^{y}=\sigma_{w}^{2}\Big{\langle}\phi^{2}\left(h_ {\alpha}^{(0)}\right)\Big{\rangle}_{\mathbf{h}^{(0)}}+\sigma_{\text{reg}}^{2}\,, \tag{54}\] \[b^{\text{(Non-lin)}} =K_{\alpha\beta}^{y}=\sigma_{w}^{2}\Big{\langle}\phi\left(h_{ \alpha}^{(0)}\right)\phi\left(h_{\beta}^{(0)}\right)\Big{\rangle}_{\mathbf{h} ^{(0)}}\,\,. \tag{55}\]

where the averages over \(h^{(0)}\) are evaluated with regard to the Gaussian (51) for \(\phi(x)=\text{erf}(x)\). We further require in 55 that \(\alpha\neq\beta,c(\alpha)=c(\beta)\).

To evaluate the corrections in (40), we also need to understand how the presence of the non-linearity \(\phi(x)\) shapes the parameters \(K,v\) that control the variability. Under the assumption of small covariance \(\Sigma_{(\alpha\beta)(\gamma\delta)}\) one can use (53) to compute \(C_{(\alpha,\beta)(\gamma,\delta)}\) using linear response theory. As \(K_{\alpha\beta}^{0}\) is stochastic and provided by (20), we decompose \(K_{\alpha\beta}^{0}\) into a deterministic kernel \(\mu_{\alpha\beta}\) and a stochastic perturbation \(\eta_{\alpha\beta}\sim\mathcal{N}\left(0,\Sigma_{(\alpha\beta)(\gamma\delta)}\right)\). Linearizing (55) around \(\mu_{\alpha\beta}\) via Price's theorem [31], the stochasticity in the read-out layer yields

\[C_{(\alpha\beta)(\gamma\delta)} =\sigma_{w}^{4}\,K_{\alpha\beta}^{(\phi^{\prime})}\,K_{\gamma \delta}^{(\phi^{\prime})}\,\Sigma_{(\alpha\beta)(\gamma\delta)}\,, \tag{56}\] \[K_{\alpha\beta}^{(\phi^{\prime})} =\left\langle\phi^{\prime}\left(h_{\alpha}^{(0)}\right)\phi^{ \prime}\left(h_{\beta}^{(0)}\right)\right\rangle, \tag{57}\]

where \(h^{(0)}\) is distributed as in (51). This clearly shows that the variability simply transforms with a prefactor

\[K^{\text{(Non-lin)}} =\sigma_{w}^{4}\,K_{\alpha\beta}^{(\phi^{\prime})}\,K_{\alpha \beta}^{(\phi^{\prime})}\,\kappa\,,\] \[v^{\text{(Non-lin)}} =\sigma_{w}^{4}K_{\alpha\beta}^{(\phi^{\prime})}\,K_{\alpha \delta}^{(\phi^{\prime})}\,v\,, \tag{58}\]

with \(\kappa,\nu\) defined as in (15). Evaluating the integral in \(\left\langle\phi^{\prime}\left(h_{\alpha}^{(0)}\right)\phi^{\prime}\left(h_{ \beta}^{(0)}\right)\right\rangle\) is hard in general. In fact, the integral which occurs is equivalent to the one in [24] for the Lyapunov exponent and, equivalently, in [34, 30] for the susceptibility in the propagation of information in deep feed-forward neural networks. This is consistent with
the assumption that our treatment of the non-linearity follows a linear response approach as in [24]. For the erf-activation we can evaluate the kernel \(K^{(\phi^{\prime})}_{\alpha\beta}\) as

\[K^{(\phi^{\prime})}_{\alpha\beta} =\frac{4}{\pi\left(1+2a^{(0)}\right)}\left(1-\left(\frac{2b^{(0)}} {1+2a^{(0)}}\right)^{2}\right)^{-\frac{1}{2}}\,, \tag{59}\] \[a^{(0)} =\sigma_{v}^{2}\quad,\quad b^{(0)}=\sigma_{v}^{2}u\,,\] (60) \[u =4p(p-1)+1\,, \tag{61}\]

which allows us to evaluate (58). Already in the one hidden-layer setting we can see that the behavior is qualitatively different from a linear setting: \(K^{(\mathrm{Non-lin})}\) and \(v^{(\mathrm{Non-lin})}\) scale with a linear factor which now also involves the parameter \(\sigma_{v}^{2}\) in a non-linear manner.

#### iii.2.3 Multilayer-Kernel

So far we considered single-layer networks. However, in practice the application of multi-layer networks is often necessary. One can straightforwardly extend the results from the non-linear case (III.2.2) to the deep non-linear case. We consider the architecture introduced in (1) in Section II.1.1 where the variable \(L\) denotes the number of hidden layers, and \(1\leq l\leq L\) is the layer index. Similar to the computations in Section III.2.2 one can derive a set of relations to obtain \(K^{y}_{\alpha\beta}\)

\[K^{0}_{\alpha\beta} =\frac{\sigma_{v}^{2}}{N_{\mathrm{dim}}}\;K^{x}_{\alpha\beta}\,,\] \[K^{(\phi)l}_{\alpha\beta} =\sigma_{w}^{2}\;\left\langle\phi\left(h_{\alpha}^{(l-1)}\right) \phi\left(h_{\beta}^{(l-1)}\right)\right\rangle\,,\] \[K^{y}_{\alpha\beta} =\sigma_{u}^{2}\;\left\langle\phi\left(h_{\alpha}^{(L)}\right) \phi\left(h_{\beta}^{(L)}\right)\right\rangle+\delta_{\alpha\beta}\,\sigma_{ \mathrm{reg}}^{2}\,. \tag{62}\]

As [30; 34; 42] showed for feed-forward networks, deep non-linear networks strongly alter both the variance and the covariance. So we expect them to influence the generalization properties. In order to understand how the fluctuations \(C_{(\alpha\beta)(\gamma\delta)}\) transform through propagation, one can employ the chain rule to linearize (62) and obtain

\[C^{yy}_{(\alpha,\beta)(\gamma,\delta)}=\sigma_{u}^{4}\,\prod_{l=1}^{L}\left[K ^{(\phi^{\prime})l}_{\alpha\beta}K^{(\phi^{\prime})l}_{\gamma\delta}\right]\, \Sigma_{(\alpha\beta)(\gamma\delta)}\,. \tag{63}\]

A systematic derivation of this result as the leading order fluctuation correction in \(N_{h}^{-1}\) is found in the appendix of [35].

Equation (62) and (63) show that the kernel performance will depend on the non-linearity \(\phi\), the variances \(\sigma_{v}^{2}\), \(\sigma_{w}^{2}\), \(\sigma_{u}^{2}\), and the network depth \(L\).

Figure 4 (a) shows the comparison of the mean inferred network output \(\langle y^{*}\rangle\) for the true test label \(y=-1\) between empirical results and the first order corrections. The regime (\(\sigma_{w}^{2}<1\)) in which the kernel vanishes, leads to a poor performance. The marginal regime (\(\sigma_{w}^{2}\simeq 1\)) provides a better choice for the overall network performance. Equation (4) (b) shows that the maximum absolute value for the predictive mean is achieved slightly in the supercritical regime \(\sigma_{w}^{2}>1\). With larger number of layers, the optimum becomes more and more pronounced and approaches the critical value \(\sigma_{w}^{2}=1\) from above. The optimum for the predictive mean to occur slightly in the supercritical regime may be surprising with regard to the expectation that network trainability peaks precisely at \(\sigma_{w}^{2}=1\)[30]. In particular at shallow depths, the optimum becomes very wide and shifts to \(\sigma_{w}>1\). For few layers, even at \(\sigma_{w}^{2}>1\) the increase of variance \(K^{y}_{\alpha\alpha}\) per layer remains moderate and stays within the dynamical range of the activation function. Thus differences in covariance are faithfully transmitted by the kernel and hence allow for accurate predictions. The theory including corrections to linear order throughout matches the empirical results and hence provides good estimates for choosing the kernel architecture.

Figure 4: **Predictive mean in a deep non-linear feed forward network with heterogeneous kernel.****(a)** Comparison of mean inferred network output for non-linear network with \(\phi(x)=\mathrm{erf(x)}\), five layers for different values of the gain \(\sigma_{w}\). The figure displays numerical results (bars), zeroth-order approximation (dashed) and first-order corrections (solid). **(b)** Similar comparison as in (a) for different network depths \(L=5,10,20,50\). In all settings we used \(N_{\mathrm{dim}}=50\) for \(D=100\), \(p=0.8\), \(\sigma_{v}^{2}=1\), \(\sigma_{\mathrm{reg}}^{2}=1\). Empirical results display mean and standard deviation over 1000 trials with 1000 test points per trial.
#### iii.3.4 Experiments on Non-Symmetric Task Settings and MNIST

In contrast to the symmetric setting in the previous subsections, real data-sets such as MNIST exhibit asymmetric statistics so that the different blocks in \(m_{\alpha\beta}\) and \(C_{(\alpha\beta)(\gamma\delta)}\) assume different values in general. All theoretical results from Section III.1 still hold. However, as the tensor elements of \(m_{\alpha,\beta}\) and \(C_{(\alpha,\beta)(\gamma,\delta)}\) change, one needs to reconsider the evaluation in Section III.1.5 in the most general form which yields a more general version of the result.

Finite MNIST datasetFirst we consider a setting, where we work with the pure MNIST dataset for two distinct labels \(0\) and \(4\). In this setting we estimate the class-dependent tensor elements \(m_{\alpha\beta}\) and \(C_{(\alpha\beta)(\gamma,\delta)}\) directly from the data. We define the data-set size per class, from which we sample the theory as \(D_{\text{base}}\). The training points are also drawn from a subset of these \(D_{\text{base}}\) data points. To compare the analytical learning curve for \(\langle y_{\ast}\rangle\) at \(D\) training data-points to the empirical results, we need to draw multiple samples of training datasets of size \(D<D_{\text{base}}\). As the amount of data in MNIST is limited, these samples will generally not be independent and therefore violate our assumption. Nevertheless we can see in Figure 5 that if \(D\) is sufficiently small compared to \(D_{\text{base}}\), the empirical results and theoretical results match well.

Gaussianized MNIST datasetTo test whether deviations in Figure 5 at large \(D\) stem from correlations in the samples of the dataset we construct a generative scheme for MNIST data. This allows for the generation of infinitely many training points and hence the assumption that the training data is i.i.d. is fulfilled. We construct a pixel-wise Gaussian distribution for MNIST images from the samples. We use this model to sample as many MNIST images as necessary for the evaluation of the empirical learning curves. Based on the class-dependent statistics for the pixel means and the pixel covariances in the input data one can directly compute the elements of the mean \(\mu_{\alpha\beta}\) and the covariance \(\Sigma_{(\alpha\beta)(\gamma\delta)}\) for the distribution of the input kernel matrix \(K^{0}_{\alpha\beta}\). We see in Figure 6 that our theory describes the results well for this data-set also for large numbers of training samples.

Furthermore we can see that in the case of an asymmetric data-set the learning curves depend on the balance ratio of training data \(\rho=D_{c_{1}}/D_{c_{2}}\). The bias towards class one in Figure 6 b) is evident from the curves with \(\rho>0.5\) predicting a lower mean inferred network output, closer to the target label \(y=-1\) of class \(1\).

## IV Discussion

In this work we investigate the influence of data variability on the performance of Bayesian inference. The probabilistic nature of the data manifests itself in a heterogeneity of the entries in the block-structured kernel matrix of the corresponding Gaussian process. We show

Figure 5: **Predictive mean for a linear network with MNIST data**: Comparison of mean inferred network output for a linear network with 1 layer for different training set sizes \(D\). The figure displays numerical results (bars), zeroth-order prediction (dashed) and first-order corrections (solid). Settings \(N_{\text{dim}}=784\), \(\sigma^{2}_{\text{reg}}=2\), \(D_{\text{base}}=4000\). MNIST classes \(c_{1}=0\), \(c_{2}=4\), \(y_{c_{1}}=-1\), \(y_{c_{2}}=1\); balanced data-set in \(D_{\text{base}}\) and at each \(D\). Empirical results display mean and standard deviation over 1000 trials with 1000 test points per trial.

Figure 6: **Predictive mean for a erf-network with Gaussianized MNIST data**: **a)** Mean inferred network output for MNIST classification with \(\phi(x)=\text{erf}(x)\). Figure shows zeroth-order (dashed line), first-order (solid line), and empirical results (bars). **b)** Mean inferred network output in first order approximation (solid lines) and empirical results (bars) for MNIST classification with different ratios \(\rho=D_{c_{1}}/D_{c_{2}}\) between numbers of training samples per class \(D_{c_{1}}\) and \(D_{c_{2}}\), respectively; \(\rho=0.5\) (yellow), \(\rho=0.6\) (blue), \(\rho=0.7\) (red). Empirical results display mean and standard deviation over 50 trials with 1000 test points per trial.
that this heterogeneity for a sufficiently large number of \(D\) of data samples can be treated as an effective non-Gaussian theory. By employing a time-dependent formulation for the mean of the predictive distribution, this heterogeneity can be treated as a disorder average that circumvents the use of the replica trick. A perturbative treatment of the variability yields first-order corrections to the mean in the variance of the heterogeneity that always push the mean of the predictive distribution towards zero. In particular, we obtain limiting expressions that accurately describe the mean in the limit of infinite training data, qualitatively correcting the zeroth-order approximation corresponding to homogeneous kernel matrices, is overconfident in predicting the mean to perfectly match the training data in this limit. This finding shows how variability fundamentally limits predictive performance and provides not only a quantitative but also a qualitative difference. Moreover at a finite number of training data the theory explains the empirically observed performance accurately. We show that our framework captures predictions in linear, non-linear shallow and deep networks. In non-linear networks, we show that the optimal value for the variance of the prior weight distribution is achieved in the super-critical regime. The optimal range for this parameter is broad in shallow networks and becomes progressively more narrow in deep networks. These findings support that the optimal initialization is not at the critical point where the variance is unity, as previously thought [30], but that super-critical initialization may have an advantage when considering input variability. An artificial dataset illustrates the origin and the typical statistical structure that arises in heterogeneous kernels, while the application of the formalism to MNIST [17] demonstrates potential use to predict the expected performance in real world applications.

The field theoretical formalism can be combined with approaches that study the effect of fluctuations due to the finite width of the layers [28; 43; 45; 35]. In fact, in the large \(N_{h}\) limit the NNGP kernel is inert to the training data, the so called lazy regime. At finite network width, the kernel itself receives corrections which are commonly associated with the adaptation of the network to the training data, thus representing what is known as feature learning. The interplay of heterogeneity of the kernel with such finite-size adaptations is a fruitful future direction.

Another approach to study learning in the limit of large width is offered by the neural tangent kernel (NTK) [15], which considers the effect of gradient descent on the network output up to linear order in the change of the weights. A combination of the approach presented here with the NTK instead of the NNGP kernel seems possible and would provide insights into how data heterogeneity affects training dynamics.

The analytical results presented here are based on the assumption that the variability of the data is small and can hence be treated perturbatively. In the regime of large data variability, it is conceivable to employ self-consistent methods instead, which would technically correspond to the computation of saddle points of certain order parameters, which typically leads to an infinite resummation of the perturbative terms that dominate in the large \(N_{h}\) limit. Such approaches may be useful to study and predict the performance of kernel methods for data that show little or no linear separability and are thus dominated by variability. Another direction of extension is the computation of the variance of the Bayesian predictor, which in principle can be treated with the same set of methods as presented here. Finally, since the large width limit as well as finite-size corrections, which in particular yield the kernel response function that we employed here, can be obtained for recurrent and deep networks in the same formalism [35] as well as for residual networks (ResNets) [9], the theory of generalization presented here can straight forwardly be extended to recurrent networks and to ResNets.

###### Acknowledgements.

 We thank Claudia Merger, Bastian Epping, Kai Segadlo and Alexander van Meegen for helpful discussions. This work was partly supported by the German Federal Ministry for Education and Research (BMBF Grant 01IS19077A to Julich and BMBF Grant 01IS19077B to Aachen) and funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 368482240/GRK2416, the Excellence Initiative of the German federal and state governments (ERS PFJARA-SDS005), and the Helmholtz Association Initiative and Networking Fund under project number SO-092 (Advanced Computing Architectures, ACA). Open access publication funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 491111487.
