Sharp Calibrated Gaussian Processes
Authors: Alexandre Capone
Sandra Hirche
Geoff Pleiss
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian Processes
Frequentist Statistics
Kernel Methods
Model Selection and Structure Learning
Regression
TL;DR: We show how to obtain sharp calibrated Gaussian process models by exploiting kernel properties.

Abstract: While Gaussian processes are a mainstay for various engineering and scientific applications, the uncertainty estimates don't satisfy frequentist guarantees and can be miscalibrated in practice. State-of-the-art approaches for designing calibrated models rely on inflating the Gaussian process posterior variance, which yields confidence intervals that are potentially too coarse. To remedy this, we present a calibration approach that generates predictive quantiles using a computation inspired by the vanilla Gaussian process posterior variance but using a different set of hyperparameters chosen to satisfy an empirical calibration constraint. This results in a calibration approach that is considerably more flexible than existing approaches, which we optimize to yield tight predictive quantiles. Our approach is shown to yield a calibrated model under reasonable assumptions. Furthermore, it outperforms existing approaches in sharpness when employed for calibrated regression.

Sharp Calibrated Gaussian Processes
Alexandre Capone Technical University of Munich ale alexandre capone@tum... de
Sandra Hirche
Technical University of Munich
hirc@cit tum de
Geoff Pleiss
University of British Columbia
Vector Institute
geoff. wwwffplleissssttttuuee ubc.
Abstract
While Gaussian processes are mainstay for various engineering and scientific applications, the uncertainty estimates don't satisfy frequentist guarantees and can miscalibrated practice. State-of-the-art approaches for designing calibrated models rely inflating the Gaussian process posterior variance, which yields confidence intervals that are potentially too coarse. To remedy this, we present calibration approach that generates predictive quantiles using computation inspired by the vanilla Gaussian process postcrior variance but using a different set of hyperparameters chosen satisfy empirical calibration constraint. This results in a calibration approach that is considerably more flexible than existing approaches, which optimize yield tight predictive quantiles. Our approach is shown to yield a calibrated model under reasonable assumptions. Furthermore, it outperforms existing approaches in sharpness when employed for calibrated regression.
Introduction
Gaussian process (GP) regression offers an ambitious proposition: by conditioning model on measurement data, we are provided with Gaussian probability distribution for the unseen data. Assuming that the posterior probability distribution holds, we can then directly calibrate our model using the inverse eтror function. Though the distribution of unseen data seldom follows the Gaussian prior distribution, and the GP generally does not adapt adequately to the observed distributions after being conditioned the data, GPs have become one of the most powerful and established regression techniques. Besides having found widespread use in machine learning (Deisenroth et al., 2015; Srinivas al., 2012), their good generalization properties have motivated applications in the fields of control (Kocijan, 2016), astrophysics (Roberts et al., 2013) and chemistry (Deringer et al., 2021), 12 to name few. Furthermore, the Bayesian paradigm offers powerful tool to analyze the theoretical properties of different regression techniques (Srinivas et al., 2012; Capone et al., 2022).
In this paper, present novel approach to obtaining sharp calibrated Gaussian processes, i.c., Gaussian processes that provide concentrated predictive distributions that accurately match the observed data. Instead of computing confidence intervals by inflating the Gaussian process posterior variance, our approach discards it and computes a new quantity inspired by the computation the posterior variance, where all hyperparameters sharp calibration. other words, we train two separate Gaussian processes: one for the predictive mean and one for obtaining predictive quantiles, which is exclusively used for calibration purposes. By doing so, we reach considerably more flexibility than existing calibration approaches, which enables us additionally optimize the sharpness of the calibration. Our approach outperforms several
chosen in way that results in both accurate and
37th Conference on Neural Information Processing Systems (NeurIPS 2023).

---

state-of-the-an calibration approaches in terms of sharpness while still yielding similar calibration performance. Furthermore, is competitive compared to neural network-based method sharpness without sacrificing calibration performance.
Notation. We use R- to denote the non-negative real numbers. Boldface lowercase/uppercase characters denote vectors/matrices For two vectors and a' in R we employ the notation a s a' to denote componentwise inequality, ie.. a, S de i 1,.... ,, d. For square matrix K, we use | to denote its determinant, and [K] to denote the entry corresponding to the i-th row and j-th column.
2 Related Work
Calibration of Classification Models. There has been extensive work on obtaining calibrated models in the domain of classification. While there are many methods that do not employ post- processing. we only focus here on methods that employ some form of post-processingg Most forms of post-processing- based calibration for classification fall into the category of conformal methods (Vovk et al., 2005), which, given an input, aim to produce sets of labels that contain the true label with a pre-specified probability. Arguably the two most common forms of calibration are isotonic regression (Niculescu-Mizii & Caruana, 2005) and Platt scaling (Platt et al., 1999). In Niculescu-Mizil & Caruana (2005), Platt scaling and isotonic regression are analyzed extensively for different types of predictive models. In Guo et al. (2017), modified form of Platt scaling for modern classification neural networks is proposed.
Calibration of Regression Models. Though initially developed for classification, conformal calibration has been extended to regression settings. In Lakshminarayanan et (2017), calibration approach was proposed for deep ensembles. Gal et al. (2017) propose dropout- based technique for calibrating deep neural networks. However, these approaches require changing the regressor, potentially deteriorating its predictive performance. It should be noted that Bayesian neural networks (MacKay, 1995), while being able to provide credible sets for the output, fully trust the posterior, resulting in a naively calibrated model that seldom reflects the data's distribution. As remedy for this, Kuleshov et al. (2018) present a recalibration approach that scales a model's predictive quantiles to satisfy the observed data's distribution. In Vovk et al. (2020), similar approach is presented, where interpolation between scaling factors randomized, and a thcoretical analysis is provided. An extension of both Kuleshov et al. (2018) and Vovk et al. (2020) and other recalibration is proposed in Marx et al. (2022), along with corresponding theoretical guarantees. While these methods have been shown to yield well-calibrated models, the resulting predictive quantiles are potentially much too crude, resulting in predictions that perform poorly in terms of sharpness, i.e., the corresponding confidence intervals will overestimate the model error by very large margin. To remedy this, Song et al. (2019) and Kuleshov & Deshpande (2022) propose optimizing the parameters of a recalibration model by obtaining calibration on distribution level. However, while Song et al. (2019) relies on complex approximations and provides no theoretical guarantees, Kuleshov & Deshpande (2022) does not allow optimize for sharpness directly, and calibration only guarantced asymptotically as the number of data grows.
3 Problem Statement
Consider compact input space X' Rdd, and output space y R. and unknown data distribution II on X' x y. Consider a model conditioned on training data Drr ~ II that, for every x E X' and con- fidence level 8 E [0, returns base prediction J'D, (x:) and an additive cut-point term (5 (6, æz), where DD (6, x) 2 0 and 33 E IR is potentially negative, such that μ.D (x) B5 (6, x) corresponds predictive d-quantile. The model is then said to be calibrated if
Fæ (y (x) < (,0
(1)
holds for every 8 E [0, 1). Furthermore, the calibrated model is also said to be sharp if the correspond- ing predictive distributions are concentrated (Gneiting et al., 2007), i.e., if the centered confidence intervals induced by the predictive quantiles
1B30D, (ō, x) -1-30p. 8,z)
(2)
2

---

are small as possible for every 5 E [0, Our goal is to find a sharply calibrated model based on GP regression.
4 Gaussian Process Regression
In this section, we briefly review GP regression, with particular focus on the choice and influence of hyperparameters.
A GP is formally defined as a collection of random variables, any finite subset of which is jointly Gaussian (Rasmussen Williams, 2006). fully specified by a prior mean function, which we set to zero without loss of generality, and a ppperparrameeerddeeeddeee kernel k θ x X' X + R, where O denotes the hyperparameter space. The core concept behind GP regression lies in assuming that any finite number of measurements of an unknown function aN E X' are jointly Gaussian with mean zero and covariance
covariance function, called
f: + R at arbitrary inputs 21,....
(0), where
K ()),, = k(0, x; x;,
consists of kernel evaluations the test inputs.
Given set of noisy observations (a;, 1/ ) where y; == fxx; + and E, ~ N of) iid Gaussian measurement noise, can condition the GP on them to obtain the posterior distribution. The posterior distribution for a new observation y* at an arbitrary input x" is again Gaussian distributed, with mean and variance
PD. (D,(0,)) 00) (K()) + o( y.
(3a)
k()) (K) + TēT) k())) + o.. with k(0) == ((00,"".11)....... k(0, x", (0,))",
2 (k (
(3b)
WN). (y1,.... , and I denoting the identity
matrix. In this paper, we restrict ourselves to kernels k(:,:s) that yield a posterior variance that is monotonically increasing with respect to the hyperparameters 0, as specified in the following assumption.
Assumption 4.1. The posterior variance o,, x²) ) is continuous function of 0, Furthermore, for all hyperparameters 0, 0' E with 0 < 0', holds that o (0,x*) ) S ob (0',
Assumption 4.1 holds trivially for the signal variance of kernel. Moreover, it holds for any hyperparameters that lead to a monotonous increase the Fourier transform, which the case, e.g., for the inverse lengtschale of stationary kernels up to multiplicative factor corresponding to the ratio of the lengthscales (Capone et al., 2022). Furthermore, several results that employ the so-called fill-distance indicate that Assumption 4. holds for the inverse lengthscale of a broad class of stationary kernels as the lengthscale becomes very large or very small (Wendland, 2004). In our experiments, observed that Assumption 4.1 was never violated for the inverse lengthscale of the squared-exponential kernel. Assumption 4.1 will be leveraged to define a cumulative density function by also changing the hyperparameters corresponding to the posterior covariance, as opposed to simply scaling it.
Arguably one of the most challenging aspects of GP regression lies in the choice of hyperparameters 0, as they ultimately determine various characteristics the posterior, c.g., smoothness and amplitude. In practice, the most common way of choosing 0 is by maximizing the log marginal likelihood
N log(2x) " (K(0) + o²T) y.
log p 0) -2 log [((0) of
(4)
In terms of posterior mean quality, i.e., predictive performance of μDD, (0, x: choosing the hyper- parameters in this manner is often the most promising option, since it seeks a trade-off between model complexity and data fit, and has repeatedly been shown yield satisfactory mean square error when applied to test data (Rasmussen & Williams, 2006). However, even when employing log-likelihood maximization, the posterior GP distribution seldom well calibrated, and, in practice, the data often has a significantly different distribution. As a result, quantiles obtained with a purely Bayesian approach are either too confident or too conservative in practice (Capone al., 2022; Fong & Holmes, 2020). Furthermore, the restrictions imposed by the GP prior often produce posterior variance o² 22 x²) that grossly conservative. See Figure 1 for illustration.
3

---

Training data
Calibration (holdout) data
Test data
Target of
Input N
Input
(a) Bayesian approach.
(b) approach.
Figure 1: Confidence interval of 95% (shaded regions) obtained with purely Bayesian approach (a), where the inverse error function employed to compute 95, and our approach (b). Solid lines represent the predictive mean. The confidence interval obtained with the Bayesian approach is not only grossly overconfident (contains less than 80% of the total data) but also partially extremely loose, exhibiting unnecessarily large confidence intervals far away from the data. By contrast, approach is both accurate and tight.
Proposed Approach
In this section, we present our approach to obtaining calibrated GPs by discarding the posterior variance and obtaining alternative predictive quantiles using quantity inspired by the posterior variance and new hyperparmetors. As we will demonstrate, our approach has numerous advantages. First, by varying the hyperparameterss used obtain predictive quantiles, the resulting quantiles are sharper than what can be obtained simply by multiplying the posterior variance with an appropriate constant. Secondly, by exploiting the monotonicity of the hyperparameters,, to obtain intervals for multiple confidence levels ō in a very efficient manner. Finally, our method is backed by tight theoretical guarantees, obtained by exploiting its connection conformal prediction.
approach can be used
5.1 Sharp Calibrated GP for Single Confidence Level 5
In the following, describe how obtain sharply calibrated GP for fixed desired calibration level d. Instead of scaling the GP posterior variance to meet the desired calibration level we propose computing predictive quantiles by training new quantity similar the posterior variance but with new hyperparameters. In other words, we discard the posterior variance of the first GP and replace by a quantity that corresponds to the posterior variance different GP, which we train separately. This way, we allow for more degrees of freedom during calibration. We then leverage this additional freedom minimize the distance of the quantiles the predictive mean, which yields a sharply calibrated model.
We assume to have data set DD which split into training data Drr and calibration (holdout) data Deat, with D = Drr J Deat- The training data Drr will be used to compute the posterior, whereas Deat will be used to calibrate the model. Note that while splitting the data might be reasonable when using other types of regression models, e.g., Bayesian or ensemble neural networks, splitting the data in the case of GPs beneficial for providing accurate quantiles for data out of distribution. This is because, for many commonly used kernels, the GP posterior distribution is considerably more concentrated for test points close to the training data Dr than for those far away from Dr. Since we wish to obtain model that is calibrated for data both close and far away from the training data, we take this into account during training by splitting the data. We also assume have a predictive mean uDD R corresponding to GP posterior mean function, where the regressor hyperparameters R were obtained, e.g., via log-likelihood maximization, as discussed in Section 4. However, note that any other way of choosing the posterior mean hyperparameters 0 is permitted.
Given HDD. (0 and D UDeal we then follow the convention of other recalibration approaches (Kuleshov et al., 2018; Marx et al., 2022) and aim to obtain, for an arbitrary < 5 < 1, scalar 33 and vector of hyperparameters ls, such that the corresponding predictive quantile (D,(0 (0 +

---

DA (0s.)) contains 5 times the total amount of data points. For this reason, we henceforth refer to Og as calibration hyperparameters.
In order to obtain a sharply calibrated GP model, ideally we would like to choose 33 and 0g such that they minimize the expected length of the centered intervals (2) subject to calibration. However, this optimization problem is hard to solve. Hence, we instead attempt to improve model sharpness by concentrating the predictive distribution around the predictive mean jrD, (D, .). i.c. by minimizing the deviation of the quantiles from j.D, (ORR.).
(OR, :). This corresponds solving the optimization problem
Nai Σ BB (0s; alaa))
min P&ER
O.EO
(5)
Nai
I20 (Ayesl - Br (0s, reaa)) 6
S.t.
Neat +1
where { Veut } E Deal samples from the calibration data set, Avia == Vew J.D. (0 x corresponds to the difference between predicted and measured output, and Neal |Deal to the number of data points used for calibration.
Remark 5.1. Note that we employ Nea + in the denominator in (5) instead Neal- Though this choice makes little difference practice, require it for theoretical guarantees.
The equality constraint in (5) is generally infeasible, e.g., if Neat 2 and 8 0.5, and is discontinu- ous, making it hard to solve. As turns out, this can be casily remedied without any detriment sharpness calibration, and the problem can be rendered considerably casier solve by substituting the equality constraint in (5) with
3B - Qu1 Ep, Ayau),
(6)
where Qum (8, Ep, Aya is monotonically increasing piecewise linear function that maps 8
j/ (Nael + 1) to the j-th smallest entry of Ep. Ayaal, and the entries of the vector
Nat
Yeal
A Yeal Drp
ODD (0s, x
TD, (05, DD, Feal
Not
correspond to the z-scores of the data under the calibration standard deviation dD, (0s, The original problem (5) then becomes
min Σ
[Tin (6, Ep. D, (0s, z(a))
which is considerably easier to solve due to the lack of constraints, and enables us to use gradient- based approaches, since Qin (6, Σn! yeal) is differentiable with respect to oD, (0,, aaa))
Remark 5.2. The choice of interpolant (6) is due to its simplicity. However, other forms of monotone interpolation are also possible. For high data sizes, the choice of interpolant becomes of little relevance, since we only perform small interpolation steps.
5.2 Calibrated GP for Arbitrary Confidence Level 5
While (7) is useful for obtaining a sharply calibrated model for single confidence level 6, solving (7) multiple times whenever we want to obtain sharply calibrated models for different confidence levels can be time-consuming. Furthermore, interpolating between any two arbitrary solutions of (7) won't 't necessarily yield result close to the desired calibration. Fortunately, we can leverage Assumption
show that interpolating between two solutions of (7) will yield a result close to the
For a E Na permutation 41, ".....,AN, ANa E [1p..y Nea) where S aift2 for all j, and e 1/ + 1), (1 1)/(Nas 1)),
(6, a) ai; (6(N (a441 au)
5

---

Algorithm Training Calibration Hyperparamekers for Arbitrary Confidence Level
Input: kernel k(-).)), predictor HDD, (0", :), calibration data Peal confidence levels ô1,...ON
for i = 1 to M do Compute Ba,, , 0,,,... NōNe
5 Nea
by solving (6) and (7) subject to
0g. S 08, 6 < 6j and 3s. 2 0,
0.. 05, if 6 6g and Bs, S
end for
Fit a continuous, monotonically increasing interpolation model (6) and a continuous model (())
using the training data {6,, 9b. Output: (()), 0(8)
desired calibration, provided that we interpolate between two strictly increasing or decreasing sets of hyperparameters. Formally, this is achieved by solving (7) Neat + times obtain fo. Ba,..... Bäy R and Oo .s....., O& O subject to two additional constraints. First, the calibration scaling parameters 35 must be monotonically increasing with i.c.
Ss, S Ss, if 6, < 6/,
and the calibration hyperparameters O, must be decreasing with ò Bg is negative, and increasing if 38 is positive, i.c.,
0, 0ö,, if 6, 6, and 38_ 2 0,
0, 0s,, of < 8, and 3s, S 0.
In other words, the entries of 03 are strictly decreasing with As up until the sign of Bs switches, after which they are increasing. The reason why impose these restrictions that we then confidently interpolate between any values of B3 and 0,, since Assumption 4.1 implies that the quantile stipulated by D (0, -) is monotonically increasing with ó. We then train simple piecewise linear interpolation models [0, 1] > and 0 [o, 1] - O, such that ((,)) = 38, and (())) Os., with additional constraint that (()) reaches a minimum whenever 3(5) = O, which can be potentially achieved by adding an artificial vector of training hyperparameters Os for computing 0(5). Note, however, that any other form of monotone interpolation is also acceptable for obtaining ((6) and (()). The procedure is summarized in Algorithm 1.
Remark 5.3. While lengthscale constraints of the form 0,, < 08 can be casily enforced when solving (7) by substituting with 0s, + AO and minimizing over the logarithm of ÅO, the constraint 35, < Bs,. is not enforced in (7). However, in practice we were often able to find local minima of (7) that satisfy this requirement.
We can then easily show that our approach achieves an arbitrary calibration level as the amount of data grows, provided that we choose the confidence levels ó accordingly.
(D,,( Theorem 5.4. Let y absolutely continuous conditioned on æ, μu, (QR,-) be a posterior GP mean, and let JD (.,.).))) be GP posterior variance conditioned Drr Then, for any calibration data set Deat = { Vial ), choose
2
Neal
61
82
NN...
N cal
Neat
N +
and let B() and (-)) be interpolation models obtained with Algorithm I and confidence levels bi.....c &N.. Then
Pay-ll (y PDd (0(6), (1000
0(0)
5 Neat 1 Neal +
1
Proof.
The proof can be found in the supplementary material.
Note that Theorem 5.4 also implies that a single set of calibration parameters |3 and Og obtained by solving (7), since we can substitute ((8) = .3 and ((8) 0g into (8).

---

σ (a) Our approach.
T
I
(b) Kuleshov et al. (2018).
(c) Vovk et al. (2020).
(d) Mcan constant (Marx al., Neural network-based recalibra- 2022).
Base model (vanilla GP).
tor (Kuleshov & Deshpande, 2022).
Figure 2: Centered 99 % confidence intervals (shaded regions) obtained with our method, vanilla GPs recalibrated using the approaches of Kuleshov et al. (2018) and Vovk et al. (2020), naive vanilla GP, the point-predictor (variance-free) approach proposed in Marx al. (2022), and naive fully Bayesian GP. Solid lines represent the predictive mean, crosses represent data used to train the base model, and squares represent calibration (holdout) data. Our approach yields model that is both well-calibrated and sharp. This is because of the added flexibility that comes from being able to also change the lengthscale design calibrated model.
Discussion
Computational Complexity. Much like hyperparameter optimization standard GPs (Rasmussen & Williams, 2006), the major driver behind the computational complexity in our approach stems from need invert the covariance matrix, an operation that scales cubically with the amount data. In order to alleviate the computational cost of our approach, we resort to different tools that improve scalability (Liu et al., 2020). One approach consists of employing only a subset of the training data D to choose the calibration hyperparameters 0g, and then the full data set choose B3. While this potentially leads to a loss in sharpness compared to when using the full data set, it still guarantees calibrated model. Our technique is also readily applicable sparse GPs (Snelson & Ghahramani, 2005; Titsias, 2009), as Assumption 4.1 typically still holds. This option is also explored in numerical experiments, in Section 7. Moreover, calibration is often required, as opposed to several different ones, e.g., in stochastic model predictive control, where chance constraints corresponding a fixed risk have be satisfied (Mesbah, 2016). In such scttings, we polentially only have to train single vector of calibration hyperparameters 03. which reduces computational cost.
many settings specific level of
The choice of initial hyperparameters can affect the optimization results considerably, and choosing good hyperparameter initialization can be challenging, is true when choosing the hyperparameters for the predictive mean. While this can be partially addressed by employing random restarts, can also reuse trained models for similar calibration levels, since is reasonable to expect that only small changes to the calibration hyperparameters are required also simplify the
Initialization and Solution.
to achieve slight increase or decrease in confidence level. Furthermore, problem by considering only scaled version the regression hyperparameters 8 to compute 0g, which would reduce the optimization problem to line search.
7 Experiments
In this section, we apply and analyze our approach using a toy data set and different regression benchmark data from the UCI repository. In the supplementary material, we also compare our

---

approach to that of Capone et al. (2022) when used to obtain uniform error bounds and apply our method to two different Bayesian optimization problems.
The goal is for our approach to obtain a sharp calibrated regression model for each data set in the calibrated regression experiments. We test our approach on various data sets and compare it to the state-of-the-ars recalibration approaches by Kulcshov et al. (2018) and Vovk et al. (2020), the point predictor (posterior variance-free) approach proposed in Marx et al. (2022), as well as the check- score-based approach of Kuleshov & Deshpande (2022). The technique proposed by Kuleshov et al. (2018) essentially multiplies the vanilla posterior standard deviation dDD (@ with the recalibrated z-score, such that the confidence level observed on the calibration data matches that of the desired confidence level. Vovk et al. (2020) employ a similar approach, except that random interpolation is employed to compute new scaling values. The point predictor-based method proposed in Marx et al. (2022) discards the posterior standard deviation dDD DD, .) and computes a constant scalar that is added to the predictive mean and used to compute quantiles everywhere within the input space. The method of Kuleshov & Deshpande (2022) trains a neural network using a quantile loss, which takes base quantiles as inputs and returns new, recalibrated quantiles. As recalibrator for the method of Kuleshov & Deshpande (2022), we employ the same neural network architecture suggested in their paper, trained over 200 epochs, with additional pretraining over 2000 epochs using single dataset for the UCI experiments. In all experiments except kin8nm and Facebook comment volume 2, we employ standard GPs with automatic relevance determination squared-exponential (ARD-SE) kernels and zero prior mean as base models, trained using log-likelihood maximization. For the kin8mm and Facebook comment volume 2 datasets, we employ sparse GPs (Titsias, 2009) with zero prior mean, ARD-SE kernels, and 300 inducing points.
7.1 Toy Data Set
The first regression data set corresponds to a one-dimensional synthetic data set, where the results can be easily displayed visually. The main purpose of this section is to give an intuition as to how our approach computes confidence intervals compared to other techniques. We investigate the performance of our approach and compare it to other methods when employed to compute centered 99% confidence intervals. We observe that the confidence intervals obtained with our approach peak less strongly far away from the data while being tight near the data compared to all other approaches except the one of Kuleshov & Deshpande (2022). This is because we allow the lengthscale to change to obtain calibrated model. In contrast, all other methods except that of Kulcshov & Deshpande (2022) scale the standard GP posterior variance without changing hyperparameters. The method of Kuleshov& & Deshpande (2022), which uses neural network as a recalibrator, offers additional flexibility, resulting in sharper confidence intervals. However, calibration is not explicitly enforced during training, resulting in poor calibration. The results are depicted in Figure 2.
7.2 Benchmark Data Sets
We now experiment with seven different regression data sets from the UCI repository, two containing over eight thousand data points and requiring sparse GP approximations. The training/clliraaiiontees split is 0.6, 0.2, and 0.2 for all data sets except the Facebook comment volume 2 data set, which contains over 80000 data points, and where the split 0.08, 0.02, and 0.9. For the approach of Kuleshov & Deshpande (2022), we follow the steps in their paper and limit the calibration data size to 500.
We assess performance by employing diagnostic tools commonly used to assess calibration and sharpness (Kuleshov al., 2018; Marx al., 2022; Gneiting et al., 2007)). The score used quantify calibration is the calibration error (Kuleshov el al., 2018), given by
m
cal (00,
Σ (p) - -11
(8)
where p5 corresponds to the j-th desired confidence level, chosen, e.g., evenly spaced between 0 and 1, and βú is the observed confidence level, i.e.,
Ayi < B, oD, (0,, æi), t= I,..., T,.
(9)
T

---

Table 1: Expected calibration and sharpness of different methods 100 repetitions per experiment. We report the expected calibration error (ECE), the average predictive standard deviation (STD), negative log-likelihood (NLL) and 95% confidence interval width (95% CT) obtained with approach, vanilla GPs recalibrated using the methods of Kuleshov et al. (2018) (RK) and Vovk et al. (2020) (RV), the variance- free approach proposed in Marx et al. (2022) (RM), and the neural network- based recalibrator of Kuleshov & Deshpande (2022) (NN). We additionally report performance for the base model (B), which corresponds all metrics. In all experiments except the Facebook2 dataset, our method is sharpest compared to all other methods except that of Kuleshov Deshpande (2022). However, Kuleshov & Deshpande (2022) performs more poorly in terms of expected calibration error.
vanilla GP without the holdout data. Lower is better for
DATA SET
METRIC
OURS
RK
RV
RM
NN
B
ECE
0.003
0.0029
0.0029
0.0029
0.0056
0.041
STD
0.16
0.31
0.33
0.22
1.9
BOSTON
NLL
0.21
0.39
0.4
0.42
-0.24
1.6
95% CI
0.76
1.4
1.4
1.4
0.73
7.4
ECE
0.0044
0.0043
0.0044
0.0043
O. 0.0081
0.039
STD
0.16
0.5
0.47
0.5
0.14
2.8
YACHT
NLL
0.26
0.68
0.69
0.68
-2
2
95%
0.76
2.3
2.3
11
ECE 0.0036
0.0035
0.0035
0.0035
0.0053
0.044
STD
0.13
0.38
0.38
0.38
0.29
2.8
MPG
NLL
0.032
0.63
0.64
0.63
0.02
2
95% CI
0.6
1.7
1.8
1.7
0.96
11
ECE
0.00047
0.00047
0.00047
0.00047
0.0067
0.0058
STD
0.54
1
0.88
0.72
1.4
WINE
NLL
1.3
1.3
-0.36
1.4
95% CI
2.1
3.8
3.8
3.9
2.8
5.4
ECE
0.00071
0.00064
0.00064
0.00064
0.00076
0.032
STD
0.25
0.64
0.64
0.64
0.57
2.8
CONCRETE
NLL 0.72
1.1
1.1
1.1
0.85
2
95% CI
0.93
2.5
2.5
2.5
2.1
11
ECE
0.00016
0.00016
0.00016
0.00016
0.00053
0.028
STD
0.074
0.12
0.12
0.12
0.098
0.4
KINSNM
-0.65
-0.65
-0.63
-0.76
0.1
95% CI
0.26
0.47
0.47
0.48
0.44
1.6
ECE 0.00044
0.00043
0.00043
0.00045
0,0089
0.044
STD
0.068
0.18
0.18
0.18
0.18
1.2
FACEBOOK2
NLL
3.6
-1.3
-1.3
-1.2
-2.3
1.2
95%
0.6
1.7
1.7
3.4
4.6
Here the superscript denotes test inputs and outputs, T denotes the total number of test points, and Ay; == L'D, (0 xx;) yr- We employ m = 21 evenly spaced values between 0 and for Pz- To measure sharpness, we employ the average length of the 95% confidence interval, the average standard deviation of the predictive distribution, and the average negative log-likelihood of the predictions (Gneiting et al., 2007; Marx al., 2022). Note that since every model outputs a quantile for any desired calibration level, the corresponding negative log-likelihood and average standard deviation are well specified. These are computed by employing the cumulative distribution function, obtained by inverting the quantile function specified by each model.
We carried out each experiment 100 times and report the resulting average expected calibration error, standard deviation, negative log-likelihood, and length of the centered 95% confidence intervals in Table 1. Our approach performs best marginally worse than all other calibration approaches regarding expected calibration eиrоr. This is to be expected from Theorem 5.4. Furthermore, it oulperforms all approaches except that of Kuleshov & Deshpande (2022) in sharpness. However, the improved sharpness of the method of Kuleshov & Deshpande (2022) comes at the expense of calibration.

---

8 Conclusion
We have presented calibration method for Gaussian process regression that leverages the monotonic- ity properties of the kernel hyperparameterss to obtain sharp calibrated models. We show that, under reasonable assumptions, our method yields an accurately calibrated model as the size of data used for calibration increases. When applicd to different regression benchmark data sets, approach was shown to be competitive in sharpness compared to state-of state-of-the of-the-art recalibration methods without sacrificing calibration performance. It is worth stressing that, though the tools presented here emerge naturally from a Gaussian process setting, we do not require our predictor to be a Gaussian process to obtain theoretical guarantees. In future work, aim to leverage similar monotonicity characteristics to get sharply calibrated models using tools different from Gaussian processes. Furthermore, we aim experiment with inducing variables hyperparameters when optimizing the models for sharpness.
Acknowledgements
This work supported in part by the European Research Council Consolidator Grant Safe data- driven control for human-centric systems (CO-MAN) under grant agreement number 864686.
References
Berkenkamp. F., Schoellig. A. P.. and Krause, A. No-regret baycsian optimization with unknown hyperparameters. Journal of Machine Learning Research, 20(50):1-24, 2019. URL http: //jmlr. wwtpg/apprr/vve satparss/v//18/ 1110018-2111 213 html.
Capone, A., Lederer, A., and Hirche, S. Gaussian process uniform error bounds with unknown hyper- parameters for safcty-critical applications. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 2609-2624. PMLR, 2022. URL https: httpss////ccceco.g...e.. mlr. httpss//vvc/cccoo22222 html.
Chowdhury, S. R. and Gopalan, A. On kernelized multi-armed bandits. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 844-853, International Convention Centre, Sydncy, Australia, 2017. PMLR. URL attp://proceedings mlr, press/v70/ https//////hwwwurr/////// chowdhury176 html.
Deisenroth, M. P., Fox, D., and Rasmussen, C. E. Gaussian processes for data-efficient learning in robotics and control. JEEE Transactions on Pattern Analysis and Machine Intelligence, 37(2): 408-423, 2015.
Deringer, L., Bartók, A. P., Bernstein, N., Wilkins, D. M., Ceriotti, M., and Csányi, G. Gaussian pro- regression for materials and molecules. Chemical Reviews, 121(3)))10033 10141,, 2021. doi: wwwwwee.eee...ceemmm 1000022. URL https: http::://doioor//1.... wttp//...////////...... 1021/acs..... chemrev sesemrevv.c0022. 1c00022. PMID: 34398616.
Fong, E. and Holmes, C. C. On the marginal likelihood and cross-validation. Biometrika, 107(2):489- 496, 2020. doi: 10.1093/biomet/asz077 URL attpp:://oi.. 1ttp:///...9////ommettte.e /tsp///iimmmaa////..... asz077.
Gal, Y., Hron, J., and Kendall, A. Concrete dropout. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL https: h/tpttpep////ooeeiinggsss www.eeeiinnnnniiiinn.oeeee newrips..s.... cc/ paper/2011 http://iieeeeaaaaa.4...884666663 Pape capper.ddt pdf..
Gneiting, T., Balabdaoui, F., and Raftery, A. E. Probabilistic forecasts, calibration and sharpness. Journal of the Royal Statistical Society: Series (Statistical Methodology), 69(2):243-268, 2007.
Guo, C., Pleiss, G., Sun, Y., and Weinberger, K. Q. On calibration of modern neural networks. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pp. 1321-1330. PMLR, 2017. URL https: https://proceed mlr. www.ss/vs00.. /tto////o.a.....// html.
Kocijan, J. Modelling and control of dynamic systems using Gaussian process models. Springer, 2016.
10

---

Marx, C., Zhao, S., Neiswanger, W., and Ermon, S. Modular conformal calibration. In Proceedings of the 39th International Conference on Machine Learning. volume 162 of Proceedings of Machine mlr press/ Kuleshov, V. and Deshpande, S. Calibrated and sharp uncertainties in deep learning via density estimation. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 11683-11693. PMLR, 2022. URL mlr. wttps///1///ulllssopooeee
html.
Kuleshov, V., Fenner, N., and Ermon, Accurate uncertainties for deep learning using calibrated regression. In Proceedings of the 35th International Conference Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 2796-2804. PMLR, 2018. URL https: //proceedings.co mlr prress/v80/kullsso/// html.
Lakshminarayanan. B.. Pritzel, A., and Blundell, C. Simple and scalable predictive uncertainty estimation using deep ensembles. In Advances in Neural Information Processing Systems, vol- ume 30. Curran Associates, Inc., 2017. URL https: /ttppp/:/s...oeeing..... meg.ouuednnggn..ur......... aatereeffaaaabbee33-32a .ddf.
wwwrips cctccc/cppp//ac..c paper/
Liu, H., Ong, Y.-S., Shen, X., and Cai, J. When Gaussian process meets big data: A review of scalable GPs. IEEE Transactions on Neural Networks and Learning Systems, 31(11):44054422) 2020. doi: 10.1109/TNN.S220111 2019..99711111171111 109.
MacKay, D. J. Bayesian neural networks and density networks. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 354(1):73-80, 1995.
Learning Research, pp. 15180 15195. PMLR, 2022. URL v162/mare m2t6mmax222...... .html.
Mesbah, A. Stochastic model predictive control: An overview and perspectives for future research. IEEE Control Systems Magazine, 36(6):30-44, 2016.
MaateescuMMzzil A. and Caruana, R. Predicting good probabilities with supervised learning. In Proceedings of the 22nd International Conference on Machine Learning, ICML '05, pp. 625 632, New York. NY. USA, 2005. Association for Computing Machinery. ISBN 159593180s doi:
10.1145/102335.111 102351.1102430 URL https: //doi. http:///oooooooo/.///o//// 110...1111111115/1111 1102430.
Platt, J. et al. Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Advances in large margin classifiers, 10(3):61-74, 1999.
Rasmussen, C. E. and Williams, C, K. Gaussian processes for machine learning. 2006. The MIT Press, Cambridge, MA, USA, 2006.
Roberts, S., Osborne, M., Ebden, M., Reece, S., Gibson, N., and Aigrain, S. Gaussian processes for time-series modelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, (1184))))201115555 2013.
Snelson, E. and Ghahramani, Z. Sparse gaussian processes using pseudo-inputs. In Advances in Neu- ral Information Processing Systems, volume 18. MIT Press, 2005. URL https::/proceediii neur ips. cc/paper///eeecpppcccecceee..eat /4t.0////../44777771222222226660
Paper. Sapper.pd.. .pdf.
Song, H., Diethe, T., Kull, M., and Flach, P. Distribution calibration for regression. In Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 5897-5906. PMLR, 2019. URL htRp::tts:://rooeeeeee wmw.ediinnnn.........mmmmm.mmr.m www.mmrreess.............mmmmmm. v97/song" 9a. html.
Srinivas, N., Krause, A., Kakade, M., and Seeger, M. W. Inforeition-theerrtie regret bounds for Gaussian process optimization in the bandit setting. IEEE Transactions on Information Theory, 58 (5):3250-322 3265, 2012. doi: 10. 2019/IT.2011.218220
Titsias, M. Variational learning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics, pp. 567 574, 2009.
Vovk, V., Gammerman, A., and Shafer, G. Algorithmic learning in a random world, volume 29. Springer, 2005. doi: 10. 10.10//116661.
11

---

Vovk, V., Petej, I., Toccaceli, P., Gammerman, A., Ahlberg, E., and Carlsson, L. Conformal calibrators. In Proceedings of the Ninth Symposium on Conformal and Probubilistic Prediction and Applications, volume 128 of Proceedings of Machine Learning Research, pp. 84-99. PMLR, 2020. URL https://proceednng. mlr press/ Setess/v128/vo/// vovk20a. html.
Wendland, H. Scattered Data Approximation. Cambridge Monographs on Applicd and Computational Mathematics. Cambridge University Press, 2004. doi: 10. 10.10177CB0998051111759
12


Reviews:
[{"review_id": "eJBVkrd1si", "reviewer": "Reviewer_jTuD", "summary": "The paper tackles calibration of Gaussian processes in regressions. The authors argue that while maximizing the evidence is a good way to choose hyperparameters to obtain an accurate posterior mean, it generally does not produces an accurate posterior variance. For this purpose, they propose a different way to obtain hyperparameters for the posterior variance which minimizes the length of centered confidence intervals, while constraining the parameters to provide accurate empirical quantiles. By replacing the constraint with a more treatable piece-wise expressions, they transform the problem into an unconstrained one. The further a practical algorithm and a theoretical result to obtain calibrated confidence intervals simultaneously over different confidence thresholds. Experimentation on toy and UCI datasets show the benefits of the method.", "questions": "- Not immediately clear to me what is the effect of replacing equation (5) with equation (6). I understand that the problem becomes easier to solve. Is the replacement, however, detrimental in some way? It is also not straightforward why the expression is (6) is the right one to use. I would encourage the authors to expand on this section.\n\n- It is not super clear how easily the inequality in Assumption 4.1 holds. Perhaps making some examples, and discussing limitations of this assumption, would make this part clearer.\n", "limitations": "The work handles calibration for regression only.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "I think the paper is well written and clear. The solution proposed is reasonable and it tackles an important problem. The algorithms derived seem to be practical, and they are motivated by theoretical arguments that appear valid.", "weaknesses": "The authors bring up themselves that the core methodology they propose shares the same scalability issues as standard GPs, as it involves they inversion of a very large matrix. Hence, approximation methods are needed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your review and positive remarks. Based on your and other reviews, we have made several changes to improve the paper. Below you can find those that address your questions and comments specifically.\n\nReplacing (5) with (6) is not detrimental, as it still allows us to obtain a model that is sharp and calibrated. The choice of expression (6) is the simplest and most straightforward to compute. However, other forms of monotone interpolation are also possible without any loss of theoretical guarantees. In practice, the choice of interpolant eventually becomes of little relevance for high data sizes, since we only require small interpolation steps. We will include this discussion in the revised paper.\n\nAssumption 4.1 holds for any hyperparameter of a stationary kernel such that the Fourier transform of the kernel is monotonic with respect to that hyperparameter. This can be used to show, e.g., that the lengthscale of a stationary kernel satisfies Assumption 4.1 up to a scaling factor depending on the lengthscale (see Capone et al., (2022), Lemma 3.3). We will discuss this in more detail in the revision."}, {"review_id": "fXtsGkXmiI", "reviewer": "Reviewer_XLKC", "summary": "This paper addresses the issue that the posterior variance of Gaussian processes are often poorly calibrated, typically underestimating quantile estimation. They propose a new method to calibrated uncertainty bounds, by training a quantity related to the posterior variance with new hyper parameters. This method further leverages optimization of the distance between the quantiles and the predictive mean, enabling sharp calibration. The method is tested on synthetic toy examples and standard UCI benchmark datasets, and appears to perform well against existing methods.", "questions": "How does the proposed method compare with those proposed in Song et al, Kuleshov & Deshpande, as well as Capone, Lederer & Hirche?\n", "limitations": " There is little view for potential negative societal impact due to improving the calibration of uncertainty bounds for Gaussian processes.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well written and presented and addresses a relevant problem, offering a principled solution.", "weaknesses": "The definition of a Gaussian process should specify that any **finite** collection of points are jointly Gaussian.\n\nThe most concerning issue is that the paper appears to be dismissive of two recent competing methods (Song et al.; Kuleshov & Deshpande) in the introduction, and does not mention them again. Indeed, the current work has advantages over each of these methods in terms of its construction and theoretical guarantees, however comparisons to these methods in terms of performance, as well as some further discussion, should help demonstrate the benefits of the proposed method. For example, the current method still relies upon a reasonably large dataset size, so dismissing a competing method because it only has asymptotic guarantees is maybe premature.\n\nI would also question the presentation of the numerical results, in particular the use of both highlighting and boldface to denote best performing model in the tables. Since confidence intervals are provided, it may be advisable to bold the best performing ones, and using green to denote statistical significance that a single model is the best performing (i.e. multiple models can be bolded but none or one should be green).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your review. We have made changes with your comments in mind as described below. If you feel that the changes sufficiently address your comments, we would be very thankful if you would consider raising your score.\n\n**Additional comparisons.** As you suggested, we compared our approach to that of Song et al. (2019), Capone et al. (2022) and Kuleshov and Deshpande (2022). We have included the results in the supplementary PDF. Note that the method of Capone et al. (2022) only computes uniform error bounds, i.e., 100 percent credible intervals, so we only considered this setting in the comparison. \n\nWe noticed that training takes longer with the approaches of Kuleshov and Deshpande (2022) and Song et al. (2019) than our approach. This is because the optimization problem in our approach is fairly straightforward to solve, whereas Kuleshov and Deshpande (2022) and Song et al. (2019) require several optimization steps. Moreover, the approach of Song et al. (2019) is fairly involved and more difficult to implement than ours.\n\nWe implemented the method of Kuleshov and Deshpande (2022) with the same neural network architecture suggested in their paper, trained over 2000 epochs. Our results suggest that the method of Kuleshov and Deshpande (2022) favors sharpness over calibration. While it achieves sharper intervals, it yielded significantly worse expected calibration errors than ours for most settings. Note that this is partially to be expected since their approach aims to achieve calibration in distribution, which only corresponds to quantile calibration if it is exact. Furthermore, it does not explicitly enforce calibration during training. We also noticed that it does not return valid intervals in some pathological cases (e.g., due to insufficient training steps).\n\nThe method of Song et al. (2019) seemingly performs similarly to our approach in calibration, whereas it performed better in sharpness in some cases. However, we could only implement it on small data sets due to slow training. Furthermore, their method is fairly involved, and we relied heavily on the code kindly provided by the authors. Furthermore, although we employed the same test and training splits for the data as in our case, we noticed differences in the GP model used in their starting setup, potentially leading to significantly different starting models. We will ensure an identical setup for the revision.\n\nThe method of Capone et al. (2022) is purely Bayesian and thus heavily dependent on the prior. We employed their approach as is from the code available online. The resulting credible intervals are well-calibrated, consistent with the results presented in their paper. However, our approach is much better regarding sharpness. This is because Capone et al. (2022) require symmetric intervals ($\\text{mean} \\pm b \\times \\text{standard deviation}$), whereas our approach allows for asymmetric credible intervals. Furthermore, our method does not rely on the prior distribution.\n\n**Presentation of numerical results.** As you suggested, we also removed the bolding everywhere except when methods perform significantly better, i.e., whenever the statistics indicate it is indeed best. In the revised version, we present the mean plus minus standard error and clearly state this to avoid confusion.\n"}, {"review_id": "SuROXSCmOc", "reviewer": "Reviewer_XbYQ", "summary": "The authors propose a novel method for calibrating Gaussian process posterior variances using held-out data. In particular, they train a new, separate, GP using the held out data for the variance, using the GP trained on the original dataset for the mean. They do this in a way which approximately maximises the sharpness (roughly speaking, minimises the variance). The calibration properties of the method are supported by theory. This is supported by an empirical evaluation of the calibration error on synthetic data and of the sharpness on real-world data.\n\nPost-discussion: The authors have provided considerable improvement in their presentation of the results, and have added some extra results which further improve the clarity.", "questions": "1. It's not clear to me how to interpret figure 2. I guess that the crosses are training data, what about the squares? Are they the calibration data?\n2. I guess Table 1 is for the synthetic data plotted in figure 2, but I don't think this is stated anywhere.\n2. When you train the vanilla GP, do you add the calibration data to the training data?\n3. The sharpness of the predictives is claimed to be an advantage, but really it looks like you cover less of the data as a result in the toy example?\n4. I appreciate that the main purpose of Table 2 is to compare the sharpness amongst the well-calibrated methods,  but I think it would be highly worthwhile to see how a vanilla GP compares on these same metrics.\n5. With respect to line 103, it is mentioned that the assumption holds for isotropic kernels. I think that it also works out if you transform an isotropic kernel to have different lengthscales in each dimension, and I think you've implied this elsewhere, is that right?\n\nSome suggestions and typos:\n* For the hyperparameters, could you use $θ$?\n* It would be helpful to have a bit more explanation in section 3. For example, instead of saying 'are small for every $δ$' perhaps you could say that the goal is to minimise that quantity for every $δ$.\n* line 127: the z-score -> z-score\n* Summations in equations 5, 7: provide the starting value (i=1 I guess)", "limitations": "The main limitation of the method appears to be its applicability to different covariance functions, and the added computational cost, both of which are noted clearly.", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 3, "strengths": "1. (major) the idea is novel,  and appears to be theoretically well supported, although I did not review the proofs in the supplementary material.\n2. (major) The main claimed advantage of the proposed method is that is produces sharper predictives than comparable methods, which is well supported by the real-world experimental results, not withstanding the points raise in questions.\n3. (major) Code is provided, which should improve reproducibility. It is based on a widely used framework, which increases potential for impact.", "weaknesses": "1. (major) The paper is not sufficiently clear. The method is fairly well explained, but the evaluation is extremely hard to follow, details of which are in the questions section. One particular area for straightforward correction is with Table 2: the values with the lowest mean are highlighted but there are several rows where the estimates are overlapping, so this should be clarified by highlighting.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review. We have made several improvements to the paper, particularly in the experimental section, to address your concerns about clarity. Below you can find our answers to your questions and comments. Corresponding modifications can also be seen in the PDF submitted with this rebuttal. Should you feel that the steps undertaken to clarify the paper address your concerns, we would be very grateful if you would consider raising your final score.\n\n**Presentation of numerical results.** We have removed bolding in Tables 1 and 2, except for cases that are statistically best-performing. The new tables show the mean plus minus standard error and we clearly indicate this in the revised version.\n\n**Responses to questions.**\n\n1. We have modified Figure 2 and the toy example. The main purpose of the toy example was to illustrate differences amongst how each method generated confidence intervals, as opposed to showing how well-calibrated they are. To this end, we deliberately picked a setting with little data and training inputs that are far apart, as this is where the difference between approaches is most obvious. However, we do not expect our method to be calibrated in this low-data region, as our method is only guaranteed to work in regions near training data (i.e. we assume that the test data and train data are drawn i.i.d. from the same distribution).\n\nWe realize that this setting did not clearly illustrate the calibration properties of each method, which made Figure 1 confusing. In the submitted PDF, we include a revised Figure 1 that uses a different toy dataset to better illustrated the calibration differences between each approach. Please refer to the submitted PDF for details. \n\n2. Table 1 only summarizes the average calibration error for all methods and datasets. A full table with detailed experimental results can be found in the supplementary material (it did not fit in the main text). We now state this explicitly in the table caption. Furthermore, we have moved Figure 2 to a different page from that of Table 1 to underscore that Table 1 is unrelated to Figure 2.\n\n3. The training data is not added to the vanilla GP. This is now stated explicitly in the experimental section.\n\n4. Please refer to point 1.\n\n5. The GP posterior standard deviation is considerably small in some examples, which leads to low STD values. However, this of course comes at the cost of a very poor calibration score, as can be seen in Table 1. We felt that including this data in the tables in the main paper would likely cause some confusion, since the models in this case are simply overconfident, not sharp AND calibrated. We will include the scores for the vanilla GP and fully Bayesian GP in the supplementary material. We have also included this information in the PDF submitted with the rebuttal.\n\n6. Thank you for this observation. The statement indeed holds for the more general class of stationary kernels, i.e., kernels where k(x,y) = K(x-y) for some function K. We have rewritten the statement accordingly in the revised version.\n"}, {"review_id": "hC4ZPEv3QH", "reviewer": "Reviewer_k7Z4", "summary": "Motivated by the observation that the posterior variance of a Gaussian process is often poorly calibrated, the authors propose an alternative approach of attaching predictive quantiles to the posterior mean. In essence, their approach minimises the width of the predictive quantiles under an empirical calibration constraint computed on held-out validation data. By satisfying the empirical calibration constraint, the authors prove that the predictive quantiles are indeed approximately the right quantiles (Theorem 5.3). The approach is tested and compared against other calibration approaches in a toy example and on seven data sets from the UCI data set repository. The results show that it outperforms other techniques in terms of sharpness.", "questions": "## Conclusion\n\nAlthough I think the proposed approach is very interesting, I am mainly worried about the validity of Theorem 5.3, the soundness the result in Section 7.1, the use of boldness to present the experimental results, and the STD and NLL metrics in Table 2.1. Therefore, at this point, I cannot accept this submission and must unfortunately recommend rejection. However, if the authors are able to address the following points, then I am willing to change my reject to an accept:\n\n* Please see if Theorem 5.3 is really flawed, and, if it is, whether it can be fixed.\n\n* Please correct the use of boldness throughout the main paper and the supplement.\n\n* Please investigate whether the predictions in Figure 2.(a) really are calibrated.\n\n* Please explain what the STD and NLL metrics in Table 2.1 are.\n\nEDIT\n\nThe authors have largely addressed my criticisms, so I have increased my score to a borderline accept. I believe that the submission might deserve a higher score, but I am not comfortable recommending any higher without seeing a revised version of the PDF.", "limitations": "See above.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "To begin with, I would like to thank the authors for their submission. \n\n## Strengths\n\n* The paper is easy to follow and generally well written. I found only a few typos.\n\n* The problem that the posterior variance of a GP may be poorly calibrated is highly relevant, so methods that attempt to attack this problem, like the proposed approach by the authors, are certainly important.\n\n* Accompanying the method with a theoretical result that guarantees correctness (Theorem 5.3) reassures the practitioner.\n\n* I really like Section 5.2, where the authors spend some effort on reducing the computational cost in practice. They could have just left it at (7) and state that the optimisation problem would have to be computed for every confidence level $\\delta$ of interest, but of course the approach in Section 5.2 is much more elegant.\n\n* The experiments appear to show that the proposed approach is sharper than alternatives. I, however, have some doubts about the experimental section. Please see below.\n\n## General Comments and Suggestions\n\n* Should the title of the paper be \"Sharply Calibrated GPs\" instead of \"Sharp Calibrated GPs\"? I think \"Sharp(ly)\" should be an adverb, since it modifies \"calibrated\".\n\n* On line 75, you use $\\mathcal{D}\\_{\\text{tr}}$ without introducing the symbol first. This is confusing, because you have just introduced $\\mathcal{D}$, not $\\mathcal{D}\\_{\\text{tr}}$, as the training data.\n\n* On line 75, you introduce the cut-point term as $\\beta_\\delta \\sigma_{\\mathcal{D}}(\\delta, x)$. At this point, you should really better explain this term by answering the following questions: Should $\\sigma_{\\mathcal{D}}(\\delta, x)$ be interpreted as a standard deviation? If so, how do you make sense of that it depends on $\\delta$? If not, then what is $\\sigma$? If you're learning $\\sigma$, why do you also need $\\beta_\\delta$? Can you not absorb $\\beta_\\delta$ into $\\sigma$? You should also mention that $\\beta_\\delta$ may be negative. Without that information, (1) doesn't make much sense for small $\\delta$.\n\n* On line 81: Do you mean \"are small\" or \"are as small as possible\"? This nuance is very important and changes the meaning substantially.\n\n* On line 119, you say that the log-marginal likelihood does not account for calibration. Depending on what precisely you mean by calibration, I think that this is false: the log-marginal likelihood is an empirical estimate of the KL divergence, and the KL divergence certainly accounts for the \"whole distribution\" and therefore the calibration.\n\n* In Section 5.1, is $\\sigma\\_{\\mathcal{D}\\_{\\text{tr}}}(\\theta, x)$ the posterior variance of the GP where the kernel parameters are now the parameters that we optimise over? From line 146 onwards, you don't actually explain this! Since this is a crucial part of your construction, I think the exposition would be better if you were to clearly explain this somewhere around line 158.\n\n* On line 175, you state that (5) can be replaced by (6) with further explanation. Since this is an important step in the derivation of your algorithm, I think that it deserves a careful explanation. Moreover, it is not true that (5) and (6) are equivalent, since $q\\_{\\text{lin}}$ is linearly interpolated. I think the exposition would benefit from a little more care here.\n\n* Line 268: \"calibration\" -> \"calibrated\"\n\n* Could you add to the legend of Figure 2 what the squares and crosses are? Why does it look like there are two \"lines of squares/crossed\"?", "weaknesses": "## Weaknesses\n\n### Assumption 4.1 Not Obviously Satisfied\n\nI agree that Assumption 4.1 is obviously satisfied for the variance of a kernel. However, for the inverse length scale, I can believe that Assumption 4.1 might be satisfied, but this is not obvious at all. Would the authors be able to produce a proof that the posterior variance of a GP is monotonic in the inverse length scale?\n\n### Theorem 5.3 Might Not Be Valid\n\nThe proof of Theorem 5.3 crucially relies on Theorem 1 by Marx et al. (2022). This Theorem 1, however, operates in the setting where all pairs $(x_i, y_i)$ are sampled i.i.d. (see Section 2 of Marx et al., 2022). But in the setting of GPs, which is the setting of the submission, the pairs aren't independent, because they are correlated by the sample from the underlying GP! This means that Theorem 1 might not apply, which means that Theorem 5.3 might not be valid. Could the authors comment on this?\n\n### Result of Section 7.1 Looks Questionable\n\nIn Figure 2, you state that your approach produces a 99% confidence interval. However, if you look at the right side of Figure 2.(a), then the confidence interval completely misses the data! It hence looks like the shaded region is not actually a 99% confidence interval, which makes me wonder whether the predictions are actually well calibrated.\n\n### Bolded Results Are Not Significantly Best Results\n\nThroughout the main paper and the supplement, you bold the score with the best average. However, if $x_1 \\pm e_1$ and and $x_2 \\pm e_2$ are such that $x_1 < x_2$, but $x_1 + e_1 \\ge x_2 - e_2$, then you cannot actually conclude with confidence that $x_1 \\pm e_1$ is really lower than $x_2 \\pm e_2$, because the difference might be explained by random fluctuations. In other words, you should really only bold results that are the best results at reasonable statistical significance. Currently, because of this issue, I think that the results throughout are misrepresented.\n\n### What Are STD and NLL in Table 2.1?\n\nThe premise of the paper is that you discard the predictive variance and instead produce calibrated predictive quantiles at one or multiple given confidence levels. This means that the predictions now consist of a mean and associated intervals. Therefore, the predictions are no longer probabilistic, so I really don't understand what the STD and NLL in Table 2 are! (For given a mean and an interval, how can you compute a probability?)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your review and your helpful comments. Please find our answers to your questions below. If you feel that we have adequately addressed your concerns and questions, we would appreciate if you would consider updating your score.\n\n**Assumption 4.1.** Capone et al. (2022) have shown that Assumption 4.1 holds for the lengthscale of a class of stationary kernels (Lemma 3.3). In fact, the proof presented therein can be employed to show that Assumption 4.1 holds for any hyperparameters that lead to a monotonous increase in the Fourier transform. We will include these details in the revision, and give thorough examples of hyperparameters that satisfy these properties. Note also that several results that employ the so-called \"fill distance\" show that Assumption 4.1 holds for a broad class of kernels as the lengthscale becomes very large or very small (see, e.g., Chapter 11 in *Scattered Data Approximation* by Holger Wendland, 2011). \n\n**Theorem 5.3.** We believe there may be some confusion regarding our assumptions as well as the requirements of Theorem 5.3. In a standard Bayesian setting, a GP prior would, as you point out, result in joint dependence amongst the observations. However, we note that our method is inherently frequentist in nature: we use GPs simply as mechanisms to produce kernelized mean functions and confidence intervals, rather than as a true Bayesian prior. Because of this frequentist setting, we make the standard frequentist assumptions about data: there is some unobserved (but fixed) function $f$, and $y \\mid f(\\boldsymbol x)$ is an i.i.d. distribution, meeting the requirements of Theorem 1 from Marx et al. (2022). Again, we emphasize that our model's mean and confidence intervals---which follow the same functional form as a Bayesian GP posterior---are purely for computation and do not arise from any assumed prior distribution.\n\nWe will make this more clear in the revision, and we will explicitly specify in Section 2 that we make the assumption of i.i.d. input/observation pairs.\n\n**Use of bold font.** We have rewritten the tables and changed the exposition of the results; see the included PDF. In particular, we now only use bold font whenever the results are best statistically. Furthermore, we employ mean plus minus standard error, which best illustrates statistical significance in this case, and state this clearly in the revised version.\n\n**Section 7.1/Figure 2.** The main purpose of the toy example was to illustrate differences amongst how each method generated confidence intervals, as opposed to showing how well-calibrated they are. To this end, we deliberately picked a setting with little data and training inputs that are far apart, as this is where the difference between approaches is most obvious. Note that, in this setting, our method has no calibration guarantees, as our guarantees only hold in regions sufficiently near training data (i.e. in regions where test data and train data are i.i.d. draws from the same data distribution). This has caused some confusion, so we have revised the toy example. The toy example is now calibrated and also highlights the differences between approaches. Please refer to the submitted PDF for details.\n\n**STD and NLL in Table 2.1.** Thank you for this comment. Following the method presented in Section 5.2, we can use our method to compute predictive quantiles for arbitrary $\\delta$. Hence, our model implicitly specifies a cumulative distribution function, obtained by inverting the quantile function. This allows us to compute the standard deviation and the negative log likelihood of the predictions."}]

Author rebuttal:
We kindly thank all the reviewers for their helpful comments. They have helped significantly toward improving our paper. Below you will find a summary of the changes made to the paper. We have also attached a PDF with additional experiments and the revised toy problem.

**Presentation of numerical results.** 
We have changed the presentation of the numerical results. In particular, we have removed the bolding whenever statistical evidence is insufficient to deem one method superior. We now present the numerical results as mean+- standard error, as opposed to mean+-standard deviation, as this better illustrates the expected deviation from the presented mean.

We have also modified the presentation of the toy example, as this caused some confusion. The toy example had been handpicked to illustrate the differences between methods concerning sharpness since all methods guarantee calibration. The new example presents calibrated results while clearly illustrating the differences between methods.

**Comparison to new methods.**
We have compared our approach to that of Song et al. (2019), Kuleshov and Deshpande (2022), and Capone et al. (2022). The results are presented in the PDF and detailed in the response to reviewer XLKC.

**Assumption 4.1.**
Assumption 4.1 can be shown to hold for kernels with Fourier transforms that are monotonically increasing with the hyperparameters. The revised paper shows how this can be leveraged to construct kernels that satisfy Assumption 4.1.



Meta-review: This paper proposes a novel method for calibrating the predictive probabilities of a Gaussian process. The idea is to minimize the width of confidence intervals, while satisfying a calibration constraint on holdout data. The authors argue this produces sharper, yet calibrated estimates of uncertainty.

This paper addresses an important problem, and a new and useful idea is presented. Reviewers mostly complain about clarity, but the author feedback is helping a lot. While the experiments in the original submission have some limitations, new results posted during the rebuttal clarified concerns of reviewers. If this paper gets accepted, it is strongly recommended to improve the presentation of results (see reviewer comments), and to improve the overall clarity of presentation.
Decision: Accept (poster)