Deep Kernel Posterior Learning under Infinite Variance Prior Weights

## Abstract

Abstract Neal ( 1996 ) proved that infinitely wide shallow Bayesian neural networks (BNN) converge to Gaussian processes (GP), when the network weights have bounded prior variance. Cho & Saul ( 2009 ) provided a useful recursive formula for deep kernel processes for relating the covariance kernel of each layer to the layer immediately below. Moreover, they worked out the form of the layer-wise covariance kernel in an explicit manner for several common activation functions. Recent works, including Aitchison et al. ( 2021 ) , have highlighted that the covariance kernels obtained in this manner are deterministic and hence, precludes any possibility of representation learning, which amounts to learning a non-degenerate posterior of a random kernel given the data. To address this, they propose adding artificial noise to the kernel to retain stochasticity, and develop deep kernel inverse Wishart processes. Nonetheless, this artificial noise injection could be critiqued in that it would not naturally emerge in a classic BNN architecture under an infinite-width limit. To address this, we show that a Bayesian deep neural network, where each layer width approaches infinity, and all network weights are elliptically distributed with infinite variance, converges to a process with α 𝛼 \alpha italic_α -stable marginals in each layer that has a conditionally Gaussian representation. These conditional random covariance kernels could be recursively linked in the manner of Cho & Saul ( 2009 ) , even though marginally the process exhibits stable behavior, and hence covariances are not even necessarily defined. We also provide useful generalizations of the recent results of Loría & Bhadra ( 2024 ) on shallow networks to multi-layer networks, and remedy the computational burden of their approach. The computational and statistical benefits over competing approaches stand out in simulations and in demonstrations on benchmark data sets.

## 1 Introduction

The study of deep kernel processes as the infinite-width limit of deep Bayesian neural networks can be traced back to the work of Cho & Saul (), and further developed by Damianou & Lawrence (); Wilson et al. (); Aitchison et al. (), among others. This connection is built upon the foundational work of Neal (), who proved under the condition of bounded variance network weights, the infinite-width limit of a shallow (one hidden layer) Bayesian neural network (BNN) is a Gaussian process (GP), whose covariance kernel depends on the choice of the nonlinear activation function. Neal () studied two bounded activations explicitly: the sign function and $\tanh$, and worked out the covariance kernels in these cases: respectively the exponential and the squared exponential covariances. Subsequently, Williams () obtained explicit expressions for the kernel under an error function activation. Through a standard use of what is now known as the *kernel trick*, these expressions then allow one to work on with the covariance matrix induced by the kernel at the observation points for the purpose of nonparametric estimation and prediction under a GP model, rather than having to work with a possibly infinite-dimensional and nonlinear feature space.

The importance of the work of Cho & Saul () in the deep GP literature over the early works on shallow networks is that they extend the kernel processes to *deep* architectures, which emerge as infinite width limits of deep BNNs under bounded variance weights in various architectures, such as simple deep feedforward BNN (Lee et al., ; de G. Matthews et al., ), or deep convolutional BNNs (Garriga-Alonso et al., ). In fact, the tensor program of Yang () shows that a limiting GP would appear under an arbitrary network architecture under bounded variance priors. The work of Cho & Saul () provides an explicit recursive formula relating the kernels at each layer to the layer below, so that one can evaluate the covariance matrix of the features in layer $\ell+1$, using only the covariance matrix of the features in layer $\ell$. However, in the case of deep Gaussian processes, the covariance kernel in each layer is deterministic (non-stochastic) and the downside is that any possibility of learning the posterior distribution of the features disappears (Aitchison, ), as the posterior distribution of a degenerate point mass kernel is again a degenerate point mass.

To address the deterministic limit, Aitchison et al. () propose deep inverse Wishart processes (DIWP), which they argue approximates neural network GPs and have convenient properties for learning a variational posterior. Before introducing DIWP they also propose deep Wishart processes (DWP), which work with a noisy version of the sample covariance matrices for finite-width DNNs, the downside being a deterministic infinite-width limit. While these approaches succeed in making the covariance kernel random through an artificial noise injection, so that the posterior is non-degenerate, these processes do not result naturally as an infinite-width limit as in Neal () or Lee et al. (). Consequently, the choice of the artificial noise distribution could itself be considered a somewhat subjective design parameter, even though some recommendations and desiderata for such choices do exist, such as the deterministic limiting kernel is a *centrality parameter* of the distribution assigned to the random kernel, which is ensured in both DIWP and DWP.

Although deep GPs emerging as the infinite-width limit of deep BNNs are definitely interesting in their own right, there are some limitations of a Gaussian scaling limit. One is the lack of feature learning mentioned above. Second is the mean-square continuity of GPs, which means a GP prior cannot model jump discontinuities or even just localized smoothness. The third limitation is one noted by Neal (), that in the case of multiple outputs, GPs cannot learn the covariance structure, as the covariance is always zero under a Gaussian limit. To address these limitations we develop $\alpha$-stable kernel processes as the infinite-width limit of deep BNNs with infinite prior variance at each layer. The realizations of $\alpha$-stable random variables are able to model jumps (Kyprianou, , Section. 1.2), and their ability to learn features via a stochastic kernel is developed in this paper.

### 1.1 Prior works on infinite-width limits of BNNs under infinite variance priors and our contributions in context

The literature on the scaling limits of infinitely wide BNNs with weights that have infinite variance priors has traditionally focused on the properties of the limiting process and less on posterior inference under such prior processes (analogous to kriging under GP priors). The first non-Gaussian limit result of its kind was given by Der & Lee (), who proved that in a single hidden layer neural network the limiting process is sub-Gaussian with stable margins, and derived the characteristic function of the process. Recent extensions to deep feedforward networks are by Peluchetti et al. (), Favaro et al. () and Lee et al. (); and to the convolutional setting by Bracale et al. (). To the best of our knowledge, the first computationally viable method for posterior inference under these prior processes was explored by Loría & Bhadra (). However, their method has two major drawbacks. The first is that the method is limited to shallow networks. The second is the computational complexity that scales exponentially in the input dimension. Specifically, with $n$ training points of dimension $I$, the computational complexity of their method is $\mathcal{O}(n^{I+2})$. Although Loría & Bhadra () provided considerable evidence in one and two dimensions that a stable process prior is beneficial when the true underlying data-generating function contains jump-type discontinuities (which arise with zero probability under a GP prior), the second drawback, i.e., the exponential complexity in $I$, is especially limiting, since it makes it infeasible to apply the method in all but one or two dimensions, as also evidenced by our numerical experiments presented later in this paper. An intuitive explanation for this exponential complexity is that Loría & Bhadra () work in the *feature space* with sign activation function in the hidden layer, for which the resulting separating hyperplanes require enumeration over each of the $I$ dimensions on whether the output of the activation function is marked as 1 or 0. Enumerating over all $n$ points results in the exponential in $I$ complexity. Avoiding this combinatorial enumeration is of critical importance for large $I$.

To address these serious limitations, the current work eschews working in the feature space altogether, and focuses on kernel methods. Setting aside the issue of computational complexity for a moment, a key observation of Loría & Bhadra () is that the limiting stable process has a conditionally Gaussian representation, i.e., it is a Gaussian process, conditioned on some positive $\alpha$-stable variables. To see the main idea, we introduce the necessary notations on $\alpha$-stable random variables at this point. Define $X\sim S(\alpha,\beta)$ an $\alpha$-stable random variable with index parameter $\alpha\in(0,2]$ and skewness parameter $\beta\in[-1,1]$, by its characteristic function: $\phi_{X}(t)=\mathbb{E}[\exp(\mathrm{i}tX)]=\exp\{-\lvert t\rvert^{\alpha}[1-
\mathrm{i}\beta\omega(\alpha,t)]\}$, with $\omega(\alpha,t)=\tan(\alpha\pi/2)\mathrm{sign}(t)$ for $\alpha\neq 1$, and $\omega(1,t)=-(2/\pi)\log(\lvert t\rvert)\mathrm{sign}(t)$ (Samorodnitsky & Taqqu, , p. 5), with no closed form to the density of $X$ in general, apart from specific $\alpha$. A remarkable property of $\alpha$-stable variables is that the $\alpha=2$ case is Gaussian, with the usual Gaussian limit for the scaled sum established by the classical central limit theorem (CLT); but they possess infinite variance for $\alpha<2$, which does not admit a Gaussian scaling limit due to an inapplicability of the classical CLT. The generalized CLT (Gnedenko & Kolmogorov, ), that still applies, gives the non-Gaussian scaling limit in this case, and is stated in Appendix A. Two cases of interest are: symmetric $\alpha$-stable random variables, which have $\beta=0$; and positive $\alpha$-stable, which requires $\alpha<1$ and $\beta=1$ and which we denote by $S_{\alpha}^{+}$. Equation 5.4.6 of Uchaikin & Zolotarev () states for $\alpha_{0}\in(0,1)$ and for all positive $\lambda$ one has: ${\exp(-\lambda^{\alpha_{0}})=\int_{0}^{\infty}\exp(-\lambda s)p_{S^{+}_{\alpha
_{0}}}(s)ds},$ where $p_{S^{+}_{\alpha_{0}}}(\cdot)$ is the density function of a positive $\alpha_{0}$-stable random variable, with no closed form in general. Using $\lambda=t^{2},\alpha_{0}=\alpha/2$ and the fact that $t^{2}=\lvert t\rvert^{2}$, we obtain for $\alpha\in(0,2)$, a Gaussian scale mixture representation over a random scale parameter $s$ for the symmetric $\alpha$-stable characteristic function as:

$$ $\exp(-\lvert t\rvert^{\alpha})=\int_{0}^{\infty}\exp(-t^{2}s)p_{S^{+}_{\alpha/ 2}}(s)ds.$ (1) $$

This observation is exploited in Theorem 1 of Loría & Bhadra (), to represent the posterior of the stable process as a Gaussian process posterior, conditional on mixing positive $\alpha/2$-stable variables, which results in posterior sampling analogous to the GP case. However, in computing the covariance kernel of this GP, they suffer from the aforementioned combinatorial bottleneck.

Elliptical $\alpha$-stable random vectors (Samorodnitsky & Taqqu, , Definition 2.5.1) centered at zero have a characteristic function defined similarly to the scalar case as: $\phi_{\mathbf{Z}}(\mathbf{t})=\mathbb{E}[\exp(\mathrm{i}\mathbf{t}^{T}\mathbf{
Z})]=\exp\{-(\mathbf{t}^{T}\Sigma\mathbf{t})^{\alpha/2}\}$, where $\Sigma$ is a positive definite matrix, termed the shape parameter. The main advantage of zero-centered elliptical $\alpha$-stable vectors is that they admit the Gaussian mixture representation: $\mathbf{Z}\overset{d}{=}S^{1/2}\mathbf{G}$, where $\overset{d}{=}$ denotes equality in distribution, $S\sim S^{+}_{\alpha/2}$ and $\mathbf{G}\sim\mathcal{N}(0,\Sigma)$, with $S$, $\mathbf{G}$ independent, which can be shown using a multivariate extension of Equation (1).

What the current work does is to find a recursive formula for deep kernel processes exploiting the technique of Cho & Saul () to relate the conditional *stochastic covariance kernels* in a recursive manner, resulting in a far more computationally practical procedure in the kernel space rather than one in the feature space. One must note that although the limit processes are conditionally Gaussian for each layer, they possess stable margins. Nevertheless, the conditionally Gaussian representation naturally results in a stochastic covariance kernel in each layer, whose posterior can be learned given the data, to enable data-dependent representation learning. In this way, the proposed method bypasses the need for artificial noise injection in the manner of Aitchison et al. ().

### 1.2 Summary of main contributions

- 1.
The development of a novel deep $\alpha$-stable kernel process, arising as the infinite-width limit of a deep BNN under elliptical infinite variance priors on the weights at each layer. We further present a *conditionally Gaussian* representation of the resulting process at each layer, with a recursive formula relating these conditional covariance kernels, even when covariances for the marginal processes do not exist at any layer.
- 2.
Theoretical and numerical demonstrations that the covariance kernel in the conditionally Gaussian representation for each layer is stochastic, allowing for learning a non-degenerate posterior of the kernels and identifying a clear demarcation with the deep GP literature.
- 3.
A Markov chain Monte Carlo method for posterior sampling, which allows out-of-sample prediction at new inputs as well as uncertainty quantification of the predictions via the full posterior predictive distribution.
- 4.
Numerical demonstrations in simulations and on benchmark UCI data sets that our method performs better in prediction than the competing methods, offers full predictive uncertainty estimates, and is considerably less time-consuming than the method of Loría & Bhadra () that is limited to shallow (one hidden layer) BNNs.

## 2 Deep α 𝛼 \alpha italic_α -stable kernel processes as infinite-width limits of deep Bayesian neural networks under heavy-tailed weights

Consider a deterministic input $\mathbf{x}\in\mathbb{R}^{I}$ with a corresponding one-dimensional output $\psi(\mathbf{x})$. We define an $L$ layer feedforward deep neural network (DNN) with $L-1$ hidden layers by the recursion:

$$ $\displaystyle f_{j}^{(\ell+1)}(\mathbf{x})$ $\displaystyle=g\left(b_{j}^{(\ell)}+\sum_{i=1}^{M_{\ell}}w_{ij}^{(\ell)}f_{i}^ {(\ell)}(\mathbf{x})\right),\quad j=1,\dots,M_{\ell+1};\;\ell=1,\ldots,L-1,$ (2) $\displaystyle\psi(\mathbf{x})$ $\displaystyle=\sum_{j=1}^{M_{L}}w_{j}^{(L)}f_{j}^{(L)}(\mathbf{x}),$ (3) $$

where $f^{(1)}\equiv\mathbf{x},\;g(\cdot)$ is a nonlinear activation function, and $M_{\ell}$ is the width of the $\ell$-th layer. Neal () proved under mild conditions, with $L=2$ and the weights $w_{j}^{(2)}\stackrel{{\scriptstyle ind}}{{\sim}}M_{2}^{-1/2}\mathcal{N}(0,1)$, by taking the infinite-width limit $M_{2}\to\infty$ for the last layer of this shallow network, the resulting stochastic process is a Gaussian process. Cho & Saul () provided a recursive formula for relating the layer-wise covariance kernels when the activation function $g(\cdot)$ is of the form $g_{\delta}(\zeta)=\zeta^{\delta}\mathbf{1}_{\{\zeta>0\}}$, for $\delta$ a non-negative integer, and for any number of hidden layers while taking the limits $M_{1}\to\infty,\dots,M_{L}\to\infty$.
Unlike these works, we use *infinite variance* prior weights at each layer and study the resulting infinite-width limit posterior. For the first layer we take the weights $w_{mj}^{(1)}\overset{indep.}{\sim}\mathcal{N}(0,s^{(1)}_{m,+})$, where $s^{(1)}_{m,+}\overset{i.i.d.}{\sim}S_{\alpha/2}^{+}$, and the bias of the first layer is simply taken to be $b_{j}^{(1)}\overset{i.i.d.}{\sim}\mathcal{N}(0,1)$. For the subsequent layers, we set $b_{j}^{(\ell)}$ as i.i.d. $\mathcal{N}(0,1)$ and the weights ($w_{ij}^{(\ell)}$, for $i=1,\dots,M_{\ell}$, and $j=1,\dots,M_{\ell+1}$) of each layer $\ell>1$ are given a distribution with infinite variance, as specified below in Theorem 1. This theorem establishes the resulting infinite-width limit as a non-Gaussian process with stable marginals. Further, it provides a *conditionally Gaussian* representation, conditioned on mixing positive $\alpha/2$ stable variables at each layer. The covariance kernel of this conditional GP can be derived explicitly under several commonly used activation functions, similarly to Cho & Saul (). Once this connection is established, posterior inference and prediction becomes feasible even in this infinite-variance regime.

###### Theorem 1 .

For $j=1,\ldots,M_{\ell+1}$ and $k=1,\ldots,n$, define: $z_{j}^{(\ell)}(\mathbf{x}_{k})=\frac{1}{M_{\ell}^{1/2}}\sum_{i=1}^{M_{\ell}}w^
{(\ell)}_{ij}f_{i}^{(\ell)}(\mathbf{x}_{k})$, where $w^{(\ell)}_{ij}=(s_{+}^{(\ell)})^{1/2}{\tilde{w}}^{(\ell)}_{ij}$, with $s_{+}^{(\ell)}$ i.i.d. positive $\alpha/2$ stable random variables for $\ell=2,\ldots,L$, and without loss of generality, ${\tilde{w}}_{ij}^{(\ell)}$ are i.i.d. with zero mean and unit variance and independent of $s_{+}^{(\ell)}$. Then, the weights $w^{(\ell)}_{ij}$ have infinite variance and as $M_{\ell}\to\infty$, the limiting distribution of $\mathbf{z}_{j}^{(\ell)}=(z_{j}^{(\ell)}(\mathbf{x}_{1}),\dots,z_{j}^{(\ell)}(
\mathbf{x}_{n})),$ in the $\ell$-th layer is elliptical $\alpha$-stable, with characteristic function:

$$ $\displaystyle\phi_{\mathbf{z}_{j}^{(\ell)}\mid\Sigma^{(\ell)}}(\mathbf{t})$ $\displaystyle=\exp\left\{-(\mathbf{t}^{T}\Sigma^{(\ell)}\mathbf{t})^{\alpha/2} \right\},$ $$

where $\mathbf{t}=(t_{1},\ldots,t_{n})$ and the $n\times n$ random matrix $\Sigma^{(\ell)}$ is positive definite with probability 1. As such, the limiting $\mathbf{z}^{(\ell)}_{j}$ admits the representation:
$\mathbf{z}^{(\ell)}_{j}\mid s_{+}^{(\ell)},\Sigma^{(\ell)}\sim\mathcal{N}(0,s_
{+}^{(\ell)}\Sigma^{(\ell)}),$ where $s^{(\ell)}_{+}\sim S_{\alpha/2}^{+}$.
For activation functions $g_{\delta}(\zeta)=\zeta^{\delta}\mathbf{1}_{\{\zeta>0\}}$, and ${\tilde{w}}_{ij}^{(\ell)}$ i.i.d. standard Gaussian, the matrix $\Sigma^{(\ell)}$ is recursively given as a function of $\Sigma^{(\ell-1)}$ and $s^{(\ell-1)}_{+}$ by:

$$ $\displaystyle\Sigma_{k,h}^{(\ell)}$ $\displaystyle=\pi^{-1}\left[\left(1+s^{(\ell-1)}_{+}\Sigma_{k,k}^{(\ell-1)} \right)\left(1+s^{(\ell-1)}_{+}\Sigma_{h,h}^{(\ell-1)}\right)\right]^{\delta/2 }J_{\delta}(\theta_{k,h}^{(\ell)}),$ $\displaystyle\theta_{k,h}^{(\ell)}$ $\displaystyle=\cos^{-1}\left\{\left[1+s^{(\ell-1)}_{+}\Sigma_{k,h}^{(\ell-1)} \right]\left[1+s^{(\ell-1)}_{+}\Sigma_{k,k}^{(\ell-1)}\right]^{-1/2}\left[1+s^ {(\ell-1)}_{+}\Sigma_{h,h}^{(\ell-1)}\right]^{-1/2}\right\},$ $$

for $k,h=1,\dots,n$, $\ell=2,\dots,L$; where $J_{\delta}(\theta)$ can be computed explicitly for $\delta\in\mathbb{N}$, and $\Sigma^{(1)}=\mathbf{X}\mathbf{S}^{(1)}_{+}\mathbf{X}^{T}$, where $\mathbf{X}=[\mathbf{x}_{1},\dots,\mathbf{x}_{n}]^{T},\;\mathbf{S}^{(1)}_{+}=
\mathrm{diag}(s^{(1)}_{1,+},\dots,s^{(1)}_{I,+})$, with $s^{(1)}_{m,+}\overset{i.i.d.}{\sim}S^{+}_{\alpha/2}$ for $m=1,\dots,I$.

A proof is in Appendix B. For completeness, we note that $g_{0}(\zeta)$ and $g_{1}(\zeta)$ correspond to respectively the step function and ReLU activations, for which $J_{0}(\theta)=\pi-\theta$ and $J_{1}(\theta)=\sin(\theta)+(\pi-\theta)\cos(\theta)$, as derived by Cho & Saul ().
This result has similarities to those previously obtained by both Lee et al. () and de G. Matthews et al. (). However, we emphasize the key difference that in our case the kernels $s_{+}^{(\ell)}\Sigma^{(\ell)}$ at all layers are *random* for $\alpha<2$, conditional on scales with a positive $\alpha/2$-stable distribution, although they are positive definite with probability one. The limiting kernel for the Gaussian $(\alpha\to 2)$ case reduces to a deterministic one, as in earlier works, since $S^{+}_{\alpha/2\to 1}$ converges to a degenerate point mass at 1 with probability one. This permits an explicit characterization of a deep kernel process with marginally infinite variance at each layer via a conditional deep Gaussian process. Two significant gains, when compared to Loría & Bhadra () are obtained: (1) there is no need for exponential complexity computations in the kernel space as opposed to in the feature space, and (2) our derivations work with multi-layer BNNs, in contrast to the results of Loría & Bhadra () that are limited to BNNs with a *single* hidden layer.

The next proposition describes the posterior predictive density implied by the probabilistic structure in Theorem 1. The posterior predictive density we present combines the simplicity of predicting in GPs while ameliorating the challenges posed by variables that have infinite variance.

###### Proposition 2 .

Consider a vector $\mathbf{y}$ of $n$ real-valued observations with corresponding inputs $\mathbf{X}=[\mathbf{x}_{1},\dots,\mathbf{x}_{n}]^{T}$, where $\mathbf{x}_{k}\in\mathbb{R}^{I}$ for all $n$. Consider the model of Equations (2) and (3) and suppose we observe data with an additive error: $y_{k}=\psi(\mathbf{x}_{k})+\varepsilon_{k}$ where $\varepsilon_{k}\overset{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2})$. Then, under the setting of Theorem 1, and taking $M_{\ell}\to\infty$ for all $\ell$, at $m$ new inputs $\mathbf{X}^{*}=[\mathbf{x}_{n+1},\dots,\mathbf{x}_{n+m}]^{T}$, the posterior predictive distribution of the observations $\mathbf{y}^{*}$ is given by:

$$ $\displaystyle\mathbf{y}^{*}$ $\displaystyle\mid\mathbf{y},\mathbf{X},\mathbf{X}^{*},\{s^{(\ell)}_{+}\}_{\ell =2}^{L},\mathbf{S}^{(1)}_{+}\sim\mathcal{N}_{m}(\bm{\mu}^{*},\bm{\Lambda}^{*}) \text{, where, }$ $\displaystyle\bm{\mu}^{*}$ $\displaystyle=\bm{\Lambda}_{(n+1):(n+m),1:n}\bm{\Lambda}_{1:n,1:n}^{-1}\mathbf {y},$ $\displaystyle\bm{\Lambda}^{*}$ $\displaystyle=\bm{\Lambda}_{(n+1):(n+m),(n+1):(n+m)}-\bm{\Lambda}_{(n+1):(n+m) ,1:n}\bm{\Lambda}_{1:n,1:n}^{-1}\bm{\Lambda}_{1:n,(n+1):(n+m)},$ $\displaystyle s^{(\ell)}_{+}$ $\displaystyle\overset{i.i.d.}{\sim}S_{\alpha/2}^{+};s^{(1)}_{m,+}\overset{i.i. d.}{\sim}S_{\alpha/2}^{+},\text{ for, }\ell=2,\dots,L;\;m=1,\dots,I,$ $$

and $\bm{\Lambda}=\Sigma^{(L)}+\sigma^{2}\mathbf{I}_{n+m}$, where $\Sigma^{(L)}$ is the matrix obtained in the last layer in Theorem 1 with the inputs, in order, $[\mathbf{X},\mathbf{X}^{*}]$ and $\mathbf{I}_{n+m}$ denotes the identity matrix of size $(n+m)$.

Proof of Proposition 2 is in Appendix C. This proposition enables computation of the kernel posterior and prediction, and is key to the numerical demonstrations presented later. We also remark that the purpose of Proposition 2 is to tackle regression tasks rather than classification. However, similar results can also be obtained for classification.
Note that for $\alpha<2$, the conditional kernel $\bm{\Lambda}^{*}$ is a stochastic matrix, which we verify numerically in Sections 3 and 4.

Making use of Proposition 2, the next algorithm provides a Markov chain Monte Carlo (MCMC) procedure for sampling from the posterior predictive density $p(\mathbf{y}^{*},\{s^{(\ell)}_{+}\}_{\ell=2}^{L},\mathbf{S}^{(1)}_{+}\mid
\mathbf{y},\mathbf{X},\mathbf{X}^{*})$. Algorithm 1 obtains a valid sample of the posterior density $p(\mathbf{y}^{*}\mid\mathbf{y},\mathbf{X},\mathbf{X}^{*})$ by classic results in the MCMC literature (Chib & Greenberg, ). Our algorithm also provides valid samples of the posterior stochastic covariance matrix induced by the kernel at the design points. We remark that simulating in previously unseen locations can be done in an off-line manner, by simply utilizing the scales and the collection of input points, which permits an efficient prediction algorithm without a need to resample the scales. We also make use of a half-Cauchy prior (Gelman, ) on the variance of the errors $\sigma^{2}$, to specify a fully-Bayesian model.

Figure: Algorithm 1 A Metropolis–Hastings sampler for the posterior predictive distribution of the deep $\alpha$-kernel process.

Following the construction of Theorem 1 ensures we have defined a valid kernel process that is consistent under marginalization, in the manner described by Aitchison et al. (). Specifically, the process $\{\Sigma^{(\ell)}\}_{\ell=1}^{L}$ is a deep kernel process with probability 1, meaning that they induce a distribution on positive definite matrices of different sizes which are consistent under marginalization. The proof of this follows from a similar argument as in Aitchison et al. (), and reduces to their result under fixed scales. For completeness we provide further elaboration in Appendix E.

Additionally, the process $\{\Sigma^{(\ell)}\}_{\ell=1}^{L}$ has an implicit dependence on $\alpha$, the index of the $\alpha$-stable random variables. For $\alpha=2$, the characteristic function of $\mathbf{z}_{j}^{(\ell)}$ becomes that of a Gaussian random variable, and $\Sigma^{(\ell)}$ is deterministic. Remarkably this kernel process is stochastic when $\alpha<2$. The importance of having a stochastic deep kernel process is that it allows us to learn features, which deep Gaussian processes are not able to do, as highlighted by Yang et al. (). We formalize this notion in the following proposition, with proof in Appendix F.

###### Proposition 3 .

For $\alpha<2$ the posterior distribution of the features $\mathbf{z}_{j}^{(\ell)}$ is given by

$$ $\displaystyle p(\mathbf{z}_{j}^{(\ell)}\mid\mathbf{X},\mathbf{y})=\int$ $\displaystyle p(\mathbf{z}_{j}^{(\ell)}\mid\{s^{(\ell)}_{+}\}_{\ell=2}^{L}, \mathbf{S}^{(1)}_{+})p(\{s^{(\ell)}_{+}\}_{\ell=2}^{L},\mathbf{S}^{(1)}_{+} \mid\mathbf{X},\mathbf{y})\prod_{\ell=2}^{L}ds^{(\ell)}_{+}\prod_{m=1}^{I}ds^{ (1)}_{m,+}.$ $$

This implies that the features also depend on the observations. For $\alpha=2$, the posterior distribution of the features are independent of the observations, and hence, cannot be learned from the data.

The importance of this proposition is that we are able to verify that in the deep kernel process we have defined, the features *do* depend on the observations. As such, there is feature learning only when $\alpha<2$, but when $\alpha=2$ we return to a deterministic kernel and are not able to learn the features. As an additional benefit, the proposition provides a straightforward procedure for sampling from the posterior of the features. Namely, we can sample the scales from $p(\{s_{+}^{(\ell)}\}_{\ell=2}^{L},\mathbf{S}^{(1)}_{+}\mid\mathbf{X},\mathbf{y})$, and then sample from $\mathbf{z}_{j}^{(\ell)}\mid\{s^{(\ell)}_{+}\}_{\ell=2}^{L},\mathbf{S}^{(1)}_{+
}\sim\mathcal{N}(0,s^{(\ell)}_{+}\Sigma^{(\ell)})$ using the posterior samples of $s_{+}^{(\ell)}$.

###### Theorem 1 .

###### Proposition 2 .

###### Proposition 3 .

## 3 Numerical experiments

### 3.1 Conditional mutual information as a function of α 𝛼 \alpha italic_α

Figure 1 displays the conditional mutual information for different values of $\alpha$ at varying distances over a uniform grid in one-dimension with $L=2,\delta=1$, for the proposed deep $\alpha$-kernel process, i.e., a process with one hidden layer. For a GP, the mutual information is simply a function of the correlation, but this quantity is not well defined in our case, since there is no well-defined covariance. As such, we need to refer to the conditional mutual information (Cover & Thomas, , p. 23), which can be computed using the conditionally-Gaussian representation and numerical integration via Monte Carlo methods. The expressions of conditional mutual information and details of computation are provided in Appendix G. The key finding is that the conditional mutual information decays at a slower rate for smaller $\alpha$. The long memory behavior of the conditional mutual information suggests the developed deep $\alpha$-kernel processes to be especially adept at picking up distant relationships, which would be modeled as nearly independent under a GP.

Figure: Figure 1: Decay of the mutual information for the deep $\alpha$-kernel process as a function of the distance between the inputs with $L=2,\delta=1$. The limiting Gaussian case ($\alpha=2$) is also included.
Refer to caption: x1.png

### 3.2 Posterior prediction and uncertainty quantification under discontinuity

We display results on functions with jumps in one, two, and ten dimensions. The reason we focus specifically on discontinuous true functions is that they are misspecified under a GP prior under *any* covariance function, which can model varying levels of smoothness, but not a lack of mean square continuity. In contrast, heavy-tailed priors enable posterior inference on a broader function class, e.g., those belonging to Besov space; see recent theoretical works by Agapiou & Castillo () and extensive numerical works in Loría & Bhadra (). The comparisons include a Bayesian GP method (Gramacy & Taddy, ) and a frequentist approach on GPs (Dancik & Dorman, ), as well as the DIWP and NNGP method as implemented by Aitchison et al. (), and the Stable method by Loría & Bhadra (). We refer to our proposed method as D$\alpha$-KP as a shorthand for deep $\alpha$-stable kernel process. For the deep kernel methods of Aitchison et al. () we train all the models using the same hyperparameters they use, with 8000 total steps with $10^{-2}$ as the step-size for the first half, and $10^{-3}$ for the second half. For the D$\alpha$-KP method we run 3000 MCMC simulations with $\alpha=1,\delta=1$. For a fair comparison, the three kernel methods all use $L=2$, meaning they are shallow with one hidden layer. The following simulation settings are considered, respectively in one, two, and ten dimensions, all under discontinuous truth.

- 1.
The true function in one dimension is $f(\xi)=5\times\mathbf{1}_{\{\xi>0\}}$ and we generate observations as $y(\xi)=f(\xi)+\varepsilon;\;\varepsilon\sim\mathcal{N}(0,0.5^{2})$. For training we consider 40 equally-spaced input points on $[-1,1]$ and predict on 100 out-of-sample points.
- 2.
In two dimensions we use the function $f(\xi_{1},\xi_{2})=5\times\mathbf{1}_{\{\xi_{1}>0\}}+5\times\mathbf{1}_{\{\xi_
{1}>0\}}$ and generate $y(\xi_{1},\xi_{2})=f(\xi_{1},\xi_{2})+\varepsilon;\;\varepsilon\sim\mathcal{N}
(0,0.5^{2})$, using a $7\times 7$ uniform grid on $[-1,1]^{2}$ for training, and a similar $9\times 9$ grid for testing.
- 3.
In the ten-dimensional setting we use the function $f(\bm{\xi})=6\mathrm{sign}(\xi_{1})+8\mathrm{sign}(\xi_{2}+\xi_{3})+6\mathrm{
sign}(\xi_{4}+\xi_{5})+6\mathrm{sign}(\xi_{6}+\xi_{7})+6\mathrm{sign}(\xi_{8}+
\xi_{9})+6\mathrm{sign}(\xi_{10})$ and generate $y(\bm{\xi})=f(\bm{\xi})+\varepsilon$, with $\varepsilon\sim\mathcal{N}(0,0.5^{2})$. The design points are generated by $\xi_{i}\overset{i.i.d.}{\sim}\mathrm{Unif}(-0.5,0.5)$. We generate 20 splits, each with 300 training and 300 testing observations.

Under these simulation settings, we provide in Table 1 the prediction errors for the competing methods by presenting the root mean squared error (RMSE) and mean absolute error (MAE).
Our method obtains smaller out-of-sample prediction errors compared to the other methods, and is generally a close second to the Stable method of Loría & Bhadra (). However, the running time of our method is a fraction of the time of Loría & Bhadra (), which we highlight in Supplementary Section H.1. As a result of this high run time owing to the exponential complexity computations, the results for the Stable method are only available in one and two dimensions. Figure 2(a) compares the function fit for the one-dimensional function. Figure 2(b) displays the uncertainty quantification (UQ) by presenting 90% posterior predictive intervals, where the benefits of the methods using $\alpha$-stable priors over GP priors are clear under discontinuity, as is the inability of GP to model such functions. We display similar figures for the examples in two and ten dimensions in Appendix H. Furthermore, the out-of-sample prediction errors for this case demonstrate the superior performance of methods that incorporate $\alpha$-stable random variables. Appendix H.2 shows the percent of observations that are covered in the 90% posterior predictive intervals for the simulation settings under all methods. In Appendix H.3 we display graphical summaries for the two and ten dimensional settings.

Figure: (a) Function fit for the different methods.
Refer to caption: x2.png

**Table 1: Out-of-sample errors in numerical examples, in twenty different splits. Best in bold. Stable not available for more than 2 dimensions.**
|  | One Dimension | Two Dimensions | Ten Dimensions |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Method | RMSE (SD) | MAE (SD) | RMSE (SD) | MAE (SD) | RMSE (SD) | MAE (SD) |
| D$\alpha$-KP | 0.57 (0.05) | 0.45 (0.04) | 0.86 (0.09) | 0.67 (0.07) | 8.08 (0.38) | 6.48 (0.33) |
| DIWP | 1.08 (0.04) | 0.84 (0.03) | 1.69 (0.05) | 1.36 (0.05) | 8.92 (0.38) | 7.21 (0.32) |
| GP Bayes | 0.69 (0.05) | 0.52 (0.04) | 0.90 (0.06) | 0.70 (0.06) | 10.39 (0.63) | 8.32 (0.52) |
| GP MLE | 0.77 (0.06) | 0.57 (0.04) | 1.19 (0.08) | 0.92 (0.07) | 8.32 (0.38) | 6.68 (0.34) |
| NNGP | 1.08 (0.04) | 0.84 (0.03) | 1.69 (0.05) | 1.36 (0.04) | 8.92 (0.39) | 7.21 (0.34) |
| Stable | 0.52 (0.03) | 0.42 (0.03) | 0.57 (0.08) | 0.45 (0.04) | – | – |

### 3.3 The effect of depth on prediction for α 𝛼 \alpha italic_α -kernel processes

We numerically investigate the effect of depth on predictive performance in the synthetic data sets of Section 3.2, with $L=2$ results given there. As can be seen from Table 2, there is no meaningful difference between deep and shallow $\alpha$ kernel processes for these simulation settings. Recent works have theoretically investigated fundamental barriers to *trainable depth* of deep neural networks, exceeding which actually results in degraded performance (Schoenholz et al., ). A similar theoretical investigation is beyond the scope of the current work, but we conjecture that the resultant stable process can capture a function class *rich enough* with even one hidden layer, so that the effect of depth becomes a secondary issue for these test functions. This is a very different situation from a deep NNGP, which can only capture functions that are mean square continuous, with the effect of depth primarily manifesting in a richer covariance kernel, but not in a non-Gaussian process.

**Table 2: Average (SD) out-of-sample errors for deep $\alpha$-KP, with $\delta=1,\alpha=1$, in the three simulation scenarios with varying layer depth, over twenty different splits.**
|  | One Dimension | Two Dimensions | Ten Dimensions |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| L | RMSE (SD) | MAE (SD) | RMSE (SD) | MAE (SD) | RMSE (SD) | MAE (SD) |
| 3 | 0.56 (0.05) | 0.45 (0.04) | 0.82 (0.07) | 0.65 (0.05) | 8.09 (0.38) | 6.49 (0.34) |
| 6 | 0.57 (0.05) | 0.45 (0.04) | 0.81 (0.07) | 0.64 (0.05) | 8.10 (0.39) | 6.49 (0.34) |
| 11 | 0.58 (0.05) | 0.45 (0.04) | 0.85 (0.06) | 0.66 (0.05) | 8.11 (0.40) | 6.51 (0.35) |
| 16 | 0.58 (0.05) | 0.46 (0.04) | 0.88 (0.06) | 0.69 (0.05) | 8.11 (0.39) | 6.50 (0.34) |

## 4 Applications to UCI data sets

### 4.1 Out-of-sample predictive performance

We apply the five methods (D$\alpha$-KP, DIWP, NNGP, GP Bayes, and GP MLE) to the three well-known data sets from the UCI repository: Boston, Yacht, and Energy. Due to the size of these data sets and the $\mathcal{O}(n^{I+2})$ computational complexity, the Stable procedure was not feasible to run. To this end, we split each of the data sets in 20 different folds, training in 19 and testing in the remaining fold; and repeat the process for each of the folds. We display the average RMSE and average MAE in Table 3, with the respective standard deviations. We find that in two out of three cases (Energy and Yacht), $\alpha$-KP has the best predictive performance, while being a close second on the Boston data set. In Appendix I.1 we display the running times for the considered methods, and in Appendix I.2 we display the coverage of the credible intervals.

**Table 3: Comparison of out-of-sample errors in twenty splits of the three UCI datasets, using the D$\alpha$-KP ($\alpha=1,\delta=1,L=2$), DIWP, NNGP (Aitchison et al., ), GP Bayes, and GP MLE. Number of observations denoted by $n$ and number of inputs by $I$. Stable method not available for $I>2$. Best in bold.**
|  | Boston ($n=506,I=13$) | Energy ($n=769,I=8$) | Yacht ($n=308,I=6$) |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Method | RMSE (SD) | MAE (SD) | RMSE (SD) | MAE (SD) | RMSE (SD) | MAE (SD) |
| D$\alpha$-KP | 2.59 (0.73) | 1.78 (0.35) | 0.46 (0.07) | 0.32 (0.05) | 0.31 (0.12) | 0.16 (0.05) |
| DIWP | 2.85 (0.89) | 2.01 (0.41) | 0.48 (0.06) | 0.34 (0.04) | 0.60 (0.20) | 0.30 (0.10) |
| NNGP | 3.00 (0.87) | 2.04 (0.43) | 2.18 (0.23) | 1.57 (0.18) | 3.88 (0.87) | 2.58 (0.49) |
| GP Bayes | 2.58 (0.75) | 1.76 (0.35) | 0.68 (0.06) | 0.51 (0.04) | 0.48 (0.23) | 0.22 (0.07) |
| GP MLE | 3.93 (1.02) | 2.58 (0.51) | 0.49 (0.06) | 0.34 (0.04) | 0.52 (0.34) | 0.25 (0.11) |
| Stable | – | – | – | – | – | – |

### 4.2 Non-Gaussian feature learning

In this section we numerically investigate whether D$\alpha$-KP indeed learns a non-degenerate and non-Gaussian posterior of the features. Specifically, we display in Figure 3 the first two features in the hidden layer of a two layer D$\alpha$-KP for the Boston data set. The posterior of features in the $\ell$-th layer is given by: $\mathbf{z}^{(\ell)}_{j}\mid\Sigma^{(\ell)},\mathbf{y}\sim\mathcal{N}(0,\Sigma^
{(\ell)})$. Note that after marginalizing over $\Sigma^{(\ell)}$, which has stable dependence, Proposition 3 indicates we would obtain a non-Gaussian distribution. Figure 3 confirms the non-Gaussianity. In this figure we include the 2-d heat map of posterior samples with 1-d marginal box plots, which show far more values in the tails than what can be expected under Gaussianity. The normal q–q plots also confirm the heavy tails.
We include similar figures for the Yacht and Energy datasets in Appendix J. The non-Gaussianity of $\mathbf{z}_{j}^{(\ell)}$ also serves as confirmation that $\Sigma^{(\ell)}$ is indeed stochastic. Since otherwise we would obtain Gaussian features.

Figure: (a) Scatterplot of the first two features, with box plots for the marginals.
Refer to caption: x4.png

## 5 Concluding remarks

We present a viable posterior inference procedure for deep $\alpha$-stable kernel processes arising as the infinite-width limit of a BNN under infinite variance prior weights at each layer. Through a convenient representation of elliptical $\alpha$-stable random vectors as a Gaussian mixture with respect to positive $\alpha/2$ stable variables acting as random scales, we obtain an explicit relationship with a deep conditionally Gaussian process, upon which the result of Cho & Saul () can be leveraged to link the conditional covariance kernels in a recursive manner, resulting in a computationally viable method. We also present extensive evidence of benefits in prediction and uncertainty quantification under a true data generating function that contains jump-type discontinuities. Compared to the approach of Loría & Bhadra (), the current work overcomes the serious limitations of being restricted to one hidden layer and exponential complexity computations.

Several future directions could naturally follow from the current work. Although MCMC is gold standard for full posterior uncertainty quantification, for computational scalability, one may look for a variational approximation to the posterior predictive density we derived. Applying techniques such as inducing points (Snelson & Ghahramani, ), would also help scale the algorithm to data sets with larger sample sizes.

Our results indicate successful feature learning under non-Gaussian stochastic processes that is not possible under a GP (shallow or deep) and achieves this without any artificial noise injection in the manner of Aitchison et al. (). In the current work, we have not made any efforts on feature selection, which appears to be a promising direction under the infinite-width limits of suitable sparsity inducing architectures, such as dropout (Srivastava et al., ), that should now be possible to study following our work and leveraging its main ideas.

## Appendix A The generalized central limit theorem

###### Theorem 4 .

(Uchaikin & Zolotarev, , p. 62) Let $X_{1},\dots,X_{n}$ be independent and identically distributed random variables with the distribution function $F_{X}(x)$ satisfying the conditions:

$$ $\begin{split}1-F_{X}(x)\sim cx^{-\mu},x\to\infty,\\ F_{X}(x)\sim dx^{-\mu},x\to-\infty,\end{split}$ $$

with $\mu>0$, $c\geq 0,d\geq 0$. Then, there exists sequences $a_{n}$ and $b_{n}$ such that the distribution of the centered and normalized sum $Z_{n}=b_{n}^{-1}(\sum_{i=1}^{n}X_{i}-a_{n})$, converges in distribution to the stable distribution with parameters:

$$ $\alpha=\min(\mu,2),\;\beta=\frac{c-d}{c+d},$ $$

as $n\to\infty$, i.e., $F_{Z_{n}}(x)\to G(x;\alpha,\beta)$ for all $x$ where $G$ is continuous, where $F_{Z_{n}}$ is the c.d.f of $Z_{n}$, and $G$ is the c.d.f. of an $\alpha$-stable random variable with symmetry parameter $\beta$. The coefficients $a_{n}$ and $b_{n}$ are given in Table 4.

**Table 4: Centering and normalizing coefficients $a_{n}$ and $b_{n}$ for the generalized central limit theorem.**
| $\mu$ | $\alpha$ | $a_{n}$ | $b_{n}$ |
| --- | --- | --- | --- |
| $0<\mu<1$ | $\mu$ | 0 | $\pi^{1/\alpha}(c+d)^{1/\alpha}[2\Gamma(\alpha)\sin(\alpha\pi/2)]^{-1/\alpha}n^<br>{1/\alpha}$ |
| $\mu=1$ | $\mu$ | $(c-d)n\log(n)$ | $\pi(c+d)n/2$ |
| $1<\mu<2$ | $\mu$ | $n\mathbb{E}[X]$ | $(c+d)^{1/\alpha}[2\Gamma(\alpha)\sin(\alpha\pi/2)]^{-1/\alpha}n^{1/\alpha}$ |
| $\mu=2$ | $2$ | $n\mathbb{E}[X]$ | $(c+d)^{1/2}[n\log(n)]^{1/2}$ |
| $\mu>2$ | $2$ | $n\mathbb{E}[X]$ | $[(1/2)\mathrm{Var}(X)]^{1/2}n^{1/2}$ |

###### Theorem 4 .

## Appendix B Proof of Theorem 1

First, the marginal variance of $w^{(\ell)}_{ij}$ is unbounded:

$$ $\displaystyle\mathbb{V}(w^{(\ell)}_{ij})$ $\displaystyle=\mathbb{V}[(s_{+}^{(\ell)})^{1/2}\mathbb{E}[{\tilde{w}_{ij}}^{( \ell)}\mid s_{+}^{(\ell)}]]+\mathbb{E}[s_{+}^{(\ell)}\mathbb{V}[{\tilde{w}_{ij }}^{(\ell)}\mid s_{+}^{(\ell)}]]=\infty,$ $$

since for Stable $(\alpha/2)$ random variables, the expectation is unbounded. Next, we have:

$$ $\displaystyle z_{j}^{(\ell)}(\mathbf{x}_{k})$ $\displaystyle=\frac{1}{M_{\ell}^{1/2}}\sum_{i=1}^{M_{\ell}}w^{(\ell)}_{ij}f_{i }^{(\ell)}(\mathbf{x}_{k})=(s_{+}^{(\ell)})^{1/2}\frac{1}{M_{\ell}^{1/2}}\sum_ {i=1}^{M_{\ell}}{\tilde{w}}^{(\ell)}_{ij}f_{i}^{(\ell)}(\mathbf{x}_{k}).$ $$

Since ${\tilde{w}}^{(\ell)}_{ij}$ are i.i.d. with unit variance by assumption, we have, as $M_{\ell}\to\infty$:

$$ $\displaystyle\frac{1}{M_{\ell}^{1/2}}\sum_{i=1}^{M_{\ell}}{\tilde{w}}^{(\ell)} _{ij}f_{i}^{(\ell)}(\mathbf{x}_{k})$ $\displaystyle\stackrel{{\scriptstyle D}}{{\to}}\mathcal{N}(0,\lambda_{k}^{( \ell)}),$ $$

by an application of the classical central limit theorem, where $\stackrel{{\scriptstyle D}}{{\to}}$ denotes convergence in distribution. Thus,

$$ $\displaystyle z_{j}^{(\ell)}(\mathbf{x}_{k})$ $\displaystyle\stackrel{{\scriptstyle D}}{{\to}}(s_{+}^{(\ell)})^{1/2}\mathcal{ N}(0,\lambda_{k}^{(\ell)}),$ $$

which is an $\alpha$-stable random variable for each $j,k$, and $\lambda_{k}^{(\ell)}=\mathbb{V}(f_{i}^{(\ell)}(\mathbf{x}_{k})\mid\{s_{+}^{(l)
}\}_{l=1}^{\ell-1},\mathbf{S}^{(1)}_{+})$, which does not depend on $i$ since the $f_{i}^{(\ell)}(\mathbf{x}_{k})$ are i.i.d. across $i$.

The conditional covariance between two features in layer one is given by:

$$ $\displaystyle\mathrm{Cov}(z_{j}^{(1)}(\mathbf{x}_{k}),z_{j}^{(1)}(\mathbf{x}_{ k^{\prime}})\mid\mathbf{S}^{(1)}_{+})$ $\displaystyle=\mathrm{Cov}\left(\sum_{m=1}^{I}(s^{(1)}_{m,+})^{1/2}\tilde{w}_{ mj}^{(1)}x_{m,k},\sum_{m^{\prime}=1}^{I}(s^{(1)}_{m^{\prime},+})^{1/2}\tilde{w }_{m^{\prime}j}^{(1)}x_{m^{\prime},k^{\prime}}\mid\mathbf{S}^{(1)}_{+}\right)$ $\displaystyle=\sum_{m^{\prime}=1}^{I}(s^{(1)}_{m^{\prime},+})^{1/2}x_{m^{ \prime},k^{\prime}}\sum_{m=1}^{I}(s^{(1)}_{m,+})^{1/2}x_{m,k}\mathrm{Cov}( \tilde{w}_{ij}^{(1)},\tilde{w}_{m^{\prime}j}^{(1)}\mid\mathbf{S}^{(1)}_{+})$ $\displaystyle=\sum_{m=1}^{I}s^{(1)}_{m,+}x_{m,k}x_{m,k^{\prime}}.$ $$

As such the conditional covariance matrix of $\mathbf{z}_{j}^{(1)}\mid\mathbf{S}^{(1)}_{+}$ is $\Sigma^{(1)}=\mathbf{X}\mathbf{S}^{(1)}_{+}\mathbf{X}^{T}$.
Now that we have the first covariance matrix, we proceed by induction over $\ell$, with the note that in the case $\ell=2$ what we denote by $s^{(1)}_{+}$ corresponds to $\mathbf{S}^{(1)}_{+}$. Next we obtain the characteristic function of $\mathbf{z}_{j}^{(\ell)}$ for a finite $M_{\ell}$, assuming that $\Sigma^{(\ell-1)},s^{(\ell)}_{+},s^{(\ell-1)}_{+}$ are given. To this end, for $\mathbf{t}=(t_{1},\dots,t_{n})\in\mathbb{R}^{n}$:

$$ $\displaystyle\phi_{\mathbf{z}_{j}^{(\ell)}\mid\Sigma^{(\ell-1)},s^{(\ell)}_{+} ,s^{(\ell-1)}_{+}}(\mathbf{t})$ $\displaystyle=\mathbb{E}\left[\exp\left\{\mathrm{i}\mathbf{t}^{T}\mathbf{z}_{j }^{(\ell)}\right\}\mid\Sigma^{(\ell-1)},s^{(\ell)}_{+},s^{(\ell-1)}_{+}\right]$ $\displaystyle=\mathbb{E}\left[\exp\left\{\mathrm{i}\sum_{k=1}^{n}\sum_{i=1}^{M _{\ell}}t_{k}M_{\ell}^{-1/2}w_{ij}^{(\ell)}f_{i}^{(\ell)}(\mathbf{x}_{k}) \right\}\mid\Sigma^{(\ell-1)},s^{(\ell)}_{+},s^{(\ell-1)}_{+}\right]$ $\displaystyle=\mathbb{E}\left[\exp\left\{\mathrm{i}\sum_{k=1}^{n}\sum_{i=1}^{M _{\ell}}t_{k}M_{\ell}^{-1/2}(s_{+}^{(\ell)})^{1/2}\tilde{w}_{ij}^{(\ell)}f_{i} ^{(\ell)}(\mathbf{x}_{k})\right\}\mid\Sigma^{(\ell-1)},s^{(\ell)}_{+},s^{(\ell -1)}_{+}\right]$ $\displaystyle=\mathbb{E}\left[\exp\left\{\mathrm{i}(s_{+}^{(\ell)})^{1/2}M_{ \ell}^{-1/2}\sum_{i=1}^{M_{\ell}}\tilde{w}_{ij}^{(\ell)}\sum_{k=1}^{n}t_{k}f_{ i}^{(\ell)}(\mathbf{x}_{k})\right\}\mid\Sigma^{(\ell-1)},s^{(\ell)}_{+},s^{( \ell-1)}_{+}\right]$ $\displaystyle=\mathbb{E}\left[\exp\left\{\mathrm{i}(s_{+}^{(\ell)})^{1/2}M_{ \ell}^{-1/2}\sum_{i=1}^{M_{\ell}}\tilde{w}_{ij}^{(\ell)}A_{i}^{(\ell)}\right\} \mid\Sigma^{(\ell-1)},s^{(\ell)}_{+},s^{(\ell-1)}_{+}\right],$ (4) $$

where we define $A_{i}^{(\ell)}=\sum_{k=1}^{n}t_{k}f_{i}^{(\ell)}(\mathbf{x}_{k})$, which are independent and identically distributed over $i$, and the expectation $\mathbb{E}(\cdot)$ is over the distribution of ${\tilde{w}}^{(\ell)}_{ij}$.
Next, the expectation of $\tilde{w}_{ij}^{(\ell)}A_{i}^{(\ell)}$ is zero since both factors of the product are independent and the mean of one of them is zero, and the variance is $\mathbb{E}[(\sum_{k=1}^{n}t_{k}f_{i}^{(\ell)}(\mathbf{x}_{k}))^{2}\mid\Sigma^{
(\ell-1)},s^{(\ell-1)}_{+}]$.

Taking the limit $M_{\ell}\to\infty$ of the expression inside the exponential in Equation (4), we obtain:

$$ $\displaystyle\lim_{M_{\ell}\to\infty}(s_{+}^{(\ell)})^{1/2}M_{\ell}^{-1/2}\sum _{i=1}^{M_{\ell}}\tilde{w}_{ij}^{(\ell)}A_{i}^{(\ell)}$ $\displaystyle\overset{D}{\to}(s_{+}^{(\ell)})^{1/2}\mathcal{N}\left(0,\mathbb{ E}\left[\left(\sum_{k=1}^{n}t_{k}f_{i}^{(\ell)}(\mathbf{x}_{k})\right)^{2}\mid \Sigma^{(\ell-1)},s^{(\ell-1)}_{+}\right]\right)$ $\displaystyle=(s_{+}^{(\ell)})^{1/2}\times\eta_{j}^{(\ell)},$ $$

since $s^{(\ell)}$ and $\tilde{w}_{ij}^{(\ell)}$ are independent and $\eta_{j}^{(\ell)}$ denotes the normal random variable on the RHS of first line of the above display, which can be seen to be i.i.d. over $j$. Define the matrix $\Sigma^{(\ell)}$ with $\Sigma^{(\ell)}_{k,h}=\mathbb{E}[f_{i}^{(\ell)}(\mathbf{x}_{k})f_{i}^{(\ell)}(
\mathbf{x}_{h})\mid\Sigma^{(\ell-1)},s^{(\ell-1)}_{+}]$, which is a valid covariance matrix and does not depend on $i$.
Thus, since $\Sigma^{(\ell)}$ is a deterministic function of $\Sigma^{(\ell-1)}$ and $s_{+}^{(\ell-1)}$, we have:

$$ $\displaystyle\lim_{M_{\ell}\to\infty}\phi_{\mathbf{z}_{j}^{(\ell)}\mid\Sigma^{ (\ell)}}(\mathbf{t})$ $\displaystyle=\exp\left\{-(1/2)(\mathbf{t}^{T}\Sigma^{(\ell)}\mathbf{t})^{ \alpha/2}\right\},$ $$

by an application of Lévy’s continuity theorem (Billingsley, , Th. 26.3). As such, $\mathbf{z}_{j}^{(\ell)}\mid\Sigma^{(\ell)}$ is an elliptical $\alpha$-stable random vector.

Next, we proceed with the calculation of $\Sigma^{(\ell)}\mid\Sigma^{(\ell-1)},s^{(\ell-1)}_{+}$ in the case of $\tilde{w}_{ij}^{(\ell-1)}$ being standard Gaussian weights and when $g_{\delta}(\zeta)=\zeta^{\delta}\mathbf{1}_{\{\zeta>0\}}$. For this, note that $(z_{j}^{(\ell-1)}(\mathbf{x}_{k}),z_{j}^{(\ell-1)}(\mathbf{x}_{h}))^{T}\mid
\Sigma^{(\ell-1)},s^{(\ell-1)}_{+}\sim\mathcal{N}(0,s^{(\ell-1)}_{+}\Sigma^{(
\ell-1)}_{*})$, where $\Sigma^{(\ell-1)}_{*}$ is a two-by-two symmetric matrix with diagonal $(\Sigma^{(\ell-1)}_{k,k},\Sigma^{(\ell-1)}_{h,h})$, and off-diagonal entries $\Sigma^{(\ell-1)}_{h,k}$. It follows that:

$$ $\displaystyle\left(b_{j}^{(\ell-1)}+z_{j}^{(\ell-1)}(\mathbf{x}_{k}),b_{j}^{( \ell-1)}+z_{j}^{(\ell-1)}(\mathbf{x}_{h})\right)^{T}\mid\Sigma^{(\ell-1)},s^{( \ell-1)}_{+}\overset{i.i.d.}{\sim}\mathcal{N}\left(0,\mathbf{U}+s^{(\ell-1)}_{ +}\Sigma^{(\ell-1)}_{*}\right),$ (5) $$

where $\mathbf{U}$ is a two-by-two matrix with all entries equal to one, with the independence over the $j$s and conditional on $(\Sigma^{(\ell-1)},s^{(\ell-1)}_{+})$.
Then for $k,h=1,\dots,n$:

$$ $\displaystyle\Sigma^{(\ell)}_{k,h}=\mathbb{E}[f_{j}^{(\ell)}(\mathbf{x}_{k})f_ {j}^{(\ell)}(\mathbf{x}_{h})\mid\Sigma^{(\ell-1)},s^{(\ell-1)}_{+}]=\int_{ \mathbb{R}^{2}}g_{\delta}(\zeta_{1})g_{\delta}(\zeta_{2})p(\zeta_{1},\zeta_{2} )d\zeta_{2}d\zeta_{1},$ $$

where in the second line $p(\zeta_{1},\zeta_{2})$ is the density of the random vector in Equation (5). It follows that:

$$ $\displaystyle\Sigma^{(\ell)}_{k,h}$ $\displaystyle=\int_{\mathbb{R}}\int_{\mathbb{R}}\zeta_{1}^{\delta}\mathbf{1}_{ \{\zeta_{1}>0\}}\zeta_{2}^{\delta}\mathbf{1}_{\{\zeta_{2}>0\}}p(\zeta_{2}, \zeta_{1})d\zeta_{2}d\zeta_{1}$ $\displaystyle=\frac{1}{2\pi v_{1}v_{2}\sqrt{1-\rho^{2}}}\int_{0}^{\infty}\int_ {0}^{\infty}\zeta_{1}^{\delta}\zeta_{2}^{\delta}\exp\left(-\frac{1}{2(1-\rho^{ 2})}\left[\zeta_{1}^{2}v_{1}^{-2}-2\rho\frac{\zeta_{1}\zeta_{2}}{v_{1}v_{2}}+ \zeta_{2}^{2}v_{2}^{-2}\right]\right)d\zeta_{1}d\zeta_{2}$ $\displaystyle=v_{1}^{\delta}v_{2}^{\delta}\frac{1}{2\pi\sqrt{1-\rho^{2}}}\int_ {0}^{\infty}\int_{0}^{\infty}\eta_{1}^{\delta}\eta_{2}^{\delta}\exp\left(- \frac{1}{2(1-\rho^{2})}\left[\eta_{1}^{2}-2\rho\eta_{1}\eta_{2}+\eta_{2}^{2} \right]\right)d\eta_{1}d\eta_{2},$ $$

where the second line follows by the bivariate normal p.d.f. where $v_{i}$ is the standard deviation for the $i$-th entry implied by Equation (5), $\rho$ is the respective correlation and by evaluating the indicator function. In the last line we use the change of variables $\eta_{i}=\zeta_{i}/v_{i}$. This last expression corresponds to a multiple of Equation (15) of Cho & Saul (), where $\rho^{2}=\cos^{2}(\theta)$. Explicitly:

$$ $\displaystyle v_{1}^{2}$ $\displaystyle=1+s^{(\ell)}_{+}\Sigma^{(\ell-1)}_{k,k},\,\,v_{2}^{2}=1+s^{(\ell )}_{+}\Sigma^{(\ell-1)}_{h,h},\,\,\rho=\frac{1+s^{(\ell)}_{+}\Sigma^{(\ell-1)} _{k,h}}{v_{1}v_{2}}.$ $$

Finally:

$$ $\displaystyle\Sigma_{k,h}^{(\ell)}=\pi^{-1}\left[v_{1}^{2}v_{2}^{2}\right]^{ \delta/2}J_{\delta}(\theta),\,\theta=\cos^{-1}(\rho),$ $$

as required.

## Appendix C Proof of Proposition 2

The expressions for the conditional expectation and variance, $\bm{\mu}^{*}$ and $\bm{\Lambda}^{*}$ are a direct consequence of the posterior predictive distribution of Gaussian random vectors. We omit the details.

## Appendix D Supplementary algorithms

The following algorithms are used in Algorithm 1. They compute the covariance matrix $\Sigma^{(L)}$ from the random scales. The algorithm used is an independent sample Metropolis–Hastings sampler where we propose from the prior, and accept with probability corresponding to the ratio of the likelihoods.

Figure: Algorithm 2 A Metropolis–Hastings sampler for $s^{(\ell)}_{+}\mid s^{(\ell+1)}_{+},\mathbf{y}$

Figure: Algorithm 3 A Metropolis–Hastings sampler for $s^{(L)}_{+}\mid\mathbf{y}$

## Appendix E { Σ ( ℓ ) } ℓ = 1 L superscript subscript superscript Σ ℓ ℓ 1 𝐿 \{\Sigma^{(\ell)}\}_{\ell=1}^{L} { roman_Σ start_POSTSUPERSCRIPT ( roman_ℓ ) end_POSTSUPERSCRIPT } start_POSTSUBSCRIPT roman_ℓ = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_L end_POSTSUPERSCRIPT is consistent under marginalization

Consistency under marginalization means that for $K_{1}\sim G(K_{0})$, where both $K_{1},K_{0}$ are p.d. $n\times n$ matrices, then a principal $P\times P$ submatrix $K_{0}^{*}$ of $K_{0}$ also induces a distribution $K_{1}^{*}\sim G(K_{0}^{*})$, and those distributions are consistent, i.e., if we consider the full p.d.f. of $K_{1}$ and integrate out the entries that do not belong to $K_{1}^{*}$ we obtain the same p.d.f. This is a property that is known for Wishart and inverse Wishart matrices (Dawid, ).

First, consider the set of positive random scales in the first layer $\mathbf{S}^{(1)}_{+}$ fixed. Note that $\Sigma^{(1)}$ is a p.d. matrix for any number of observations $n$, since the obtained cross product in Appendix B is a valid covariance function. Similarly, for $P<n$ any $P\times P$ principal submatrix $(\Sigma^{(1)})^{*}_{1:P,1:P}$ (wlog) is a consistent submatrix, since the random scales are fixed.

Now, consider a principal submatrix (wlog) in the $\ell$ layer: $(\Sigma^{(\ell)})_{1:P,1:P}$,
and consider $s^{(\ell)}_{+}$ fixed. Then $(\Sigma^{(\ell+1)})_{1:P,1:P}$ is a deterministic function of $(\Sigma^{(\ell)})_{1:P,1:P}$, since $s^{(\ell)}_{+}$ is fixed. This implies that $(\Sigma^{(\ell+1)})_{1:P,1:P}$ is consistent under marginalization, as the only randomness comes from $s^{(\ell)}_{+}$ over which we can integrate.

## Appendix F Proof of Proposition 3

First consider $\alpha<2$. The posterior of the features is given by:

$$ $\displaystyle p(\mathbf{z}_{j}^{(\ell)}\mid\mathbf{X},\mathbf{y})=$ $\displaystyle\int p(\mathbf{z}_{j}^{(\ell)},\{s^{(\ell)}_{+}\}_{\ell=2}^{L}, \mathbf{S}^{(1)}_{+}\mid\mathbf{X},\mathbf{y})\prod_{\ell=2}^{L}ds^{(\ell)}_{+ }\prod_{m=1}^{I}ds^{(1)}_{m,+}$ $\displaystyle=$ $\displaystyle\int p(\mathbf{z}_{j}^{(\ell)}\mid\{s^{(\ell)}_{+}\}_{\ell=2}^{L} ,\mathbf{S}^{(1)}_{+},\mathbf{X},\mathbf{y})$ $\displaystyle\times p(\{s^{(\ell)}_{+}\}_{\ell=2}^{L},\mathbf{S}^{(1)}_{+}\mid \mathbf{X},\mathbf{y})\prod_{\ell=2}^{L}ds^{(\ell)}_{+}\prod_{m=1}^{I}ds^{(1)} _{m,+}.$ (6) $$

Employing the elliptical $\alpha$-stable characterization of Theorem 1, the distribution of $\mathbf{z}_{j}^{(\ell)}$ is given by:

$$ $\displaystyle\mathbf{z}_{j}^{(\ell)}\mid\{s^{(\ell)}_{+}\}_{\ell=2}^{L}, \mathbf{S}^{(1)}_{+}$ $\displaystyle\sim\mathcal{N}(0,s_{+}^{(\ell)}\Sigma^{(\ell)}).$ $$

As such, we have that the posterior of $\mathbf{z}_{j}^{(\ell)}$ can depend on $\mathbf{y}$ only through the posterior of $\{s^{(\ell)}\mid\mathbf{y}\}$ for $\ell=2,\ldots,L$, since the $\Sigma^{(\ell)}$ is a deterministic function of a fixed design matrix $\mathbf{X}$ and conditional on the random scales $s^{(2)}_{+},\ldots,s^{(\ell)}_{+},\mathbf{S}^{(1)}_{+}$.

Now, when $\alpha=2$, all $s^{(\ell)}_{+}$ are equal to one with probability 1 (i.e., they are degenerate Dirac point mass at one) and therefore the posterior $s^{(\ell)}_{+}\mid\mathbf{y}$ is also a degenerate point mass at 1. Thus, the posterior of the features cannot depend on the observations when $\alpha=2$, i.e., in the Gaussian case.

## Appendix G Conditional mutual information

For a bivariate Gaussian vector the mutual information is easily computed by: $-(1/2)\log(1-\rho^{2})$, where $\rho$ is the correlation coefficient between the two entries of the vector (Cover & Thomas, , p. 252). However, recall that in our non-Gaussian setting, marginal moments or correlations are not defined and we need a conditional construction. Consider three random variables $Y_{1},Y_{2},S$, the conditional mutual information (Cover & Thomas, , p. 23) is given by

$$ $\displaystyle I(Y_{1};Y_{2}\mid S)$ $\displaystyle=\int_{\mathcal{S}}\int_{\mathcal{Y}_{1}\times\mathcal{Y}_{2}} \log\left(\frac{p(y_{1},y_{2}\mid s)}{p(y_{1}\mid s)p(y_{2}\mid s)}\right)p(y_ {1},y_{2}\mid s)p(s)dy_{1}dy_{2}ds$ $\displaystyle=\int_{\mathcal{S}}D_{KL}[p(Y_{1},Y_{2}\mid s)\;||\;p(Y_{1}\mid s )p(Y_{2}\mid s)]p(s)ds,$ $$

where $D_{KL}$ is the Kullback–Leibler divergence. In our case, we consider $Y_{1},Y_{2}$ as two observations from the deep $\alpha$-kernel process, and $S$ corresponds to the conditional scales. Since given $S$ the observations are conditionally-Gaussian, we can apply the formula that employs the correlation and integrate over the scales to obtain the mutual information. Specifically, we use:

$$ $\displaystyle I(Y_{1};Y_{2}\mid S)$ $\displaystyle=-(1/2)\int_{\mathcal{S}}\log(1-\rho^{2}_{Y_{1},Y_{2}}(s))p(s)ds.$ $$

To perform the integration we use Monte Carlo samples of $\{s_{+}^{(\ell)}\}_{\ell=2}^{L},\mathbf{S}^{(1)}_{+}$, and use the correlation induced by the $\Sigma^{(L)}$.

## Appendix H Additional numerical results for simulated data

### H.1 Timing of simulations

In Table 5 we display the running time for the six methods employed in the simulations described in Section 3. All the times are in seconds and display CPU times. We avoid any type of parallelization throughout our experiments.

**Table 5: Average (SD) time (s) for the six methods in the numerical examples for one, two, and ten dimensions.**
| Method | 1D | 2D | 10D |
| --- | --- | --- | --- |
| D$\alpha$-KP | 181.23 (15.3) | 200.19 (31.16) | 9713.51 (212.10) |
| DIWP | 206.56 (10.24) | 312.29 (21.17) | 14509.31 (2783.10) |
| GP Bayes | 17.26 (0.44) | 21.35 (0.24) | 729.78 (26.96) |
| GP MLE | 0.08 (0.00) | 0.18 (0.01) | 49.89 (2.51) |
| NNGP | 131.03 (11.99) | 167.85 (6.36) | 11044.99 (1927.03) |
| Stable | 33.14 (5.64) | 1456.72 (85.45) | – |

### H.2 Coverage of credible intervals

In Table 6 we display the percent of observations in the out-of-sample observations captured by the 90% credible intervals produced by the Bayesian methods. The GP MLE method is a frequentist method and as such does not produce valid credible intervals. The Stable method is not available for 10 dimensions. The methods that use stable random variables are the only ones that consistently manage to be close to the nominal coverage of 90%. The other methods typically have coverage much below 90%, especially in high dimensions.

**Table 6: Average (SD) percent of out-of-sample observations captured by the 90% credible intervals in the twenty splits of the simulation examples.**
| Method | 1D | 2D | 10D |
| --- | --- | --- | --- |
| D$\alpha$-KP | 90.85% (5.12) | 92.47% (5.63) | 89.70% (2.42) |
| DIWP | 92.70% (1.89) | 95.86% (1.22) | 60.47% (3.08) |
| GP Bayes | 94.85% (2.43) | 94.14% (3.02) | 65.17% (5.82) |
| GP MLE | – | – | – |
| NNGP | 92.50% (1.79) | 95.99% (1.12) | 60.65% (2.82) |
| Stable | 91.65% (4.40) | 91.11% (4.63) | – |

### H.3 Visualizations of simulation results

#### H.3.1 Plots for two dimensions

Figure 4 displays the function fit for the example in two dimensions and the residuals with the accompanying credible intervals. We do this for the first split of the dataset. Recall that the true function is $f(\xi_{1},\xi_{2})=5\times\mathbf{1}_{\{\xi_{1}>0\}}+5\times\mathbf{1}_{\{\xi_
{2}>0\}}$. The D$\alpha$-KP and the Stable method are the only ones that are able to capture the jump without oversmoothing, and with narrower credible intervals.

Figure: (a) Function fit for the different methods.
Refer to caption: x6.png

#### H.3.2 Predictions and uncertainty quantification in high dimensional simulation

In Figure 5 we display the a scatter plot of residuals (observation - prediction) with the prediction on the first split of the simulation example for 10 dimensions. Note that the D$\alpha$-KP and GP MLE have the smallest range of residuals, while the other methods have a wider range of residuals.

Figure: Figure 5: Residual plot for first split of the high dimensional example with 90% credible intervals.
Refer to caption: x8.png

## Appendix I Additional numerical results for UCI data

Figures 6, 7 and 8 displays a scatter plot of the predictions and residuals (observation - prediction) of the test set in the first split of the UCI datasets, as well as the credible intervals for said observations.
The predictions of the D$\alpha$-KP are much closer to the observations, while the DIWP and NNGP present worse predictions relative to the observations. Moreover, the width of the credible intervals of the D$\alpha$-KP vary, while the DIWP and NNGP are much more uniform. The results indicate that the credible intervals of the D$\alpha$-KP capture more observations than the competing methods. Figure 9 displays the trace plot for the log-likelihood of the fit of the deep $\alpha$-kernel process in the first split of the Boston dataset, indicating a fast mixing.

Figure: Figure 6: Comparison of predictions and observations in the test set for the first split of the Boston dataset, with $90\%$ credible interval.
Refer to caption: x9.png

Figure: Figure 7: Comparison of predictions and observations in the test set for the first split of the Energy dataset, with $90\%$ credible interval.
Refer to caption: x10.png

Figure: Figure 8: Comparison of predictions and observations in the test set for the first split of the Yacht dataset, with $90\%$ credible interval.
Refer to caption: x11.png

Figure: Figure 9: Traceplot of the log-likelihood for the first split of the Boston dataset.
Refer to caption: x12.png

### I.1 Timing of UCI experiments

In Table 7 we display the average running time for each of the methods, except the Stable method as the input dimensions of the data sets are all greater than 2. Overall the GP MLE and GP Bayes are faster than the kernel process methods. The kernel process methods have a similar running time in each dataset.

**Table 7: Average (SD) time (s) for the methods in the three UCI datasets over twenty different splits. Times for the Stable method are not available.**
| Method | Boston | Energy | Yacht |
| --- | --- | --- | --- |
| D$\alpha$-KP | 22320.89 (690.23) | 36154.63 (1812.99) | 4590.63 (149.31) |
| DIWP | 20104.97 (3877.92) | 41291.76 (15973.11) | 22177.2 (10017.53) |
| GP Bayes | 1047.13 (6.51) | 2151.64 (12.37) | 271.82 (2.21) |
| GP MLE | 142.06 (6.68) | 390.52 (9.36) | 22.18 (1.20) |
| NNGP | 18242.54 (8474.21) | 36872.04 (21713.89) | 9177.82 (3803.24) |
| Stable | – | – | – |

### I.2 Coverage of credible intervals for UCI datasets

In Table 8 we display the average coverage of the 90% credible intervals for the UCI datasets in the same splits used in Section 4. The D$\alpha$-KP consistently is close to the nominal $90\%$ coverage, while the other methods do not have a consistent coverage close to the nominal 90% percent.

**Table 8: Coverage (SD) of the 90% credible intervals for the methods in the UCI datasets over twenty different splits. Stable method is not available, and GP MLE does not produce credible intervals.**
| Method | Boston | Energy | Yacht |
| --- | --- | --- | --- |
| D$\alpha$-KP | 91.86% (5.14) | 87.40% (3.91) | 95.32% (2.86) |
| DIWP | 57.35% (7.50) | 48.77% (6.05) | 90.00% (6.61) |
| GP Bayes | 93.24% (4.43) | 97.01% (1.79) | 95.00% (3.22) |
| GP MLE | – | – | – |
| NNGP | 68.53% (8.26) | 68.77% (5.45) | 67.58% (7.65) |
| Stable | – | – | – |

## Appendix J Feature learning in additional datasets

In Figures 10 and 11 we provide the scatter plot of the features of the Energy and Yacht data sets under D$\alpha$-KP, along with the corresponding box plots and q-q plots for each of the features. Similar to the Boston dataset, the Energy features display a heavy-tailed behavior, while those for Yacht are closer to normality.

Figure: (a) Scatterplot of the first two features, with box plots for the marginals.
Refer to caption: x13.png

Figure: (a) Scatterplot of the first two features, with box plots for the marginals.
Refer to caption: x15.png

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 5
presentation: 2
questions: Q1. Can we summarize the key points as follows: Under a specific infinite variance prior, the covariance kernel becomes a stochastic process at each layer. Even if the marginal does not exist, sequential sampling is still possible. Furthermore, conditioning on a covariance sample at a given layer yields a Gaussian representation, while marginalizing over the covariance leads the representation to follow a heavy-tailed distribution centered at 0.

Q2. You could perhaps include a 1-dimensional regression example to illustrate the model’s behavior for varying $\alpha$. Additionally, it might be helpful to show the covariance matrix to demonstrate its properties and the correlation among input data.

Q3. Could you clarify the technical differences between the current paper and Loria & Bhadra (2024), which makes the approach feasible?
soundness: 3
strengths: - The technical aspects of the paper appear solid, although I have not checked the results in detail.
- The experiments seem to effectively support the theoretical claims.
summary: The authors consider Bayesian neural networks in the infinite-width limit under infinite variance priors on the weights. They show that the resulting covariance kernel becomes a random process that depends on the given data, in contrast to the finite variance prior case, where the kernel is deterministic. Additionally, they demonstrate that the representation in each layer is Gaussian when conditioned on the associated covariance realization. In the experiments, they show that the data-dependent stochastic covariance enables feature learning, and that the resulting predictive Gaussian process achieves improved performance on benchmarks.
weaknesses: - The paper is somewhat challenging to approach, given its niche topic and highly technical content. It also feels quite text-heavy.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: **Questions:**
- I don't understand this sentence, could you please clarify? "The third limitation is one noted by Neal (1996), that in the case of multiple outputs, GPs cannot learn the covariance structure, as the covariance is always zero under a Gaussian limit." Are you saying that for a vector-valued neural network output, any two outputs are independent? This would be false, even in the finite variance setting, as long as one considers non-isotropic Gaussian weights (which is relatively straight-forward). E.g. see weights W in section 1 of "Avoiding Kernel Fixed Points: Computing with ELU and GELU infinite networks, AAAI 2021".
- Another related work might be "Richer priors for infinitely wide multi-layer perceptrons". In particular, Theorem 1 and Proposition 2 seem to show that the resulting posterior predictive distribution also marginalises out an $\alpha$-stable distributed scale matrix in the covariance matrix. This seems similar to equation (19), where a scale and a location parameter are integrated out in the posterior predictive, and the resulting model is a kind of infinite mixture of GPs. Although I think the motivation for considering such infinite mixtures is different: this paper focusses on heavy tailed whereas the other one focusses on exchangeability, and considers more "mild" mixing. Still, it could be nice to motivate the current approach both from the perspective of heavy-tailed priors and from the perspective of exchangeability in the hidden layers.
- Is it possible to construct deep $\alpha$-stable processes for gradient flow (extending NTK), rather than BNNs?
- As far as I can tell from the proof of Theorem 1 in Appendix B, rectified monomial activations are not strictly required for the proof to work. All that is required is that a recursive equation for the conditional covariance is available, which is true beyond this activation for other activations (see e.g. the works I mentioned above). It might be nice to extend the theoretical statement to more broad activation functions.
soundness: 3
strengths: **Strengths:**
- The proposed model is capable of representation learning. Specifically, Proposition 3 nicely shows that the feature at layer l depends on the data X, y observed in training, which is unlike many finite width models. It is nice to have this notion of representation learning formalised in this simple way.
- The idea of the paper is relatively straightforward: everything is conditionally Gaussian given this scale parameter, which induces the heavy-tailed behaviour. I see this simplicity as a strength, as it doesn't require complicated constructs to convey a nice and useful result.
- The method is demonstrated, firstly on some artifical but enlightening discontinuous functions (where it either almost matches the shallow stable process or is able to compute a posterior where the shallow stable process fails), and also on some UCI benchmarks, where it usually outperforms other methods and the previous shallow stable process is intractable. There are also some nice graphics comparing the posterior predictives of other models in 1D and a QQ plot showing non-Gaussian behaviour of the features. More extreme non-Gaussian behaviour is shown in the appendix (in particular, Figure 10 - I am wondering why the authors did not choose to include this as their main non-Gaussian demonstration in the main paper?)
- Appendix B and D, which contain the proofs, appear to not have any major technical errors.
summary: **Summary**
Neal's infinite width neural networks, even when extended to the deep setting, do not account for representation learning. That is, they result in a GP with a fixed covariance function. This is unlike typical but intractable Bayesian neural networks, which are capable of representation learning. The authors consider heavy tailed priors, elliptically distributed with infinite variance, and show that the network converges to a process with $\alpha$-stable marginals. Each layer is conditionally Gaussian, even though the covariances are not necessarily even defined.  The main technical tool is the CLT for alpha-stable RVs (Theorem 4 from another paper, with alpha=2 being the classical CLT), and a (conditionally Gaussian) Gaussian mixture representation of the characteristic function of the alpha-stable RV (equation (1)). Furthermore, a computationally viable MCMC method for doing inference is also provided. The derived method is demonstrated on some UCI benchmarks, showing better prediction than other methods, and being faster than a previous heavy-tailed method.
weaknesses: **Weaknesses:**
- Unless I am mistaken, "The key finding is that the conditional mutual information decays at a slower rate for smaller α" should be "The key finding is that the mutual information [which itself is computed via MCMC as an expected conditional mutual information] decays at a slower rate for smaller α." Figure 1 shows a mutual information, not conditional mutual information.


**Minor:**
- Theorem 1. $J_\delta(\theta)$ can be computed explicitly, and this is obvious to people familiar with the result of Cho & Saul (as mentioned in the following paragraph). But some readers might find it alarming to see $J_\delta(\theta)$ before it is defined. Perhaps the result of Cho & Saul could be made more explicit before theorem 1 (i.e. move $J$ up in the page)
- Figure 1 and 2 labels are too small. Sans-serif font clashes with the text in the main paper.
- last sentence of first paragraph of section 4.2 seems to have odd grammar.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 2
questions: * Could the authors please offer more comments on the motivation for considering this particular model? 
* Could the authors please offer comments on what the concrete and significant technical advances are?
soundness: 2
strengths: The article 
* provides a discussion of references, theoretical results, and experiments. 
* discusses potential benefits in prediction and uncertainty quantification. 
* suggests feature learning that is not possible under a Gaussian process.
summary: The article presents a deep kernel process with infinite variance priors including a recursive formula and experiments illustrating dependence on the inputs and behavior for discontinuous targets.
weaknesses: 1. My main concern with the article is the writing and presentation, which I think need to be improved. 
The abstract gives a long discussion of prior works but ideally it should instead give a crisp description of the main points in the article. The lengthy discussion in the introduction comments on prior works and perceived limitations, but does not provide a sufficiently concise and clear description of the objective, motivation, and contributions of the present work. Terminology could be introduced more clearly in preliminaries, rather than in prior works. The writing is in parts repetitive and the thread is not sufficiently clear.  \
&nbsp;
\
Suggestions proposed by the ICLR Associate Program Chair: 
    1. Revise the abstract to focus on succinctly stating the key contributions and results of the paper, rather than discussing prior work.
    2. Restructure the introduction to clearly state the objectives, motivation and main contributions upfront before discussing related work.
    3. Add a preliminaries section to introduce key terminology and concepts.
    4. Review the paper to remove repetition and improve the logical flow of ideas between sections.  \
&nbsp;

2. My second concern is the motivation and innovations, which I do not find sufficiently clearly articulated in the article. 
For Theorem 1 it is mentioned that the result permits characterization of the deep kernel process in an infinite variance setting and for Proposition 3 that it results in a posterior distribution of the features that depends on the observations. I can see that there are differences to prior works and that feature learning is an interesting topic, but in my opinion it is not sufficiently clearly explained what the motivation is for formulating this particular deep kernel process. Another point that is not described in sufficient clarity is what the concrete and significant conceptual and technical innovations are that go into obtaining the proofs and how they could be useful for obtaining further relevant advances.  \
&nbsp;
\
Suggestions proposed by the ICLR Associate Program Chair: 
    1. Provide a clearer explanation of why an infinite variance setting is important or beneficial for deep kernel processes.
    2. Discuss more explicitly how the ability to learn features from observations addresses limitations of previous approaches.
    3. Highlight the key technical challenges in proving Theorem 1 and Proposition 3, and how overcoming these challenges advances the field.
    4. Explain potential applications or extensions of their technical innovations beyond this specific model.
    5. Provide a clearer explanation of why this particular deep kernel process formulation was chosen and what advantages it offers.  \
&nbsp;

3. Description of the experiments. 
In my opinion the motivation and description of the experiments are not sufficiently clearly presented.  \
&nbsp;
\
Suggestions proposed by the ICLR Associate Program Chair: 
    1. Provide clearer motivation for each experiment, explaining how it relates to the theoretical results.
    2. Improve the description of experimental setups, including details on datasets, baselines, and evaluation metrics.
    3. Discuss the implications of the experimental results in more depth, relating them back to the theoretical contributions.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: N/A
soundness: 3
strengths: - Due to my limited familiarity with α-stable processes, I didn’t review all of the theoretical details. However, the paper reads as well-structured and complete, and the writing clearly highlights the contributions in the context of prior work with a solid mix of theoretical grounding, ablation studies, and performance benchmarks on public datasets.
- The problem statement is otherwise well-described, and I think this work should be of interest to future work on BNNs by addressing core limitations of smoothness and finite variance in existing models. I also wonder whether the non-Gaussian priors explored in this work could also be considered in other deep models?
- I appreciated the extensive numerical experiments that showcased the model’s performance in controlled settings that properly validate it meets the intended criteria of modelling non-Gaussian behaviours. I do wish the section real-world datasets was more compelling however.

**Post-rebuttal update**

Raised score from 6 to 8.
summary: Although deep kernel processes allow for more flexibility than standard GPs, they still assume smoothness and finite variance, which is unsuitable for modelling data with heavy tails or discontinuities. This is particularly relevant to representation learning, where deterministic covariance kernels limit the ability to learn non-degenerate, data-dependent representations. In contrast, $\alpha$-stable distributions are well-suited to problems were large fluctuations or discontinuities can be expected. Taking inspiration from these strengths, in this paper the authors assign $\alpha$-stable distributions with infinite variance to the weights in each layer of a BNN, such that the BNN now converges to a deep $\alpha$-stable kernel process rather than GP in the infinite-width limit.

They also develop conditional Gaussian representations of $\alpha$-stable processes that allow for recursive computation of the stochastic kernels at each layer. Through a combination of numerical experiments on synthetic discontinuous functions and evaluations on widely-used datasets from the UCI repository, the authors demonstrate that the deep $\alpha$-stable kernel process outperforms competing approaches when modelling discontinuous or heavy-tailed data. However, simpler models or GPs with traditional Gaussian priors might be preferable for smoother data due to their lower computational cost.
weaknesses: Although the paper is well-written overall, the dense presentation of different concepts makes it challenging to intuitively grasp without visual aids or simplified examples. I would recommend the inclusion of more figures illustrating how the α-stable kernel process differs from traditional GPs, ideally showing its behavior on simple synthetic data to highlight the flexibility and stochasticity introduced by their method. Some detail can likely also be moved to the appendix instead.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: N/A
soundness: 3
strengths: The theoretical contributions are extremely strong.  They consider an extremely interesting regime, and the question of how to avoid deterministic kernels in infinite-width networks is a critical one for Bayesian neural networks.

They have an extensive and excellent Appendix.
summary: This paper considers an infinite-width Bayesian neural network, with infinite-variance weights.  Interestingly, in this setting, representations/kernels do not become deterministic as usual in infinite-width limits, but remain stochastic, and hence learnable.
weaknesses: The experiments are weaker, but this can be expected in a paper with a mainly theoretical focus and contribution, such as this one.

Meta-review: This paper explores infinite-width Bayesian neural networks with elliptically-distributed weights with infinite variance, demonstrating convergence to a process with α-stable marginals and conditionally Gaussian representations. This results in a stochastic, data-dependent covariance kernel capable of representation learning, unlike the deterministic kernels arising from finite-variance priors. The authors also provide a computationally efficient MCMC method for inference and demonstrate performance benefits on benchmark datasets, particularly in modeling discontinuous or heavy-tailed data.

The reviewers noted that the paper's presentation is dense and technical, making it challenging for a broader audience, particularly in the abstract and introduction. They also found that the motivation behind the specific model formulation and its advantages over previous approaches could be articulated more clearly, and that the experimental descriptions lack depth and connection to theoretical findings. Nevertheless, the novel perspectives, strong theoretical results, and interesting conclusions suggest that this paper should be published and I recommend acceptance.

Decision: Accept (Poster)