Posterior Inference on Shallow Infinitely Wide Bayesian Neural Networks
  under Weights with Unbounded Variance
From the classical and influential works of Neal (1996), it is known that the
infinite width scaling limit of a Bayesian neural network with one hidden layer
is a Gaussian process, when the network weights have bounded prior variance.
Neal's result has been extended to networks with multiple hidden layers and to
convolutional neural networks, also with Gaussian process scaling limits. The
tractable properties of Gaussian processes then allow straightforward posterior
inference and uncertainty quantification, considerably simplifying the study of
the limit process compared to a network of finite width. Neural network weights
with unbounded variance, however, pose unique challenges. In this case, the
classical central limit theorem breaks down and it is well known that the
scaling limit is an $\alpha$-stable process under suitable conditions. However,
current literature is primarily limited to forward simulations under these
processes and the problem of posterior inference under such a scaling limit
remains largely unaddressed, unlike in the Gaussian process case. To this end,
our contribution is an interpretable and computationally efficient procedure
for posterior inference, using a conditionally Gaussian representation, that
then allows full use of the Gaussian process machinery for tractable posterior
inference and uncertainty quantification in the non-Gaussian regime.
Jorge Loría, Anindya Bhadra
Posterior Inference on Infinitely Wide Bayesian Neural Networks under Weights with Unbounded Variance

###### Abstract

From the classical and influential works of Neal (1996), it is known that the infinite width scaling limit of a Bayesian neural network with one hidden layer is a Gaussian process, _when the network weights have bounded prior variance_. Neal's result has been extended to networks with multiple hidden layers and to convolutional neural networks, also with Gaussian process scaling limits. The tractable properties of Gaussian processes then allow straightforward posterior inference and uncertainty quantification, considerably simplifying the study of the limit process compared to a network of finite width. Neural network weights with unbounded variance, however, pose unique challenges. In this case, the classical central limit theorem breaks down and it is well known that the scaling limit is an \(\alpha\)-stable process under suitable conditions. However, current literature is primarily limited to forward simulations under these processes and the problem of posterior inference under such a scaling limit remains largely unaddressed, unlike in the Gaussian process case. To this end, our contribution is an interpretable and computationally efficient procedure for posterior inference, using a _conditionally Gaussian_ representation, that then allows full use of the Gaussian process machinery for tractable posterior inference and uncertainty quantification in the non-Gaussian regime.

## 1 Introduction

Gaussian processes (GPs) have been studied as the infinite width limit of Bayesian neural networks with priors on network weights that have finite variance (Neal, 1996). This presents some key advantages over Bayesian neural networks with finite widths that usually require computation intensive Markov chain Monte Carlo (MCMC) posterior calculations (Neal, 1996) or variational approximations (Goodfellow et al., 2016, Chapter 19); in contrast to straightforward posterior inference and probabilistic uncertainty quantification afforded by the GP machinery (Williams and Rasmussen, 2006). In this sense, the work of Neal (1996) is foundational. The technical reason for this convergence to a GP is due to an application of the central limit theorem under the bounded second moment condition. More specifically, given an \(I\) dimensional input \(\mathbf{x}\) and a one-dimensional output \(y(\mathbf{x})\), a \(K\) layer feedforward deep neural network (DNN) with \(K-1\) hidden layers is defined by the recursion:

\[z_{j}^{(l+1)}(\mathbf{x}) =g\left(b_{j}^{(l)}+\sum_{i=1}^{p_{l}}w_{ij}^{(l)}z_{i}^{(l)}( \mathbf{x})\right),\quad l=1,\ldots,K-1, \tag{1}\] \[y(\mathbf{x}) =\sum_{j=1}^{p_{K}}w_{j}^{(K)}z_{j}^{(K)}(\mathbf{x}), \tag{2}\]
where \(z^{(1)}\equiv\mathbf{x},\ p_{1}=I,\ p_{K}=D\) and \(g(\cdot)\) is a nonlinear activation function. Thus, the network repeatedly applies a linear transformation to the inputs at each layer, before passing it through a nonlinear activation function. Sometimes a nonlinear transformation is also applied to the final hidden layer to the output layer, but in this paper it is assumed the output is a linear function of the last hidden layer. Neal (1996) considers the case of a Bayesian neural network with a single hidden layer, i.e., \(K=2\). So long as the hidden to output weights \(w^{(2)}\) are independent and identically distributed Gaussian, or at least, have a common bounded variance given by \(c/p_{2}\) for some \(c>0\), and \(g(\cdot)\) is bounded, an application of the classical central limit theorem shows the network converges to a GP as the number of hidden nodes \(p_{2}\rightarrow\infty\). Recently, Neal's result has been extended to prove fully connected multi-layer feedforward networks (Lee et al., 2018) and convolutional neural networks (Garriga-Alonso et al., 2019) also converge to GPs. This is useful for uncertainty quantification by designing emulators for deep neural networks (DNNs) based on GPs, since the behavior of finite-dimensional DNNs for direct uncertainty quantification is much harder to characterize. In contrast, once a convergence to GP can be ensured, well established tools from the GP literature (see, e.g., Williams and Rasmussen, 2006) can be brought to the fore to allow straightforward posterior inference. The induced covariance function depends on the choice of the nonlinear activation function \(g(\cdot)\), and is in general anisotropic. However, it can be worked out in explicit form under a variety of activation functions for both shallow (Neal, 1996; Cho and Saul, 2009) and deep (Lee et al., 2018) feedforward neural networks, where for deep networks usually a recursive formula is available that expresses the covariance function of a given layer conditional on the layer just below it. The benefit of depth is that it allows a potentially very rich covariance function at the level of the observed data, even if the covariances in each layer conditional on the layer below are simple. Viewing a GP as a prior on the function space, this allows for a rich class of prior structures. However, the process is still Gaussian in all these cases and our intention in this paper is a departure from the Gaussian world.

### Challenges posed by network weights with unbounded prior variance

Although the GP literature has been immensely influential for uncertainty quantification in DNNs, it is obvious that a DNN does not converge to a GP if the final hidden to output layer weights are allowed to have unbounded variance, e.g., belonging to \(t\) or others in the stable family, such that the scaling limit distribution is non-Gaussian (Gnedenko and Kolmogorov, 1954). This was already observed by Neal (1996) who admits: _"in contrast to the situation for Gaussian process priors, whose properties are captured by their covariance functions, I know of no simple way to characterize the distributions over functions produced by the priors based on non-Gaussian stable distributions."_ Faced with this difficulty, Neal (1996) confines himself to forward simulations from DNNs with \(t\) weights, and yet, observes that the network realizations under these weights demonstrate very different behavior (e.g., large jumps) compared to normal priors on the weights. This is not surprising, since Gaussian processes, with their almost surely continuous sample paths, are not necessarily good candidate models for functions containing sharp jumps, perhaps explaining their lack of popularity in certain application domains, e.g., finance, where jumps and changepoints need to be modeled (see, e.g., Chapter 7 of Cont and Tankov, 2004). Another key benefit of priors with polynomial tails, pointed out by Neal (1996), is that it allows a few hidden nodes to make a large contribution to the output, while drowning out the others, akin to feature selection. In contrast, in the GP limit, the contributions of individual nodes are averaged out. Thus, there are clear motivations for developing computationally feasible posterior inference machinery under these non-Gaussian limits.

Neal (1996) further hints that it may be possible to prove an analogous result using priors that do not have finite variance. Specifically, suppose the network weights are given symmetric \(\alpha\)-stable priors. If \(X\) is an \(\alpha\)-stable random variable, the density does not in general have a closed form, but the characteristic function is \(\phi_{X}(t)=\exp[it\mu-\nu^{\alpha}|t|^{\alpha}\{1-i\beta\text{sign}(t)\omega(t; \alpha)\}]\), where \(\omega(t;\alpha)=\tan(\alpha\pi/2)\), for \(\alpha\neq 1\) and \(\omega(t;\alpha)=-(2/\pi)\log[t]\), for \(\alpha=1\). Here \(\mu\in\mathbb{R}\) is called the shift parameter, \(\alpha\in(0,2]\) is the index parameter, \(\beta\in[-1,1]\) is the symmetry parameter, and \(\nu>0\) is the scale parameter (Samorodnitsky and Taqu, 1994, p. 5). Throughout, we use a zero shift (\(\mu=0\)) stable variable, and denote it by \(X\sim S(\alpha,\beta,\nu)\). Here \(\beta=0\) corresponds to the symmetric case, and when \(\beta=1,\alpha<1,\nu=1\), the random variable is strictly positive, which we denote by \(S^{+}(\alpha)\). Der and Lee (2005) confirm Neal's conjecture by establishing that the scaling limit of a shallow neural network under \(\alpha\)-stable priors on the weights is an \(\alpha\)-stable process. Proceeding further, Peluchetti et al. (2020) show that the limit process for infinitely wide DNNs with infinite-variance priors is also an \(\alpha\)-stable processes. However, both Der and Lee (2005) and Peluchetti et al. (2020) only consider
the forward process and neither considers posterior inference. Inference using \(\alpha\)-stable densities is not straightforward, and some relevant studies are by Samorodnitsky and Taqqu (1994), Lemke et al. (2015), and more recently by Nolan (2020). The main challenge is that a covariance function is not necessarily defined, precluding posterior inference analogous to the GP case, for example, using the _kriging_(Stein, 1999) machinery. To this end, our contribution lies in using a representation of the characteristic function of symmetric \(\alpha\)-stable variables as a normal scale mixture, that then allows a _conditionally Gaussian_ representation. This makes it possible to develop posterior inference and prediction techniques under stable priors on network weights using a latent Gaussian framework.

### Summary of main contributions

Our main contributions consist of:

1. An explicit characterization of the posterior predictive density function under infinite width scaling limits for shallow (one hidden layer) Bayesian neural networks under stable priors on the network weights, using a conditionally Gaussian representation.
2. An MCMC algorithm for posterior inference and prediction, with publicly available code.
3. Numerical experiments in one and two dimensions that validate our procedure by obtaining better posterior predictive properties for functions with jumps and discontinuities, compared to both Gaussian processes and Bayesian neural networks of finite width.
4. A real world application to the air quality data in the city of York, where our results show a smaller prediction error in out-of-sample observations.

## 2 Infinite width limits of Bayesian neural networks under weights with unbounded variance

Consider the case of a shallow, one hidden layer network, with the weights of the last layer being independent and identically distributed with symmetric \(\alpha\)-stable priors. Our results are derived under this setting using the following proposition of Der and Lee (2005).

**Proposition 1**.: _(_Der and Lee, 2005_)__. Let the network specified by Equations (1) and (2), with a single hidden layer \((K=2)\), have i.i.d. hidden-to-output weights \(w_{j}^{(2)}\) distributed as a symmetric \(\alpha\)-stable with scale parameter \((\nu/2)^{1/2}p_{2}^{-1/\alpha}\). Then \(y(\mathbf{x})\) converges in distribution to a symmetric \(\alpha-\)stable process \(f(\mathbf{x})\) as \(p_{2}\to\infty\). The finite dimensional distribution of \(f(\mathbf{x})\), denoted as \((f(\mathbf{x}_{1}),\ldots,f(\mathbf{x}_{n}))\) for all \(n\), where \(\mathbf{x}_{i}\in\mathbb{R}^{I}\), is multivariate stable and admits a characteristic function:_

\[\phi(\mathbf{t})=\exp\left\{-(\nu/2)^{\alpha/2}\mathbb{E}[|\langle\mathbf{t}, \mathbf{g}\rangle|^{\alpha}]\right\}, \tag{3}\]

_where angle brackets denote the inner product, \(\mathbf{g}=\left(g(\mathbf{x}_{1}),\ldots,g(\mathbf{x}_{n})\right)\), and \(g(\mathbf{x})\) is a random variable with the common distribution (across \(j\)) of \(z_{j}^{(2)}(\mathbf{x})\)._

Following Neal (1996), assume for the rest of the paper that the activation function \(g(\cdot)\) corresponds to the sign function: \(\mathrm{sign}(z)=1\), if \(z>0\); \(\mathrm{sign}(z)=-1\), if \(z<0\); and \(\mathrm{sign}(0)=0\). For \(\mathbf{z}\in\mathbb{R}^{I}\) we define \(g(\mathbf{z})=\mathrm{sign}\left(b_{0}+\sum_{i=1}^{I}w_{i}z_{i}\right)\), where \(b_{0}\) and \(w_{i}\) are i.i.d. standard Gaussian variables. The next challenge is to compute the expectation within the exponential in Equation (3). To resolve this, we break it into simpler cases. Define \(\Lambda\) as the set of all possible functions \(\tau:\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}\to\{-1,+1\}\). Noting that each \(\mathbf{x}_{j}\) can be mapped to two possible options: \(+1\) and \(-1\), indicates that there are \(2^{n}\) elements in \(\Lambda\). For each \(\ell=1,\ldots,2^{n}\), consider \(\tau_{\ell}\in\Lambda\), the event \(A_{\ell}=\{\tau_{\ell}(\mathbf{x}_{j})=g(\mathbf{x}_{j})\}_{j=1}^{n}\), and the probability \(q_{\ell}=\mathbb{P}(A_{\ell})\). By definition \(\{A_{\ell}\}_{\ell=1}^{2^{n}}\) is a set of disjoint events. Next, using the definition of the expectation of discrete disjoint events we obtain:

\[\mathbb{E}[|\langle\mathbf{t},\mathbf{g}\rangle|^{\alpha}]=\sum_{\ell=1}^{2^ {n}}q_{\ell}\left|\sum_{j=1}^{n}t_{j}\tau_{\ell}(\mathbf{x}_{j})\right|^{ \alpha}. \tag{4}\]
### Computation of \(q_{\ell}\) and \(\tau_{\ell}\)

Since the size of \(\Lambda\) is \(2^{n}\), indicating an exponential complexity of naive enumeration, we further simplify Equation (4). When the input points are arranged in a way that \(\tau_{\ell}\) is not possible, then \(q_{\ell}\) must be zero. We can identify the elements in \(\Lambda\) with positive probabilities by considering arbitrary values of \(b_{0},w_{1},\ldots,w_{I}\). The corresponding \(\tau\) function is determined by: \(\tau(\mathbf{x}_{j})=\mathrm{sign}(b_{0}+w_{1}x_{j1}+\cdots+w_{I}x_{j1})\), for each \(j\). This corresponds to labeling with \(+1\) the points above a hyperplane, and with \(-1\) the points that lie below the hyperplane. When \(I=1\), without loss of generality, let \(x_{1}<\cdots<x_{n}\). In this case \(\Lambda\) corresponds to the possible changes in sign that can occur between the input variables, which is equal to \(n\). Specifically, the sign change can occur before \(x_{1}\), between \(x_{1}\) and \(x_{2}\),..., between \(x_{n-1}\) and \(x_{n}\), and after \(x_{n}\). A similar argument is made in Example 2.1.1 of Der and Lee (2005), suggesting the possibility of considering more than one dimensions.

For \(I>1\), Harding (1967) studies the possible partitions of \(n\) points in \(\mathbb{R}^{I}\) by an \((I-1)\)-dimensional hyperplane--which corresponds to our problem, and determines that for points in general configuration there are \(O(n^{I})\) partitions. Goodman and Pollack (1983) give an explicit algorithm for finding the elements of \(\Lambda\) that have non-zero probabilities in any possible configuration. We summarize their algorithm for \(I=2\) as Algorithm S.1 in Supplementary Section S.2. This algorithm runs in a computational time of order \(n^{I}\log(n)\), which is reasonable for moderate \(I\). For the rest of the article, we denote the cardinality of elements in \(\Lambda\) that have positive probability by \(L\), with the understanding that \(L\) will depend on the input vectors \(\mathbf{x}_{i}\) that are used and their dimension. This solves the issue of computing deterministically the values of \(\tau_{\ell}\) that have positive probability.

Next we compute the probability \(q_{\ell}\) for the determined \(\tau_{\ell}\). For \(I=1\), the \(q_{\ell}\)s correspond to probabilities obtained from a Cauchy cumulative density function, which we state explicitly in Supplementary Section S.1. For general dimension of the input \(I>1\), the value of \(q_{\ell}\) is given by \(\mathbb{P}(\mathbf{Z}^{(\tau_{\ell})}>0)\), where the \(n\)-dimensional Gaussian vector \(\mathbf{Z}^{(\tau_{\ell})}\) has \(i\)-th entry given by \(\tau_{\ell}(\mathbf{x}_{i})(b_{0}+\sum_{j=1}^{I}w_{j}x_{ij})\). This implies that \(\mathbf{Z}^{(\tau_{\ell})}\sim\mathcal{N}_{n}(0,\mathbf{\Sigma}^{(\tau_{\ell} )})\), where the (possibly singular) variance matrix is \(\Sigma^{(\tau_{i,j})}_{i,j}=\tau_{\ell}(\mathbf{x}_{i})\tau_{\ell}(\mathbf{x}_ {j})(1+\sum_{k=1}^{n}x_{ik}x_{jk})\). This means we can compute \(q_{\ell}=\mathbb{P}(\mathbf{Z}^{(\tau_{\ell})}>0)\) using for example the R package mvtnorm, which implements the method of Genz and Bretz (2002) for evaluating multivariate Gaussian probabilities. Since the \(q_{\ell}\)s require independent procedures to be computed, this is easily parallelized after obtaining the partitions.

A characterization of the posterior predictive density under stable network weights using a conditionally Gaussian representation

While the previous section demonstrated the characteristic function of Equation (4) can be computed, the resulting density, obtained via its inverse Fourier transform, does not necessarily have a closed form, apart from specific values of \(\alpha\), such as \(\alpha=2\) (Gaussian), \(\alpha=1\) (Cauchy) or \(\alpha=0.5\) (inverse Gaussian). In this section we show that a _conditionally Gaussian_ characterization of the density function is still possible for the entire domain of \(\alpha\in(0,2]\), facilitating posterior inference. First, note that the result of Der and Lee (2005) is obtained assuming that there is no intrinsic error in the observation model, i.e., they assume the observations are obtained as \(y_{i}=f(\mathbf{x}_{i})\), and the only source of randomness is the network weights. We generalize this to more realistic scenarios and consider an additive error term, known as the nugget effect in spatial statistics (Stein, 1999, p. 94). That is, we consider the observation model \(y_{i}=f(\mathbf{x}_{i})+\varepsilon_{i}\), where the error terms \(\varepsilon_{i}\) are independent identically distributed normal random variables with constant variance \(\sigma^{2}\). Using the expression for the expectation in the characteristic function from Proposition 1, we derive the full probability density function, as specified in the following theorem, with proof in Appendix A.1.

**Theorem 1**.: _For real-valued observations \(\mathbf{y}=(y_{1},\ldots,y_{n})\) under the model \(y_{i}=f(\mathbf{x}_{i})+\varepsilon_{i};\) where \(\varepsilon_{i}\overset{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2})\) and \(f(\cdot)\) is as specified in Proposition 1, denote the matrix \(\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]^{T}\). The probability density function of \((\mathbf{y}\mid\mathbf{X})\) is:_

\[p(\mathbf{y}\mid\mathbf{X})=(2\pi)^{-n/2}\int_{(\mathbb{R}^{+})^{L}}\exp \left(-\frac{1}{2}\mathbf{y}^{T}\mathbf{Q}^{-1}\mathbf{y}\right)\det(\mathbf{ Q})^{-1/2}\prod_{\ell=1}^{L}p_{S^{+}}(s_{\ell})ds_{\ell},\]
_where \(p_{S^{+}}\) is the density for a positive \(\alpha/2\)-stable random variable, and \(\mathbf{Q}\) is a positive definite matrix with probability one that depends on \(\{s_{\ell}\}_{\ell=1}^{L}\). Specifically, \(\mathbf{Q}\) has diagonal entries \(\sigma^{2}+\nu\sum_{\ell=1}^{L}q_{\ell}^{2/\alpha}s_{\ell}\), and the entry \((i,j)\) is given by \(\nu\sum_{\ell=1}^{L}q_{\ell}^{2/\alpha}s_{\ell}\tau_{\ell}(\mathbf{x}_{i})\tau _{\ell}(\mathbf{x}_{j})\), for \(i\neq j\)._

Theorem 1 is the main machinery we need for posterior inference. We emphasize that \(\mathbf{Q}\) is a matrix with random entries, conditional on \(\{s_{\ell}\}_{\ell=1}^{L}\); and \(\{q_{\ell}\}_{\ell=1}^{L}\) and \(\{\tau_{\ell}\}_{\ell=1}^{L}\) are deterministic. Further, the input variables \(\mathbf{X}\) are also deterministic. Theorem 1 implies the hierarchical Gaussian model:

\[\mathbf{y}\mid\mathbf{X},\{s_{\ell}\}_{\ell=1}^{L}\sim\mathcal{N}_{n}(0, \mathbf{Q}),\ s_{\ell}\stackrel{{ i.i.d.}}{{\sim}}S^{+}(\alpha/2).\]

Each of the \(\tau\)s defines a level set for the points that lie in the \(+1\) side in contraposition to those that lie in the \(-1\) side. A forward simulation of this model is a weighted sum of the \(\tau\)s, with corresponding positive weight for those that lie closer and a negative weight for the points that lie further. In a similar vein, we present a corollary to interpret the distribution of the covariance matrix.

**Corollary 1**.: _The marginal distribution of the diagonal entries of the matrix \(\mathbf{Q}\) is \(\sigma^{2}+\nu S^{+}(\alpha/2)\), where the \(\sigma^{2}\) acts as a shift parameter, and the marginal distribution of the entry \(i,j\) in the \(\mathbf{Q}\) matrix is \(S(\alpha/2,\nu,2p_{ij}-1)\), where \(p_{ij}=\sum_{\ell:\tau_{\ell}(\mathbf{x}_{i})=\tau_{\ell}(\mathbf{x}_{j})}q_{\ell}\), the probability that \(\mathbf{x}_{i}\) and \(\mathbf{x}_{j}\) lie on the same side of the hyperplane partition. Further, the entries of \(\mathbf{Q}\) are not independent._

Proof.: For \(\{\mathbf{Q}\}_{ii}-\sigma^{2}\), we apply Property 1.2.1 of Samorodnitsky and Taqqu (1994), which we refer as the closure property, to obtain \(\nu S^{+}(\alpha/2)\). Next for \(\{\mathbf{Q}\}_{ij}\), we split the summation in two cases: \(\tau_{\ell}(\mathbf{x}_{i})=\tau_{\ell}(\mathbf{x}_{j})\), and \(\tau_{\ell}(\mathbf{x}_{i})\neq\tau_{\ell}(\mathbf{x}_{j})\). Using the closure property in the separate splits, we obtain \(\{\mathbf{Q}\}_{ij}\sim S(\alpha/2,\nu p_{ij}^{2/\alpha},1)-S(\alpha/2,\nu(1- p_{ij})^{2/\alpha},1)\). The result follows by applying the closure property once more. The entries of \(\mathbf{Q}\) are not independent, as they are obtained from a linear combination of the independent variables \(\{s_{\ell}\}_{\ell=1}^{L}\). 

The value of Corollary 1 does not lie in a numerical or computational speed-up, since the obtained marginals are not independent, but rather in the interpretability that it lends to the model. Explicitly, it indicates that when the points lie closer their conditional covariance is more likely to be positive and when they lie farther apart the covariance is more likely to be negative (see Proposition 1.2.14 of Samorodnitsky and Taqqu, 1994).

Now that the probability model is clear, we proceed to the next problem of interest: prediction. To this end, we present the following proposition, characterizing the posterior predictive density.

**Proposition 2**.: _Consider a vector of \(n\) real-valued observations \(\mathbf{y}=(y_{1},\ldots,y_{n})\), each with respective input variables \(\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\in\mathbb{R}^{I}\), under the model \(y_{i}=f(\mathbf{x}_{i})+\varepsilon_{i};\varepsilon_{i}\stackrel{{ i.i.d.}}{{\sim}} \mathcal{N}(0,\sigma^{2}),\) and \(m\) new input variable locations: \(\mathbf{x}_{1}^{*},\ldots,\mathbf{x}_{m}^{*}\in\mathbb{R}^{I}\), with future observations at those locations denoted by \(\mathbf{y}^{*}=(\mathbf{y}_{1}^{*},\ldots,\mathbf{y}_{m}^{*})\). Denote the matrices \(\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]^{T}\) and \(\mathbf{X}^{*}=[\mathbf{x}_{1}^{*},\ldots,\mathbf{x}_{m}^{*}]^{T}\). The posterior distribution at these new input variables satisfies the following properties:_

1. _The conditional posterior of_ \(\mathbf{y}^{*}\mid\mathbf{y},\mathbf{X},\mathbf{X}^{*},\mathbf{Q}\) _is an_ \(m\)_-dimensional Gaussian. Specifically:_ \[\mathbf{y}^{*}\mid\mathbf{y},\mathbf{X},\mathbf{X}^{*},\mathbf{Q}\sim\mathcal{ N}_{m}(\boldsymbol{\mu}^{*},\boldsymbol{\Sigma}^{*}),\] (5) _where_ \(\boldsymbol{\mu}^{*}=\mathbf{Q}_{*,1:n}\mathbf{Q}_{1:n,1:n}^{-1}\mathbf{y}\)_, and_ \(\boldsymbol{\Sigma}^{*}=\mathbf{Q}_{*,*}-\mathbf{Q}_{*,1:n}\mathbf{Q}_{1:n,1:n} ^{-1}\mathbf{Q}_{1:n,*}\)_, using_ \(\mathbf{Q}\) _as previously defined for the_ \(n+m\) _input variables, and denoting by_ \({}^{*}\)_' the entries_ \((n+1):(n+m)\)_._
2. _The posterior predictive density at the_ \(m\) _new locations conditional on the observations_ \(\mathbf{y}\)_, is given by:_ \[p(\mathbf{y}^{*}\mid\mathbf{y},\mathbf{X},\mathbf{X}^{*})=\int_{(\mathbb{R}^{ +})^{L}}p(\mathbf{y}^{*}\mid\mathbf{y},\mathbf{X},\mathbf{X}^{*},\mathbf{Q})p( \mathbf{Q}\mid\mathbf{y},\mathbf{X})d\mathbf{Q},\] (6) _where_ \(p(\mathbf{y}^{*}\mid\mathbf{y},\mathbf{X},\mathbf{X}^{*},\mathbf{Q})\) _is the conditional posterior density of_ \(\mathbf{y}^{*}\)_,_ \(\mathbf{Q}\) _is as previously described for the_ \(n+m\) _input variables, and the integral is over the values determined by_ \(\{s_{\ell}\}_{\ell=1}^{L}\)_._

Proof.: The first is an immediate application of Theorem 1 and the conditional density of a multivariate Gaussian. The second part follows from a standard application of marginal probabilities.
An MCMC algorithm to sample from the posterior predictive distribution

Dealing with \(\alpha\)-stable random variables includes the difficulty that the moments of the variables are only finite up to an \(\alpha\) power. Specifically, for \(\alpha<2\), if \(X\sim S(\alpha,\nu,\beta)\), then \(\mathbb{E}[|X|^{r}]=\infty\) if \(r\geq\alpha\), and is finite otherwise (Samorodnitsky and Taqqu, 1994, Property 1.2.16). To circumvent dealing with potentially ill-defined moments, we propose to sample from the full posterior. For fully Bayesian inference, we assign \(\sigma^{2}\) a half-Cauchy prior (Gelman, 2006) and iteratively sample from the posterior predictive distribution by cycling through \((\mathbf{y}^{*},\mathbf{Q},\sigma^{2})\) in an MCMC scheme, as described in Algorithm 1, which has computational complexity of the order of \(\mathcal{O}(T[(n+m)^{I}n^{2}+m^{3}])\), where \(T\) is the number of MCMC simulations used. The method of Chambers et al. (1976) is used for simulating the stable variables. An implementation of our algorithms is freely available online.

```
0: Vector of \(n\) real-valued observations \(\mathbf{y}\) with \(I\)-dimensional input variables \(\mathbf{X}\in\mathbb{R}^{n\times I}\), new input variables \(\mathbf{X}^{*}\in\mathbb{R}^{m\times I}\), and number of MCMC iterations \(T\).
0: Posterior predictive samples \(\{\mathbf{y}_{t}^{*}\}_{t=1}^{T}\)  Obtain \(\Lambda\) for \((\mathbf{X},\mathbf{X}^{*})\) using Algorithm S.1.  Compute \(\{q_{\ell}\}_{\ell=1}^{L}\) as described in Section 2.1.  Initialize \(\mathbf{Q}_{0}\) using independent samples of \(s_{\ell}\) from the prior distributions. for\(t=1,\ldots,T\)do  Simulate \(\mathbf{Q}_{t}\mid\mathbf{y},\mathbf{Q}_{t-1}\) using Algorithm S.2.  Compute \(\boldsymbol{\mu}_{t}^{*}\) and \(\boldsymbol{\Sigma}_{t}^{*}\) using \(\mathbf{Q}_{t}\) in Part 1 of Proposition 2.  Simulate \(\mathbf{y}_{t}^{*}\mid(\mathbf{y},\mathbf{Q}_{t})\sim\mathcal{N}_{m}(\boldsymbol {\mu}_{t}^{*},\boldsymbol{\Sigma}_{t}^{*})\). endfor return\(\{\mathbf{y}_{t}^{*}\}_{t=1}^{T}\).
```

**Algorithm 1** A Metropolis-Hastings sampler for the posterior predictive distribution

There are two hyper-parameters in our model: \(\alpha,\nu\). We propose to select them by cross validation on a grid of \((\alpha,\nu)\), and selecting the result with smallest mean absolute error (MAE). Another possible way to select \(\alpha\) is by assigning a prior. The natural choice for \(\alpha\) is a uniform prior on \((0,2)\), however the update rule would need to consider the densities of the \(L\) such \(\alpha/2\)-stable densities \(p_{S^{+}}\), which would be computationally intensive as there is no closed form expression for this density apart from specific values of \(\alpha\). A prior for \(\nu\) could be included but a potential issue of identifiability emerges, similar to what has been identified for the Matern kernel (Zhang, 2004). Our proposed method avoids dealing with these issues, which we leave open for future research.

## 4 Numerical experiments in one and two dimensions

To validate our method, we compare it against the predictions obtained from three methods. The first two are methods for Gaussian processes that correspond to the two main approaches in GP inference: maximum likelihood with a Gaussian covariance kernel (Dancik and Dorman, 2008), and an MCMC based Bayesian procedure using the Matern kernel (Gramacy and Taddy, 2010); the third method is a two-layer Bayesian neural network (BNN) using a single hidden layer of 100 nodes with Gaussian priors, implemented in pytorch(Paszke et al., 2019), and fitted via a variational approach. The choice of a modest number of hidden nodes is intentional, so that we are away from the inifinite width GP limit, and the finite-dimensional behavior can be visualized. The respective implementation are in the R packages mlegp, tgp, and the python libraries pytorch and torchbnn. The estimates used from these methods are respectively the kriging estimate, posterior median and posterior mean. We tune our method by cross-validation over a grid of \((\alpha,\nu)\). We use point-wise posterior median as the estimate, and report the values with smallest mean absolute error (MAE) and the optimal parameters. Results on uncertainty quantification, timing, and additional simulations are in Supplementary Section S.3. We use a data generating mechanism of the form \(y=f(x)+\varepsilon\), where the \(f\) is the true function. The overall summary is that when \(f\) has at least one discontinuity, our method performs better at prediction than the competing methods, and when \(f\) is continuous the proposed method performs just as well as the other methods. This provides empirical support that the assumption of _continuity_ of the true function cannot be disregarded in the _universal approximation_ property of neural networks (Hornik et al., 1989), and the adoption of infinite variance prior weights might be a crucial missing ingredient for successful posterior prediction when the truth is discontinuous.
### Experiments in one dimension

We consider a function with three jumps: \(f(x)=5\times\mathbf{1}_{\{x\geq 1\}}+5\times\mathbf{1}_{\{-1\leq x<0\}}\), to which we add a Gaussian noise with \(\sigma=0.5\). We consider \(x\in[-2,2]\), with 40 equally spaced points as the training set, and 100 equally spaced points in the testing set. We display the two-panel Figure 1, showing the comparisons between the four methods. The boxplots in the left panel show that the proposed method - which we term _"Stable,"_ has the smallest prediction error. The right panel shows that the BNN, the GP based fully Bayesian, and maximum likelihood methods have much smoother predictions than the Stable method. This indicates the inability of these methods to capture sharp jumps as well as the Stable method, which very clearly captures them. The Stable method obtains the smallest cross validation error for this case with \(\alpha^{*}=1.1\) and \(\nu^{*}=1\).

### Experiments in two dimensions

We consider the function: \(f(x_{1},x_{2})=5\times\mathbf{1}_{\{x_{1}>0\}}+5\times\mathbf{1}_{\{x_{2}>0\}}\), with additive Gaussian noise with \(\sigma=0.5\), and observations on an equally-spaced grid of 49 points in the square \([-1,1]^{2}\). In Figure 2 we display the boxplots and contour plots for all methods for out-of-sample prediction on an equally spaced grid of \(9\times 9\) points in the same square. The methods that employ Gaussian processes (GP Bayes and GP MLE) and BNN seem to have smoother transitions between the different quadrants, whereas the Stable method captures the sharp jumps better. This is reinforced by the prediction errors displayed in the left panel. For this example, the Stable method obtains the smallest cross validation error with \(\alpha^{*}=1.1,\nu^{*}=1\).

Figure 1: _Left:_ Boxplots of mean absolute error of out-of-sample prediction over test points, and _Right:_ predicted values over 100 points on a regular grid on \([-2,2]\). Training points in black dots.

Figure 2: _Left:_ Boxplots of mean absolute error (MAE) of out-of-sample prediction over test points, and _Right:_ predicted values over a \(9\times 9\) grid on \([-1,1]^{2}\).
## 5 A demonstration of out of sample prediction using York air quality data

Measurements in parts per million (ppm) on the level of nitrogen dioxide (\(\mathrm{NO}_{2}\)), a known environmental pollutant, are available in the city of York, UK, 1 at different locations. We use the data from July 2020, to obtain the mean absolute error in out-of-sample prediction. We split the outcomes randomly into two parts: validation and training, comprising of 60 and 128 observations respectively. To select the hyper-parameters in the Stable method we use random splits of the training set, and obtain the optimal point \((\alpha^{*},\nu^{*})=(1.9,1)\) by cross validation. We compare it to the three methods described in Section 4, and display the results in Table 1, with the Stable method performing the best. We include visualizations in Supplementary Section S.4.

Footnote 1: [https://data.yorkopendata.org/dataset/diffusion-tubes-locations](https://data.yorkopendata.org/dataset/diffusion-tubes-locations)

## 6 Conclusions and future work

We develop a novel method for posterior inference and prediction for infinite width limits of shallow (one hidden layer) BNNs under weights with infinite prior variance. While the \(\alpha\)-stable forward scaling limit in this case has been known in the literature (Der and Lee, 2005; Peluchetti et al., 2020), the lack of a covariance function precludes the inverse problem of feasible posterior inference and prediction, which we overcome using a conditionally Gaussian representation. There is a wealth of literature on the universal approximation property of both shallow and deep neural networks, following the pioneering work of Hornik et al. (1989), but they work under the assumption of a _continuous_ true function. Our numerical results demonstrate that when the truth has jump discontinuities, it is possible to obtain much better results with a BNN using weights with unbounded prior variance. The fully Bayesian posterior also allows straightforward probabilistic uncertainty quantification for the infinite width scaling limit under \(\alpha\)-stable priors on network weights.

Several future directions could naturally follow from the current work. The most immediate is perhaps an extension to posterior inference for deep networks under stable priors, where the width of each layer simultaneously approaches infinity, and we strongly suspect this should be possible. Developing analogous results under non-i.i.d. or tied weights to perform posterior inference under the scaling limits for Bayesian convolutional neural networks should also be of interest. Finally, one may of course investigate alternative activation functions, such as the hyperbolic tangent, which will lead to a different characteristic function for the scaling limit.

## Supplementary Material

The supplementary material contains technical details, numerical results and computer code to reproduce the simulation results on Github at: [https://github.com/loria/alphastableNNet](https://github.com/loria/alphastableNNet).

## Appendix A Appendix

### Proof of Theorem 1

Our derivations rely on Equation 5.4.6 of Uchaikin and Zolotarev (1999), which states that for \(\alpha_{0}\in(0,1)\) and for all positive \(\lambda\) one has: \(\exp(-\lambda^{\alpha_{0}})=\int_{0}^{\infty}\exp(-\lambda t)p_{S^{+}}(t)dt,\) where \(p_{S^{+}}\) is the density function of a positive \(\alpha_{0}\)-stable random variable. Using \(\lambda=\nu z^{2}\) and \(\alpha_{0}=\alpha/2\) and the fact that \(z^{2}=|z|^{2}\), we obtain for \(\alpha\in(0,2)\) that:

\[\exp(-\nu^{\alpha/2}|z|^{\alpha})=\int_{0}^{\infty}\exp(-\nu z^{2}t)p_{S^{+}}( t)dt, \tag{7}\]

\begin{table}
\begin{tabular}{l c c c c} \hline  & Stable & GP MLE & GP Bayes & Bayes NNet \\ \hline MAE & 0.57 & 0.59 & 0.58 & 0.59 \\ \hline \end{tabular}
\end{table}
Table 1: Mean absolute error of predictions in the York \(\mathrm{NO}_{2}\) air quality data set, by method.
where \(p_{S^{+}}\) is the density function of a positive \(\alpha/2\)-stable random variable and \(z\in\mathbb{R}\). By Equation (4) of Der and Lee (2005), and the characteristic function of independent normally distributed error terms with a common variance \(\sigma^{2}\), we have that:

\[\phi_{\mathbf{y}}(\mathbf{t}) =\exp\left(-\sum_{\ell=1}^{L}2^{-\alpha/2}\nu^{\alpha/2}q_{\ell} \left|\sum_{j=1}^{n}t_{j}\tau_{\ell}(\mathbf{x}_{j})\right|^{\alpha}-\frac{1} {2}\sum_{j=1}^{n}\sigma^{2}t_{j}^{2}\right)\] \[=\prod_{\ell=1}^{L}\exp\left(-2^{-\alpha/2}\nu^{\alpha/2}q_{\ell} \left|\sum_{j=1}^{n}t_{j}\tau_{\ell}(\mathbf{x}_{j})\right|^{\alpha}\right) \times\exp\left(-\frac{1}{2}\sum_{j=1}^{n}\sigma^{2}t_{j}^{2}\right)\] \[=\prod_{\ell=1}^{L}\int_{0}^{\infty}\exp\left(-\frac{1}{2}\nu s_{ \ell}q_{\ell}^{2/\alpha}\left(\sum_{j=1}^{n}t_{j}\tau_{\ell}(\mathbf{x}_{j}) \right)^{2}\right)p_{S^{+}}(s_{\ell})ds_{\ell}\times\exp\left(-\frac{1}{2} \sum_{j=1}^{n}\sigma^{2}t_{j}^{2}\right)\] \[=\prod_{\ell=1}^{L}\int_{0}^{\infty}\exp\left(-\frac{1}{2}\nu s_{ \ell}q_{\ell}^{2/\alpha}\mathbf{t}^{T}\mathbf{M}_{\ell}\mathbf{t}\right)p_{S^ {+}}(s_{\ell})ds_{\ell}\times\exp\left(-\frac{1}{2}\sum_{j=1}^{n}\sigma^{2}t_ {j}^{2}\right),\]

where the third line follows by using Equation (7), and in the last line we define \(\mathbf{M}_{\ell}\) as the matrix with ones in the diagonal and with the \((i,j)\)th entry given by \(\tau_{\ell}(\mathbf{x}_{i})\tau_{\ell}(\mathbf{x}_{j}),\ i\neq j\). Next, using the fact that the densities are over the independent variables \(\{s_{\ell}\}_{\ell=1}^{L}\) we bring the product inside the integrals and employ the property of the exponential to obtain:

\[\phi_{\mathbf{y}}(\mathbf{t}) =\int_{(\mathbb{R}^{+})^{L}}\exp\left(-\frac{1}{2}\sum_{j=1}^{n}t _{j}^{2}\sigma^{2}-\frac{1}{2}\nu\sum_{\ell=1}^{L}s_{\ell}q_{\ell}^{2/\alpha} \mathbf{t}^{T}\mathbf{M}_{\ell}\mathbf{t}\right)\prod_{\ell=1}^{L}p_{S^{+}}(s _{\ell})ds_{\ell}\] \[=\int_{(\mathbb{R}^{+})^{L}}\exp\left(-\frac{1}{2}\mathbf{t}^{T} \mathbf{Q}\mathbf{t}\right)\prod_{\ell=1}^{L}p_{S^{+}}(s_{\ell})ds_{\ell},\]

using on the second line the definition of \(\mathbf{Q}\). The required density is now obtained by the use of the inverse Fourier transform on the characteristic function:

\[p(\mathbf{y}\mid\mathbf{X}) =\int_{\mathbb{R}^{n}}\phi_{\mathbf{y}}(\mathbf{t})\exp(i( \mathbf{t},\mathbf{y}))\prod_{j=1}^{n}dt_{j}\] \[=\int_{\mathbb{R}^{n}}\left\{\int_{(\mathbb{R}^{+})^{L}}\exp \left(-\frac{1}{2}\mathbf{t}^{T}\mathbf{Q}\mathbf{t}\right)\prod_{\ell=1}^{L}p _{S^{+}}(s_{\ell})ds_{\ell}\right\}\exp(i(\mathbf{t},\mathbf{y}))\prod_{j=1}^{ n}dt_{j}\] \[=\int_{(\mathbb{R}^{+})^{L}}\left\{\int_{\mathbb{R}^{n}}\exp \left(-\frac{1}{2}\mathbf{t}^{T}\mathbf{Q}\mathbf{t}\right)\exp(i(\mathbf{t}, \mathbf{y}))\prod_{j=1}^{n}dt_{j}\right\}\prod_{\ell=1}^{L}p_{S^{+}}(s_{\ell}) ds_{\ell},\]

where the second line follows by the derived expression for the characteristic function, and the third line follows by Fubini's theorem since all the integrals are real and finite. We recognize that the term \(\exp(-(1/2)\mathbf{t}^{T}\mathbf{Q}\mathbf{t})\) corresponds to a multivariate Gaussian density with covariance matrix \(\mathbf{Q}^{-1}\), though it is lacking the usual determinant term. After completing it, we obtain the density and can use the characteristic function of Gaussian variables to finally obtain the result:

\[p(\mathbf{y}\mid\mathbf{X})=(2\pi)^{-n/2}\int_{(\mathbb{R}^{+})^{L}}\exp \left(-\frac{1}{2}\mathbf{y}^{T}\mathbf{Q}^{-1}\mathbf{y}\right)\det( \mathbf{Q})^{-1/2}\prod_{\ell=1}^{L}p_{S^{+}}(s_{\ell})ds_{\ell}.\]

In using \(\mathbf{Q}^{-1}\) freely, we assumed through the previous steps that \(\mathbf{Q}\) is positive-definite. We proceed to prove this fact. Note that \(\mathbf{Q}\) is obtained from the sum of \(L\) rank-one matrices and a diagonal matrix, where each of the rank-one matrices is \(q_{\ell}^{2/\alpha}s_{\ell}\nu\tau_{\ell}\boldsymbol{\tau}_{\ell}^{T}\), denoting by \(\boldsymbol{\tau}_{\ell}\) the vector with entries \((\tau_{\ell}(\mathbf{x}_{1}),\ldots,\tau_{\ell}(\mathbf{x}_{n}))\). Consider a non-zero vector \(\mathbf{w}\in\mathbb{R}^{n}\). Then:

\[\mathbf{w}^{T}\mathbf{Q}\mathbf{w}=\mathbf{w}^{T}\left(\sigma^{2}\mathbf{I}+ \nu\sum_{\ell=1}^{L}s_{\ell}q_{\ell}^{2/\alpha}\boldsymbol{\tau}_{\ell}^{T} \boldsymbol{\tau}_{\ell}\right)\mathbf{w}\]
\[=\sigma^{2}\mathbf{w}^{T}\mathbf{w}+\nu\sum_{\ell=1}^{L}s_{\ell}q_{ \ell}^{2/\alpha}\mathbf{w}^{T}\boldsymbol{\tau}_{\ell}\boldsymbol{\tau}_{\ell}^ {T}\mathbf{w}\] \[=\sigma^{2}\sum_{j=1}^{n}w_{j}^{2}+\nu\sum_{\ell=1}^{L}s_{\ell}q_ {\ell}^{2/\alpha}\left(\sum_{j=1}^{n}w_{j}\tau_{\ell}(x_{j})\right)^{2}\] \[>0,\]

which implies that \(\mathbf{Q}\) is positive-definite with probability 1.
## Acknowledgments

Bhadra is supported by U.S. National Science Foundation grant DMS-2014371.
