Marginalised Gaussian Processes with Nested Sampling
Authors: Fergus Simpson
Vidhi Lalchand
Carl Edward Rasmussen
Venue: NeurIPS 2021 Poster (2021)
Keywords: Gaussian Processes
nested sampling
Bayesian inference

Abstract: Gaussian Process models are a rich distribution over functions with inductive biases controlled by a kernel function. Learning occurs through optimisation of the kernel hyperparameters using the marginal likelihood as the objective. This work proposes nested sampling as a means of marginalising kernel hyperparameters,  because it is a technique that is well-suited to exploring complex, multi-modal distributions. We benchmark against Hamiltonian Monte Carlo on time-series and two-dimensional regression tasks, finding that a principled approach to quantifying hyperparameter uncertainty substantially improves the quality of prediction intervals.


Marginalised Gaussian Processes with Nested Sampling
Fergus Simpson" Secondmind Cambridge, UK fergus@secondmind.. ai
Vidhi Lalchand University of Cambridge, UK vr3080cam. ac uk
Carl E. Rasmussen
University of Cambridge, UK
cer540cam ac. uk
Abstract
Gaussian Process models are rich distribution over functions with inductive biases controlled by kernel function. Learning occurs through optimisation of the kernel hyperparameters using the marginal likelihood as objective. This work proposes nested sampling as means of marginalising kernel hyperparameters, because it is technique that is well-suited to exploring complex, multi-modal distributions. We benchmark against Hamiltonian Monte Carlo time-series and two-dimensional regression tasks, finding that a principled approach to quantifying hyperparameter uncertainty substantially improves the quality of prediction intervals.
Introduction
Gaussian processes (GPs) represent a powerful non-parametric and probabilistic framework for performing regression and classification. An important step in using GP models the specification of covariance function and in turn setting the parameters of the chosen covariance function. Both these steps are jointly referred to as the model selection problem [Rasmussen and Williams 2006). The parameters of the covariance function are called hyperparameters the latent function values take the place of parameters a GP.
The bulk of the GP literature addresses the model selection problem through maximisation of the GP marginal likelihood p(y (). This approach called MIL-IIyypically involves using gradient based opti- misation methods to yield point estimates of hyperparameters. The posterior predictive distribution is then evaluated these optimised hyperparameter point estimates. Despite its popularity this classical approach to training GPs suffers from two issues: 1) Using point estimates of hyperparameters yield overconfident predictions, by failing to account for hyperparameter uncertainty, and 2) Non-convexity of the marginal likelihood surface can lead to poor estimates located at local minima. Further, the presence of multiple modes can affect the interpretability of kernel hyperparameters. This work proposes principled treatment of model hyperparameters and assesses its impact the quality of prediction uncertainty. Marginalising over hyperparameters can also be seen as a robust approach to the model selection question in GPs.
The form of the kernel function influences the geometry of the marginal likelihood surface. For instance, periodic kernels give rise to multiple local minima functions with different periodicities can be compatible with the data. Expressive kernels which are derived by adding/multiplying together primitive kernels to encode different types of inductive biases typically have many hyperparameters, exacerbating the local minima problem.
"Equal contribution
"where 'MI. stands for marginal likelihond and 'TT' denotes the second level of the hicrarchy pertaining hyperparameters
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

The spectral mixture (SM) kernel proposed in Wilson and Adams [2013] is an expressive class of kernels derived from the spectral density reparameterisation of the kemnel using Bochner's Therorem [Bochner 1959]. The SM kernel has prior support over all stationary kernels which means it can recover sophisticated structure provided sufficient spectral components are used. Several previous works [Kom Samo and Roberts 2015 Remes et al. 2017 2018 |eettnn et al. (019) have attempted to further enhance the flexibility of spectral mixture kernels, such as the introduction of a time-dependent spectrum. However, we postulate that the key limitation in the SM kernel's performance lies not in its stationarity or expressivity, but the optimisation procedure. The form of the SM kernel gives rise to many modes in the marginal likelihood, making optimisation especially challenging. therefore presents an excellent opportunity to test out nested sampling's capabilities, as it is an inference technique which is highly effective at navigating multimodal distributions. We note that aside from being successfully applied to graphical models by Murray et al. 20066. nested sampling has been largely overlooked in the machine learning literature.
Below we provide summary of our main contributions:
Highlight some of the failure modes of MI.III training. We provide insights into the effectiveness of MI.III training in weak and strong data regimes.
Present viable set of priors for the hyperparameters of the spectral mixture kernel.
Propose the relevance of Nested Sampling (NS) as an effective means of sampling from the hyperparameter posterior [see also Faria et al., 2016 Aksulu al., 2020], particularly the
presence of a multimodal likelihood.
We demonstrate that several time series modelling tasks (where we predict the future given past observations), incorporating hyperparameter uncertainty yields superior prediction intervals.
2 Background
This section provides a brief account of marginalised Gaussian process regression, followed by an outline of the spectral mixture kernel [Wilson and Adams 2013 which we shall use in the numerical experiments.
Marginalised Gaussian Processes. Given some input-output pairs (X,y) = (x,, Bib where y
are noisy realizations of latent function values fo with Gaussian noise, = fi,,, ~ N o..). we seek to infer some as-yet unseen values y². Let kg xz), denote a positive definite kernel function parameterized with hyperparameters 0. Following the prescription of Lalchand and Rasmussen
2020), the marginalised GP framework is given by,
Hyperprior: p(0), GP Prior: f| 0 N/ Ko), Likelihood: y/f ~ (!)) (1)
where Ko denotes the N N covariance matrix, (Ko)is (x xx), The predictive distribution for unknown test inputs X* integrates over the joint posterio
p(f = /f ())(((
(2)
For Gaussian noise, the integral over f can be treated analytically as follows
p((*ly)= = p(rly 0)p(Oly)d0
p( (
E Nig, Ej). M
M
where 0 is dealt with numerically and O3 correspond to draws from the hyperparameter posterior p(fly). The distribution inside the summation represents the standard GP posterior predictive for fixed hyperparameter setting 0j. Thus, the final form of the posterior predictive in marginalised GPs approximated by a mixture of Gaussians (see Rasmussen and Williams) [2006) for a review of GPs).
Throughout this work we shall adopt a Gaussian likelihood, hence the only intractable integrand we need to consider is the hyperparameter posterior p(0ly),
p(0ly) ply (0)
(3)
'Impliitlly conditioning over inputs X, X" for compactness.

---

For non- Non-Gaussian likelihoods, one would have to approximate the joint posterior over f and 6.
In related work, Filipponee and Girolami 2014) compared range of MCMC algorithms for marginalised Gaussian processes. Murray and Adams [2010) used a slice sampling scheme for general likelihoods specifically addressing the coupling between and 0. Titsias and Låzaro-Grredilla
[2014] presented an approximation scheme using variational inference. Jang et al. [0018] considered a Lévy process over the spectral density. Most recently, Lalchand and Rasmussen marginalisation with HMC and variational inference.
2020 demonstrated
Spectral Mixture Kernels. The spectral mixture kernel is of particular interest in this work, due
to its complex likelihood surface, thereby offering a difficult challenge for sampling methods and optimisation methods alike. The kernel is characterised by spectral density S(v), defined mixture of Q Gaussian pairs [Wilson and Adams 2013]:
S(v) Σ 10 [G(,, + 2
(4)
Here the weight w; specifies the variance contributed by the ith component while G(v, [4i, oz) denotes Gaussian function with mean [4; and standard deviation o,- To avoid confusion with other standard deviations, and convey its physical significance, we shall refer to o; as the bandwidth.
The formalism of 4 is readily expanded to higher dimensions. If the power spectrum is represented by a sum of multivariate Gaussians, each with a diagonal covariance, with entries collected in the D dimensional vector og, the corresponding kernel k(r = x - x') takes the form,
(T) wy c2TT (4) II eee(p- 2 270
(5)
(d) where Tax Mi of (d)
the dth elements of the D dimensional vectors 7, My and o; respectively. The (w,, M4, ) is typically unknown, we account for this
vector of kernel hyperparameters uncertainty by treating them as random variables.
Sampling and Symmetries. Fig. 1 shows the marginal likelihood surface for 2-component SM kernel, given two datasets of different size. One of the striking features of these surfaces lies in their symmetry: this is due to the kernel's invariance to the ordering of its components. The marginal likelihood of SM kernel with spectral components possesses Q! identical regions of parameter space. A naive attempt to explore the full posterior distribution of a spectral mixture kernel would try to quantify probability mass across these degenerate regions, a much more computationally intensive task than is necessary. One solution is to only sample from one region, and ignore its symmetric counterparts. To achieve this, we adopt approach known as forced identifiability [Buscicchio ct al. 2019] to ensure that the components are defined sequential order with respect to their frequency J.
While sampling from degenerate (symmetric) modes will not improve performance, in low data regimes the modes are not always well-identified the LML surface favours several local optima. Hence, a sampling scheme that is capable sampling from multiple modes will give rise to a diversity of functions compatible with the data.
3 Applying Nested Sampling to Gaussian Processes
This section summarises our proposal for how Nested Sampling (NS) can be used to marginalise the hyperparameters of Gaussian Processes, while simultaneously yielding estimate of the model evidence.
Nested Sampling, The nested sampling algorithm was developed by Skilling 2004 (see also Skilling et al. (2006)) as a means of estimating the model evidence Z = (+)) ( ( O)d0. Here ा (0) denotes the prior, and v is the likelihood, which for our purposes is given by the GP marginal likelihood p(y|8 ()). Irrespective of the dimensionality of 6, the integral can be collapsed to a one-
dimensional integral over the unit interval, Z fo (X)dX. Here X is the quantile function
3

---

Otam 10, on 0.1
Pbrin 100,00 00
True
True
20
20
10
1.25459 x
5
5
1.23977 x10
Visualising the negative log marginal likelihood surface as function of mean frequencies in 2 component spectral mixture kernel GP. The training data was generated from latent functions evaluated on the input domain [-1.1] and refers the intrinsic data noise level. The identical peaks correspond different re-ordersngs the 2 component mean frequency vector with value [2,5].
Figure 1:
Algorithm 1: Nested Sampling for hyperparameter inference
Initialisation: Draw n. 'live' points (0)-1 from the prior 0, (0), model evidence 0. while convergence criterion unmer
Compute ; min(()(0))
(0m,)) the lowest marginal likelihood from the current
set of live points.
Sample a new live point subject to (00)) V
Remove the point corresponding to lowest marginal likelihood dyy moving set of 'saved' points
Assign estimated prior mass at this step X, = mz Assign weight for the saved point, V, X X3 Accumulate evidence, 2 1/y V
Evaluate convergence criterion, if triggered then break:
end
Add final n.. points the 'saved' list of K samples:
Each of these final points is assigned weight p: 8K/nL, K,..., ,L K // final slab of enclosed prior
Final evidence is given by, 2 = Σ +K V.VV
Importance weights for each sample are given by, p: V.V.2
return set of sumples [0.)24+* aalon with importance weights {p: 1
- and evidence estimate Z.
associated with the likelihood: it describes the prior mass lying below the likelihood value s/. The extreme values of the integrand, =0)=0) and (xx=1), therefore correspond to the minimum and maximum likelihood values found under the support of the prior π (0).
Sampling proceeds in accordance with Algorithm 1. begins with a random selection set of live points drawn from the prior. The number of live points, n,, dictates the resolution at which the likelihood surface will be explored. At each iteration, a new sample drawn from the prior and replaces the live point with the lowest likelihood, under the strict conditions that it possesses higher likelihood. By following this iterative procedure the set of live points gradually converge upon regions of higher likelihood. The prior mass enclosed by these points shrinks exponentially. By the jth iteration, the expectation of X has reduced to (X,) [ (nz + +)]} Xg & -j/N Xo- Higher values n therefore come at the cost of a slower convergence rate. While the sequence of samples provides estimate of Z, invaluable quantity in the context of Bayesian model selection, they also represent importance weighted samples the posterior.
If we wished to perform inference in low-dimensional setting, then uniform sampling strategy would suffice. However since we wish to explore higher dimensions, we employ the PolyChord algorithm HHndley et al. 2015 Handley al. 2015). which performs slice sampling [Neal 2003

---

Frequences
Bandwidths
HMC
HMC
15.0
NS
12.5
X
ML-8
Truth
10.0
1.5
7.5
5.0
10
0.5
0.0
0.0
0.0
25
5.0
7.5
10.0
12.5
15.0
Figure 2: comparison of methods for estimating the mean (left hand panel) and bandwidth (right hand panel) of two spectral components. In each case the ground truth is indicated by the magenta star, and the grey shading indicates the prior. We note that the nested sampling algorithm successfully samples from both modes, whereas HMC and ML-II can only identify single mode.
at each iteration. Unless otherwise stated, we use 100 live points, which are bounded in a set of ellipsoids (Feroz et al. 2009). These bounding surfaces allow the macroscopic structure of the likelihood contours to be traced, enabling much more efficient sampling process. This approach has proven particularly adept at navigating multi-modal likelihood surfaces [Allison and Dunkley 2014). These attractive properties have motivated numerous scientific applications, including the detection of gravitational waves (Veitch al. 2015). the categorisation of cosmic rays [Cholis al. 2015].
and the imaging of a supermassive black hole [Akiyama et al. 2019).
We implemented [De G. Matthews et al., 2017] software packages. For the former, we make use the "rslice" sampling option provided by DYNESTY, along with the default number of five slices, and adopt the default criterion for convergence. This is defined as the point at which the estimated posterior mass contained within the set of live points falls below 1% of the total, specifically log(2 (2 + 2 log(2) (2) < 0.01,
the the above algorithm via combination of the DYNESTY [Speagle 2020) and GPflow
where AZ Z = Vimax
Hyperpriors for the Spectral Mixture Kernel
In this section we outline a physically motivated set of priors for the hyperparameters which govern the spectral mixture kernel. As defined in (4). the three fundamental parameters of each spectral component are the mean frequency μ. bandwidth a, and weight w. For most of these hyperparameters, we find it is sufficient to impose a weakly informative prior of the form (a, w, a LogNormal(0. 2) where of is the data noise variance. However the behaviour of the spectral component's 's charactcristic frequency μ merits closer attention.
Two properties of the data strongly influence our perspective on which frequencies we can expect to observe: the fundamental frequency and the highest observable frequency. The fundamental frequency VF is the lowest frequency observable within the data, and is given by the inverse of the interval spanned by the observed x locations. Meanwhile the maximum frequency 1/x represents the highest observable frequency. For gridded data, this naturally determined by the Nyquist frequency, which is half the sampling frequency.
is crucial bear in mind that the spectral density we wish model is that of the underlying process, and not the spectral density of the data. These two quantities are often very different, due to the limited scope of the observations. For example, the change in stock prices over the period of single day cannot exhibit frequencies above 10 HHz. Yet the process will have received contributions from long term fluctuations, such as those due macroeconomic factors, whose periods can span many years. If we make no assumption regarding the relationship between the process we wish to model and the finite range over which it observed, then a priori, some as-y as-yet undiscovered frequency within the process ought to be considered equally likely to lie above below the fundamental frequency.
'Code available https:: htops:://gittee //githubb..ttttt com/ www.///gggsmssooo/nsaaoaaa
5

---

train
train
(X260)
7.5
HMC
Nested
NLPD
NLPD
5.0
2.5
-2
0.0
0.01
0.1
0.5
0.01
0.1
0.5
0.01
0.1
0.5
Naise level
Naise level
Noise level
Figure 3: The impact of marginalising hyperparameters on the predictive performance (NLPD) a spectral mixture kernel. The error bars at cach noise level correspond to the mean NLPD and standard error of the mean across diversity of latent functions. Nested sampling significantly outperforms ML-II in all nine experiments. The subplots contain the performance summaries under the different configurations. Left: Sparse distribution training data (n frain = 10). Middle: Moderate sized distribution training data (n_fraate = 30), Right: Dense distribution of training data (r_train = 100) all on a fixed domain [-1,+1].
Furthermore, given the large initial uncertainty in frequency 14, it is appropriate to adopt a prior which spans many orders of magnitude.
Towards very low frequencies, μ < vr, sinusoid contributes very little variance to the observations an annual modulation makes tiny contribution to minute's worth of data. As we consider frequencies much lower than the fundamental frequency, therefore becomes less likely that we will detect their contributions. We model this suppressed probability of observation with broad Gaussian log frequency for the regime lr < up. Meanwhile, at frequencies above the Nyquist frequency, lr UF, we encounter degenerate behaviour: these sinusoids are indistinguishable from their counterparts at lower frequencies: they are said to be aliases of each other. As result this aliasing behaviour, the likelihood surface is littered with degenerate modes with identical likelihood values. From a computational perspective, it is advantageous to restrict our frequency parameter to a much narrower range than is permitted by our prior, while maintaining the same probability mass. As illustrated in the supplementary, mapping these higher frequencies down to their corresponding alias μ vN yields uniform prior μ.
for < VF Uniform(1, for uF P VN.
(0g
n/ve
(6)
5 Experimental Results
In this section we present results from series of experiments on synthetic data one and two dimensions, as well as real world time series data.
Baselines. For all of experiments, the key benchmark we compare nested sampling is the conventional MI-II method. This involves maximizing L(0) = log p(y(()) over the kernel hyperparameters, 0, argmax, L(0). In addition ML-II, it is also instructive to compare nested sampling against an MCMC method for marginalising the hyperparameters. Hamiltonian Monte Carlo [Duane al. 1987 is fundamental tool for inference intractable Bayesian models. It relies upon gradient information to suppress random walk behaviour inherent to samplers like Metropolis- Hastings and its variants [for a detailed tutorial, see Neal al. 2011]. In the experiments we use a self-tuning variant of HMC called the No-U-Turn Sampler [NUTS, Hoffman and Gelman 2014] where the path length is adapted for every iteration. NUTS is frequently shown to work as well as hand-tuned HMC; hence in this way we avoid the compute overhead in tuning for good values of the step-size (t) and path length (I.). We use the version of NUTS available in the python package pymc3.
Synthetic data. As a first demonstration of how the different methods perform, we draw four samples from a two-component SM kernel with known weights, frequencies and bandwidths. For each latent function sample we construct noisy training data across three noise levels [0.01, 0.1, 0.5) and three training sizes [10, 30, 100). the quality of predictions and prediction uncertainty under weakly identified regimes characterised by very few data points in fixed domain to strongly identified regimes characterized by dense
fixed input domain [-1, 1]. In this way we seek to quantify

---

Ground
7,5
5.
5.0
3
2.5
225
-2.
5500
-5.0
-75
-7.5
-10.0
-10
-10
-5
--
7.5
7.5
5.0
5.0
-
225
-25
5500
-5.0
-55
-7.5
Figure 4: Learning a two-dimensional pattern with 10 component spectral mixture kernel. Top Left: The ground truth, and the black dots denote the locations of the 50 training points. Top Right: Posterior predictive under ML-II. Bottom Left: Posterior predictive under HMC. Bottom Right: Posterior predictive mean under nested sampling.
distribution of data in fixed domain. Further, the intrinsic noise level of of the data also impact inference in weakly and strongly identified data regimes. In order to analyse the impact of o and training size calculate the average performance across each of three different noise levels for each training set size.
We train under each of the three inference methods (ML-II, HMC, Nested) for each of the 9 x 4 data sets created and report prediction performance terms of the average negative log predictive density (NLPD) in Figure 3 ML-II uses five random restarts with initialisation protocol tied to the training data. Following protocols from Wilson and Adams [2013), the SM weights (wy) were initialised the standard deviation of the targets y scaled by the number of components (Q = 2).
The SM bandwidths (a;) were initialised to points randomly drawn from a truncated Gaussian (W(0, max d(æ., x²)? where max d(xx, e') is the maximum distance between training points and mean frequencies (14:) were drawn from Uni vn) ) to bound against degenerate frequencies. HMC used LogNormal 2) priors for all the hyperparameters. The ML-II experiments used gpytorch while the HMC experiments used the NUTS sampler in pymc3.
The advantage of synethetic data is we know the values ought to recover. Figure compares the three inference methods in their ability to recover the true hyperparameters for a synthetic dataset. This involved 100 datapoints, a noise amplitude of o, 0.1 and signal-to sgnal-to-noise ratio of x 3.2 on a fixed domain [-1,1]. The magenta star * denotes the true value and the red cross x denotes the ML-II estimate. The HMC and NS sampling schemes are both better at recovering the ground truth hyperparameters than the point estimates. Of particular significance, the nested sampling scheme is able to simultaneously sample from both modes inherent the marginal likelihood.
Figure we look the impact on predictive performance, and notice that ML-II struggles with small training data sets, catastrophically underestimating the noise level (in the left hand panel the bars are rescaled for visibility). The two sampling methods, HMC NS, perform comparably well when there are only 10 datapoints, but NS consistently outperforms HMC when there are 100 datapoints,

---

200
HIAT
35.0
230
12.5
10.0
750
7.5
200
150
2.5
2.3
50
1400
1405
7440
14000 1447 1985
1035
1640
1907
1450
7475
1825
1967
3445
1923
5005
systematic underestimation of predictive uncertainty when adopting point estimates of kernel hyperparameters. Here we show 95% test confidence intervals derived from the spectral mixture kernel. Black lines denote training data and the black dots denote test data. Left: 'BIRTHS' data set the ML-II intervals capture 84% of the test data where HMC and NS capture 93% and 99% of the test points respectively. Right: 'RADIO' data set (also showing performance however, the average predictive density under test higher under sampling methods.
Figure An illustration
prediction) the three methods show similar
across all noise levels. This may be due to its inherent ability to sample from multiple modes, we saw in Figure
Two-dimensional Pattern Extrapolation. To demonstrate how the inference methods adapt to higher dimensional problems, we revisit the challenge presented in Sun et al [26]. The two- dimensional ground truth function is given by y = (cos 2r1 x cos 2rz / 2 This experiment marks significant increase in the dimensionality of parameters space, since in two dimensions each spectral component has five degrees of freedom. We train with 50 points chosen at random across [-6, +6] in the zy-dimain. The test points are defined 20 20 grid.
In Figure 4 we visualise the mean values of the posterior predictive distribution from three different inference methods: ML-II, HMC and nested sampling (NS). Visually, the reconstruction under the marginalised GP schemes (HMC/ NS) appears be superior to ML-II. Further, the 95% confidence intervals (not visualised) differ markedly,
is evident from the NLPD values. These given by 216, 2.56, and 2.62 for ML-II, HMC and NS respectively (lower better). For reference we also trained the neural kernel network (Sun et al [26]), which achieved NLPD of 3.8. Marginalised Gaussian processes comfortably outperform competing methods. ML-II was trained with Adam (learning rate=0.05) with 10 restarts and 2,000 iterations.
Time series benchmarks. To test of our inference methods on a selection of real world data, we evaluate their predictive performance against the thirteen 'industry' benchmark time series, as used by Lloyd et al. 22014] The time series are of variable length, with up to 1, 000 data points in each.
"The raw data available
http:///hwc..ttooooooooooo/oo//a renamee
Table NLPD values for various GP methods across a range time series tasks.
KERNEL
RBF
RBF
RBF
SPECTRAL
SPECTRAL
SPECTRAL
NKN
INFERENCE
ML-II
HMC
NS
ML-II
HMC
NS
ML-II
AIRLINE
13
6.37
11.1
7.25
5.83
5.22
5.6
BIRTHS
5.27
5.26
5.26
5.17
5.23
4.96
5.42
CALL CENTRE
9.21
8.44
8.67
7.32
7.29
7.76
PRODUCTION
12.39
18.8
11
12.4
INTERNET
11.2
11.35
11.4
11.3
11.4
13.2
12.6
MAUNA
2.53
2.64
2.53
1.5
3.32
1.8
3,4
RADIO
2.47
2.45
2.46
2.19
2.12
2.07
4.12
SOLAR
1.69
1.03
1.74
1.4
0.82
0.58
2.38
SULPHURIC
5.17
5.15
5.16
5.13
5.15
5.11
2.79
2.81
2.8
2.8
2.49
2.52
4.2
MEET
11.9
11.03
11.3
12.8
11.2
10.6
8.5
WAGES
57.1
30.20
33.4
159
8.71
6.59
4.28
WHEAT
9.8
8.61
9.73
8.47
6.44
6.58
6.24
MEAN
11.31
8.29
10.2
18.7
6.34
5.96
6.40
11.2
10.03
10.03
8

---

Time (minutess)
10a-
HMC
30
Nested
20
10.
10
10
airline
radio
call
wheat
temp
mauua
salar
supp
untmp.
births
wages
intemet
mean
NKN
Figure 6: Left: Training time for each of the time series benchmarks, Right: Average training time across all data sets for the three inference methods we consider. Perhaps surprisingly, nested sampling is not substantially slower than NKN which is trained using MI.-II.
All models are trained normalised data, but the quoted test performance is evaluated in the original space. The output is scaled to zero mean and unit variance, while the imput space is scaled to span the unit interval, thereby ensuring uF 1. Given the complex nature of the task, we utilise more spectral components than for the synthetic case. Our fiducial kernel has seven components (Q =7) 7), yielding 22-dimensional hyperparameter space to explore. For reference, we also include results from the Neural Kernel Network (Sun et al. 2018) where a flexible kernel is learnt incrementally through a weighted composition primitive base kernels. This trained with the Adam optimiser for 100,000 iterations and a learning rate of 10-3.
In Table I we report the negative log predictive density (NLPD). The evaluation conducted with a 60/40 train/test split. Each set of experiments was repeated three times with a different random seed (and used to generate the quoted uncertainties. For ML-II, each run included initialisation protocol involving evaluating the marginal likelihood at one thousand points in the hyperparameter space, choosing the initial setting the be the point with the highest marginal likelihood value.
We find that the spectral mixture kernel exhibits significant performance gains when using a sampling- based inference method compared to the conventional ML-II approach. The high average NLPD score for the 'WAGES' data-set underscores the lack robustness associated with the ML-II approximation. We note that excluding this outlying score does not change the ranking of the methods. NS outperforms HMC on 11 of the 13 tasks, and also carries an advantage of faster evaluation times, and it provides estimate of the model evidence (which could potentially be used to advise how many spectral components are necessary). While NS excels with the spectral mixture kernel, it cannot offer similar gains with the simpler RBF model. We again interpret this as being due to its ability to traverse the multi-modal likelihood surface better than other methods. Examples of the posterior distributions are illustrated in Figure while more detailed results for various configurations of the nested sampler can be found in the Appendix.
While Nested Sampling is lesser known inference technique, it is one which has gained dramatic increase in usage within the physics community over the past decade. We note that much the same could have been said for HMC in the early 1990s. The key advantage of NS is that the location of the next sample is not governed by the single previous point (as in MCMC), but by whole clouds of 'live points'. Hence it readily samples across multiple modes in the marginal likelihood surface, as can be seen in Figure 2 The hyperparameter samples obtained from NS induce a greater diversity of functions than those from HMC. This also sheds light on the question about the uncertainty intervals, NS samples give wider intervals in prediction tasks as their samples are more diverse, and these more rigorous uncertainty estimates led to improved test performance.
Computation. Figure depicts the wall clock time required for training during our time-series experiments, all of which utilised single Nvidia GTX1070 GPU. For HMC this corresponds to single chain with 500 warm-up iterations and 500 samples. Across these benchmarks, HMC was found to be slower by approximately factor of two compared to nested sampling. While sampling hyperparameters will add some computational overhead, there are many means by which the marginal likelihood could be computed more rapidly. This includes inducing point methods (sparse GPs) and structure exploiting Kronecker methods when the data lie a grid. Overall, the nested sampling approach was not found to be significantly slower than sophisticated kernel learning schemes like the NKN kernel learning method with ML-II inference.

---

What are the potential societal impacts? By developing inference methods with more reliable quantification uncertainty. this may lead to safer decision making. On the other hand, greater computational overheads can be accompanied with greater carbon footprint.
Consequences.
6 Discussion
We present a first application of nested sampling to Gaussian processes, previously only performed in the context of astronomy. It offers a promising alternative to established GP inference methods, given its inherent capability of exploring multimodal surfaces, especially since modern incarnations such as PolyChord (Handley et al. 2015 enable higher dimensional spaces (D S 50) to be explored.
The spectral mixture kernel possesses a challenging likelihood surface, which made it an excellent target to showcase the versatility of nested sampling. We have performed the first demonstration of marginalised spectral mixture kernels, finding significant performance gains are readily available when using either HMC or NS, relative to the conventional approach of using point estimates for hyperparameters. The benefits of performing full GP marginalisation are three-fold. Firstly, defends against overfitting of hyperparameters, which may occur when their uncertainty is substantial. Secondly, there is lower susceptibility to getting trapped in sub-optimal local minima. Thirdly, they provide more robust prediction intervals by accounting for hyperparameter uncertainty. These three advantages are most pronounced when many hyperparameters are present, or the data is noisy.
While a pathological marginal likelihood geometry can pose problems for both gradient based optimisation and sampling; sampling schemes are able to quantify them better if the practical difficulties of deploying them are overcome. The nested sampling scheme, which does not require gradient information, provides an opportunity to sample the multimodal posterior at fraction of the cost of running HMC. Further, it is crucial to ask if we are compromising the full expressivity of kernels by resorting to point estimates in the learning exercise? A marginalised GP with a standard SM kernel was found to improve considerably upon its conventional counterpart. We also note that the methods discussed here could be generalised to operate in sparse variational GP framework when dealing with larger datasets. Future work will focus in this direction.
Acknowledgments and Disclosure of Funding
The authors would like to thank João Faria, Joshua Albert, Elvijs Sarkans, and the anonymous reviewers for helpful comments. VL acknowledges funding from the Alan Turing Institute and Qualcomm Innovation Fellowship (Europe).
References
Carl Edward Rasmussen and Christopher KI Williams. Gaussian processes in machine learning. Springer, 2006.
Andrew Wilson and Ryan Adams. Gaussian process kernels for pattern discovery and extrapolation. In International Conference on Machine Learning, pages 1067-1075, 2013.
Salomon Bochner. Lectures on Fourier integrals, volume 42. Princeton University Press, 1959.
Yves-Laurent Kom Samo and Stephen Roberts. Generalized spectral kernels. arXiv preprint arXiv: 120..0236..2236 2015.
Sami Remes, Markus Heinonen, and Samuel Kaski. Non-stationary spectral kernels. In Advances in Neural Information Processing Systems, pages 4642-4651, 2017.
Sami Remes, Markus Heinonen, and Samuel Kaski. Neural non-stationary spectral kernel. arXiv preprint 1r11v:1811..000 2018.
Gregory W Benton, Wesley J Maddox, Jayson P Salkey, Julio Albinati, and Andrew Gordon Wilson. Function-space distributions over kernels. arXiv preprint arXiv:1910.13:: 19101.990...55. 2019.
lain Murray, David JC. MacKay, Zoubin Ghahramani, and John Skilling. Nested sampling for potts models. Advances in Neural Information Processing Systems, 18:947, 2006.
10

---

JP Faria, RD Haywood, BJ Brewer, P Figueira, M Oshagh, A Santerne, and NC Santos. Uncovering the planets and stellar activity of corot-7 using only radial velocities. Astronomy Astrrphysics. A88:A31, 2016.
Mehmet Deniz Aksulu, RAMJ Wijers, HJ van Eerten, and AJ van der Horst. A new approach to modelling gamma-ray burst afterglows: using gaussian processes to account for the systematics. Monthly Notices of the Royal Astronomical Society, 497(4):4672 461444662,266,,,,, 2020.
Vidhi Lalchand and Carl Edward Rasmussen. Approximate inference for fully bayesian gaussian process regression. volume 1 118 of Proceedings of Machine Learning Research. PMLR, 08 Dec 2020.
Maurizio Filippone and Mark Girolami. Pseudo-marginal Bayesian inference for Gaussian processes. IEEE Transactions on Pattern Analysis and Machine Intelligence, 36(11):2214-2226, 2014.
lain Murray and Ryan P Adams. Slice sampling covariance hyperparameters of latent gaussian models. In Advances in neural information processing systems, pages 1732-1740, 2010.
Michalis Titsias and Miguel Lázaro-Gredilla. Doubly stochastic variational bayes for non-conjugate inference. In International conference on machine learning, pages 1971 2014.
Phillip A Jang, Andrew E Loeb, Matthew B Davidow, and Andrew Gordon Wilson. Scalable 1\ 'evv 2018.
process priors for spectral kernel learning. arXiv preprint
Riccardo Buscicchio, Elinore Roebber, Janna M. Goldstein, and Christopher Moore. Label switching problem in bayesian analysis for gravitational wave astronomy. Physical Review D, 100(8), Oct 2019. ISSN 2470-0029, doi: 10.1103/physeevd. wttp:seeeed..0......00 100.084041. URL http: doi
100. 084041
John Skilling. Nested Sampling. In Rainer Fischer, Roland Preuss, and Udo Von Toussaint, editors, American Institute of Physics Conference Series, volume 735 of American Institute of Physics Conference Series, pages 395-405, November 2004. doi: 10.
John Skilling et al. Nested sampling for general bayesian computation. Bayesian analysis, 1(4): 833-859, 2006.
W. J. Handley, M. Hobson, and A. N. Lasenby. polychord: nested sampling for cosmology. MNRAS, 450:L61-L65, June 2015. doi: 10.1093/mmras///sh0044.
WJ Handley, MP Hobson, and AN Lasenby. Polychord: next-generation nested sampling. Monthly Notices of the Royal Astronomical Society,453(4)))))) 453(4):4384-4398, 2015.
Radford M. Neal. Slice sampling. Ann. Statist., 31(3):705-767, 06 2003. doi: 10.1214/aos/ 1056562461. Utt.httss:://o/ wwtoooo////o 10. 1101111111//////0555555555
F. Feroz, M. P. Hobson, and M. Bridges. MULTINEST: an efficient and robust Bayesian inference tool for cosmology and particle physics. MNRAS, 398(4):1601-1114 October 2009. doi: 10.1111/
j.1365-2966.209) 14548.x.
Rupert Allison and Joanna Dunkley. Comparison of sampling techniques for bayesian parameter estimation. Monthly Notices of the Royal Astronomical Society, 437(4):3918 3928, 2014.
John Veitch, Vivien Raymond, Benjamin Farr, W Farr, Philip Graff, Salvatore Vitale, Ben Aylott, Kent Blackburn, Nelson Christensen, Michael Coughlin, et al. Parameter estimation for compact binaries with ground- based gravitational-war wave observations using the lalinference software library. Physical Review D, (44:00420,,, 2015.
Ilias Cholis, Carmelo Evoli, Francesca Calore, Tim Linden, Christoph Weniger, and Dan Hooper. The galactic center gev excess from series of leptonic cosmic-ray outbursts. Journal of Cosmology and Astroparticle Physics, 2015(12) 2015.
Kazunori Akiyama, Antxon Alberdi, Walter Alef, Keiichi Asada, Rebecca Azulay, Anne-Kathrin Baczko, David Ball, Mislav Baloković, John Barrett, Dan Bintley, et al. First m87 event horizon telescope results. iv. imaging the central supermassive black hole. The Astrophysical Journal Letters, L7, 2019.
11

---

Joshua S. Speagle. DYNESTY: a dynamic nested sampling package for estimating Bayesian posteriors and evidences. MNRAS, 49333)3132-33188, April 2020. doi: 10. 10.1093/mmraa//tuu2275
Alexander De G. Matthews, Mark Van Der Wilk, Tom Nickson, Keisuke Fujii, Alexis Boukouvalas, Pablo León- Villagrá, Zoubin Ghahramani, and James Hensman. Gpflow: A gaussian process library using tensorflow. The Journal of Machine Learning Research, 18(1):1299- 1304, 2017.
Simon Duane, Anthony D Kennedy, Brian Pendleton, and Duncan Roweth. Hybrid monte carlo. Physics letters B, 195(2)211-2222, 1987.
Radford M Neal al. MCMC using hamiltonian dynamies. Handbook of Markov chain Monte Carlo, 2(11):2, 2011.
Matthew D Hoffman and Andrew Gelman. The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1):1593-162, 2014.
James Robert Lloyd, David Duvenaud, Roger Grosse, Joshua Tenenbaum, and Zoubin Ghahramani. Automatic construction and natural- language description of nonparametrie regression models. In Tienty-eighth AAAI conference on artificial intelligence. 2014.
Shengyang Sun, Guodong Zhang, Chaoqi Wang, Wenyuan Zeng, Jiaman Li, and Roger Grosse. Dif- ferentiable compositional kernel learning for gaussian processes. arXiv preprint arXiv:1806.04326. 2018.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the pa- per's contributions and scope? [Yes] We list our key contributions at the end of the Introduction.
(b) Did you deseribe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [Yes] Final paragraph of Section 5.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of theoretical results? [N/A] Our work mostly empirical in nature.
(b) Did you include complete proofs of all theoretical results? [N/A] Our work is mostly empirical in nature.
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL)? [No] The code we used to run the time series benchmarks proprietary.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] Details can be found in Sections 4 and 5.
(c) Did you report error bars (e.g., with respect the random seed after running experi- ments multiple times)? [Yes] Table 1 includes uncertainties based upon three different runs.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See the final paragraph of Section
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] We cite papers associated with GPflow and Dynesty in Section 3.
(b) Did you mention the license of the assets? [No]
12

---

(c) Did you include any new assets either in the supplemental material or as URL? [No]
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [No]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [No]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A] We did not make use of human subjects.
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "znLeJ4qdXDC", "reviewer": "Reviewer_XcBt", "summary": "The paper suggests using the nested sampling algorithm to marginalize hyperparameters of Gaussian processes with spectral mixture kernels. The motivation is that the posterior over hyperparameters of such kernels is highly multimodal and nested sampling is known to work well with multimodal distributions. Furthermore, practical priors for the parameters of spectral mixture kernels are suggested. The proposed nested sampling driven marginalised Gaussian process is compared to NUTS driven marginalised Gaussian process and the Gaussian process with hyperparameters estimated by maximum likelihood on a set of low-dimensional benchmark problems. The comparison is in terms of negative log predictive density. It is shown that the proposed method usually outperforms other methods while also performing faster than NUTS.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper highlights the method of nested sampling, popular in physics, as a method for marginalizing over hyperparameters of spectral mixture kernels. This is a valuable contribution. The conducted empirical study is convincing in showing that the method is indeed promising in this new setting. The method is only benchmarked in low dimensions but may be useful in the time series modeling problems which are quite broad and important themselves.\n\n# Originality\nAs far as I can tell, the contribution is novel.\n\n# Quality\nThe claims are supported by an extensive empirical evaluation. On the downside, the implementation is not provided and said to be proprietary.\n\n# Clarity\nThe paper is very well written. The most serious problem of the sort spotted by me in the paper is fonts being too small in some figures (especially in Figure 2).\n\n# Significance\nThe method suggested by the paper and the corresponding empirical study will be valuable for machine learning practitioners who use low-dimensional Gaussian processes. Mostly, this work is relevant for time series modeling. Regrettably, this positive impact is somewhat limited by the fact that the code will not be openly released for people to be able to quickly and easily apply the method to their problems.\n\n# Remarks\n- Line 17. I suggest giving a reference or briefly discussing the abbreviation ML-II, it might not be a universally known one.\n- Line 39. “lies not its”->”lies not in its”.\n- Line 44. “some the”-> “some”.\n- Lines 67-68. Phrase “the final form of the posterior predictive in a marginalised GOs is a mixture of Gaussians” is not rigorous because it is only approximately so, I suggest reflecting this in text.\n- Line 75. “consider”->“considered”.\n- Line 76. “demonstrate”->”demonstrated”.\n- Line 79. “complet”->”complex”.\n- All figures except Figure 3 need larger fonts.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations and societal impact are adequately addressed.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "7"}, "extra_text": {}}, {"review_id": "gO4uFCPVOwB", "reviewer": "Reviewer_whBU", "summary": "The paper proposes using nested sampling to have a fully Bayesian GP. The idea is illustrated with spectral mixture kernel, and supported (to some extend) with experiments.", "questions": "", "limitations": "", "rating": 6, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This is a paper to bring the technical of nested sampling that is well-used in physics and cosmology to the machine learning GP community. The paper is not very original and significant, given that nested sampling has been introduced previously in N(eur)IPS [Murray, 2005], and the authors of that paper have very deep knowledge of GP and are experts in this area. It is also widely accepted that margalised (or Bayesian) GP is *the correct* thing to to, though it can be somewhat challenging to execute in practice.\n\nNevertheless, this is a worthwhile paper, (re-)illustrating the significance of  Bayesian GP, esp in the context of the recently introduced spectral mixture kernel (SKM) which has multiple modes as illustrated by the authors. Also useful is to bring to attention the DYNESTY package. If the authors were to release their code, I can see that the code will be used frequently.\n\nSection 4 is a particularly helpful advice for anyone attempting to use SKM. The paper is also well-written.\n\n# Detailed comments\n1. I prefer the term \"Bayesian Gaussian Processes\" rather than \"Marginalised Gaussian Processes\" in the title, but this is a matter of personal preference.\n2. Given that the study is limited to SKM, this should be in the title.\n3. Reorganised the paper to put SMK material in section 2 together with section 4.\n4. It is paramount to cite and discuss [Murray, 2005].\n5. I beg to differ from the sentence from line 232  to 233. The pattern for NS is too clear on along the axis in Figure 4, and then it blurs out along the edges. This is not what the generating function gives. It is nevertheless reassuring that all these are better than NKN.\n6. Lines 262 to 263 tries to give something negative a positive light. What I read is that for simpler kernels/covariance-functions like the RBF/squared-exponential that are typically used in practice, there is no requirement for NS.\n7. Lines  284 to 286 requires a correct heading.\n\nIain Murray, David MacKay, Zoubin Ghahramani, John Skilling. Nested sampling for Potts models. Part of Advances in Neural Information Processing Systems 18 (NIPS 2005)\n\nChanges after author response\n=======================\nRemoved the comment on NLPD from point 6.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "See comment 6 above. The paper should say that while they strongly advice to use NS for SMK, the conclusion for simpler kernels like RBF/squared-exponential with limited number of hyperparameters required further study.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "_sX2v5REA2O", "reviewer": "Reviewer_1r2q", "summary": "The paper proposes nesting sampling techique for sampling from GP hyperparameter posteriors. Specifically, obtaining hyperparameters of the spectral mixture kernel is  discussed. Empirical results show that nesting sampling benefits GP performance while being computationally cheaper than HMC.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper addresses a problem of marginalizing GP hyperparameters. The proposed technique is nested sampling of the hyperparameters. Nested sampling (NS henceforth) was developed previuosly, however the application of the technique to GPs is original work. The authors thoroughly cite related work.\n\nThe work claims that NS improves GP behavior compared to a point estimate of the hyperparameters. On the other hand, NS takes significantly less computation time than exploring the posterior with HMC. Performance comparison between NS and HMC seems to depend largely on the kernel and the dataset. There is also a comparison with Neural Kernel Network. \nThe claims are supported by experiments. \n\nThe paper is very well-written and organized. The supplementary material provides settings to reproduce the results. Here are a some issues:\np. 3 l. 79: complet -- typo\nIt could also benefit to have larger pictures or larger fonts in the pictures, because they have to be ridiculously zoomed-in to be seen.\nAlso, please provide a clearer explanation of the choice of the uniform prior for frequencies larger than the fundamental.\n\nThe results are important because they provide a relatively low-cost (compared to HMC) technique to marginalize SM kernel hyperparameters which improves the performance of the GP. The SM kernel is indeed a very expressive family. \n\nOverall, very good paper.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors are aware of the limitations of their work and possible societal impact.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "72jbfYte5tG", "reviewer": "Reviewer_q4es", "summary": "The paper addresses, the problem of simultaneously evaluating the marginal likelihood for Bayesian model selection and the inference of the posterior of hyperparameters in Gaussian process (GP) regression. To this end the authors utilize the nested sampling (NS). They show empirically, that when little data is present or the dimensionality of the problem is high NS outperforms the more common methods, such as marginal-likelihood (ML-II) maximization and Hamiltonian Monte-Carlo (HMC) sampling.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The problem the authors consider is an important one in the Bayesian community. Apart from some minors (see below), I found the paper also well written. The experiments seem rigorously designed.\n\nWhile I find the paper interesting and the experiments insightful, my main concern is the novelty of the proposed method. While it is true, that as the authors state NS isn’t well known in the ML community, it has been discussed to evaluate the marginal likelihood integral (see e.g. the thesis of Ian Murray [https://homepages.inf.ed.ac.uk/imurray2/pub/07thesis/murray_thesis_2007.pdf], which btw contains some nice illustrations to explain NS intuitively) and was also used in several works (of whose to be fair the manuscript is citing several). The main contribution of the current work comes from more rigorously comparing NS to the dominant solutions (ML-II, HMC), to solve the integral of interest, and using NS for  complicated priors (priors with spectral mixture kernels) in the GP regression setting. \n\nOverall, I think the paper is and edge case, where on the pro-side stands a rigour empirical study to showcase the advantages of NS as a widely used model in the community, and on the contra limited novelty. In my opinion, the former slightly dominates, hence my rating.\n\n## Update\n\nAfter reading the authors response and the other reviews, I think, I was potentially underestimating the novelty of the work slightly. Hence, I will increase my score to 7.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Minors:\n+ I would appreciate if runtimes could be reported for the experiments.\n+ Notation in Algorithm1:\n++ N -> n_L in the exponent\n++ I think, the index i is a bit overloaded in meaning (running index in the initial set of live points, then a counter of samplings steps). Maybe change i to j in the loop and make clear, that it is a counter.\n+ I’d prefer to write the convergence criterion additionally as a formula, because I needed to read the sentence in l.140 several times, and was still struggling.\n+ L.79 completE\n+ I found the paragraph from l.137 onwards not helpful. Especially, when it came to the rslice option. It would be nice, if the authors could add some more details, at least in the supplement.\n\n", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5h"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes the use of nested sampling (NS) for inference in Gaussian process (GP) models with Gaussian likelihoods and shows the benefits of this fully Bayesian approach over competing methods such as HMC and type-II marginal likelihood hyper-parameter estimation. Although NS has been studied previously and the evaluation of the proposed method is mainly focused on the spectral mixture kernel and low-dimensional settings, all the reviewers recommend acceptance and I agree, as the paper further strengthens the machinery to carry out inference in GP models. I note and appreciate the authors’ commitment to release the corresponding code. 
Decision: Accept (Poster)