[a]Bastian B. Brandt

# QCD thermodynamics at non-zero isospin asymmetry

Francesca Cuteri  Gergely Endrődi

###### Abstract

We study the thermodynamic properties of QCD at nonzero isospin chemical
potential using improved staggered quarks at physical quark masses. In
particular, we discuss the determination of the equation of state at zero and
nonzero temperatures and show results. Using the results for the isospin
density $n_{I}$, we also determine the phase diagram in the $(n_{I},T)$-plane.

##  1 Introduction

The QCD phase diagram and the equation of state (EoS) at non-zero quark
densities have been the subject of intense experimental and theoretical
research in the past two decades. Physical systems typically feature a non-
zero charge and strangeness chemical potential on top of the baryon chemical
potential. The effect of the latter is usually dominant, but the other
chemical potentials can potentially also contribute a sizeable and important
contribution to the dynamics of the system. Furthermore, in some cases the
charge chemical potential can play the major role. This is the case for an
early Universe starting with large lepton flavour asymmetries [1, 2, 3]. Non-
zero charge chemical potential directly translates into a non-zero isospin
chemical potential $\mu_{I}$. While lattice simulations with a generic
combination of quark chemical potentials suffer from the infamous complex
action (or sign) problem, QCD at pure isospin chemical potential, i.e. at
vanishing other chemical potential components, has a real action and is
amenable to direct Monte-Carlo simulations.

Figure 1: Continuum phase diagram of QCD at pure isospin chemical potential
[4, 5] (taken from Ref. [6]).

In the past few years we have performed extensive studies of QCD at non-zero
isospin asymmetry. In particular, we have computed the continuum phase diagram
at physical pion mass [4, 5], first discussed in chiral perturbation theory in
Ref. [7] 111First pioneering studies on the phase structure have been reported
in Refs. [8, 9, 10, 11, 12, 13]., featuring the standard hadronic and quark-
gluon plasma (QGP) phases and a phase with a Bose-Einstein condensate of
charged pions (BEC phase). The phase diagram is shown in Fig. 1. Furthermore,
perturbation theory predicts the existence of a superconducting BCS phase at
asymptotically large $\mu_{I}$ [7] and we are currently investigating whether
this phase could already be present for the intermediate chemical potentials
at our disposal, see Ref. [14] and the contribution to these proceedings Ref.
[15].

In this proceedings article we discuss the extraction of the EoS at non-zero
$\mu_{I}$ and present results for different temperatures and $N_{t}=8,\,10$
and 12 lattices. First accounts on our study of the EoS have been given in
Refs. [3, 16, 6]. Furthermore, we present the phase diagram in the plane of
isospin density and temperature, $(n_{I},T)$-plane, for which we perform a
tentative continuum extrapolation.

##  2 The equation of state at non-zero isospin asymmetry

The extraction of the EoS at $T=0$ has so far been done on a single lattice
spacing only. The details have already been published in Ref. [17]. To extract
the EoS at $\mu_{I}\neq 0$ and $T\neq 0$, we use the ensembles of Ref. [4],
which have already been used to map out the phase diagram up to
$\mu_{I}/m_{\pi}\lesssim 0.9$. We refer the interested reader to this
reference for the simulation details. The ensembles we use entail three
different temporal extents, $N_{t}=8,\,10$ and 12, with an aspect ratio
$N_{s}/N_{t}\approx 3$, and are set up with physical mass $u$, $d$ and $s$
quarks ($u$ and $d$ are mass-degenerate). To simulate this setup, we use
improved rooted staggered fermions with two levels of stout smearing.
Simulations are done at non-vanishing regulator $\lambda$ (pionic source) and
results are extrapolated to $\lambda=0$ using the improvement program
described in Ref. [4]. The most important observable for this proceedings
article is the isospin density $n_{I}$, for which the definition and
extrapolation details are discussed in Ref. [17]. From now on we will only use
quantities which have already been extrapolated to $\lambda=0$.

###  2.1 Extracting the EoS from the isospin density

To compute the EoS from the lattice, the main task is the computation of the
isospin density, pressure $p$ and the interaction measure $I$, which determine
the other related quantities. For energy and entropy density, $\epsilon$ and
$s$, for instance, we have

$\epsilon=I+3p\quad\text{and}\quad s=\frac{\epsilon+p-\mu_{I}n_{I}}{T}\,.$ (1)

The isospin density can be computed directly from the simulations, see Ref.
[17]. What remains is the computation of $p$ and $I$. It is convenient to
rewrite these quantities as

$p(T,\mu_{I})=p(T,0)+\Delta p(T,\mu_{I})\quad\text{and}\quad
I(T,\mu_{I})=I(T,0)+\Delta I(T,\mu_{I})\,.$ (2)

The $\mu_{I}=0$ contributions are already available in the literature, e.g.
[18, 19]. For our action results are also available in Ref. [20]. This leaves
the computation of the modifications of the EoS due to the isospin chemical
potential, $\Delta p$ and $\Delta I$. The computation of the modification of
the pressure has already been discussed in Refs. [3, 6, 16]. It can be
computed using

$\Delta p(T,\mu_{I})=\int_{0}^{\mu_{I}}d\mu\,n_{I}(T,\mu)\,.$ (3)

The computation of the interaction measure uses the relation

$\frac{I(T,\mu)}{T^{4}}=T\frac{\partial}{\partial
T}\left(\frac{p(T,\mu_{I})}{T^{4}}\right)+\frac{\mu_{I}n_{I}(T,\mu_{I})}{T^{4}}\,.$
(4)

There are two basic possibilities to compute $I$: (a) rewrite it as a
derivative of the partition function with respect to the lattice scale at non-
zero chemical potentials (see, e.g., Ref. [21]) and evaluate the relevant
observables using $\mu_{I}=0$ subtraction; (b) use a two-dimensional
interpolation of the results for $n_{I}$ to obtain the function
$n_{I}(T,\mu_{I})$ and use Eq. (4). Both have different systematics and
challenges and eventually one wants to compare the results from the two
methods. Overall we found a careful implementation of method (b) to lead to
more accurate results, which we will present in the following.

To be able to evaluate $I$, we can insert the pressure from Eq. (3), which
results in (see also Ref. [3])

$\Delta
I(T,\mu_{I})=\mu_{I}n_{I}(T,\mu_{I})+\int_{0}^{\mu_{I}}d\mu_{I}^{\prime}\Big{[}T\frac{\partial}{\partial
T}-4\Big{]}n_{I}(T,\mu_{I}^{\prime})\,.$ (5)

Thus, given a suitable differentiable and integrable interpolation for
$n_{I}(T,\mu_{I})$, we can compute the interaction measure using this
expression.

###  2.2 Interpolation of the isospin density

The remaining task is performing the interpolation, ideally in a way which is
as model-independent as possible. Note, that an interpolation of an unknown
function based on a discrete set of points with uncertainties is an ill-posed
inverse problem and, as such, does not have a unique solution. So one of the
tasks is to find a solution which is close to the actual physical solution.
This is already true for the interpolations of the EoS at $\mu_{I}=0$ and will
remain true through the necessary interpolation for $n_{I}$ at $\mu_{I}\neq
0$. Here we will try to remain as model-independent as possible by using all
possible two-dimensional spline interpolations which provide a “good”
description of the data. These spline interpolations are generated by a spline
Monte-Carlo discussed in Ref. [22]. The included weight function, as well as
the imposed boundary conditions are chosen carefully to reduce model
dependence as much as possible. In addition, the spline boundary conditions
can be used to include additional physical information on the function at
hand.

Figure 2: Two-dimensional interpolation of the isospin density on the
$N_{t}=8$ ensembles as explained in the text.

For the interpolation of $n_{I}$, we use two-dimensional cubic splines for
which we set the outer spline grid point at small $\mu_{I}$ to $\mu_{I}=0$ and
impose the physical conditions that $n_{I}$ and its second derivative are zero
on the $\mu_{I}=0$ axis ($n_{I}$ is an odd function of $\mu_{I}$). The
position of the other outer grid point in $\mu_{I}$-direction and the outer
grid points in $T$-direction is variable and we keep the second derivatives at
the outer points as free parameters (note, that all other boundary conditions
are determined by the interpolation). In the weight for the Monte-Carlo
average we use the Akaike information criterion as action (see Ref. [22] and
references therein for details) and include another term which is designed to
suppress oscillatory solutions (see section 4 in Ref. [23]).222In practice, we
have used $\mathcal{D}$ as defined in eq. (15) in Ref. [23], replacing $f$ in
the denominators by its statistical uncertainty and with $\epsilon$ set
according to a small fraction of the width between the datapoints closest to
the varied grid point. The associated free parameter for a relative weighting
of the two terms is kept as small as possible, while not leading to solutions
which oscillate strongly.

The resulting interpolation on a $N_{t}=8$ lattice is shown for a set of
temperatures in Fig. 2. The uncertainties include the uncertainty due to the
individual data points (computed using the bootstrap procedure with 1000
samples) and from the Monte-Carlo over spline interpolations.

###  2.3 Results for the EoS

Figure 3:  Results for the pressure (top-left), the interaction measure (top-
right), the energy density (bottom-left) and and the entropy density (bottom-
right). The results have been obtained on our lattice with $N_{t}=8$ from the
spline interpolation of the isospin density discussed in the main text. The
dark gray curves are the $T=0$ results from Ref. [17].

Given the spline interpolation discussed in the previous section, we can now
compute $\Delta p$ and $\Delta I$ from Eqs. (3) and (4), evaluating the
derivatives and integrals analytically. Using Eq. (1) we further evaluate the
energy and entropy densities. The results are shown for different temperatures
and the $N_{t}=8$ ensemble in Fig. 3. We also show the $T=0$ results from Ref.
[17] (dark gray curves). Evidence for the presence of the pion condensate in
the EoS is given by the characteristic behavior of the interaction measure
(top-right). At $T=0$ it initially increases starting at the BEC phase
boundary, $\mu_{I}/m_{\pi}=0.5$, before reaching a maximum at around
$\mu_{I}/m_{\pi}\approx 0.66$ from where it decreases and becomes negative at
around $\mu_{I}/m_{\pi}\approx 0.84$. Note, that for $T=0$ the interaction
measure and the other quantities vanish at $\mu_{I}=0$. These findings are in
agreement with the results from chiral perturbation at next-to-leading order
[24]. The same behavior can still be seen for $T=120$ MeV and to some extent
for $T=145$ MeV, where, however, the maximum (if it still exists) is shifted
towards larger values of $\mu_{I}$. For $T=165$ MeV the system hardly enters
the BEC phase and no sign of pion condensation is present in $\Delta I$. The
BEC phase also leads to non-monotonous behavior of the entropy density with
$\mu_{I}$ at low temperatures.

Figure 4: Comparison of the results for pressure (left) and interaction
measure (right) for different values of $N_{t}$. The bands with the solid
bounding lines are the ones for $N_{t}=8$ as shown already in Fig. 3, the ones
with the shaded bounding lines the ones for $N_{t}=10$ and the ones with
dotted bounding lines the ones for $N_{t}=12$.

Similar results are available for the $N_{t}=10$ and 12 lattices and we show a
comparison between the results for pressure and interaction measure in Fig. 4.
While for the interaction measure the bands overlap in most regions of
parameter space, the pressure shows clear signs for sizeable lattice artifacts
at $N_{t}=8$. This is not unexpected given the large lattice artifacts for the
leading order Taylor coefficient at $N_{t}=8$ observed in Ref. [25] together
with the fact that the Taylor expansion provides a good description for the
data at small $\mu_{I}$ [5] and the accumulation of differences in $n_{I}$ in
the integration to obtain $\Delta p$, see Eq. (3). Consequently, $N_{t}=8$ is
not yet within the scaling region for a leading order continuum limit, for
which $N_{t}=16$ lattices are needed for a reliable extrapolation.

##  3 The phase diagram in the $(n_{I},T)$-plane

Figure 5: Tentative continuum extrapolation for the phase diagram in the
$(n_{I},T)$-plane, as explained in the main text. The red line corresponds to
the next-to-leading chiral perturbation theory prediction at small
temperatures [26].

With the interpolation for $n_{I}$ we can also determine the phase diagram in
the $(n_{I},T)$-plane. A tentative continuum extrapolation of the phase
diagram using the $N_{t}=10$ and 12 results only (for the same reasoning as
above we do not expect $N_{t}=8$ to be in the scaling region) is shown in Fig.
5. Note, that lattice results are not available below $T=120$ MeV, but we have
matched the continuum extrapolation to next-to-leading order chiral
perturbation theory at a temperature of 30 MeV [26], shown in Fig. 5 as the
red line. We note, that the $N_{t}=10$ and 12 phase boundaries are fully
included in the uncertainties for the continuum phase boundaries. A final
continuum extrapolation will be obtained once the $N_{t}=16$ results become
available.

##  4 Conclusions

We have shown first results for the full equation of state at non-zero isospin
chemical potential from $N_{t}=8,\,10$ and 12 ensembles. The results have been
obtained from a two-dimensional spline interpolation of the isospin density
with reduced model dependence due to a spline Monte-Carlo analysis. The
interaction measure shows a characteristic behavior with $\mu_{I}$ due to the
presence of pion condensation with an initial rise before reaching a miximum
and decreasing. For small temperatures it eventually becomes negative. This
behavior is also seen in chiral perturbation at next-to-leading order [24].
Using the interpolation for the isospin density, we have also determined the
phase diagram in the $(n_{I},T)$-plane and performed a tentative continuum
extrapolation matching to next-to-leading order chiral perturbation theory at
small temperatures. For the isospin density, the $N_{t}=8$ results seem to lie
outside of the scaling region. Consequently results for $N_{t}=16$ are
required for reliable continuum extrapolations and we are currently extending
our set of ensembles accordingly.

Acknowledgements:
We thank Prabal Adhikari, Jens Oluf Andersen and Martin Mojahed for
discussions and for providing the chiral perturbation theory data from Ref.
[26]. This work has been supported by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) via the Emmy Noether Programme EN 1064/2-1
and TRR 211 – project number 315477589. The authors gratefully acknowledge the
Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this
project by providing computing time on the GCS Supercomputer SuperMUC-NG at
Leibniz Supercomputing Centre (www.lrz.de).

## References

  * [1] M. M. Wygas, I. M. Oldengott, D. Bödeker and D. J. Schwarz, _Cosmic QCD Epoch at Nonvanishing Lepton Asymmetry_ , _Phys. Rev. Lett._ 121 (2018) 201302 [1807.10815].
  * [2] M. M. Middeldorf-Wygas, I. M. Oldengott, D. Bödeker and D. J. Schwarz, _The cosmic QCD transition for large lepton flavour asymmetries_ , 2009.00036.
  * [3] V. Vovchenko, B. B. Brandt, F. Cuteri, G. Endrődi, F. Hajkarim and J. Schaffner-Bielich, _Pion Condensation in the Early Universe at Nonvanishing Lepton Flavor Asymmetry and Its Gravitational Wave Signatures_ , _Phys. Rev. Lett._ 126 (2021) 012701 [2009.02309].
  * [4] B. B. Brandt, G. Endrődi and S. Schmalzbauer, _QCD phase diagram for nonzero isospin-asymmetry_ , _Phys. Rev. D_ 97 (2018) 054514 [1712.08190].
  * [5] B. B. Brandt and G. Endrődi, _Reliability of Taylor expansions in QCD_ , _Phys. Rev. D_ 99 (2019) 014518 [1810.11045].
  * [6] B. B. Brandt, G. Endrődi and S. Schmalzbauer, _QCD at nonzero isospin asymmetry_ , _PoS_ Confinement2018 (2018) 260 [1811.06004].
  * [7] D. T. Son and M. A. Stephanov, _QCD at finite isospin density_ , _Phys. Rev. Lett._ 86 (2001) 592 [hep-ph/0005225].
  * [8] J. B. Kogut and D. K. Sinclair, _Quenched lattice QCD at finite isospin density and related theories_ , _Phys. Rev. D_ 66 (2002) 014508 [hep-lat/0201017].
  * [9] J. B. Kogut and D. K. Sinclair, _Lattice QCD at finite isospin density at zero and finite temperature_ , _Phys. Rev. D_ 66 (2002) 034505 [hep-lat/0202028].
  * [10] J. B. Kogut and D. K. Sinclair, _The Finite temperature transition for 2-flavor lattice QCD at finite isospin density_ , _Phys. Rev. D_ 70 (2004) 094501 [hep-lat/0407027].
  * [11] P. de Forcrand, M. A. Stephanov and U. Wenger, _On the phase diagram of QCD at finite isospin density_ , _PoS_ LATTICE2007 (2007) 237 [0711.0023].
  * [12] W. Detmold, K. Orginos and Z. Shi, _Lattice QCD at non-zero isospin chemical potential_ , _Phys. Rev. D_ 86 (2012) 054507 [1205.4224].
  * [13] P. Cea, L. Cosmai, M. D’Elia, A. Papa and F. Sanfilippo, _The critical line of two-flavor QCD at finite isospin or baryon densities from imaginary chemical potentials_ , _Phys. Rev. D_ 85 (2012) 094512 [1202.5700].
  * [14] B. B. Brandt, F. Cuteri, G. Endrődi and S. Schmalzbauer, _The Dirac spectrum and the BEC-BCS crossover in QCD at nonzero isospin asymmetry_ , _Particles_ 3 (2020) 80 [1912.07451].
  * [15] B. B. Brandt, F. Cuteri and G. Endrődi, _Searching for the BCS phase at nonzero isospin asymmetry_ , _PoS_ LATTICE2021 (2021) 232\.
  * [16] B. B. Brandt, G. Endrődi and S. Schmalzbauer, _QCD at finite isospin chemical potential_ , _EPJ Web Conf._ 175 (2018) 07020 [1709.10487].
  * [17] B. B. Brandt, G. Endrődi, E. S. Fraga, M. Hippert, J. Schaffner-Bielich and S. Schmalzbauer, _New class of compact stars: Pion stars_ , _Phys. Rev. D_ 98 (2018) 094510 [1802.06685].
  * [18] S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg and K. K. Szabó, _Full result for the QCD equation of state with 2+1 flavors_ , _Phys. Lett. B_ 730 (2014) 99 [1309.5258].
  * [19] HotQCD collaboration, _Equation of state in ( 2+1 )-flavor QCD_ , _Phys. Rev. D_ 90 (2014) 094503 [1407.6387].
  * [20] S. Borsányi, G. Endrődi, Z. Fodor, A. Jakovác, S. D. Katz, S. Krieg et al., _The QCD equation of state with dynamical quarks_ , _JHEP_ 11 (2010) 077 [1007.2580].
  * [21] C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann et al., _The Equation of state for two flavor QCD at nonzero chemical potential_ , _Phys. Rev. D_ 68 (2003) 014507 [hep-lat/0305007].
  * [22] B. B. Brandt and G. Endrődi, _QCD phase diagram with isospin chemical potential_ , _PoS_ LATTICE2016 (2016) 039 [1611.06758].
  * [23] G. Endrődi, _Multidimensional spline integration of scattered data_ , _Comput. Phys. Commun._ 182 (2011) 1307 [1010.2952].
  * [24] P. Adhikari and J. O. Andersen, _QCD at finite isospin density: chiral perturbation theory confronts lattice data_ , _Phys. Lett. B_ 804 (2020) 135352 [1909.01131].
  * [25] S. Borsányi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti and K. Szabó, _Fluctuations of conserved charges at finite temperature from lattice QCD_ , _JHEP_ 01 (2012) 138 [1112.4416].
  * [26] P. Adhikari, J. O. Andersen and M. A. Mojahed, _Condensates and pressure of two-flavor chiral perturbation theory at nonzero isospin and temperature_ , _Eur. Phys. J. C_ 81 (2021) 173 [2010.13655].


# Posterior Inference on Infinitely Wide Bayesian Neural Networks under
Weights with Unbounded Variance

Jorge Loría
Department of Statistics
Purdue University
West Lafayette, IN 47907, USA
<EMAIL_ADDRESS>
&Anindya Bhadra
Department of Statistics
Purdue University
West Lafayette, IN 47907, USA
<EMAIL_ADDRESS>

###### 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:

$\displaystyle z_{j}^{(l+1)}(\mathbf{x})$
$\displaystyle=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,$ (1) $\displaystyle y(\mathbf{x})$
$\displaystyle=\sum_{j=1}^{p_{K}}w_{j}^{(K)}z_{j}^{(K)}(\mathbf{x}),$ (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}\to\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.

###  1.1 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}\lvert
t\rvert^{\alpha}\\{1-i\beta\mathrm{sign}(t)\omega(t;\alpha)\\}]$, where
$\omega(t;\alpha)=\tan(\alpha\pi/2),$ for $\alpha\neq 1$ and
$\omega(t;\alpha)=-(2/\pi)\log\lvert t\rvert$, 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 Taqqu, 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.

###  1.2 Summary of main contributions

Our main contributions consist of:

  1. 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. 2.

An MCMC algorithm for posterior inference and prediction, with publicly
available code.

  3. 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. 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}),\dots,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}[\lvert\langle\mathbf{t},\mathbf{g}\rangle\rvert^{\alpha}]\right\\},$
(3)

where angle brackets denote the inner product,
$\mathbf{g}=(g(\mathbf{x}_{1}),\dots,g(\mathbf{x}_{n}))$, 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},\dots,\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,\dots,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}[\lvert\langle\mathbf{t},\mathbf{g}\rangle\rvert^{\alpha}]=\sum_{\ell=1}^{2^{n}}q_{\ell}\left\lvert\sum_{j=1}^{n}t_{j}\tau_{\ell}(\mathbf{x}_{j})\right\rvert^{\alpha}.$
(4)

###  2.1 Computation of $q_{\ell}$ and $\tau_{\ell}$

Since the size of $\Lambda$ is $2^{n}$, indicating an exponential complexity
of naïve 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},\dots,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_{jI})}$,
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,\boldsymbol{\Sigma}^{(\tau_{\ell})})$,
where the (possibly singular) variance matrix is
${\Sigma^{(\tau_{\ell})}_{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.

###  2.2 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},\dots,y_{n})$ under the model
$y_{i}=f(\mathbf{x}_{i})+\varepsilon_{i};$ where
$\varepsilon_{i}\stackrel{{\scriptstyle
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},\dots,\mathbf{x}_{n}]^{T}$. The probability
density function of $(\mathbf{y}\mid\mathbf{X})$ is:

$\displaystyle p(\mathbf{y}\mid\mathbf{X})$
$\displaystyle=(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:

$\displaystyle\mathbf{y}\mid\mathbf{X},\\{s_{\ell}\\}_{\ell=1}^{L}$
$\displaystyle\sim\mathcal{N}_{n}(0,\mathbf{Q}),\;s_{\ell}\overset{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},\dots,y_{n})$, each with respective input variables
${\mathbf{x}_{1},\dots,\mathbf{x}_{n}\in\mathbb{R}^{I}}$, under the model
$y_{i}=f(\mathbf{x}_{i})+\varepsilon_{i};\varepsilon_{i}\overset{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2}),$
and $m$ new input variable locations:
$\mathbf{x}_{1}^{*},\dots,\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},\dots,\mathbf{x}_{n}]^{T}$ and
$\mathbf{X}^{*}=[\mathbf{x}_{1}^{*},\dots,\mathbf{x}_{m}^{*}]^{T}$. The
posterior distribution at these new input variables satisfies the following
properties:

  1. 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}^{*},\mathbf{\Sigma}^{*}),$
(5)

where
$\boldsymbol{\mu}^{*}=\mathbf{Q}_{*,1:n}\mathbf{Q}_{1:n,1:n}^{-1}\mathbf{y}$,
and
$\mathbf{\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. 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. ∎

##  3 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}[\lvert X\rvert^{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.

Algorithm 1 A Metropolis–Hastings sampler for the posterior predictive
distribution

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$.

Output: 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,\dots,T$ do

Simulate $\mathbf{Q}_{t}\mid\mathbf{y},\mathbf{Q}_{t-1}$ using Algorithm S.2.

Compute $\boldsymbol{\mu}^{*}_{t}$ and $\mathbf{\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})$.

end for

return $\\{\mathbf{y}_{t}^{*}\\}_{t=1}^{T}$.

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 Matérn 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 Matérn 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.

###  4.1 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$.

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.

###  4.2 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 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, 111https://data.yorkopendata.org/dataset/diffusion-tubes-
locations 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.

Table 1: Mean absolute error of predictions in the York $\mathrm{NO}_{2}$ air quality data set, by method. | Stable | GP MLE | GP Bayes | Bayes NNet
---|---|---|---|---
MAE | 0.57 | 0.59 | 0.58 | 0.59

##  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/loriaJ/alphastableNNet.

##  Appendix A Appendix

###  A.1 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}=\lvert z\rvert^{2}$, we obtain
for $\alpha\in(0,2)$ that:

$\exp(-\nu^{\alpha/2}\lvert z\rvert^{\alpha})=\int_{0}^{\infty}\exp(-\nu
z^{2}t)p_{S^{+}}(t)dt,$ (7)

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:

$\displaystyle\phi_{\mathbf{y}}(\mathbf{t})$
$\displaystyle=\exp\left(-\sum_{\ell=1}^{L}2^{-\alpha/2}\nu^{\alpha/2}q_{\ell}\left\lvert\sum_{j=1}^{n}t_{j}\tau_{\ell}(\mathbf{x}_{j})\right\rvert^{\alpha}-\frac{1}{2}\sum_{j=1}^{n}\sigma^{2}t_{j}^{2}\right)$
$\displaystyle=\prod_{\ell=1}^{L}\exp\left(-2^{-\alpha/2}\nu^{\alpha/2}q_{\ell}\left\lvert\sum_{j=1}^{n}t_{j}\tau_{\ell}(\mathbf{x}_{j})\right\rvert^{\alpha}\right)\times\exp\left(-\frac{1}{2}\sum_{j=1}^{n}\sigma^{2}t_{j}^{2}\right)$
$\displaystyle=\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)$
$\displaystyle=\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:

$\displaystyle\phi_{\mathbf{y}}(\mathbf{t})$
$\displaystyle=\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}$
$\displaystyle=\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:

$\displaystyle p(\mathbf{y}\mid\mathbf{X})$
$\displaystyle=\int_{\mathbb{R}^{n}}\phi_{\mathbf{y}}(\mathbf{t})\exp(i\langle\mathbf{t},\mathbf{y}\rangle)\prod_{j=1}^{n}dt_{j}$
$\displaystyle=\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\langle\mathbf{t},\mathbf{y}\rangle)\prod_{j=1}^{n}dt_{j}$
$\displaystyle=\int_{(\mathbb{R}^{+})^{L}}\left\\{\int_{\mathbb{R}^{n}}\exp\left(-\frac{1}{2}\mathbf{t}^{T}\mathbf{Q}\mathbf{t}\right)\exp(i\langle\mathbf{t},\mathbf{y}\rangle)\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:

$\displaystyle p(\mathbf{y}\mid\mathbf{X})$
$\displaystyle=(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\boldsymbol{\tau}_{\ell}\boldsymbol{\tau}^{T}_{\ell}$,
denoting by $\boldsymbol{\tau}_{\ell}$ the vector with entries
$(\tau_{\ell}(\mathbf{x}_{1}),\dots,\tau_{\ell}(\mathbf{x}_{n}))$. Consider a
non-zero vector $\mathbf{w}\in\mathbb{R}^{n}$. Then:

$\displaystyle\mathbf{w}^{T}\mathbf{Q}\mathbf{w}$
$\displaystyle=\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}$
$\displaystyle=\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}^{T}_{\ell}\mathbf{w}$
$\displaystyle=\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}$
$\displaystyle>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.

## References

  * Chambers et al. (1976) J. M. Chambers, C. L. Mallows, and B. Stuck.  A method for simulating stable random variables.  _Journal of the American Statistical Association_ , 71(354):340–344, 1976.
  * Cho and Saul (2009) Y. Cho and L. Saul.  Kernel methods for deep learning.  _Advances in Neural Information Processing Systems_ , 22, 2009.
  * Cont and Tankov (2004) R. Cont and P. Tankov.  _Financial modelling with jump processes_.  Chapman & Hall/CRC financial mathematics series. Chapman & Hall/CRC, Boca Raton, Fla, 2004.  ISBN 1584884134.
  * Dancik and Dorman (2008) G. M. Dancik and K. S. Dorman.  mlegp: statistical analysis for computer models of biological systems using R.  _Bioinformatics_ , 24(17):1966, 2008.
  * Der and Lee (2005) R. Der and D. Lee.  Beyond Gaussian processes: On the distributions of infinite networks.  In _Proceedings of the 18th International Conference on Neural Information Processing Systems_ , NIPS’05, page 275–282, Cambridge, MA, USA, 2005\. MIT Press.
  * Garriga-Alonso et al. (2019) A. Garriga-Alonso, C. E. Rasmussen, and L. Aitchison.  Deep convolutional networks as shallow Gaussian Processes.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=Bklfsi0cKm.
  * Gelman (2006) A. Gelman.  Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper).  _Bayesian Analysis_ , 1:515–534, 2006.
  * Genz and Bretz (2002) A. Genz and F. Bretz.  Comparison of methods for the computation of multivariate t probabilities.  _Journal of Computational and Graphical Statistics_ , 11(4):950–971, 2002.  doi: 10.1198/106186002394.  URL https://doi.org/10.1198/106186002394.
  * Gnedenko and Kolmogorov (1954) B. V. Gnedenko and A. N. Kolmogorov.  _Limit distributions for sums of independent random variables_.  Addison-Wesley, Cambridge, Mass. Transl. from Russian by K. L. Chung, 1954\.
  * Goodfellow et al. (2016) I. Goodfellow, Y. Bengio, and A. Courville.  _Deep Learning_.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * Goodman and Pollack (1983) J. E. Goodman and R. Pollack.  Multidimensional sorting.  _SIAM Journal on Computing_ , 12(3):484–507, 1983\.  doi: 10.1137/0212032.  URL https://doi.org/10.1137/0212032.
  * Gramacy and Taddy (2010) R. B. Gramacy and M. Taddy.  Categorical inputs, sensitivity analysis, optimization and importance tempering with tgp version 2, an R package for treed Gaussian process models.  _Journal of Statistical Software_ , 33(6):1–48, 2010.  doi: 10.18637/jss.v033.i06.  URL https://www.jstatsoft.org/v33/i06/.
  * Harding (1967) E. F. Harding.  The number of partitions of a set of n points in k dimensions induced by hyperplanes.  _Proceedings of the Edinburgh Mathematical Society_ , 15(4):285–289, 1967.
  * Hornik et al. (1989) K. Hornik, M. Stinchcombe, and H. White.  Multilayer feedforward networks are universal approximators.  _Neural networks_ , 2(5):359–366, 1989.
  * Lee et al. (2018) J. Lee, J. Sohl-dickstein, J. Pennington, R. Novak, S. Schoenholz, and Y. Bahri.  Deep neural networks as Gaussian processes.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=B1EA-M-0Z.
  * Lemke et al. (2015) T. Lemke, M. Riabiz, and S. Godsill.  Fully Bayesian inference for $\alpha$-stable distributions using a poisson series representation.  _Digital Signal Processing_ , 47, 10 2015.  doi: 10.1016/j.dsp.2015.08.018.
  * Neal (1996) R. M. Neal.  Priors for infinite networks.  In _Bayesian Learning for Neural Networks_ , pages 29–53. Springer, 1996.
  * Nolan (2020) J. J. P. Nolan.  _Univariate stable distributions : models for heavy tailed data_.  Springer Series in Operations Research and Financial Engineering. Springer, Cham, Switzerland, 1st ed. 2020. edition, 2020.  ISBN 3-030-52915-0.
  * Paszke et al. (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Peluchetti et al. (2020) S. Peluchetti, S. Favaro, and S. Fortini.  Stable behaviour of infinitely wide deep neural networks.  In S. Chiappa and R. Calandra, editors, _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , volume 108 of _Proceedings of Machine Learning Research_ , pages 1137–1146. PMLR, 26–28 Aug 2020.  URL https://proceedings.mlr.press/v108/peluchetti20b.html.
  * Samorodnitsky and Taqqu (1994) G. Samorodnitsky and M. Taqqu.  _Stable non-Gaussian random processes : stochastic models with infinite variance_.  Stochastic modeling. Chapman & Hall, New York, 1994.  ISBN 0412051710.
  * Stein (1999) M. L. Stein.  _Interpolation of Spatial Data: Some Theory for Kriging_.  Springer Science & Business Media, July 1999.
  * Uchaikin and Zolotarev (1999) V. V. Uchaikin and V. M. Zolotarev.  _Chance and Stability_.  De Gruyter, Dec. 1999.  URL https://doi.org/10.1515/9783110935974.
  * Williams and Rasmussen (2006) C. K. Williams and C. E. Rasmussen.  _Gaussian processes for machine learning_ , volume 2.  MIT Press, Cambridge, MA, 2006.
  * Zhang (2004) H. Zhang.  Inconsistent Estimation and Asymptotically Equal Interpolations in Model-Based Geostatistics.  _Journal of the American Statistical Association_ , 99(465):250–261, Mar. 2004.

Supplementary Material for
_Posterior Inference on Infinitely Wide Bayesian Neural Networks under Weights
with Unbounded Variance_

##  Appendix S.1 Computation of $q_{\ell}$ and $\tau_{\ell}$ in one dimension

Assume, since we are in one dimension, that $x_{1}<\cdots<x_{n}$. In the one-
dimensional case $\Lambda$ consists of the different locations where the
change in sign can be located. This corresponds to:

  1. 1.

Before the first observation, which corresponds to $\tau(x_{k})=1$ for all
$k$, we call this $\tau_{0}$.

  2. 2.

Between $x_{j}$ and $x_{j+1}$ for some $j=1,\dots,n-1$, which corresponds to
$\tau(x_{k})=-1$ for $k<j$ and $\tau(x_{k})=+1$ otherwise, we call this
$\tau_{j}$.

  3. 3.

After $x_{n}$, $\tau(x_{k})=-1$ for all $k$, and we call this $\tau_{n}$.

Note that the first and last items correspond to linearly dependent vectors as
$\tau_{0}(x_{k})=-\tau_{n}(x_{k})$, for all $k=1,\dots,n$. Now, we compute the
probability for the first and last items:

$\displaystyle q_{\tau_{0}}+q_{\tau_{n}}$
$\displaystyle=\mathbb{P}(b_{0}+w_{1}x_{1}<0)+\mathbb{P}(b_{0}+w_{1}x_{n}>0)$
$\displaystyle=\mathbb{P}(x_{1}<-b_{0}/w_{1})+\mathbb{P}(x_{n}>-b_{0}/w_{1})$
$\displaystyle=\mathbb{P}(x_{1}<-C)+\mathbb{P}(x_{n}>-C)$
$\displaystyle=\frac{1}{2}+\frac{1}{\pi}\arctan(x_{1})+\frac{1}{2}-\frac{1}{\pi}\arctan(-x_{n})$
$\displaystyle=1+\frac{1}{\pi}(\arctan(x_{1})-\arctan(-x_{n})),$

where $C\sim\mathrm{Cauchy}(0,1)$ since $b_{0},w_{1}$ are independent standard
normal variables.

Next, the change in sign occurring between $x_{j}$ and $x_{j+1}$, yields:

$\displaystyle q_{\tau_{j}}$
$\displaystyle=\mathbb{P}(\mathrm{sign}(b_{0}+w_{1}x_{j})=-1,\;\mathrm{sign}(b_{0}+w_{1}x_{j+1})=1)$
$\displaystyle=\mathbb{P}(b_{0}+w_{1}x_{j}<0,\;b_{0}+w_{1}x_{j+1}>0)$
$\displaystyle=\mathbb{P}(x_{j}<-b_{0}/w_{1},\;x_{j+1}>-b_{0}/w_{1})$
$\displaystyle=\mathbb{P}(x_{j}<C<x_{j+1})$
$\displaystyle=\frac{\arctan(x_{j+1})-\arctan(x_{j})}{\pi},$

which corresponds to the desired probability.

##  Appendix S.2 Supporting algorithms

Algorithm S.1 is modified from Goodman and Pollack (1983) for $I=2$
dimensions, which we employ in Algorithm 1, and has a computational complexity
of $\mathcal{O}(n^{I}\log(n))$ for $n$ input points and a general input
dimension $I$. We present Algorithm S.2 to sample the latent scales
$\\{s_{\ell}\\}_{\ell=1}^{L}$ and error standard deviation $\sigma$, which
consists of a independent samples Metropolis–Hastings procedure where we
sample from the priors, and iteratively update the matrix $\mathbf{Q}$. For
computation of the density functions we use the Woodbury formula, and an
application of the matrix-determinant lemma, for an efficient update of
$\mathbf{Q}$ and to avoid computationally intensive matrix inversions.
Algorithm S.2 has computational complexity of
$\mathcal{O}(Ln^{2})=\mathcal{O}((n+m)^{I}n^{2})$.

Algorithm S.1 (Goodman and Pollack, 1983). Multidimensional sorting for $I=2$

Matrix $\mathbf{X}\in\mathbb{R}^{n\times 2}$.

Output: partition vectors $\\{\tau_{\ell}\\}_{\ell=1}^{L}$.

for $i=1,\dots,n-1$ do

for $j=i+1,\dots,n$ do

Let $u_{j}=x_{j,1}-x_{i,1}$, and $v_{j}=x_{j,2}-x_{i,2}$. If
$(u_{j},v_{j})=(0,0)$ call $j$ “good”.

Let $u_{n+j}=-u_{j},v_{n+j}=-v_{j}$, and let $m_{j}=m_{n+j}=v_{j}/u_{j}$.

end for

Sort the indices $\\{j:j\text{ is good}\\}\cup\\{n+j:j\text{ is good}\\}$ into
subsets:

  1. 1.

for those for which $u_{j}>0$, using $m_{j}$ as key

  2. 2.

for those for which $u_{j}=0$ and $v_{j}>0$

  3. 3.

for those for which $u_{j}<0$, using $m_{j}$ as key

  4. 4.

for those for which $u_{j}=0$, and $v_{j}<0$

From the sorting in Items 1 and 3 we obtain a list of subsets. Say:
$\\{J_{11},\dots,J_{1p_{1}},\dots,J_{r1},\dots,J_{rp_{r}}\\}$, where the
points with indices $J_{k1},\dots,J_{kp_{k}}$ constitute an entire subset, and
denote $J^{(k)}$ as their union, and there are $r$ subsets all together.
Denote by $k(j)$ the number of the subset within which $j$ lies.

For each $k=1,\dots,r$, let: $n_{k}=\\#\\{m:1\leq m,J_{km}\leq n\\}$, the
number of points in each ray.

For each good $j$, consider
$A_{0}^{(ij)}=\\{i\\}\cup(J^{(k(j))}-\\{n+1,\dots,2n\\})$ as the points in the
same ray as $ij$.

if $k(n+j)>k(j)$ then

Define the points in the positive side by:
$A_{+}^{(ij)}=\cup_{k=k(j)+1}^{k(n+j)-1}J^{(k)}-\\{n+1,\dots,2n\\}$, the
points in the negative side by:
$A_{-}^{(ij)}=\\{1,\dots,n\\}-A_{+}^{(ij)}-A_{0}^{(ij)}$.

else if $k(n+j)<k(j)$ then

Define the points in the positive side by:
$A_{+}^{(ij)}=\cup_{k=k(j)+1}^{r}J^{(k)}\cup\cup_{k=1}^{k(n+j)-1}J^{(k)}$, and
the points in the negative side by:
$A_{-}^{(ij)}=\\{1,\dots,n\\}-A_{+}^{(ij)}-A_{0}^{(ij)}$.

end if

For each $j=i+1,\dots,n$, if $A_{0}^{(ij)}$ has $L_{j}$ ordered items denoted
by $\\{a_{m}:m=1,\dots,L_{j}\\}$. Add $2L_{j}$ vectors:
$\boldsymbol{\tau}_{\ell}$ for $\ell=1,\dots,L_{j}$, with entry $k$ equal to
$+1$ if $k\in A_{+}^{(ij)}\cup\\{a_{m}:m=1,\dots,\ell\\}$, and $-1$ otherwise,
and the vectors $\boldsymbol{\tau}_{\ell+L_{j}}$ for $\ell=1,\dots,L_{j}$,
with entry $k$ equal to $+1$ if $k\in
A_{+}^{(ij)}\cup\\{a_{m}:m=\ell+1,\dots,L_{j}\\}$, and $-1$ otherwise.

Repeat using $A_{+}^{(ij)}$ as the negative, and $A_{-}^{(ij)}$ as the
positive.

end for

Discard repeated vectors.

return the collection of vectors $\\{\boldsymbol{\tau}_{\ell}\\}_{\ell=1}^{L}$

Algorithm S.2 A Metropolis–Hastings sampler for $\mathbf{Q}$ by simulating the
latent scales from the prior

Vector $\mathbf{y}\in\mathbb{R}^{n}$, previous latent scales
$\\{s_{\ell}\\}_{\ell=1}^{L}$, vectors
$\\{\boldsymbol{\tau}_{\ell}\\}_{\ell=1}^{L}$, partition probabilities
$\\{q_{\ell}\\}_{\ell=1}^{L}$, previous variance matrix
${\mathbf{Q}=\nu\sum_{\ell=1}^{L}q_{\ell}^{2/\alpha}s_{\ell}\boldsymbol{\tau}_{\ell}\boldsymbol{\tau}_{\ell}^{T}+\sigma^{2}\mathbf{I}}$,
and magnitude of errors $\sigma^{2}$.

Output: Updated $\mathbf{Q}$ matrix.

for $k=1,\dots,L$ do

Propose $s_{k}^{*}\sim S^{+}(\alpha/2)$.

Define $\mathbf{Q}^{(prop)}=\nu\sum_{\ell\neq
k}q_{\ell}^{2/\alpha}s_{\ell}\boldsymbol{\tau}_{\ell}\boldsymbol{\tau}_{\ell}^{T}+\nu
s_{k}^{*}q_{k}^{2/\alpha}\boldsymbol{\tau}_{k}\boldsymbol{\tau}_{k}^{T}+\sigma^{2}\mathbf{I}$
.

Accept $s_{k}^{*}$ with probability
$\min\\{p(\mathbf{y}\mid\mathbf{Q}^{(prop)}_{1:n,1:n})/p(\mathbf{y}\mid\mathbf{Q}_{1:n,1:n}),1\\}$.

If $s_{k}^{*}$ is accepted, replace $s_{k}$ by $s_{k}^{*}$.

end for

Propose $\sigma^{2}_{*}\sim\mathrm{Cauchy}^{+}(0,1)$.

Compute
${\mathbf{Q}^{(prop)}=\nu\sum_{\ell=1}^{L}q_{\ell}^{2/\alpha}s_{\ell}\boldsymbol{\tau}_{\ell}\boldsymbol{\tau}_{\ell}^{T}+\sigma^{2}_{*}\mathbf{I}}$.

Accept $\mathbf{Q}^{(prop)}$ with probability
$\min\\{p(\mathbf{y}\mid\mathbf{Q}^{(prop)}_{1:n,1:n})/p(\mathbf{y}\mid\mathbf{Q}_{1:n,1:n}),1\\}$.

return $\mathbf{Q}^{(prop)}$ if it was accepted, return $\mathbf{Q}$
otherwise.

##  Appendix S.3 Additional numerical results

We report MCMC convergence diagnostics, run times and posterior predictive
uncertainty quantification for the competing methods for the simulations in
the main paper. We also include simulation results for a variety of functions
in one and two dimensions, where we demonstrate the flexibility of the
proposed method.

###  S.3.1 MCMC diagnostics and computation times

Figure S1: Trace plots for the proposed MCMC sampler (Algorithm 1) for the simulations in Section 4, indicating good mixing in about 1000 burn-in iterations. _Left:_ one dimensional case, _Right:_ two dimensional case. Numerical results were obtained using the last 2000 iterations. Table S1: Total (in seconds) and per iteration (in milliseconds) computation times for the simulations in Section 4, for the competing methods. |  | Stable | GP MLE | GP Bayes | Bayes NNet
---|---|---|---|---|---
Total time (s) | 1-d | 24.047 | 0.067 | 18.073 | 6.91
| 2-d | 1339.543 | 0.178 | 22.185 | 7.193
Per iteration time (ms) | 1-d | 8.016 | 13.400 | 0.602 | 2.303
| 2-d | 446.514 | 35.600 | 0.740 | 2.398

###  S.3.2 Posterior uncertainty quantification

For the one-dimensional example of Section 4.1, we display the $90\%$
posterior predictive intervals for GP Bayes and Stable methods in Figure S2.
For the two dimensional example of Section 4.2, we display in Figure S3 the
$5$th, $50$th and $95$th percentiles of the posterior predictive distribution
for the same methods. In general, the credible intervals are narrower under
the Stable procedure, with better behavior at the discontinuities.

Figure S2: The point-wise $90\%$ posterior predictive intervals for GP Bayes
and Stable for the example in Section 4.1. Figure S3: Posterior predictive
quantiles for GP Bayes and Stable for the example in Section 4.2.

###  S.3.3 Additional results in one dimension

We show using a variety of one-dimensional functions that the Stable procedure
results in better performance in presence of discontinuities, while performing
similarly to GP-based methods or finite width networks for smooth functions.
Posterior uncertainty quantification results are omitted.

#### One-dimensional one-jump function.

Consider the function with a single jump given by
$f(x)=5\times\mathbf{1}_{\\{x>0\\}}$. We use forty equally-spaced observations
between $-2$ and $2$ with a Gaussian noise with standard deviation of $0.5$.
We display the obtained results on Figure S4, with optimal hyper-parameters
$\alpha^{*}=1.1$ and $\nu^{*}=1$.

Figure S4: Out-of-sample error comparison and predictions with scatter plot of
the observations for the four methods for a function with a single jump.

#### One-dimensional two-jump function.

Consider the function with two jumps given by
$f(x)=5\times\mathbf{1}_{\\{-2/3\leq x<2/3\\}}$. We use forty equally-spaced
observations between $-2$ and $2$ with a Gaussian noise with standard
deviation of $0.5$. We display the obtained results on Figure S5, with optimal
hyper-parameters $\alpha^{*}=1.3$ and $\nu^{*}=1$.

Figure S5: Out-of-sample error comparison and predictions with scatter plot of
the observations for the four methods for a function with two jumps.

#### One-dimensional piece-wise smooth.

Consider the piece-wise smooth function with a single jump, given by

$\displaystyle f(x)=\begin{cases}-2x^{2}+8,&x\geq 0,\\\
-3x+2,&x<0.\end{cases}$

We use forty equally-spaced observations between $-2$ and $2$ with a Gaussian
noise with standard deviation of $0.5$, and display the obtained results in
Figure S6, using the optimal hyper-parameters $\alpha^{*}=1$ and $\nu^{*}=1$.

Figure S6: Out-of-sample error comparison and predictions with scatter plot of
the observations for the four methods for a piece-wise smooth function.

#### One-dimensional smooth function.

Finally, consider the smooth function $f(x)=-2\cos(x)^{2}+3\tanh(x)-2x$. We
use forty equally-spaced observations between $-2$ and $2$ with a Gaussian
noise with standard deviation of $0.5$. The obtained results are shown in
Figure S7. The optimal hyper-parameters are $\alpha^{*}=1.9$ and $\nu^{*}=1$.

Figure S7: Out-of-sample error comparison and predictions with scatter plot of
the observations for the four methods for a smooth function.

###  S.3.4 Additional results in two dimensions

We show using a variety of two-dimensional functions that the Stable procedure
results in better performance in presence of discontinuities, while performing
similarly to GP-based methods or finite width networks for smooth functions.
Posterior uncertainty quantification results are available, but omitted.

#### Two-dimensional one-jump function.

Consider the function
$f(x_{1},x_{2})=5\times\mathbf{1}_{\\{x_{1}+x_{2}>0\\}}$. Using the grid of
points on $[-1,1]^{2}$ detailed in Section 4, and additive Gaussian noise with
$\sigma=0.5$, we obtain the predictions results as shown in Figure S8. Optimal
hyper-parameters are $\alpha^{*}=0.1$ and $\nu^{*}=1$.

Figure S8: Out-of-sample error comparison and predictions for the four methods
for a jump function in two-dimensions.

#### Two-dimensional smooth edge.

Consider the function
$f(x_{1},x_{2})=5\times\mathbf{1}_{\\{x_{1}^{2}+2x_{2}-0.4>0\\}}$. Note that
the jump boundary is determined by a smooth curve. Using the grid of points on
$[-1,1]^{2}$ detailed in Section 4, and additive Gaussian noise with
$\sigma=0.5$, we obtain the predictions results as shown in Figure S9. Optimal
hyper-parameters are $\alpha^{*}=1$ and $\nu^{*}=1$. The Stable method is able
to capture the smoothness of the jump boundary without losing predictive
power.

Figure S9: Out-of-sample error comparison and predictions for the four methods
for a function with a jump that is parameterized by a smooth function.

#### Two-dimensional smooth function.

Consider the smooth function $f(x_{1},x_{2})=x_{1}^{2}+x_{2}^{2}-x_{1}x_{2}$.
We use the grid of points on $[-1,1]^{2}$ detailed in Section 4, and additive
Gaussian noise with $\sigma=0.5$. Since this function is continuous, it would
be expected that the Stable method would perform similarly as the competing
methods. We obtain the predictions results as shown in Figure S10. The optimal
hyperparameters are $\alpha^{*}=1.3$ and $\nu^{*}=1$.

Figure S10: Out-of-sample error comparison and predictions for the four
methods for a smooth function in two-dimensions.

##  Appendix S.4 Additional results on York air quality data

The training and validation data sets are displayed in Figure S11.

Figure S11: Scatter plot of NO2 meeasurements in York, UK in July 2020.

Figure S12 displays the width of the posterior predictive credible intervals
under GP Bayes and Stable on the validation set. In the 60 validation data
points, the $90\%$ credible intervals from the Stable method are shorter in 50
points when compared to the GP Bayes, while 57 of the validation data points
lie inside the CIs for the Stable method. Meanwhile, GP Bayes has shorter CIs
only in 10 of the unseen points, while the CIs capture 56 of the observed
values.

Figure S12: Comparison of out of sample prediction, with $90\%$ predictive
intervals for the Stable and GP Bayes methods (centered at the median).

Finally, we display the quantiles of the posterior predictive distribution on
a regular 2-d grid over the study domain in Figure S13.

Figure S13: Posterior predictive intervals at one hundred unseen locations for
$\mathrm{NO}_{2}$ at the city of York with the Stable method. From left to
right the 5th, 50th and 95th percentiles.


Physics | 96.64 | 0.05 | 96.81 | 0.09 | 96.23 | 0.05 | 96.11 | 0.24 | 96.62 | 0.04 | 96.90 | 0.04
citeseer | 95.15 | 0.20 | 95.32 | 0.16 | 88.53 | 0.27 | 94.26 | 0.19 | 95.40 | 0.10 | 94.61 | 0.45

Table 12. The table shows micro $F_{1}$-measure in small homophilic graphs.

Large Heterophilic Graphs | GSAGE | GSAINT | LINKX | AGS-NS | AGS-GS
---|---|---|---|---|---
$\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
genius | 81.76 | 0.33 | 82.09 | 0.19 | 82.59 | 0.02 | 82.84 | 0.01 | 81.36 | 0.52
pokec | 68.91 | 0.03 | 68.18 | 0.12 | 70.57 | 0.3 | 70.08 | 0.00 | 65.77 | 0.00
arxiv-year | 47.5 | 0.15 | 40.04 | 0.32 | 49.89 | 0.18 | 50.27 | 0.10 | 38.14 | 0.00
snap-patents | 48.36 | 0.01 | 32.86 | 0.11 | 43.19 | 1.86 | 48.22 | 0.02 | 28.68 | 0.00
twitch-gamer | 61.41 | 0.00 | 61.39 | 0.32 | 59.62 | 0.18 | 61.38 | 0.05 | 60.9 | 0.00
AmazonProducts | 62.96 | 0.00 | 75.25 | 0.05 | 50.66 | 0.28 | 73.78 | 0.01 | 75.07 | 0.06
Yelp | 65.15 | 0.00 | 77.06 | 0.07 | 52.84 | 2.4 | 75.82 | 0.01 | 77.09 | 0.08

Table 13. The table shows micro $F_{1}$-measure in large heterophilic graphs.

Large Homophilic Graphs | GSAGE | GSAINT | LINKX | AGS-NS | AGS-GS
---|---|---|---|---|---
$\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
Reddit | 95.17 | 0.01 | 93.91 | 0.03 | 92.53 | 0.26 | 95.64 | 0.06 | 94.64 | 0.03
Reddit2 | 88.94 | 0.32 | 73.3 | 2.2 | 86.54 | 0.44 | 92.23 | 0.11 | 70.14 | 0.82
Reddit0.525 | 91.31 | 0.06 | 90.46 | 0.05 | 87.13 | 0.18 | 91.91 | 0.06 | 91.40 | 0.05
Reddit0.425 | 89.33 | 0.07 | 89.13 | 0.06 | 84.38 | 0.61 | 90.27 | 0.1 | 90.17 | 0.20
Reddit0.325 | 87.19 | 0.12 | 87.69 | 0.15 | 81.97 | 0.31 | 88.48 | 0.05 | 87.50 | 0.05

Table 14. The table shows micro $F_{1}$-measure in large homophilic graphs.

###  11.2. AGS-GNN vs other homophilic and heterophilic GNNs

Homophilic GNNs: Some notable homophilic GNNs from the literature are Graph
Isomorphic Network (GIN) (Xu et al., 2018a), Graph Attention Network (GAT)
(Veličković et al., 2017), Graph Convolutional Network (GCN) (Kipf and
Welling, 2016). For scaling widely used method for node sampling is GSAGE
(Hamilton et al., 2017) and for sampling based methods ClusterGCN (Chiang et
al., 2019) and GSAINT (Zeng et al., 2019).

Heterophlic GNNs: Currently, the SOTA performance in heterophilic graph is
from ACM-GCN (Luan et al., 2022). Some other notable heterophilic GNNs in the
literature are H2GCN (Zhu et al., 2020), MixHop (Abu-El-Haija et al., 2019),
GPRGNN (Chien et al., 2020), GCNII (Chen et al., 2020a), However, most of
these are not scalable for large graphs as they need to process entire graphs
for improved performance. LINKX (Lim et al., 2021) is the only scalable GNN
for heterophilic graphs. However, the input dimension of MLP depends on the
number of nodes and may cause scaling issues for very large graphs. Some other
related GNN for heterophilic graphs are jumping knowledge networks (GCNJK,
GATJK)(Xu et al., 2018b), and APPNP (Gasteiger et al., 2018). The method uses
only node features, MLP (Goodfellow et al., 2016). Uses only graph structure,
label propagation, LINK (Zheleva and Getoor, 2009), MultiLP (Zhou et al.,
2003). And simple network methods, Simple Graph Convolutional Network (SGC)
(Wu et al., 2019), C&S (Huang et al., 2020).

Table. 15 shows the performance of AGS relative to the $18$ recent algorithms
for small heterophilic graphs. We can see that ACM-GCN, AGS-NS, and LINKX are
the best-performing, with AGS-NS the best among them.

| Cornell | Texas | Wisconsin | reed98 | amherst41 | penn94 | Roman-empire | cornell5 | Squirrel
---|---|---|---|---|---|---|---|---|---
GNNs | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
GSAGE | 71.35 | 6.07 | 77.3 | 5.57 | 79.61 | 3.64 | 61.87 | 0.53 | 66.62 | 0.33 | 75.65 | 0.42 | 79.52 | 0.42 | 69.22 | 0.12 | 38.66 | 1.24
GCN | 49.19 | 4.32 | 60 | 6.26 | 55.69 | 3.64 | 60.83 | 0.53 | 57.23 | 1.91 | 61.78 | 0.35 | 42.78 | 1.02 | 57.62 | 0.11 | 28.07 | 0.79
GAT | 56.76 | 7.05 | 61.62 | 3.15 | 56.08 | 4.22 | 59.59 | 0.46 | 57.09 | 0.54 | 62.15 | 0.25 | 41.73 | 0.73 | 57.64 | 0.5 | 35.52 | 0.79
GIN | 52.43 | 6.3 | 61.08 | 4.71 | 56.08 | 6.27 | 60.1 | 0.87 | 59.82 | 1.45 | 63.18 | 0.24 | 47.34 | 0.26 | 58.74 | 0.18 | 30.55 | 0.49
GSAINT | 67.03 | 3.15 | 79.46 | 6.53 | 79.61 | 6.86 | 64.15 | 0.69 | 69.57 | 0.71 | 75.11 | 0.33 | 77.51 | 0.47 | 68.1 | 0.15 | 39.14 | 1.45
LINKX | 76.76 | 4.32 | 81.62 | 2.02 | 83.53 | 4.74 | 66.63 | 1.37 | 79.64 | 0.35 | 85.92 | 0.32 | 59.14 | 0.45 | 80.1 | 0.27 | 35.91 | 1.09
ACM-GCN | 74.59 | 1.32 | 84.32 | 7.13 | 84.31 | 3.28 | 66.11 | 1.25 | 78.12 | 0.3 | 85.38 | 0.53 | 71.42 | 0.39 | 78.43 | 0.5 | 72.06 | 2.21
LINK | 54.59 | 4.32 | 64.86 | 7.83 | 56.47 | 6.61 | 60.62 | 2.84 | 70.11 | 1.31 | 79.25 | 0.42 | 8.59 | 0.5 | 71.49 | 0.44 | 75.41 | 1.22
MLP | 70.27 | 4.83 | 67.57 | 5.92 | 80.39 | 4.47 | 41.14 | 0.9 | 51.81 | 2.75 | 71.22 | 0.57 | 65.31 | 0.32 | 61.12 | 1.45 | 31.22 | 0.88
CS | 69.73 | 5.77 | 69.73 | 5.51 | 79.22 | 6.86 | 44.66 | 3.23 | 49.71 | 3.39 | 70.86 | 0.76 | 65.42 | 0.3 | 60.12 | 0.96 | 31.01 | 0.81
SGC | 46.49 | 5.24 | 52.43 | 5.57 | 53.73 | 2.66 | 58.86 | 2.38 | 68.86 | 2.48 | 68.29 | 0.18 | 42.32 | 0.63 | 69.14 | 0.66 | 24.32 | 0.46
GPRGNN | 50.81 | 6.26 | 55.68 | 8.82 | 60.78 | 4.64 | 50.67 | 3.95 | 58.08 | 1.92 | 75.43 | 0.99 | 69.68 | 0.3 | 63.75 | 0.96 | 26.94 | 1.13
APPNP | 45.41 | 7.53 | 52.43 | 7.17 | 50.2 | 6.02 | 53.37 | 1.99 | 62.55 | 2.22 | 73.54 | 0.41 | 56.99 | 0.26 | 66.78 | 0.79 | 24.17 | 1.64
MIXHOP | 62.7 | 9.27 | 54.05 | 9.21 | 72.94 | 6.49 | 57.41 | 2.61 | 68.99 | 2.07 | 78.08 | 1.18 | 78.49 | 0.24 | 70.09 | 0.98 | 32.8 | 1.33
GCNJK | 45.95 | 7.25 | 47.03 | 6.3 | 50.2 | 2.66 | 58.34 | 3.01 | 71.41 | 2.74 | 78.33 | 1.32 | 58.32 | 0.46 | 68.85 | 0.77 | 26.13 | 0.64
GATJK | 54.59 | 3.15 | 52.97 | 7.57 | 56.86 | 4.11 | 61.35 | 4.23 | 70.29 | 1.11 | 79.47 | 0.32 | 71.54 | 0.67 | 69.7 | 0.31 | 31.2 | 0.72
LINKConcat | 76.22 | 6.26 | 74.05 | 2.76 | 81.57 | 4.04 | 59.79 | 0.96 | 76.06 | 1.47 | 83.81 | 0.28 | 12.06 | 1.61 | 76.87 | 0.33 | 67.36 | 1.57
GCNII | 60.36 | 8.92 | 54.95 | 1.27 | 71.9 | 7.22 | 57.51 | 2.77 | 70.32 | 1.48 | 80.21 | 0.4 | 74.82 | 0.26 | 71.97 | 0.28 | 32.69 | 0.84
AGS-NS | 74.59 | 2.16 | 84.86 | 6.19 | 81.96 | 4.71 | 66.74 | 1.37 | 79.19 | 0.47 | 76.06 | 0.41 | 80.49 | 0.48 | 82.84 | 0.01 | 68.24 | 0.97
| johnshopkins55 | Actor | Minesweeper | Questions | Chameleon | Tolokers | Flickr | Amazon-ratings |  |
| $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ |  |
GSAGE | 67.37 | 0.54 | 34.82 | 0.55 | 85.74 | 0.25 | 97.13 | 0.01 | 51.18 | 2.7 | 79.15 | 0.32 | 50.86 | 0.32 | 48.08 | 0.38 |  |
GCN | 61.81 | 0.95 | 29.18 | 0.78 | 80.12 | 0.06 | 97.02 | 0.00 | 41.89 | 2.15 | 78.17 | 0.01 | 49.05 | 0.06 | 43.63 | 0.43 |  |
GAT | 60.33 | 0.72 | 30.64 | 0.8 | 80.00 | 0.00 | 97.02 | 0.00 | 48.11 | 1.82 | 78.16 | 0.00 | 48.29 | 0.16 | 39.21 | 0.37 |  |
GIN | 63.4 | 0.31 | 28.09 | 0.25 | 80.74 | 0.24 | 97.1 | 0.01 | 40.13 | 2.21 | 78.29 | 0.09 | 47.2 | 0.88 | 45.85 | 0.4 |  |
GSAINT | 67.43 | 0.2 | 35.24 | 0.81 | 85.46 | 0.49 | 97.18 | 0.04 | 52.32 | 2.47 | 78.89 | 0.37 | 50.28 | 0.11 | 52.21 | 0.27 |  |
LINKX | 79.63 | 0.16 | 33.93 | 0.82 | 80.02 | 0.03 | 97.06 | 0.03 | 50.18 | 2.01 | 80.07 | 0.53 | 53.81 | 0.31 | 52.68 | 0.26 |  |
ACM-GCN | 77.37 | 0.61 | 34.42 | 1.08 | 80.33 | 0.23 | 97.02 | 0.00 | 75.81 | 1.67 | 80.45 | 0.54 | 52.19 | 0.24 | 52.94 | 0.23 |  |
LINK | 70.95 | 1.37 | 22.54 | 1.54 | 70.74 | 0.72 | 96.41 | 0.1 | 79.08 | 1.29 | 78.13 | 0.56 | 51.18 | 0.02 | 51.43 | 0.2 |  |
MLP | 52.26 | 2.31 | 33.91 | 0.7 | 80.00 | 0.00 | 97.14 | 0.06 | 40.26 | 2.96 | 78.52 | 0.14 | 46.6 | 0.07 | 41.35 | 0.46 |  |
CS | 52.61 | 2.68 | 32.41 | 0.68 | 80.00 | 0.00 | 97.14 | 0.04 | 37.87 | 1.19 | 78.34 | 0.04 | 46.6 | 0.14 | 40.61 | 0.45 |  |
SGC | 69.75 | 0.64 | 27.91 | 0.71 | 82.22 | 0.18 | 97.07 | 0.03 | 33.29 | 1.41 | 77.48 | 0.98 | 44.46 | 0.11 | 41.27 | 0.1 |  |
GPRGNN | 58.65 | 1.47 | 29.49 | 0.97 | 85.67 | 0.31 | 97.08 | 0.05 | 34.17 | 2.51 | 79.41 | 0.51 | 51.69 | 0.19 | 46.28 | 0.49 |  |
APPNP | 63.2 | 1.48 | 25.75 | 0.4 | 79.95 | 0.11 | 97.06 | 0.04 | 31.97 | 1.7 | 78.03 | 0.26 | 49.65 | 0.19 | 46.69 | 0.43 |  |
MIXHOP | 68.86 | 2.26 | 32.21 | 0.86 | 85.54 | 0.38 | 97.15 | 0.03 | 40.88 | 2.53 | 81.82 | 0.4 | 52.23 | 0.31 | 46.61 | 0.99 |  |
GCNJK | 69.21 | 0.62 | 26.14 | 0.68 | 84.86 | 0.43 | 97.1 | 0.02 | 34.43 | 0.84 | 81.51 | 0.57 | 53.17 | 0.18 | 45.18 | 0.26 |  |
GATJK | 68.32 | 1.09 | 26.42 | 0.91 | 84.66 | 0.37 | 96.96 | 0.17 | 44.25 | 1.89 | 81.29 | 0.38 | 51.35 | 0.32 | 49.41 | 0.61 |  |
LINKConcat | 74.48 | 1.22 | 28.09 | 1.65 | 62.72 | 2.5 | 96.79 | 0.17 | 69.04 | 1.97 | 77.96 | 0.91 | 46.71 | 1.73 | 51 | 0.37 |  |
GCNII | 70.53 | 0.61 | 26.29 | 0.27 | 80.92 | 0.38 | 97.07 | 0.03 | 32.68 | 1.17 | 81.8 | 0.61 | 52.74 | 0.24 | 47.74 | 0.24 |  |
AGS-NS | 78.13 | 0.42 | 36.55 | 0.93 | 85.56 | 0.28 | 97.27 | 0.04 | 73.46 | 2.29 | 80.52 | 0.41 | 51.52 | 0.13 | 53.21 | 0.46 |  |

Table 15. F1-Score of other homophilic and heterophilic GNN performance on
small heterophilic graphs (node $<100k$).

##  12\. Appendix: Ablation Study Results

###  12.1. AGS with existing GNNs

Table 16 compares the performance cof AGS with existing GNNs (GSAGE, ChebNet,
GSAINT, GIN, GAT, and GCN) evaluated on the heterophilic graphs (Reed98,
Roman-empire, Actor, Minesweeper, Tolokers). These results show that AGS with
GSAGE, ChebNet, and GSAINT performed the best. Fig. 6 in the ablation study
(§5.3) is based on the result from this table.

Methods | AGS-GSAGE | AGS-GSAINT (rw) | AGS-Chebnet | AGS-GCN | AGS-GAT | AGS-GIN
---|---|---|---|---|---|---
Dataset | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
reed98 | 63.32 | 0.76 | 55.44 | 3.04 | 64.56 | 0.84 | 56.58 | 1.20 | 55.85 | 1.33 | 59.48 | 2.00
Roman-empire | 77.35 | 0.43 | 74.39 | 0.41 | 77.44 | 0.24 | 44.36 | 1.76 | 46.57 | 0.75 | 45.91 | 0.26
Actor | 35.36 | 0.49 | 29.97 | 0.52 | 33.76 | 0.73 | 28.64 | 0.75 | 26.58 | 1.01 | 26.74 | 0.30
Minesweeper | 84.45 | 0.41 | 82.31 | 0.68 | 84.45 | 0.62 | 80.21 | 0.08 | 80.08 | 0.16 | 80.88 | 0.13
Tolokers | 78.33 | 0.14 | 78.69 | 0.39 | 78.31 | 0.35 | 78.16 | 0.00 | 78.16 | 0.00 | 78.38 | 0.11

Table 16. $F_{1}$ score of heterophilic graphs using AGS sampler with
different underlying GNNs. The best-performing models are AGS-GSAGE and AGS-
CHEB.

###  12.2. Different similarity functions for Nearest Neighbor and Submodular
optimization

Table 17 shows a performance comparison of different submodular and nearest
neighbor samplers in heterophilic graphs. Here, we generate three synthetic
graphs of average degree $200$ with strong ($0.05$), moderate ($0.25$), and
weak ($0.50$) heterophily. The underlying GNN is ChebNet, and only the
samplers are changed. Neighborhood samplers of two hops, with sizes set as
$k=[25,25]$.

| Cora1000.050.110True | Cora1000.250.110True | Cora1000.500.110True
---|---|---|---
Samplers | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
Wholegraph | 61.20 | 0.16 | 66.49 | 0.25 | 95.45 | 0.41
Random | 56.61 | 0.33 | 61.16 | 1.08 | 88.61 | 0.66
knncosine | 53.02 | 0.88 | 66.67 | 0.30 | 90.72 | 0.52
knneuclidean | 52.95 | 0.85 | 63.53 | 0.70 | 88.82 | 0.67
submodularcosine | 56.65 | 0.75 | 58.77 | 0.45 | 84.23 | 0.79
submodulareuclidean | 56.16 | 1.42 | 59.05 | 0.36 | 87.72 | 0.53
fastlink | 45.96 | 0.58 | 64.76 | 0.68 | 84.87 | 1.01
link-nn | 44.20 | 0.33 | 59.44 | 0.62 | 75.73 | 0.76
link-sub | 65.26 | 0.35 | 67.69 | 0.57 | 84.09 | 0.65
apricotfacilityeuclidean | 56.12 | 0.73 | 59.05 | 0.28 | 87.44 | 0.50
apricotfacilitycosine | 56.40 | 0.54 | 59.79 | 0.49 | 86.56 | 0.26
apricotcoverage | 55.24 | 0.71 | 59.58 | 0.81 | 85.19 | 1.02
apricotfeature | 54.99 | 0.47 | 57.07 | 0.43 | 75.49 | 0.77
apricotgraph | 56.40 | 0.54 | 59.79 | 0.49 | 86.56 | 0.26

Table 17. Performance of different sub-modular functions and nearest neighbors
on Synthetic Cora graphs.

We consider the following sampler variants.

  * •

Wholegraph: All neighbors of vertices are taken.

  * •

Random: $25$ Random neighbors are taken from each node’s first and second hop.

  * •

knncosine: Similarity sampler with $k$-Nearest Neighbor using Cosine
similarity function.

  * •

knneuclidean: $k$-Nearest Neighbor sampler with Euclidean distance function.
Note that $distance.max()-distance$ operation is performed to convert it into
similarity.

  * •

submodularcosine: submodular optimization with facility location using cosine
similarity as $\phi$ is used. The general form of the facility location
function is:

$f(X,Y)=\sum\limits_{y\in Y}\max_{x\in X}\phi(x,y)$

  * •

submodulareuclidean: submodular optimization with facility location using
Euclidean distance measures is used. The distance is also converted into
similarity.

  * •

fastlink: The learned weight from the Siamese model using regression task
directly used for sampling.

  * •

link-nn: The learned weight is used with the nearest neighbor ranking and
weight assignment.

  * •

link-sub: The learned function generates a pairwise similarity kernel for the
submodular optimization facility location function.

  * •

apricotfacilityeuclidean: This is the same as submodulareuclidean, but the
implementation and optimization operations are performed using the Apricot
library.

  * •

apricotfacilitycosine: This is the same as submodularcosine, but the
implementation and optimization operations are performed using the Apricot
library.

  * •

apricotcoverage: Maximum Coverage submodular function is used. The general
form of coverage function is:

$f(X)=\sum\limits_{d=1}^{D}\min\left(\sum\limits_{n=1}^{N}X_{i,d},1\right)$

  * •

apricotfeature: Feature-based submodular function is used. The general form of
feature-based function is,

$f(X)=\sum\limits_{d=1}^{D}\phi\left(\sum\limits_{i=1}^{N}X_{i,d}\right)$

This experiment uses the $\phi=\mathrm{sqrt}$ concave function.

  * •

apricotgraph: Graph-based submodular function is used. The general form of
Graph cut selection is:

$f(X,V)=\lambda\sum_{v\in V}\sum_{x\in X}\phi(x,v)-\sum_{x,y\in X}\phi(x,y)$