Bayesian inference with finitely wide neural networks
The analytic inference, e.g. predictive distribution being in closed form,
may be an appealing benefit for machine learning practitioners when they treat
wide neural networks as Gaussian process in Bayesian setting. The realistic
widths, however, are finite and cause weak deviation from the Gaussianity under
which partial marginalization of random variables in a model is
straightforward. On the basis of multivariate Edgeworth expansion, we propose a
non-Gaussian distribution in differential form to model a finite set of outputs
from a random neural network, and derive the corresponding marginal and
conditional properties. Thus, we are able to derive the non-Gaussian posterior
distribution in Bayesian regression task. In addition, in the bottlenecked deep
neural networks, a weight space representation of deep Gaussian process, the
non-Gaussianity is investigated through the marginal kernel.
Chi-Ken Lu
# Bayesian inference with finitely wide neural networks

###### Abstract

The analytic inference, e.g. predictive distribution being in closed form, may be an appealing benefit for machine learning practitioners when they treat wide neural networks as Gaussian process in Bayesian setting. The realistic widths, however, are finite and cause weak deviation from the Gaussianity under which partial marginalization of random variables in a model is straightforward. On the basis of multivariate Edgeworth expansion, we propose a non-Gaussian distribution in differential form to model a finite set of outputs from a random neural network, and derive the corresponding marginal and conditional properties. Thus, we are able to derive the non-Gaussian posterior distribution in Bayesian regression task. In addition, in the bottlenecked deep neural networks, a weight space representation of deep Gaussian process, the non-Gaussianity is investigated through the marginal kernel and the accompanying small parameters.

## I Introduction

Neal in his seminal work [1] pointed out that a shallow but infinitely wide random neural network is a Gaussian process (GP) [2] in statistical sense. Subsequent work [3; 4] in interpreting neural network with specific nonlinear activation units as kernel machines was also inspired by such idea. More recent reports [5; 6] further claimed the equivalence between GP and deep neural networks when each hidden layer in latter is of infinite width. Consequently, machine learning practitioners can perform Bayesian inference by treating deep and wide neural network as a GP, and exploit the analytic marginal and conditional properties of multivariate Gaussian distribution. Otherwise, one needs to employ gradient-based learning and bootstrap sampling for obtaining predictive distribution [7].

In reality, all neural networks have finite width. Therefore, the deviation from Gaussianity requires further quantitative account as practitioners may wonder the corrections to the predictive mean and variance in, for example, a regression task. Yaida [8] and colleagues [9] proposed a perturbative approach for computing the multivariate cumulants by direct application of Wick's contraction theorem. Moreover, the fourth cumulants are shown to be nonzero, scaled by the sum of reciprocal widths, \(1/N_{1}+1/N_{2}+\cdots+1/N_{L-1}\), and signaling non-Gaussian aspect of the random processes representing finite-width deep neural networks with \(L\) hidden layers [8]. A quartic energy functional for a fixed set of network outputs is formulated field-theoretically with which the corrections to the posterior mean and variance due to the _weak_ non-Gaussianity are obtained [8].

While Yaida's approach is appealing from a field theory perspective (see also [10; 11]), the loss of elegant marginal property due to the presence of fourth power term in exponent is critical for analytic inference with finite-width networks. An alternative thinking is to modify the multivariate Gaussian distribution so that the new distribution can match the higher cumulants associated with the networks. In this paper, we shall use the multivariate Edgeworth series expansion [12] to construct the non-Gaussian distribution for the network's outputs. In particular, we find that the differential representation of Edgeworth expansion greatly facilitates the derivation of marginal and conditional properties of the non-Gaussian distribution. Three main results are reported in this paper. First, the marginal property is intact, and the corrections to conditional mean and variance are derived. Second, with observed data, the non-Gaussian posterior distribution associated with an unobserved output is derived. Third, we derive the marginal covariance [13] of a bottlenecked deep neural network [14], which represents deep Gaussian Process [15] in weight space. It is worthwhile to note that some of the hidden layers in the bottlenecked network are narrow and _strong_ non-Gaussianity may be induced.

The paper has the following organization. In Sec. II we begin by reviewing the computational structure of a shallow cosine network, its equivalence to GP with a Gaussian kernel, and the emergence of non-Gaussianity due to finite width. The shallow network with random parameters in Bayesian setting is a non-Gaussian prior over function. In Sec. III the non-Gaussian prior is represented as a differential representation of Edgeworth expansion around Gaussian multivariate distribution, and its marginal and conditional properties are established. Application of the non-Gaussian prior in Bayesian regression task is discussed in Sec. IV. Finally, Sec. V is devoted to investigating the combined effect
of nonlinear activation, depth and finite width on the non-Gaussian prior for a deep bottlenecked cosine network, which is followed by a discussion in Sec. VI.

## II Wide feed forward neural network

Let us start with the discussion of a random shallow feed forward network with cosine activation. In the infinite width limit, the network is statistically equivalent to a GP with Gaussian kernel [16]. The goal is to see the emergence of nonzero fourth cumulant when the width is finite. Consider a single output network with \(N\) activation units, the real function value \(z_{\alpha}\in\mathbb{R}\) with greek subscript is an indexed random variable in association with its input \(\mathbf{x}_{\alpha}\in\mathbb{R}^{d}\). Explicitly, the output-input of the network has the following relation,

\[z_{\alpha}=\sqrt{\frac{2}{N}}\sum_{i=1}^{N}w_{i}\cos(\frac{\Omega_{i}\cdot \mathbf{x}_{\alpha}}{\sqrt{d}}+\phi_{i})\,, \tag{1}\]

where the weight variables \(w\)'s are sampled from the Gaussian distribution \(\mathcal{N}(0,1)\), the scaling variables \(\Omega\)'s sampled from \(\mathcal{N}(0,I_{d})\), and the phase variables \(\phi\)'s sampled from uniform distribution \(\mathcal{U}([0,2\pi])\). The normalization factors \(\sqrt{2/N}\) and \(1/\sqrt{d}\) above follow the parameterization used in [17]. Because \(w\) is zero-mean, the first relevant statistical moment is the covariance, \(k(\mathbf{x}_{\alpha},\mathbf{x}_{\beta})\), namely the expectation of product of two function values at two inputs,

\[\begin{split} k(\mathbf{x}_{\alpha},\mathbf{x}_{\beta})& :=\mathbb{E}[z_{\alpha}z_{\beta}]\\ &=\frac{1}{N}\sum_{i}\mathbb{E}\bigg{\{}\cos[\frac{\Omega_{i} \cdot(\mathbf{x}_{\alpha}-\mathbf{x}_{\beta})}{\sqrt{d}}]+\cos[\frac{\Omega_{ i}\cdot(\mathbf{x}_{\alpha}+\mathbf{x}_{\beta})}{\sqrt{d}}+2\phi_{i}]\bigg{\}}\\ &=e^{-\frac{1}{2d}|\mathbf{x}_{\alpha}-\mathbf{x}_{\beta}|^{2}} \,.\end{split} \tag{2}\]

In above, the independence between the random variables \(w\)'s is used to arrive at the second equality. The third equality is due to the fact that Fourier transformation of Gaussian is Gaussian and the vanishing average of cosine with uniformly random phase. The fourth moment can be computed in similar manners, but one shall notice that the product of two cosine terms containing the random phase can generate a nonzero contribution. Here, we focus on the fourth cumulant tensor as the signature of non-Gaussianity,

\[\begin{split} V_{\alpha\beta\gamma\delta}&:= \mathbb{E}[z_{\alpha}z_{\beta}z_{\gamma}z_{\delta}]-K_{\alpha\beta}K_{\gamma \delta}-K_{\alpha\gamma}K_{\beta\delta}-K_{\alpha\delta}K_{\beta\gamma}\\ &=\frac{1}{N^{2}}\sum_{ij}\mathbb{E}\big{\{}\cos[\frac{\Omega_{i} \cdot(\mathbf{x}_{\alpha}+\mathbf{x}_{\beta})}{\sqrt{d}}+2\phi_{i}]\cos[\frac {\Omega_{j}\cdot(\mathbf{x}_{\gamma}+\mathbf{x}_{\delta})}{\sqrt{d}}+2\phi_{j }]\big{\}}+\text{sym. perm.}\\ &=\frac{1}{2N}\bigg{(}e^{-\frac{1}{2d}|\mathbf{x}_{\alpha}+ \mathbf{x}_{\beta}-\mathbf{x}_{\gamma}-\mathbf{x}_{\delta}|^{2}}+e^{-\frac{1 }{2d}|\mathbf{x}_{\alpha}+\mathbf{x}_{\gamma}-\mathbf{x}_{\beta}-\mathbf{x}_ {\delta}|^{2}}+e^{-\frac{1}{2d}|\mathbf{x}_{\alpha}+\mathbf{x}_{\delta}- \mathbf{x}_{\beta}-\mathbf{x}_{\gamma}|^{2}}\bigg{)}\,.\end{split} \tag{3}\]

Hence, the fourth cumulant tensor receives a contribution proportional to the reciprocal width, \(1/N\), and it is symmetric with respect to permutation of indices.

## III Multivariate Edgeworth expansion and non-Gaussian prior

The fourth cumulant tensor in Eq. (3) signifies that the distribution over the network outputs is non-Gaussian. In other words, the random parameters \(w\)'s, \(\Omega\)'s, and \(\phi\)'s from a prior distribution induce the non-Gaussian prior function distribution due to the finite width. The first interesting consequence of the nonzero fourth cumulant is that the prior distribution over a single output, say \(z_{\alpha}\), receives a correction term deviating from Gaussian [18],

\[q(z_{\alpha})=\mathcal{N}(z_{\alpha}|0,\sigma_{0}^{2})\bigg{[}1+\frac{V_{ \alpha\alpha\alpha\alpha}}{24}(3-6\frac{z_{\alpha}^{2}}{\sigma_{0}^{2}}+\frac {z_{\alpha}^{4}}{\sigma_{0}^{4}})\bigg{]}\,. \tag{4}\]

The variance \(\sigma_{0}^{2}:=k(z_{\alpha},z_{\alpha})\) and the new distribution \(q\) has matched moments, i.e. \(\mathbb{E}_{q}[1]=1\), \(\mathbb{E}_{q}[z_{\alpha}^{2}]=\sigma_{0}^{2}\), and \(\mathbb{E}_{q}[z_{\alpha}^{4}]=3\sigma_{0}^{4}+V_{\alpha\alpha\alpha\alpha}\). Although the work [18] derived the above non-Gaussian distribution from a renormalization group perspective, such expansion with Hermite polynomials around Gaussian distribution has been known as Edgeworth expansion in statistics literature [19; 20]. When the width \(N\rightarrow\infty\), the fourth cumulant \(V\) vanishes and for notational convenience we denote the limiting distribution by \(q_{\infty}\), which is Gaussian.
### Joint distribution

As we are mainly interested in the prior distribution over a finite set of function values \(\{z_{1},z_{2},\cdots\}\), the objective here is to construct the non-Gaussian joint distribution as a prior in Bayesian learning. The work [12] suggested the construction of multivariate Edgeworth expansion by replacing the univariate Gaussian \(\mathcal{N}(z|0,\sigma_{0}^{2})\) with multivariate one \(\mathcal{N}(\mathbf{z}|0,K)\) and the Hermite polynomials with contracted tensor terms containing the inverse covariance matrix \([K^{-1}]\) and vector of function values \(\mathbf{z}\). Please see Appendix A for the expression.

An apparent difficulty with such representation is the demonstration of consistency after marginalization, i.e., \(\int dz_{1}q(z_{1},z_{2},z_{3},\cdots)=q(z_{2},z_{3},\cdots)\), which is of critical importance for deriving conditional distribution and subsequent Bayesian learning. Inspired by the fact that the Hermite polynomials are obtained by taking derivative of Gaussian, we can rewrite the multivariate Edgeworth expansion in the following differential form,

\[q(\mathbf{z})=\big{(}1+\frac{V_{ijml}}{24}\partial_{z_{i}}\partial_{z_{j}} \partial_{z_{m}}\partial_{z_{l}}\big{)}\mathcal{N}(\mathbf{z}|0,K) \tag{5}\]

where the notation of summation by repeated indices is employed. An advantage of such representation is that the adjoint property of derivative, \(\partial_{z}^{\dagger}=-\partial_{z}\), can be exploited in the integration of functions which vanish at infinity. Hence, it is easy to see the distribution \(q\) is normalized, and the second and fourth moments match, for instance,

\[\begin{split}\mathbb{E}_{q}[z_{1}z_{2}]&=K_{12}\,, \\ \mathbb{E}_{q}[z_{1}z_{2}z_{3}z_{4}]&=K_{12}K_{34}+K _{13}K_{24}+K_{14}K_{23}+V_{1234}\;.\end{split} \tag{6}\]

The above non-Gaussian prior can be viewed as a perturbative extension of multivariate Gaussian. In the setting of Bayesian regression, the marginal and conditional properties of Gaussian are appealing when conjugate likelihood functions are employed. In the following, we shall examine the effects of the perturbation in differential form on these properties.

### Marginal consistency and conditional statistics

We first show that the marginal consistency is intact with the Edgeworth expansion in differential form. Namely, without loss of generality, it suffices to show that marginalization of \(z_{1}\) of the joint distribution \(q(z_{1},z_{2},z_{3},\cdots)\),

\[\begin{split}\int dz_{1}q(z_{1},\tilde{\mathbf{z}})& =\mathcal{N}(\tilde{\mathbf{z}}|0,\tilde{K})+v_{\widetilde{ij}\widetilde{j} \widetilde{m}\widetilde{l}}\partial_{z_{l}}\partial_{\tilde{z}_{j}}\partial_ {z_{\tilde{m}}}\partial_{\tilde{z}_{l}}\int dz_{1}\mathcal{N}(z_{1},\tilde{ \mathbf{z}}|0,K)\\ &=\big{(}1+v_{\widetilde{ij}\widetilde{m}\widetilde{l}}\partial _{z_{l}}\partial_{z_{j}}\partial_{z_{\tilde{m}}}\partial_{\tilde{z}_{l}} \big{)}\mathcal{N}(\tilde{\mathbf{z}}|0,\tilde{K})\\ &=q(\tilde{\mathbf{z}})\;,\end{split} \tag{7}\]

indeed reproduces the joint distribution for \(z_{2}\), \(z_{3},\cdots\) in the original form Eq. (5) as the smaller covariance matrix \(\tilde{K}\) is the submatrix of \(K\) excluding the first row and column. For simplicity, we use \(v\) to denote \(V/24\). In deriving above, we have employed the fact that \(\int dz_{1}(\partial_{z_{1}})^{*}\mathcal{N}(z_{1},\tilde{\mathbf{z}})=0\) for any integer power \(s\geq 1\), meaning that the indices \(i,j,m,l\) shall exclude any contribution associated with \(z_{1}\). Thus, the tilded indices \(\tilde{i}\in\{2,3,4,\cdots\}\), and they can be factored out of the integral.

Here we wish to stress that the above marginal property of non-Gaussian prior may seems straightforward with the differential representation of Edgeworth expansion, but it does not seem promising to reach the same conclusion if the marginalization is carried out with the explicit representation in Appendix A. Nevertheless, one can easily write the conditional distribution \(q(z_{1}|\tilde{\mathbf{z}})=q(z_{1},\tilde{\mathbf{z}})/q(\tilde{\mathbf{z}})\). More interestingly, the conditional mean is useful for prediction with noiseless observation of data, and also shed light on the effect of non-Gaussianity due to the finite width. The details of deriving the following conditional mean

\[\begin{split}\mathbb{E}[z_{1}|\tilde{\mathbf{z}}]& =\int dz_{1}z_{1}\frac{q(z_{1},\tilde{\mathbf{z}})}{q(\tilde{ \mathbf{z}})}\\ &=\mu(\tilde{\mathbf{z}})+\frac{1}{q(\tilde{\mathbf{z}})}A_{ \widetilde{ij}\widetilde{m}}\partial_{\tilde{z}_{1}}\partial_{\tilde{z}_{ \tilde{j}}}\partial_{z_{\tilde{m}}}\mathcal{N}(\tilde{\mathbf{z}}|0,\tilde{K} )\;,\end{split} \tag{8}\]

with the finite-width correction term proportional to the third-order tensor,

\[A_{\widetilde{ij}\widetilde{m}}=\frac{V_{\widetilde{ij}\widetilde{m}\widetilde {l}}}{6}[\tilde{K}^{-1}\mathbf{k}]_{\widetilde{l}}-\frac{V_{\widetilde{1} \widetilde{ij}\widetilde{m}}}{6}\;, \tag{9}\]
can be found in Appendix B. In the GP limit, the conditional mean becomes the well-known result \(\mu(\tilde{\mathbf{z}})=\mathbf{k}^{t}\tilde{K}^{-1}\tilde{\mathbf{z}}\), as the fourth cumulant \(V\) vanishes in large \(N\) limit. We also remind the readers that the tilded symbols are associated with the conditioned outputs \(z_{2,3,\cdots}\). In a similar manner, we can also show the conditional second moment,

\[\begin{split}\mathbb{E}[z_{1}^{2}|\tilde{\mathbf{z}}]& =\int dz_{1}z_{1}^{2}\frac{q(z_{1},\tilde{\mathbf{z}})}{q(\tilde{ \mathbf{z}})}\\ &=\mu^{2}(\tilde{\mathbf{z}})+\sigma^{2}+\frac{1}{q(\tilde{ \mathbf{z}})}\bigg{[}2\mu(\tilde{\mathbf{z}})A_{\tilde{\imath}\tilde{\jmath} \tilde{m}}\partial_{z_{\tilde{\imath}}}\partial_{z_{\tilde{\jmath}}}\partial_{z _{\tilde{\imath}}}+B_{\tilde{\imath}\tilde{\jmath}}\partial_{z_{\tilde{\imath} }}\partial_{z_{\tilde{\jmath}}}\bigg{]}\mathcal{N}(\tilde{\mathbf{z}}|0, \tilde{K})\,,\end{split} \tag{10}\]

with the second-order tensors denoting

\[B_{\tilde{\imath}\tilde{\jmath}}=\frac{V_{\tilde{\imath}\tilde{\jmath}\tilde {m}\tilde{\jmath}}}{2}[\tilde{K}^{-1}\mathbf{k}]_{\tilde{m}}[\tilde{K}^{-1} \mathbf{k}]_{\tilde{l}}-V_{1\tilde{\imath}\tilde{\jmath}\tilde{m}}[\tilde{K}^ {-1}\mathbf{k}]_{\tilde{m}}+\frac{V_{11\tilde{\imath}\tilde{\jmath}}}{2}\;. \tag{11}\]

The conditional second moment coincides with that in the GP limit, \(\sigma^{2}=K_{11}-\mathbf{k}^{t}\tilde{K}^{-1}\mathbf{k}\). The details of derivation can also be found in Appendix B. It is worthwhile to note that the conditional variance, \(\mathbb{E}_{q}[z_{1}^{2}|\tilde{\mathbf{z}}]-(\mathbb{E}_{q}[z_{1}|\tilde{ \mathbf{z}}])^{2}\), can include contributions related to function values \(\tilde{\mathbf{z}}\) through the nonzero \(A\) and \(B\) terms in above expressions. The lack of such dependence in GP limit is in fact a shortcoming for modelling, which has motivated the study of non-Gaussian prior in machine learning community (for example the Student-t process in [21]).

### An example: bivariate distribution and prediction

Now let us pause to consider a simple example where a shallow network defined in Eq. (1) with width \(N\) is used to model two input-output pairs \((\mathbf{x}_{1},z_{1})\) and \((\mathbf{x}_{2},z_{2})\). In the end, we shall present the predictive mean and variance for \(z_{1}\) conditioning on \(z_{2}\). The following notations are used for simplification: the covariance \(c:=\exp(-|\mathbf{x}_{1}-\mathbf{x}_{2}|^{2}/2d)\), the 2-by-2 covariance matrix \(\Sigma=\left(\begin{smallmatrix}1&c\\ c&1\end{smallmatrix}\right)\) for the bivariate prior, and \(z_{21}:=(z_{2}-cz_{1})/(\sqrt{1-c^{2}})\). Besides, the relations of derivatives of Gaussian in terms of Hermite polynomial shall be useful,

\[\partial_{z_{1}}^{n}\partial_{z_{2}}^{m}\mathcal{N}(z_{1},z_{2}|0,\Sigma) \propto\partial_{z_{1}}^{n}\left[\mathcal{N}(z_{1}|0,1)\partial_{z_{2}^{ \prime}}^{m}\mathcal{N}(z_{21}|0,1)\right]\propto(-1)^{m}\frac{\mathcal{N}(z_{ 2}|0,1)}{\sqrt{(1-c^{2})^{m}}}\partial_{z_{1}}^{n}[H_{m}(z_{21})\mathcal{N}(z _{12}|0,1)]\,, \tag{12}\]

where we have used \(\mathcal{N}(z_{1},z_{2}|0,\Sigma)\propto\mathcal{N}(z_{1}|0,1)\mathcal{N}(z_{ 21}|0,1)\) up to some irrelevant constant. The probabilist's Hermite polynomials are defined as \(H_{n}(z)=(-1)^{n}[\mathcal{N}(z|0,1)]^{-1}d_{z}^{n}\mathcal{N}(z|0,1)\)[19]. Consequently, following the previous discussion, the non-Gaussian bivariate distribution is shown to be

\[\begin{split} q(z_{1},z_{2})&=\mathcal{N}(z_{1},z_ {2}|0,\Sigma)\bigg{\{}1+\frac{\gamma^{4}}{16N}[H_{4}(z_{21})+H_{4}(z_{12})]+ \frac{c\gamma^{4}}{4N}[3cH_{2}(z_{21})+H_{3}(z_{21})H_{1}(z_{12})+z_{12} \leftrightarrow z_{21}]\\ &+\frac{(c^{4}+2)\gamma^{4}}{8N}[2c^{2}+4cH_{1}(z_{21})H_{1}(z_{ 12})+H_{2}(z_{21})H_{2}(z_{12})]\bigg{\}}\;,\end{split} \tag{13}\]

and the symbol \(\gamma=1/\sqrt{1-c^{2}}\) is used to simplify the expression. Besides, the fourth cumulant for the cosine network, \(V_{1111}=3/2N\), \(V_{1112}=3c/2N\), and \(V_{1122}=(c^{4}+2)/2N\), have been plugged in above. As for the conditional distribution, e.g. \(q(z_{1}|z_{2})=q(z_{1},z_{2})/q(z_{2})\), for the bivariate case, one can simply divid the above with \(q(z_{2})=\mathcal{N}(z_{2}|0,1)[1+H_{4}(z_{2})/16N]\), which shall yield a conditional Gaussian \(\mathcal{N}(z_{1}|z_{2})\) multiplied by a factor representing the finite-width effect.

We can directly apply Eq. (8) to obtain the conditional mean in the simple example where the tilded indices there only apply to \(z_{2}\). Surprisingly, the conditional mean for such noiseless case is

\[\mathbb{E}[z_{1}|z_{2}]=cz_{2}\,, \tag{14}\]

coinciding with that in the GP limit because the _accidental_ cancelation in the third-order tensor \(A_{222}\propto(cV_{222}-V_{1222})\). We shall show that the cancelation does not occur in the noisy (next Section) and more general cases. Application of Eq. (10) together with the above conditional mean results in the following conditional variance,

\[\text{Var}[z_{1}|z_{2}]=(1-c^{2})\bigg{[}1+\frac{4(2-c^{2})}{16N+H_{4}(z_{2})}H_ {2}(z_{2})\bigg{]}\;, \tag{15}\]

which consists of the corresponding variance in GP limit along with a term depending on \(z_{2}^{2}\).
Bayesian regression with non-Gaussian prior

Having established the marginal and conditional properties of the non-Gaussian distribution in Eq. (5), we now continue investigating the posterior distribution over the unseen function value \(z_{*}\) associated with input \(\mathbf{x}_{*}\) when the noisy observations \(\mathcal{D}=\{(\mathbf{x}_{1},y_{1}),(\mathbf{x}_{2},y_{2}),\cdots\}\) are known. According to Bayes's rule, the objective distribution is the posterior distribution defined as,

\[q(z_{*}|\mathcal{D})=\frac{q(z_{*},\mathcal{D})}{q(\mathcal{D})}=\frac{\int d \mathbf{z}\ q(z_{*},\mathbf{z})\mathcal{N}(\mathbf{y}|\mathbf{z},\Lambda)}{ \int d\mathbf{z}\ q(\mathbf{z})\mathcal{N}(\mathbf{y}|\mathbf{z},\Lambda)}\;, \tag{16}\]

with \(\Lambda\) denoting a diagonal matrix representing the noise in Gaussian likelihood. When the \(q\)'s in the numerator and denominator are both Gaussian, which corresponds to the infinitely wide network, the limiting distribution reads

\[q_{\infty}(z_{*}|\mathcal{D}) =\int d\mathbf{z}\ \mathcal{N}(z_{*}|\mathbf{k}_{*}^{\prime}K^{-1} \mathbf{z},k_{**}-\mathbf{k}_{*}^{\prime}K^{-1}\mathbf{k}_{*})\ \mathcal{N}(\mathbf{z}|K(K+\Lambda)^{-1}\mathbf{y},\Lambda K(K+\Lambda)^{-1})\] \[=\mathcal{N}\big{(}z_{*}|\mathbf{k}_{*}^{\prime}(K+\Lambda)^{-1} \mathbf{y},k_{**}-\mathbf{k}_{*}^{\prime}(K+\Lambda)^{-1}\mathbf{k}_{*}\big{)}\;.\]

In above first equality, the first and second Gaussian distributions on the right hand side represent the conditional density \(q_{\infty}(z_{*}|\mathbf{z})\) and the posterior \(q_{\infty}(\mathbf{z}|\mathcal{D})\), respectively.

As for the case of finite-width, the involved Gaussian likelihood can be dealt with by continuing application of the properties of adjoint differential operator as well as derivative of independent Gaussians. The details of derivation of evidence \(q(\mathcal{D})\) can be seen in Appendix C, and the expression in terms of derivatives with respect to the observed \(y\)'s reads,

\[\begin{split} q(\mathcal{D})&=\int d\mathbf{z}\ \mathcal{N}(\mathbf{y}|\mathbf{z},\Lambda)q(\mathbf{z})\\ &=(1+v_{ijml}\partial_{y_{i}}\partial_{y_{j}}\partial_{y_{m}} \partial_{y_{m}})\mathcal{N}(\mathbf{y}|0,K+\Lambda)\;.\end{split} \tag{17}\]

Despite its simple form, the evidence term can not be further be manipulated like in the infinite-width case. The details of derivation of the posterior distribution for output \(z_{*}\) at the unseen input,

\[\begin{split} q(z_{*}|\mathcal{D})&=\frac{\int d \mathbf{z}\ \mathcal{N}(\mathbf{y}|\mathbf{z},\Lambda)q(z_{*},\mathbf{z})}{q( \mathcal{D})}\\ &=(1+v_{\hat{1}\hat{j}\hat{m}\hat{l}}\partial_{\hat{1}\hat{j} \hat{m}\hat{l}})\mathcal{N}\big{[}\big{(}\begin{smallmatrix}z_{*}\\ \mathbf{y}\end{smallmatrix}\big{)}\big{]}0,\big{(}\begin{smallmatrix}1& \mathbf{k}_{*}^{\prime}\\ \mathbf{k}_{*}& K+\Lambda\end{smallmatrix}\big{)}\big{]}/q(\mathcal{D})\;, \end{split} \tag{18}\]

can also be found in Appendix C. Here, the notations can be understood as follows: the hatted indices, \(\hat{i}\) for instance, additionally include the symbol \(*\) associated with test function value \(z_{*}\), and the derivative refers to \(\partial_{*}:=\partial_{z_{*}}\) (with respect to function value) and \(\partial_{i}:=\partial_{y_{i}}\) (with respect to the observed value).

We conclude this Section with considering the same simple regression example but with noisy observation in \(y_{2}\) at input \(\mathbf{x}_{2}\), and we wish to predict the unobserved function value \(z_{1}\). To stress the role of noise parameter \(\sigma_{n}\), we only show the predictive mean here,

\[\mathbb{E}[z_{1}|y_{2}]=c\alpha_{n}^{2}y_{2}-\frac{c\sigma_{n}^{2}\alpha_{n}^{ 5}H_{3}(\alpha_{n}y_{2})}{4N+\alpha_{n}^{4}H_{4}(\alpha_{n}y_{2})}\;, \tag{19}\]

with \(\alpha_{n}=1/\sqrt{1+\sigma_{n}^{2}}\). The correction term vanishes as the noise parameter does, and is odd with respect to the sign change of \(y_{2}\).

## V Deep bottlenecked network

Up to this point, we have focused on the emergence of non-Gaussianity in prior function distribution in wide and shallow network, as well as how the non-Gaussian prior affects the inference in a regression task. Our approach can be extended to the deep and wide neural networks as the corresponding prior distributions approach GP [5; 6]. However, such wide-width assumption leading to _weak_ non-Gaussianity is not valid for the bottlenecked networks in which the hidden layers are alternatively wide and narrow [14; 22].
Let us consider a two-layer bottlenecked feed forward network modeling the hierarchical mapping \(\mathbf{x}\in\mathbb{R}^{d}\mapsto\mathbf{z}^{(1)}\in\mathbb{R}^{H}\mapsto z^{(2) }\in\mathbb{R}\). The hierarchy consists of the following computations,

\[\begin{split} z^{(1)}_{i,\alpha}&=\sqrt{\frac{2}{N_{ 1}}}\sum_{j=1}^{N_{1}}w^{(1)}_{ij}\cos(\frac{\Omega^{(1)}_{j}\cdot\mathbf{x}_{ \alpha}}{\sqrt{d}}+\phi^{(1)}_{j})\,,\\ z^{(2)}_{\alpha}&=\sqrt{\frac{2}{N_{2}}}\sum_{j=1}^ {N_{2}}w^{(2)}_{j}\cos(\frac{\Omega^{(2)}_{j}\cdot\mathbf{z}^{(1)}_{\alpha}}{ \sqrt{H}}+\phi^{(2)}_{j})\;,\end{split} \tag{20}\]

where the parameters \(w\)'s, \(\Omega\)'s, and \(\phi\)'s share the same prior distribution with their counterparts in the shallow network. The widths of hidden layers, \(N_{1}\) and \(N_{2}\), are large but the hidden output dimension \(H\) is not. In the limit \(N_{1,2}\to\infty\) while \(H\) remains finite, the two-layer cosine bottlenecked network corresponds to the prior of deep Gaussian process [15; 23; 24] with Gaussian kernel, which is a flexible and expressive function prior due to its compositional nature. As its name suggests, the conditional prior \(z^{(2)}|\mathbf{z}^{(1)}\) is a GP and so is each component \(z^{(1)}_{i}|\mathbf{x}\) in the hidden output. However, the marginal distribution for \(z^{(2)}|\mathbf{x}\) is non-Gaussian as [13] showed that the fourth cumulant is positive, i.e., a heavy-tailed distribution [25].

Here, we are interested in tracking how the small parameters, \(1/N_{1}\) and \(1/N_{2}\), and the bottleneck parameter \(1/H\) enters the second moment. For a deep and finitely wide linear network, the prior is non-Gaussian [8; 26] but the second moment does not receive correction due to the finite width [8; 9]. Thus, the effects of nonlinear activation on higher statistical moments are interesting (also see a recent work in [27]). For the deep cosine network in Eq. (20), we can first consider the wide limit, \(N_{1,2}\to\infty\), but the bottleneck width \(H\) remains finite. Following the result in Eq. (2), the covariance for the deep model can be computed as followings,

\[\begin{split} k^{(2)}_{\infty}(\mathbf{x}_{\alpha},\mathbf{x}_{ \beta})&=\mathbb{E}_{q_{\infty}}\bigg{\{}\exp\big{(}-\frac{| \mathbf{z}^{(1)}_{\alpha}-\mathbf{z}^{(1)}_{\beta}|^{2}}{2H}\big{)}\bigg{\}} \\ &=\bigg{[}\int dz_{\alpha}dz_{\beta}\ e^{-\frac{(z_{\alpha}-z_{ \beta})^{2}}{2H}q_{\infty}(z_{\alpha},z_{\beta})}\bigg{]}^{H}\\ &=\big{[}1+\frac{1-\exp(-\frac{|\mathbf{x}_{\alpha}-\mathbf{x}_{ \beta}|^{2}}{2d})}{H/2}\big{]}^{-H/2}\,.\end{split} \tag{21}\]

Note that the decomposition \(\mathcal{N}(z_{1},z_{2}|0,K)=\mathcal{N}(z_{1}+z_{2}|0,\sigma^{2}_{+})\mathcal{ N}(z_{1}-z_{2}|0,\sigma^{2}_{-})\)[28] has been used in derivation with the variances \(\sigma^{2}_{\pm}=2(K_{11}\pm K_{12})\). We want to stress that the above kernel is exact for any \(H\). For a very narrow case, i.e., \(H=1\), the above result coincides with that in [13], while in the very wide bottleneck case, \(H\gg 1\), it can be shown that

\[k^{(2)}_{\infty}(\mathbf{x}_{\alpha},\mathbf{x}_{\beta})\approx\exp[k( \mathbf{x}_{\alpha},\mathbf{x}_{\beta})-1]\big{\{}1+\frac{1}{H}[k(\mathbf{x} _{\alpha},\mathbf{x}_{\beta})-1]^{2}\big{\}}\,, \tag{22}\]

where the covariance \(k\) in shallow network is given in (2). In addition, the exponential of shallow kernel in above corresponds to \(H\to\infty\), which is the same as the deep kernel in [29]. Indeed, as suggested in [9], the appearance of reciprocal width, \(1/H\), in the correction term is in association with the weak non-Gaussianity even though the other widths \(N_{1,2}\) are infinite.

Next, we shall compute the covariance for the case where all widths are finite. In fact, the outer width \(N_{2}\) does not enter the kernel and we only need large but finite \(N_{1}\). The computation is similar and one only needs to replace the limiting joint distribution \(q_{\infty}(z^{(1)}_{\alpha},z^{(1)}_{\beta})\) with the non-Gaussian \(q\). The details of derivation using the property derivative of Gaussian can be seen in Appendix D. The kernel of the two-layer bottlenecked cosine network with \(N_{1,2}<\infty\) reads,

\[k^{(2)}(\mathbf{x}_{\alpha},\mathbf{x}_{\beta})=\bigg{[}1+\frac{1-k(\mathbf{x} _{\alpha},\mathbf{x}_{\beta})}{H/2}\bigg{]}^{-H/2}(1+\epsilon)^{H}\;, \tag{23}\]

with the correction due to finite \(N_{1}\),

\[\epsilon=\frac{V_{\alpha\alpha\alpha\alpha}+V_{\beta\beta\beta\beta}-4V_{ \alpha\beta\beta\beta}-4V_{\beta\alpha\alpha\alpha}+6V_{\alpha\alpha\beta \beta}}{24}[\frac{3}{H^{2}}-\frac{6\sigma^{2}_{-}}{H^{3}(1+\sigma^{2}_{-}/H) }+\frac{3\sigma^{4}_{-}}{H^{4}(1+\sigma^{2}_{-}/H)^{2}}]\;, \tag{24}\]

with \(\sigma^{2}_{-}=2(1-k)\) is used for easing the notation. Again, the above result is for general \(H\). For the wide bottleneck limit, one can show that these small parameters enter the kernel with the correction of \(O(\frac{1}{H})\) followed by \(O(\frac{1}{N_{1}H})\).
Hence, the deep model as a function prior serves as a GP with random kernel whose mean value converges to Eq. (22) when all widths \(N_{1,2}\) and \(H\) approach infinity. However, the role of \(H\) is different from \(N_{1}\) from our perturbative analysis, and \(1/N_{1}\) does not appear alone in the small parameter. The influence of finite width on the kernel of deep and nonlinear and convolutional models is also studied in [30; 31; 32].

The fourth moment can be computed in a similar manner. The closed form for \(N_{1,2}\rightarrow\infty\) can be found in [13]. For the finite width case, the fourth cumulant \(V_{\alpha\beta\gamma\delta}^{(2)}\) becomes the \(H\)-th power of the permutational symmetrization of \(\frac{1}{N_{2}}\mathbb{E}_{q}\{\exp[-(z_{\alpha}^{(1)}+z_{\beta}^{(1)}-z_{ \gamma}^{(1)}-z_{\delta}^{(1)})^{2}/2H)\}\), the small parameters of which shall consist of \(1/N_{2}\), \(1/(N_{2}H)\), and \(1/(N_{2}HN_{1})\). The analysis of scaling with respect to the reciprocal width is more complex than the deep linear network.

## VI Discussions

In essence, the finite width in random neural networks induces nonzero variance of kernel since the fourth cumulant is not zero. From function space perspective, the shallow networks with finitely large width can be regarded as a GP but the kernel itself is a random variable too, so the learned data representation is a distribution over kernels. Therefore, neural network with very narrow layers may not have enough expressive power for learning while the network with very wide layers is not flexible enough because it is equivalent to GP learning with one fixed kernel. It was recently suggested in Ref. [33] through studying the partition function [34; 35] of finite deep network that Student-t process [21] may suitably represent finite-width network, which offers an alternative description than this paper. The reports in [24; 30] observed the degradation of performance for Bayesian deep neural network when expanding the widths while the study in [36] suggested otherwise.

Taking the wide limit with the perturbative approach is appealing because of the analytic and elegant expressions for inference, which may shed light on future investigation of alternative algorithm for classification [37], study of average test error in regression [38; 39] and classification [40] tasks using finite-width networks. In this paper, the conditional statistics and perturbed posterior distribution over the unseen output using the non-Gaussian prior are obtained with help of the differential representation of multivariate Edgeworth expansion. In parallel with the equivalence between GP prior and the random wide neural network, the investigation of neural tangent kernel [41; 17; 42] is important for understanding the learning dynamics during the optimization and the relation with Bayesian neural network learning [43; 44].

## Acknowledgement

The research is supported by the Dean Office of School of Arts and Science at Rutgers University Newark. Correspondences with Jacob Zavatone-Veth and Pietro Rotondo are acknowledged.
