Title: QUANTIFYING POINT-PREDICTION UNCERTAINTY IN NEURAL NETWORKS VIA RESIDUAL ESTIMATION WITH AN I/O KERNEL

Abstract: Neural Networks (NNs) have been extensively used for a wide spectrum of realworld regression tasks, where the goal is to predict a numerical outcome such as revenue, effectiveness, or a quantitative result. In many such tasks, the point prediction is not enough: the uncertainty (i.e. risk or confidence) of that prediction must also be estimated. Standard NNs, which are most often used in such tasks, do not provide uncertainty information. Existing approaches address this issue by combining Bayesian models with NNs, but these models are hard to implement, more expensive to train, and usually do not predict as accurately as standard NNs. In this paper, a new framework (RIO) is developed that makes it possible to estimate uncertainty in any pretrained standard NN. The behavior of the NN is captured by modeling its prediction residuals with a Gaussian Process, whose kernel includes both the NN's input and its output. The framework is evaluated in twelve realworld datasets, where it is found to (1) provide reliable estimates of uncertainty, (2) reduce the error of the point predictions, and (3) scale well to large datasets. Given that RIO can be applied to any standard NN without modifications to model architecture or training pipeline, it provides an important ingredient for building real-world NN applications.

1. INTRODUCTION:

Nowadays, Neural Networks (NNs) are arguably the most popular machine learning tool among Artificial Intelligence (AI) community. Researchers and practitioners have applied NNs to a wide variety of fields, including manufacturing (Bergmann et al., 2014) , bioinformatics (LeCun et al., 2015) , physics (Baldi et al., 2014) , finance (Niaki & Hoseinzade, 2013) , chemistry (Anjos et al., 2015) , healthcare (Shahid et al., 2019) , etc. Although standard NNs are good at making a point prediction (a single outcome) for supervised learning tasks, they are unable to provide uncertainty information about predictions. For real-world decision makers, representing prediction uncertainty is of crucial importance (Krzywinski & Altman, 2013; Ghahramani, 2015) . For example, in the case of regression, providing a 95% confidence interval around the prediction allows the decision maker to anticipate the possible outcomes with explicit probability. In contrast, simply returning a single point prediction imposes increased risks on decision making, e.g., a predictively good but actually risky medical treatment may be overconfidently interpreted without uncertainty information.

Conventional Bayesian models such as Gaussian Processes (GP) (Rasmussen & Williams, 2006) offer a mathematically grounded approach to reason about the predictive uncertainty, but often come with a prohibitive computational cost and lower prediction performance compared to NNs (Gal & Ghahramani, 2016) . As a potential solution, considerable research has been devoted to the combination of Bayesian models and NNs (see Section 2 for a detailed review of such approaches), aiming to overcome the downsides of each. However, from the classical Bayesian Neural Network (Neal, 1996) , in which a distribution of weights is learned, to the recent Neural Processes (Garnelo et al., 2018a; b; Kim et al., 2019) , in which a distribution over functions is defined, all such methods require significant modifications to the model infrastructure and training pipeline. Compared to standard (non-Bayesian) NNs, these new models are often computationally slower to train and harder to implement (Gal & Ghahramani, 2016; Lakshminarayanan et al., 2017) , creating tremendous difficulty for practical uses. Gal & Ghahramani (2016) derived a theoretical tool to extract uncertainty information from dropout training, however, the method can only be applied to dropout models, and requires changes to their internal inference pipeline. Quantifying point-prediction uncertainty in standard NNs, which are overwhelmingly popular in practical applications, still remains a challenging problem with significant potential impact.

To circumvent the above issues, this paper presents a new framework that can quantitatively estimate the point-prediction uncertainty of standard NNs without any modifications to the model structure or training pipeline. The proposed approach works as a supporting tool that can augment any pretrained NN without retraining them. The idea is to capture the behavior of the NN by estimating its prediction residuals using a modified GP, which uses a new composite (I/O) kernel that makes use of both inputs and outputs of the NNs. The framework is referred to as RIO (for Residual estimation with an I/O kernel). In addition to providing valuable uncertainty estimation, RIO has an interesting side-effect: It provides a way to reduce the error of the NN predictions. Moreover, with the help of recent sparse GP models, RIO scales well to large datasets. Since classification problems can be treated as regression on class labels (Lee et al., 2018) , this paper focuses on regression tasks. In this setting, empirical studies are conducted with twelve real-world datasets, and an initial theoretical investigation is presented that characterizes the benefits of RIO. The results show that RIO exhibits reliable uncertainty estimations, more accurate point predictions, and better scalability compared to alternative approaches.

2. RELATED WORK:

There has been significant interest in combining NNs with probabilistic Bayesian models. An early approach was Bayesian Neural Networks, in which a prior distribution is defined on the weights and biases of a NN, and a posterior distribution is then inferred from the training data (MacKay, 1992; Neal, 1996) . Notice that unlike these methods, RIO is concerned only with uncertainty over NN predictions, not over NN weights. Traditional variational inference techniques have been applied to the learning procedure of Bayesian NN, but with limited success (Hinton & van Camp, 1993; Barber & Bishop, 1998; Graves, 2011) . By using a more advanced variational inference method, new approximations for Bayesian NNs were achieved that provided similar prediction performance as dropout NNs (Blundell et al., 2015) . However, the main drawbacks of Bayesian NNs remain: prohibitive computational cost and difficult implementation procedure compared to standard NNs.

Alternatives to Bayesian NNs have been developed recently. One such approach introduces a training pipeline that incorporates ensembles of NNs and adversarial training (Lakshminarayanan et al., 2017) . Another approach, NNGP, considers a theoretical connection between NNs and GP to develop a model approximating the Bayesian inference process of wide deep neural networks (Lee et al., 2018) . Deep kernel learning (DKL) combines NNs with GP by using a deep NN embedding as the input to the GP kernel (Wilson et al., 2016) . In Iwata & Ghahramani (2017) , NNs are used for the mean functions of GPs, and parameters of both NNs and GP kernels are simultaneously estimated by stochastic gradient descent methods. Conditional Neural Processes (CNPs) combine the benefits of NNs and GP, by defining conditional distributions over functions given data, and parameterizing this dependence with a NN (Garnelo et al., 2018a) . Neural Processes (NPs) generalize deterministic CNPs by incorporating a latent variable, strengthening the connection to approximate Bayesian and latent variable approaches (Garnelo et al., 2018b) . Attentive Neural Processes (ANPs) further extends NPs by incorporating attention to overcome underfitting issues (Kim et al., 2019) . The above models all require significant modifications to the original NN model and training pipeline. Compared to standard NNs, they are also less computationally efficient and more difficult for practitioners to implement. In the approach that shares the most motivation with RIO, Monte Carlo dropout was used to estimate the predictive uncertainty of dropout NNs (Gal & Ghahramani, 2016) . However, this method is restricted to dropout NNs, and also requires modifications to the NN inference process. Different from all above-mentioned methods, which are independent alternatives to standard NNs, RIO is designed as a supporting tool that can be applied on top of any pretrained NN. RIO augments standard NNs without retraining or modifying any component of them.

3. THE RIO FRAMEWORK:

This section gives the general problem statement, develops the RIO framework, and discusses its scalability. For background introductions of NNs, GP, and its more efficient approximation, SVGP (Hensman et al., 2013; 2015) , see Appendix B.

3.1. PROBLEM STATEMENT:

Consider a training dataset D = (X , y) = {(x i , y i )} n i=1 , and a pretrained standard NN that outputs a point prediction ŷi given x i . The problem is two-fold: (1) Quantify the uncertainty in the predictions of the NN, and (2) calibrate the point predictions of the NN (i.e. make them more accurate).

3.2. FRAMEWORK OVERVIEW:

RIO solves this problem by modeling the residuals between observed outcomes y and NN predictions ŷ using GP with a composite kernel. The framework can be divided into two phases: the training phase and the deployment phase.

In the training phase, the residuals between observed outcomes and NN predictions on the training dataset are calculated as

EQUATION

Let r denote the vector of all residuals and ŷ denote the vector of all NN predictions. A GP with a composite kernel is trained assuming r ∼ N (0, K c ((X , ŷ), (X , ŷ)) + σ 2 n I), where K c ((X , ŷ), (X , ŷ)) denotes an n × n covariance matrix at all pairs of training points based on a composite kernel

k c ((x i , ŷi ), (x j , ŷj )) = k in (x i , x j ) + k out (ŷ i , ŷj ), for i, j = 1, 2, . . . , n .

(2) Suppose a radial basis function (RBF) kernel is used for both k in and k out . Then,

EQUATION

The training process of GP learns the hyperparameters σ 2 in , l in , σ 2 out , l out , and σ 2 n by maximizing the log marginal likelihood log p(r|X , ŷ) given by The neural network has discovered true complex structure in the labels, so the residuals have low variance and are easy for the GP to fit with high confidence; (right) The ineffective neural network has introduced unnecessary complexity, so the residuals are modeled with high uncertainty. In both cases, RIO matches the intuition for how uncertain the NN really is.

EQUATION

In the deployment phase, a test point x * is input to the NN to get an output ŷ * . The trained GP predicts the distribution of the residual as r * |X , ŷ, r, x * , ŷ * ∼ N ( r * , var(r * )), where

EQUATION

EQUATION

where k * denotes the vector of kernel-based covariances (i.e., k c ((x * , ŷ * ), (x i , ŷi ))) between (x * , ŷ * ) and the training points.

Interestingly, note that the predicted residuals can also be used to calibrate the point predictions of the NN, so that the final calibrated prediction with uncertainty information is given by

EQUATION

In other words, RIO not only adds uncertainty estimation to a standard NN-it also provides a way to calibrate NN predictions, without any modification to its architecture or training. Figure 1 shows the overall procedure when applying the proposed framework in real-world applications. Figure 2 shows example behavior of RIO that illustrates the intuition of the approach. An algorithmic description of RIO is provided in Appendix C.

3.3. SCALABILITY:

RIO is scalable to large datasets (No. of data points×No. of features>1M) by applying sparse GP methods, e.g., SVGP (Hensman et al., 2013; 2015) . All the conclusions in previous sections still remain valid since sparse GP is simply an approximation of the original GP. In the case of applying SVGP with a traditional optimizer, e.g., L-BFGS-B (Byrd et al., 1995; Zhu et al., 1997) , the computational complexity is O(nm 2 ), and space complexity is O(nm), where n is the number of data points and m is the number of inducing variables, compared to O(n 3 ) and O(n 2 ) for traditional GP. Experiments verify that the computational cost of RIO with SVGP is significantly cheaper than other state-of-the-art approaches.

4. EMPIRICAL EVALUATION:

Experiments in this section compare nine algorithms on 12 real-world datasets. The algorithms include standard NN, the proposed RIO framework, four ablated variants of RIO, and three state-ofthe-art models that provide predictive uncertainty: SVGP (Hensman et al., 2013) , Neural Network Gaussian Process (NNGP) (Lee et al., 2018) , and Attentive Neural Processes (ANP) (Kim et al., 2019) .

In naming the RIO variants,"R" means estimating NN residuals then correcting NN outputs,"Y" means directly estimating outcomes, "I" means only using input kernel, "O" means only using output kernel, and "IO" means using I/O kernel. For all RIO variants (including full RIO), SVGP is used as the GP component, but using the appropriate kernel and prediction target. Therefore, "Y+I" amounts to original SVGP, and it is denoted as "SVGP" in all the experimental results. All 12 datasets The symbols † and ‡ indicate that the difference between the marked entry and RIO is statistically significant at the 5% significance level using paired t-test and Wilcoxon test, respectively. The best entries that are significantly better than all the others under at least one statistical test are marked in boldface (ties are allowed).

are real-world regression problems (Dua & Graff, 2017) , and cover a wide variety of dataset sizes and feature dimensionalities. Except for the "MSD" dataset, all other datasets are tested for 100 independent runs. During each run, the dataset is randomly split into training set, validation set, and test set, and all algorithms are trained on the same split. All RIO variants that involve an output kernel or residual estimation are based on the trained NN in the same run. For "MSD", since the dataset split is strictly predefined by the provider, only 10 independent runs are conducted. NNGP and ANP are only tested on the four smallest dataset (based on the product of dataset size and feature dimensionality) because they do not scale well to larger datasets. It is notable that for all the RIO variants, no extensive hyperparameter tuning is conducted; the same default setup is used for all experiments, i.e., standard RBF kernel and 50 inducing points. See Appendix D.1 for additional details of experimental setups and a link to downloadable source code. Table 1 summarizes the numerical results from these experiments. The main conclusions in terms of point-prediction error, uncertainty estimation, computational requirements, and ablation studies are summarized below.

Point-Prediction Error:

The errors between point predictions of models and true outcomes of test points are measured using Root Mean Square Error (RMSE); the mean and standard deviation of RMSEs over multiple experimental runs are shown in Table 1 . For models that return a probabilistic distribution, the mean of the distribution is the point prediction. Although the main motivation of RIO is to enhance pretrained NN rather than construct a new state-of-the-art prediction model from scratch, RIO performs the best or equals the best method (based on statistical tests) in 10 out of 12 datasets. RIO significantly outperforms original NN in all 12 datasets, while original SVGP performs significantly worse than NN in 7 datasets. For the CT dataset, which has 386 input features, SVGP fails severely since the input kernel cannot capture the implicit correlation information. ANP is unstable on the airfoil dataset because it scales poorly with dataset size. Figure 3 Uncertainty Estimation Average negative log predictive density (NLPD) is used to measure the quality of uncertainty estimation, which favours conservative models, but also effectively penalizes both over-and under-confident predictions (Quiñonero-Candela et al., 2006) . The mean and standard deviation of NLPDs over multiple experimental runs are shown in Table 1 (lower is better). RIO performs the best or equals the best method (based on statistical tests) in 8 out of 12 datasets. NNGP always yields a high NLPD; it returns over-confident predictions, because it does not include noise estimation in its original implementation. For the yacht dataset, ANP achieves the best NLPD, but with high RMSE. This is because ANP is able to correctly return high predictive variance when its prediction error is high. For all other tested datasets, RIO variants consistently outperform ANP. Among all RIO variants, the full RIO provides the most reliable overall predictive uncertainty. The conclusion is that RIO successfully extracts useful uncertainty information from NN predictions. Appendix D.2 evaluates uncertainty estimates on an additional metric: the true coverage probabilities of estimated confidence intervals. Although RIO also provides reasonable estimates for this metric in most cases, comparison among algorithms using this metric is discouraged due to its unreliability (see Appendix D.2 for examples).

Computation Time Table 1 shows the average wall clock time of each algorithm. All algorithms are implemented using Tensorflow under the same running environment (see Appendix D.1 for implementation details). The RIO variants scale well to increasing dataset sizes and feature dimensionalities. L-BFGS-B converges especially quickly for R+I on the three highest dimensional datasets, presumably because the residuals are very well-behaved compared to the raw targets or NN output.

In contrast, ANP's computation time increases significantly with the scale of the dataset, and NNGP always needs very expensive computational budgets due to its costly grid search of hyperparameters.

Ablation Study:

The RIO variants with residual estimation generally perform better than its counterparts in both point-prediction error and uncertainty estimation. This result confirms the effectiveness of residual estimation, as suggested in Section 5.1. Another important result is that Y+IO outperforms both Y+I (SVGP) and Y+O in most cases across all performance metrics, and RIO generally provides better performance than R+I and R+O in all respects. This result, in turn, confirms that the I/O kernel provides additional robustness, as suggested by the analysis in Section 5.2. In sum, both residual estimation and the I/O kernel contribute substantially to the performance of the framework. Correlation between error and noise variance Table 2 shows the results of Spearman's Rank Correlation between RMSE and noise variance σ 2 n . For each dataset, there are six pairs of data points, each of which contains the mean of RMSE and noise variance of a RIO variant. A correlation value larger than 0.8 and p-value less than 0.05 indicate a very strong positive monotonic correlation. For 10 out of 12 datasets, very strong positive monotonic correlation between RMSE and noise variance was observed. This empirical result is in accordance with the theoretical prediction in Section 5.1.

Application to a large-scale vision problem:

To show RIO's off-the-shelf applicability to modern deep convolutional architectures, it was applied to a recent pretrained NN for age estimation (Yang et al., 2018) . The pretrained NN and all data preprocessing were taken exactly from the official code release. The model is a variant of DenseNet-121 (Huang et al., 2017) , and uses the IMDB age estimation dataset, which contains ≈ 172K RGB images (Rothe et al., 2016) . The goal is to predict the age of the individual in each image. The features for the GP input kernel were simply a global max pool of the CNN's first stage output. From Table 3 , RIO substantially improves upon the mean absolute error (MAE) of the pretrained NN, outperforms SVGP in terms of both MAE and NLPD, and yields realistic confidence intervals (CIs). SVGP effectively learns nothing, so it estimates almost all variance as noise, while RIO effectively augments the pretrained model. 

5. ANALYSIS OF RIO:

This section presents an initial theoretical investigation into what makes RIO work, including why it is useful to fit residuals, and why it is beneficial to use an I/O kernel.

5.1. BENEFITS OF FITTING NEURAL NETWORK RESIDUALS:

Given a pretrained NN, this section provides a theoretical perspective on how fitting its residuals with a GP can yield useful uncertainty information, while leading to prediction error lower than the NN or GP alone. The idea is that, due to its high expressivity, a pretrained NN may have learned complex structure that a GP would model poorly, when selecting a kernel that would capture this complexity is infeasible. Fitting a GP to the residuals of the NN is easier, since this complex structure has been removed, and takes advantage of the predictive power of the NN, while providing useful uncertainty estimates.

Given a pretrained NN, when is fitting its residuals with a GP a good idea? To produce a model with uncertainty information, one could simply discard the NN and train a GP directly on D. However, for problems with enough complexity, GP model selection (i.e., specifying an appropriate class of kernels) is challenging. Fortunately, the pretrained NN may have learned a useful representation of this complexity, which can then be exploited. Suppose instead that a GP tries to model this complexity, and this attempt leads to worse generalization. This poor performance could come from many dimensions of the kernel selection process. As an initial theoretical investigation, an ideal case is considered first: suppose the GP optimizer optimally avoids modeling this complexity incorrectly by estimating the variance of such complexity as noise. Such an estimate leads to the following decomposition:

EQUATION

where h(•) is the true signal, f (•) is the apparent signal of the GP, g(•) is its apparent (spurious) noise, and ξ i ∼ N (0, σ 2 n ) is real noise. Intuitively, due to its high expressivity, it is possible for a pretrained NN hNN to correctly model part of g, along with part of f , resulting in the residuals

EQUATION

where r f (•) is now the apparent signal of the GP, and r g (•) is now its apparent noise. In such a case, it will be easier for the GP to learn the the residuals than to learn h directly. The uncertainty estimates of the GP will also be more precise, since the NN removes spurious noise. Thus, Eq. 8 immediately predicts that the noise estimates for GP will be higher when fitting h directly than when fitting NN residuals. This prediction is confirmed in experiments (Table 1 ).

To analyze these behaviors, the above decomposition is used to identify a class of scenarios for which fitting residuals leads to improvements (see Appendix A for additional details and proofs of following theorems). For simplicity, suppose GP kernels are parameterized only by their signal variance β, i.e., kernels are of the form βk(•, •) for some kernel k. Then, the following two assumptions make Eq. 8 concrete:

(1) GP is well-suited for f , i.e., f

(•) ∼ GP(0, k(•, •));

(2) GP is ill-suited for g, i.e., g is independent of f and ε-indistinguishable for GP on D (see Lemma A.2 for an example of such a function). Intuitively, this means the predictions of GP change by no more than ε in the presence of g.

Consider the residuals y:

i -hNN (x i ) = r(x i ) + ξ i = r f (x i ) + r g (x i ) + ξ i

, where r f is the remaining GP component, and r g is the remainder of g. The following assumption ensures that r g captures any spurious complexity added by an imperfectly trained NN:

r f ∼ GP(0, αk(•, •)), for some α ∈ (0, 1].

Let hGP be the GP predictor trained directly on D, and hGP+NN = hNN + rGP be the final predictor after fitting residuals; let E h GP , E h NN , and E h GP+NN be the expected generalization errors of hGP , hNN , and hGP+NN , resp. The advantage of hGP+NN becomes clear as g becomes complex:

Theorem 5.1 (Advantage of fitting residuals).

lim ε→0 E h GP -E h GP+NN ≥ E[g 2 (x)] -E[r 2 g (x)] and lim ε→0 E h NN -E h GP+NN > 0.

The experimental results in Section 4 support this result (Table 1 ). First, in all experiments, GP+NN does indeed improve over the underlying NN. Second, notice that the improvement of GP+NN over GP is only guaranteed when

E[g 2 (x)] -E[r 2 g (x)

] is positive, i.e., when the NN successfully captures some underlying complexity. In experiments, when the NN is deficient, GP+NN often outperforms GP, but not always. For example, if the NN is overfit to zero training loss, and uses the same training data as the GP, the GP cannot provide benefit. However, when the NN is well-trained, using proper overfitting prevention, the improvements are consistent and significant.

The prediction of Theorem 5.1 that lower apparent noise corresponds to improved error is confirmed in experiments (Table 2 ). Note that an optimal model would learn h exactly, and have a predictive variance of σ 2 n . So, the improvement in error of hGP+NN over hGP also corresponds to a predictive variance that is closer to that of the optimal model, since some spurious noise has been removed by hNN . Complementing this theoretical improvement in uncertainty estimation, the above setting also leads to a key practical property, which is illustrated in Figure 2 : Theorem 5.2. The uncertainty of rGP is positively correlated with the variance of NN residuals.

This property matches the intuition that the GP's variance should generally be higher for unstable NNs, i.e., NNs with high residual variance, than for stable NNs. Such a property is crucial to measuring the confidence of NN predictions in practice.

Overall, the initial theoretical results in this section (and the appendix) show that if a problem is well-suited to be learned by an NN, and a practitioner has access to a pretrained NN for that problem, then training a GP to fit the residuals of the NN can provide potential benefits in uncertainty estimates, without sacrificing prediction performance. Further theoretical investigation regarding more complicated situations, e.g., misspecification of GP kernels, is considered as one interesting direction for future work.

5.2. ROBUSTNESS OF THE I/O KERNEL:

This section provides a justification for why a GP using the proposed I/O kernel is more robust than the standard GP, i.e., using the input kernel alone. The key assumption is that the output of an NN can contain valuable information about its behavior, and, consequently, the structure of the target function. Consider the setup in Section 5.1, but now with y i = f in (x i ) + f out (x i ) + ξ i , where f in ∼ GP (0, k in ) and f out ∼ GP (0, k out ), with non-trivial RBF kernels k in , k out (as in Equation 3). Let hNN be an NN of sufficient complexity to be nontrivially non-affine, in that there exists a positive-measure set of triples (x, x , x ) s.t. xx = xx , but hNN (x) -hNN (x ) = hNN (x) -hNN (x ) . Denote the generalization errors of the standard GP, GP with output kernel only, and GP with I/O kernel by E h I , E h O , and E h I/O , respectively. The expected result follows (proof in Appendix A).

Theorem 5.3 (Advantage of I/O kernel). E h I/O < E h I and E h I/O < E h O .

The optimizer associated with the GP simultaneously optimizes the hyperparameters of both kernels, so the less useful kernel usually receives a smaller signal variance. As a result, the I/O kernel is resilient to failures of either kernel. In particular, the GP using I/O kernel improves performance even in the case where the problem is so complex that Euclidean distance in the input space provides no useful correlation information or when the input space contains some noisy features. Conversely, when the NN is a bad predictor, and h NN is no better noise, the standard GP with input kernel alone is recovered. In other words, the I/O kernel is never worse than using the input kernel alone, and in practice it is often better. This conclusion is confirmed in the experiments in Section 4.

6. DISCUSSION AND FUTURE DIRECTIONS:

In addition to the reliable uncertainty estimation, accurate point prediction, and good scalability, demonstrated in Section 4, RIO provides other important benefits.

RIO can be directly applied to any standard NN without modification to the model architecture or training pipeline. Moreover, retraining of the NN or change of inference process are not required.

The framework simply requires the outputs of an NN; it does not need to access any internal structure. This feature makes the framework more accessible to practitioners in real-world applications, e.g., data scientists can train NNs using traditional pipelines, then directly apply RIO to the trained NNs.

RIO also provides robustness to a type of adversarial attack. Consider a worst-case scenario, in which an adversary can arbitrarily alter the output of the NN with minuscule changes to the input. It is well-known that there are NNs for which this is possible (Goodfellow et al., 2015) . In this case, with the help of the I/O kernel, the model becomes highly uncertain with respect to the output kernel. A confident prediction then requires both input and output to be reasonable. In the real world, a high degree of uncertainty may meet a threshold for disqualifying the prediction as outside the scope of the model's ability.

There are several promising future directions for extending RIO: First, applying RIO to reinforcement learning (RL) algorithms, which usually use standard NNs for reward predictions, would allow uncertainty estimation of the future rewards. Agents can then directly employ efficient exploration strategies, e.g., bandit algorithms (Thompson, 1933) , rather than traditional stochastic approaches like ε-greedy. Second, RIO applied to Bayesian optimization (BO) (Močkus, 1975) would make it possible to use standard NNs in surrogate modeling. This approach can potentially improve the expressivity of the surrogate model and the scalability of BO. Third, since RIO only requires access to the inputs and outputs of NNs, it could be directly applied to any existing prediction models, including hybrid and ensemble models. For example, experiments in Appendix D.6 show that RIO achieves good results when directly applied to random forests. This general applicability makes RIO a more general tool for real-world practitioners.

7. CONCLUSION:

The RIO framework both provides estimates of predictive uncertainty of neural networks, and reduces their point-prediction errors. The approach captures NN behavior by estimating their residuals with an I/O kernel. RIO is theoretically-grounded, performs well on several real-world problems, and, by using a sparse GP approximation, scales well to large datasets. Remarkably, it can be applied directly to any pretrained NNs without modifications to model architecture or training pipeline. Thus, RIO can be used to make NN regression practical and powerful in many real-world applications. A ADDITIONAL DETAILS AND PROOFS FOR SECTIONS 5.1 AND 5.2

Section 5.1 provides a theoretical perspective on how using GP to fit NN residuals can provide advantages over GP or NN alone. The idea is to decompose the target function into two independent parts: one that is well-suited for a GP to learn, and the other that is difficult for it to capture. The following definition provides a formalization of this notion of difficulty. Definition A.1 (ε-indistinguishable). Let fA be the predictor of a learner A trained on {(x i , y i )} n i=1 , and hA that of

A trained on D = {(x i , y i + g(x i ))} n i=1 . Then, if | fA (x) -hA (x)| < ε ∀ x, we say g is ε-indistinguishable for A on D. 1

In other words, the predictions of A change by no more than ε in the presence of g. If ε is small compared to E[g 2 (x)], then g contains substantial complexity that A cannot capture.

For notation, let ψA denote the predictor of a learner A trained on {(x i , ψ(

x i ) + ξ i )} n i=1 , with ξ i ∼ N (0, σ 2 n ). Let E ψ A = E[(ψ(x) -ψA (x)) 2 ] denote the expected generalization error of A on ψ given training locations {x i } n i=1 . Now, consider a dataset D = {(x i , y i )} n i=1 with y i = h(x i ) + ξ i = f (x i ) + g(x i ) + ξ i , where f (•) ∼ GP(0, k(•, •)), ξ i ∼ N (0, σ 2 n ), and E[g 2 (x)] = σ 2 g .

For simplicity, assume the GP learner selects optimal hyperparameters for k and σ 2 n . Suppose g is independent of f and εindistinguishable for GP on D, but a neural network hNN may have successfully learned some of g's structure. These assumptions induce a concrete class of examples of the decomposition in Eq. 8. To demonstrate the existence of such decompositions, Lemma A.2 provides a constructive example of such a function g. Lemma A.2 (Construction of independent ε-indistinguishable functions for GP). Given a dataset {(x i , f (x i ) + ξ i )} n i=1 , for any σ 2 g and ε > 0, there exists g with E[g 2 (x)] = σ 2 g , and

{x i } 2n i=n+1 , s.t. g is ε-indistinguishable on {(x i , f (x i ) + g(x i ) + ξ i )} 2n

i=1 for GP with a continuous kernel k.

Proof. Let x n+i = x i + v ∀ i = 1, . . . , n, for a fixed vector v. Consider the linear smoother view of GP prediction (Rasmussen & Williams, 2006) : et al., 1984) , and y is the target vector. Since k and w(x * ) are continuous around each

fGP (x * ) = w(x * ) y, where w(x * ) = (K + σ 2 n I) -1 k(x * ) is the weight function (Silverman

x i , we can choose v such that ∀ x * |w(x * ) i -w(x * ) n+i | < ε nσ 2 g ∀ i = 1, . . . , n.

Let fGP be the predictor of GP trained on {(x i , f

(x i ) + ξ i )} 2n

i=1 , and hGP the predictor trained on

{(x i , f (x i ) + g(x i ) + ξ i )} 2n i=1 = D. Now, choose g so that g(x i ) = σ 2 g if i ≤ n, -σ 2 g if i > n. Then, hGP (x * ) = w(x * ) (f (x i ) + g(x i ) + ξ i ) 2n i=1 = w(x * ) (f (x i ) + ξ i ) 2n i=1 + w(x * ) (g(x i )) 2n i=1 = fGP (x * ) + n i=1 σ 2 g w(x * ) i - n i=1 σ 2 g w(x * ) n+i = fGP (x * ) + σ 2 g n i=1 w(x * ) i -w(x * ) n+i .

Published as a conference paper at ICLR 2020

σ 2 g n i=1 w(x * ) i -w(x * ) n+i < σ 2 g n( ε nσ 2 g ) = ε =⇒ | fGP (x * ) -hGP (x * )| < ε. Consider the residuals y i -hNN (x i ) = r(x i ) + ξ i = r f (x i ) + r g (x i ) + ξ i . Here, r f is the remaining GP component, i.e., r f ∼ GP(0, αk(•, •)), for 0 < α ≤ 1. Similarly, r g is the remainder of g ε-indistinguishable for GP on {(x i , r(x i ) + ξ i )} n i=1 , i.e., σ 2 g -E[r 2 g (x)] = δ.

These two assumptions about r f and r g ensure that any new spurious complexity added by the NN is collected in r g , since the biggest advantage of the NN (its expressivity) is also its greatest source of risk. The final predictor after fitting residuals is then hGP+NN = hNN + rGP . The following sequence of results captures the improvement due to residual estimation.

Lemma A.3 (Generalization of GP on h). E f GP + σ 2 g -2ε(E f GP ) 1 2 -2εσ g < E h GP < E f GP + σ 2 g + 2ε(E f GP ) 1 2 + 2εσ g + ε 2 .

Proof. Let ∆ fGP (x) = hGP (x) -fGP (x) denote the change in GP prediction due to g. Then,

E h GP = E[(h(x) -hGP (x)) 2 ] = E[((f (x) + g(x)) -( fGP (x) + ∆ fGP (x))) 2 ] = E[((f (x) -fGP (x)) + (g(x) -∆ fGP (x))) 2 ] = E f GP + 2E[(f (x) -fGP (x))(g(x) -∆ fGP (x))] + E[(g(x) -∆ fGP (x)) 2 ] = E f GP -2E[(f (x) -fGP (x))∆ fGP (x))] + E[(g(x) -∆ fGP (x)) 2 ]

, where the last line makes use of the fact that f and g are independent. Now,

2E[(f (x) -fGP (x))∆ fGP (x))] < 2ε(E f GP ) 1 2 , and E[(g(x) -∆ fGP (x)) 2 ] = σ 2 g -2E[g(x)∆ fGP (x)] + E[∆ fGP (x) 2 ], where 2E[g(x)∆ fGP (x)] < 2εσ g and 0 ≤ E[∆ fGP (x) 2 ] < ε 2 . So, E f GP + σ 2 g -2ε(E f GP ) 1 2 -2εσ g < E h GP < E f GP + σ 2 g + 2ε(E f GP ) 1 2 + 2εσ g + ε 2 . Lemma A.4 (Generalization of NN). E h NN = αE[f 2 (x)] + σ 2 g -δ.

Proof. Making use of the fact that E[r f ] = 0,

EQUATION

σ 2 g -E[r 2 g (x)] = δ =⇒ E[r 2 g (x)] = σ 2 g -δ. So, E h NN = αE[f 2 (x)] + σ 2 g -δ. Lemma A.5 (Generalization of GP fitting residuals). E r f GP +σ 2 g -δ-2ε(E r f GP ) 1 2 -2ε(σ 2 g -δ) 1 2 < E h GP+NN < E r f GP +σ 2 g -δ+2ε(E r f GP ) 1 2 +2ε(σ 2 g -δ) 1 2 +ε 2 .

Proof. Let ∆r f GP (x) = rGP (x) -rfGP (x) denote the change in GP prediction due to r g . Then,

E h GP+NN = E[(h(x) -hGP+NN (x)) 2 ] = E[(h(x) -hNN (x) -rGP (x)) 2 ] = E[(r f (x) + r g (x)) -(r f GP (x) + ∆r f GP (x))) 2 ] = E[(r f (x) -rfGP (x)) + (r g (x) -∆r f GP (x))) 2 ] = E r f GP + 2E[(r f (x) -rGP (x))(r g (x) -∆r f GP (x))] + E[(r g (x) -∆r f GP (x)) 2 ] = E r f GP -2E[(r f (x) -rGP (x)∆r f GP (x)] + E[(r g (x) -∆r f GP (x)) 2 ]. Similar to the case of Lemma A.3, 2E[(r f (x) -rGP (x))(r g (x) -∆r f GP (x))] < 2ε(E r f GP ) 1 2 , and E[(r g (x) -∆r f GP (x)) 2 ] = σ 2 g -δ -2E[r g (x)∆r f GP (x)] + E[∆r f GP (x) 2 ], where 2E[r g (x)∆r f GP (x)] < 2ε(σ 2 g -δ) 1 2 and 0 ≤ E[∆r f GP (x) 2 ] < ε 2 . So, E r f GP +σ 2 g -δ-2ε(E r f GP ) 1 2 -2ε(σ 2 g -δ) 1 2 < E h GP+NN < E r f GP +σ 2 g -δ+2ε(E r f GP ) 1 2 +2ε(σ 2 g -δ) 1 2 +ε 2 .

Lemma A.6 (Generalization of GP on f and r f ). From the classic result (Opper & Vivarelli, 1999; Sollich, 1999; Rasmussen & Williams, 2006 ): Consider the eigenfunction expansion k(x, x ) = j λ j φ j (x)φ j (x ) and k(x, x )φ i (x)p(x)dx = λ i φ i (x ). Let Λ be the diagonal matrix of the eigenvalues λ j , and Φ be the design matrix, i.e., Φ ji = φ j (x i ). Then,

E f GP = tr(Λ -1 + σ -2 n ΦΦ ) -1 and E r f GP = tr(α -1 Λ -1 + σ -2 n ΦΦ ) -1 . Theorem 5.1. lim ε→0 E h GP -E h GP+NN ≥ δ and lim ε→0 E h NN -E h GP+NN > 0.

Proof. From Lemmas A.3, A.4, and A.5 we have, resp.:

lim ε→0 E h GP = E f GP + σ 2 g , lim ε→0 E h NN = αE[f 2 (x)] + σ 2 g -δ, and lim ε→0 E h GP+NN = E r f GP + σ 2 g -δ.

From Lemma A.6, we have

E f GP -E r f GP = tr(Λ -1 + σ -2 n ΦΦ ) -1 -tr(α -1 Λ -1 + σ -2 n ΦΦ ) -1 ≥ 0, and αE[f 2 (x)] -E r f GP = αE[f 2 (x)] -tr(α -1 Λ -1 + σ -2 n ΦΦ ) -1 > 0. So, lim ε→0 E h GP -E h GP+NN ≥ δ and lim ε→0 E h NN -E h GP+NN > 0.

Theorem 5.2. The variance of NN residuals is positively correlated with the uncertainty of r GP .

Proof.

Increases in E[r 2 f (x)] lead to increases in α; increases in E[r 2 g (x)

] lead to decreases in δ, and thus increases in the estimated noise level σ2 n . So, an increase in either E[r 2 f (x)] or E[r 2 g (x)] leads to an increase in αk((x * , ŷ * ), (x * , ŷ * )) -αk * (αK((X , ŷ), (X , ŷ)) + σ2 n I) -1 αk * , which is the predictive variance of r GP .

Theorem 5.3.

E h I/O < E h I and E h I/O < E h O .

Proof.

x

-x = x -x =⇒ k in (x, x ) = k in (x, x ). hNN (x) -hNN (x ) = hNN (x) -hNN (x ) =⇒ k out (x, x ) = k out (x, x ).

These are true for all hyperparameter settings of k in or k out , since both are RBF kernels. So, there is no hyperparameter setting of k in that yields

k in (x 1 , x 2 ) = k in (x 1 , x 2 ) + k out (x 1 , x 2 ) ∀ x 1 , x 2 .

Similarly, there is no hyperparameter setting of k out that yields

k out (x 1 , x 2 ) = k in (x 1 , x 2 ) + k out (x 1 , x 2 ) ∀ x 1 , x 2 .

Since, neither the input nor output kernel alone can correctly specify the kernel over all sets of positive measure, their generalization error is greater than the Bayes error, which is achieved by the correct kernel (Sollich, 2002) 

, i.e., k in (x 1 , x 2 ) + k out (x 1 , x 2 ).

B BACKGROUND:

This section reviews notation for Neural Networks, Gaussian Process, and its more efficient approximation, SVGP. The RIO method, introduced in Section 3 of the main paper, uses Gaussian Processes to estimate the uncertainty in neural network predictions and reduces their point-prediction errors. 

i : (x i , y i ) ∈ D, x i = [x 1 i , x 2 i , . . . , x din i ] | i = 1, 2, . . . , n} and Y = {y i : (x i , y i ) ∈ D, y i = [y 1 i , y 2 i , . . . , y dout i ] | i = 1, 2, .

. . , n} denote the inputs and outputs (i.e., targets). A fully-connected feed-forward neural network with L hidden layers of width N l (for layer l = 1, 2, . . . , L) performs the following computations: Let z j l denote the output value of jth node in lth hidden layer given input x i , then z j l = φ( din k=1 w j,k l x k i + b j l ), for l = 1 and z j l = φ(

N l-1 k=1 w j,k l z k l-1 + b j l )

, for l = 2, . . . , L, where w j,k l denotes the weight on the connection from kth node in previous layer to jth node in lth hidden layer, b j l denotes the bias of jth node in lth hidden layer, and φ is a nonlinear activation function. The output value of jth node in output layer is then given by ŷj

i = N L k=1 w j,k out z k L + b j out , where w j,k

out denotes the weight on the connection from kth node in last hidden layer to jth node in output layer, and b j out denotes the bias of jth node in output layer. A gradient-based optimizer is usually used to learn the weights and bias given a pre-defined loss function, e.g., a squared loss function L = 1 n n i=1 (y i -ŷi ) 2 . For a standard NN, the learned parameters are fixed, so the NN output ŷi is also a fixed point. For a Bayesian NN, a distribution of the parameters is learned, so the NN output is a distribution of ŷi . However, a pretrained standard NN needs to be augmented, e.g., with a Gaussian Process, to achieve the same result.

B.2 GAUSSIAN PROCESS:

A Gaussian Process (GP) is a collection of random variables, such that any finite collection of these variables follows a joint multivariate Gaussian distribution (Rasmussen & Williams, 2006) . Given a training dataset X = {x i | i = 1, 2, . . . , n} and y = {y i = f (x i ) + | i = 1, 2, . . . , n}, where denotes additive independent identically distributed Gaussian noise, the first step for GP is to fit itself to these training data assuming y ∼ N (0, K(X , X ) + σ 2 n I), where N denotes a multivariate Gaussian distribution with mean 0 and covariance matrix K(X , X ) + σ 2 n I. K(X , X ) denotes the kernel-based covariance matrix at all pairs of training points with each entry k i,j = k(x i , x j ), and σ 2 n denotes the noise variance of observations. One commonly used kernel is the radial basis function (RBF) kernel, which is defined as k(x i , x j ) = σ 2 f exp(-

1 2l 2 f x i -x j 2

). The signal variance σ 2 f , length scale l f and noise variance σ 2 n are trainable hyperparameters. The hyperparameters of the covariance function are optimized during the learning process to maximize the log marginal likelihood log p(y|X ).

After fitting phase, the GP is utilized to predict the distribution of label y * given a test point x * . This prediction is given by y * |X , y, x * ∼ N (ȳ * , var(y * )) with ȳ * = k * (K(X , X ) + σ 2 n I) -1 y and var(y

* ) = k(x * , x * ) -k * (K(X , X ) + σ 2 n I) -1 k * ,

where k * denotes the vector of kernel-based covariances (i.e., k(x * , x i )) between x * and all the training points, and y denotes the vector of all training labels. Unlike with NN, the uncertainty of the prediction of a GP is therefore explicitly quantified.

B.3 SVGP:

The main limitation of the standard GP, as defined above, is that it is excessively expensive in both computational and storage cost. For a dataset with n data points, the inference of standard GP has time complexity O(n 3 ) and space complexity O(n 2 ). To circumvent this issue, sparse GP methods were developed to approximate the original GP by introducing inducing variables (Csató & Opper, 2002; Seeger et al., 2003; Quiñonero Candela & Rasmussen, 2005; Titsias, 2009) . These approximation approaches lead to a computational complexity of O(nm 2 ) and space complexity of O(nm), where m is the number of inducing variables. Following this line of work, SVGP (Hensman et al., 2013; 2015) further improves the scalability of the approach by applying Stochastic Variational Inference (SVI) technique, as follows:

Consider the same training dataset and GP as in section B.2, and assume a set of inducing variables as Z = {z i | i = 1, 2, . . . , m} and U = {u i = f (z i ) + | i = 1, 2, . . . , m} (f (•) and are unknown). calculate covariance matrix K c ((X , ŷ), (X , ŷ)), where each entry is given by k c ((x i , ŷi ), (x j , ŷj )) = k in (x i , x j ) + k out (ŷ i , ŷj ), for i, j = 1, 2, . . . , n 4:

optimize GP hyperparameters by maximizing log marginal likelihood log p(r|X ,

ŷ) = -1 2 r (K c ((X , ŷ), (X , ŷ)) + σ 2 n I)r -1 2 log |K c ((X , ŷ), (X , ŷ)) + σ 2 n I| -n 2 log 2π Deployment Phase: 5: calculate residual mean r * = k * (K c ((X , ŷ), (X , ŷ)) + σ 2 n I) -1 r and residual variance var(r * ) = k c ((x * , ŷ * ), (x * , ŷ * )) -k * (K c ((X , ŷ), (X , ŷ)) + σ 2 n I) -1 k * 6: return distribution of calibrated prediction ŷ * ∼ N (ŷ * + r * , var(r * ))

SVGP learns a variational distribution q(U) by maximizing a lower bound of log p(y|X ), where log p(y|X ) = log p(y|U, X )p(U)dU and p(•) denotes the probability density under original GP. Trainable hyperparameters during the learning process include values of z i and hyperparameters of the covariance function of original GP. Given a test point x * , the predictive distribution is then given by p(y * |x * ) = p(y * |U, x * )q(U)dU, which still follows a Gaussian distribution. One advantage of SVGP is that minibatch training methods (Le et al., 2011) can be applied in case of very large dataset. Suppose the minibatch size is m and m m , then for each training step/iteration, the computational complexity is O(m m 2 ), and the space complexity is O(m m). For full details about SVGP, see Hensman et al. (2013) . Since NNs typically are based on training with relatively large datasets, SVGP makes it practical to implement uncertainty estimates on NNs.

C PROCEDURE OF RIO:

This section provides an algorithmic description of RIO (see Algorithm 1).

D EMPIRICAL STUDY D.1 EXPERIMENTAL SETUPS:

Dataset Description In total, 12 real-world regression datasets from UCI machine learning repository (Dua & Graff, 2017) are tested. Table 4 summarizes the basic information of these datasets. For all the datasets except MSD, 20% of the whole dataset is used as test dataset and 80% is used as training dataset, and this split is randomly generated in each independent run. For MSD, the first 463715 samples are used as training dataset and the last 51630 samples are used as testing dataset according to the provider's guideline. During the experiments, all the datasets except for MSD are tested for 100 independent runs, and MSD datasets are tested for 10 independent runs. For each independent run, the dataset is randomly split into training set, validation set, and test set (except for MSD, in which the dataset split is strictly predefined by the provider), and the same random dataset split are used by all the tested algorithms to ensure fair comparisons. All source codes for reproducing the experimental results are provided at: (https://github.com/leaf-ai/rio-paper).

Parametric Setup for Algorithms:

• NN: For SC dataset, a fully connected feed-forward NN with 2 hidden layers, each with 128 hidden neurons, is used. For CT dataset, a fully connected feed-forward NN with 2 hidden layers, each with 256 hidden neurons, is used. For MSD dataset, a fully connected • RIO, RIO variants and SVGP (Hensman et al., 2013) : SVGP is used as an approximator to original GP in RIO and all the RIO variants. For RIO, RIO variants and SVGP, the number of inducing points are 50 for all the experiments. RBF kernel is used for both input and output kernel. For RIO, RIO variants and SVGP, the signal variances and length scales of all the kernels plus the noise variance are the trainable hyperparameters.

The optimizer is L-BFGS-B with default parameters as in Scipy.optimize documentation (https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html), and the maximum number of iterations is set as 1000. The training process runs until the L-BFGS-B optimizer decides to stop.

• NNGP (Lee et al., 2018) : For NNGP kernel, the depth is 2, and the activation function is ReLU. n g = 101, n v = 151, and n c = 131. Following the learning process in original paper, a grid search is performed to search for the best values of σ 2 w and σ 2 b . Same as in the original paper, a grid of 30 points evenly spaced from 0.1 to 5.0 (for σ 2 w ) and 30 points evenly spaced from 0 to 2.0 (for σ 2 b ) was evaluated. The noise variance σ 2 is fixed as 0.01. The grid search process stops when Cholesky decomposition fails or all the 900 points are evaluated. The best values found during the grid search will be used in the experiments. No pre-computed lookup tables are utilized.

• ANP (Kim et al., 2019) : The parametric setups of ANP are following the recommendations in the original paper. The attention type is multihead, the hidden size is 64, the max number of context points is 50, the context ratio is 0.8, the random kernel hyperparameters option is on. The size of latent encoder is 64 × 64 × 64 × 64, the number of latents is 64, the size of deterministic encoder is 64 × 64 × 64 × 64, the size of decoder is 64 × 64 × 64 × 64 × 2, and the deterministic path option is on. Adam optimizer with learning rate 10 -4 is used, and the maximum number of training iterations is 2000.

Performance Metrics:

• To measure the point-prediction error, the Root Mean Square Error (RMSE) between the method predictions and true outcomes on test datasets are calculated for each independent experimental run. After that, the mean and standard deviations of these RMSEs are used to measure the performance of the algorithms.

• To quantitatively measure the quality of uncertainty estimation, average negative log predictive density (NLPD) (Quiñonero-Candela et al., 2006) is used to measure the quality of uncertainty estimation. NLPD is given by

L = - 1 n n i=1 log p(ŷ i = y i |x i ) (10)

where ŷi indicates the prediction results, x i is the input with true associated outcome y i , p(•) is the probability density function (PDF) of the returned distribution based on input x i . • To investigate the behaviors of RIO variants during learning, the mean of estimated noise variance σ 2 n over all the independent runs are calculated. • To compare the computation time of the algorithms, the training time (wall clock time) of NN, RIO, all the ablated RIO variants, SVGP and ANP are averaged over all the independent runs as the computation time. It is notable that the computation time of all RIO variants does not include the training time of associated NN, because the NN is considered to be pretrained. For NNGP, the wall clock time for the grid search is used. In case that the grid search stops due to Cholesky decomposition failures, the computation time of NNGP will be estimated as the average running time of all the successful evaluations × 900, which is the supposed number of evaluations. All the algorithms are implemented using Tensorflow, and tested in the exactly same python environment. All the experiments are running on a machine with 16 Intel(R) Xeon(R) CPU E5-2623 v4@2.60GHz and 128GB memory.

D.2 RESULTS ON ESTIMATED CONFIDENCE INTERVALS:

Confidence interval (CI) is a useful tool to estimate the distribution of outcomes with explicit probabilities in real-world applications. In order to provide more insights for practitioners, the percentages of testing outcomes that are within the 95%/90%/68% CIs as estimated by each algorithm are calculated. Ideally, these percentages should be as close to the estimated confidence levels as possible, e.g., a perfect uncertainty estimator would have exactly 95% of testing outcomes within its estimated 95% CIs. In practice, when two models have similar quality of CI estimations, the more conservative one is favoured. Figure 4 and 5 show the distribution of the percentages that testing outcomes are within the estimated 95%/90%/68% CIs over all the independent runs for all the datasets and algorithms. Figure 6 , 7 and 8 compare the distributions of coverage percentages for 1%-99% CIs for RIO and SVGP. The experimental data are extracted from the experiments that generate Table 1 in main paper. In Figure 4 , 5, 6, 7 and 8, the box extends from the lower to upper quartile values of the data (each data point represents an independent experimental run), with a line at the median. The whiskers extend from the box to show the range of the data. Flier points are those past the end of the whiskers, indicating the outliers.

According to the results, RIO is able to provide reasonable CI estimations in most cases. However, using this metric to compare algorithms regarding uncertainty estimation quality may be noisy and misleading. A method that simply learns the distribution of the labels would perform well in CI coverage metric, but it cannot make any meaningful point-wise prediction. More detailed examples are given as follows.

For "CT" dataset, SVGP has an extremely high RMSE of around 52 while RIO variants only have RMSEs of 1. However, SVGP still shows acceptable 95% and 90% CI coverage, and even has over-confident coverage for 68% CI. After investigation, it was found that this happened only by chance, and was not due to the accurate CI estimations of SVGP. Since SVGP is not able to extract any useful information from the high-dimensional input space, it simply treated all outcomes as noise. As a result, SVGP shows a very large RMSE compared to other algorithms, and the mean of its predicted outcome distribution is always around 0. Since SVGP treats everything as noise, the estimated noise variance is very high, and the estimated 95% CI based on this noise variance is overly high and covers all the test outcomes in most cases. When the estimated 90% CI is evaluated, the large error in mean estimation and large error in noise variance estimation cancel mostly cancel each other out by chance, i.e., the estimated 90% CI is mistakenly shifted by erroneous mean then the overly wide noise variance fortunately covers slightly more than 90% test outcomes. The case of the estimated 68% CI is similar, but now the error in noise variance cannot fully cover the error in mean, so the coverage percentages are now below 68%, indicating over-confidence.

One interesting observation is that SVGP tends to be more "conservative" for high confidence levels (>90%), even in cases where they are "optimistic" for low confidence levels. After investigation, this is because SVGP normally has an overly high noise variance estimation (also comes with a higher prediction RMSE in most cases), so it has a higher chance of covering more points when the increase in CI width (influenced by noise variance) surpasses the erroneous shift of mean (depending on prediction errors).

D.3 ADDITIONAL RESULTS ON OUTPUT CORRECTION:

This section shows additional demonstrations of the detailed behaviors of RIO output correction. The experimental data are extracted from the experiments that generate Table 1 in main paper.

Figure 9 plots the distributions of ground truth labels (outcomes), original NN predictions and predictions corrected after RIO for a randomly picked run for each dataset. Figure 10 shows the pointwise comparisons between NN outputs and RIO-corrected outputs for the same experimental runs as in Figure 9 . Based on the results, it is clear that RIO tends to calibrate each NN prediction accordingly. The distribution of outputs after RIO calibration may be a shift, or shrinkage, or expansion, or even more complex modifications of the original NN predictions, depending on how different are NN predictions from ground truths. As a result, the distribution of RIO calibrated outputs are closer to the distribution of the ground truth. One interesting behavior can be observed in Figure 9 for "protein" dataset: after applying RIO, the range of whiskers shrunk and the outliers disappeared, but the box (indicating 25 to 75 percentile of the data) expanded. This behavior shows that RIO can customize its calibration to each point. Another interesting behavior is that for "wine/r" dataset (see both Figure 9 and 10), RIO shifts all the original NN outputs to lower values, which are closer to the distribution of the ground truth.

To better study the behaviors of RIO, a new performance metric is defined, called improvement ratio (IR), which is the ratio between number of successful corrections (successfully reducing the prediction error) and total number of data points. For each run on each dataset, this IR value is calculated, and the distribution of IR values over 100 independent runs (random dataset split except for MSD, random NN initialization and training) on each dataset is plotted in Figure 11 . According to the results, the IR values for RIO are above 0.5 in most cases. For 7 datasets, IR values are above 0.5 in all 100 independent runs. For some runs in yacht, ENB/h, CT, and MSD, the IR values are above 0.8 or even above 0.9. All these observations show that RIO is making meaningful corrections instead of random perturbations. Results in Figure 11 also provides useful information for practitioners: Although not all RIO calibrations improve the result, most of them do.

Same empirical analysis is also conducted for two RIO variants, namely R+I (predicting residuals with only input kernel) and R+O (predicting residuals with only output kernel). Figure 12 , 13 and 14 show results for R+I, and Figure 15 , 16 and 17 show results for R+O. From the results, the output kernel is helpful in problems where input kernel does not work well (CT and MSD), and it also shows more robust performance in terms of improvement ratio (IR) in most datasets. However it is still generally worse than full RIO. More specifically, R+I shows an extremely low IR in CT dataset (see Figure 14 ). After investigation, this is because the input kernel itself is not able to learn anything from the complex high-dimensional input space, so it treats everything as noise. As a result, it keeps the NN output unchanged during correction in most cases. Applying output kernel instead solves the issue. After comparing Figure 10 , 13 and 16, it can be observed that the behaviors of RIO are either a mixture or selection between R+I and R+O. This means RIO with I/O kernel is able to choose the best kernel among these two or combines both if needed.

D.4 EXPERIMENTAL RESULTS WITH ARD:

The experimental results shown in the main paper is based on the setup that all RIO variants are using standard RBF kernel without Automatic Relevance Determination (ARD). To investigate the performance of RIO under more sophisticated setups, same experiments are run for all RIO variants with ARD mechanism turned on (all other experimental setups are identical with section D.1, except that the NN depth for MSD dataset is reduced from 4 hidden layers to 2, due to computation failure issue during Cholesky decomposition). Table 5 shows the summarized experimental results. From the results, RIO still performs the best or equals the best method (based on statistical tests) in 9 out of 12 datasets for both RMSE and NLPD metrics. This clearly shows the effectiveness and robustness of RIO. The symbols † and ‡ indicate that the difference between the marked entry and RIO is statistically significant at the 5% significance level using paired t-test and Wilcoxon test, respectively. The best entries that are significantly better than all the others under at least one statistical test are marked in boldface (ties are allowed). The experimental results shown in the main paper is based on the setup that all RIO variants are using L-BFGS-B optimizer with a maximun number of iterations as 1,000. To investigate the performance of RIO under larger computational budget, same experiments are run for all RIO variants with maximum number of optimizer iterations as 10,000 (all other experimental setups are identical with section D.1). Table 6 shows the summarized experimental results for 8 smallest datasets. According to the results, the rankings of the methods are very similar to those in Table 1 (in the main paper), and RIO still performs best in all metrics.

D.6 RESULTS ON RANDOM FORESTS:

In principle, RIO can be applied to any prediction model since it treats the pretrained model as a black box. To demonstrate this extensibility, RIO is applied to another classical prediction model -Random Forests (RF) (Ho, 1995; Breiman, 2001) . The same experimental setups as in section D.1 are used. For RF, the number of estimators is set as 100 for all datasets. To avoid overfitting, the minimum number of samples required to be at a leaf node is 10 for all datasets, and the max depth of the tree is 7 for MSD and 5 for all other datasets. Table 7 shows the experimental results. From Table 7 , RIO performs the best or equals the best method (based on statistical tests) in 9 out of 12 datasets in terms of both RMSE and NLPD. In addition, RIO significantly improves the performance of original RF in 11 out of 12 datasets. These results demonstrates the robustness and broad applicability of RIO. The symbols † and ‡ indicate that the difference between the marked entry and RIO is statistically significant at the 5% significance level using paired t-test and Wilcoxon test, respectively. The best entries that are significantly better than all the others under at least one statistical test are marked in boldface (ties are allowed). The symbols † and ‡ indicate that the difference between the marked entry and RIO is statistically significant at the 5% significance level using paired t-test and Wilcoxon test, respectively. The best entries that are significantly better than all the others under at least one statistical test are marked in boldface (ties are allowed).

:

Note that a related information theoretic notion of indistinguishability has recently been defined for analyzing the limits of deep learning(Abbe & Sandon,

2018).

Figures and Tables:

Figure1: Complete model-building process. Given a dataset, first a standard NN model is constructed and trained by a data scientist. The RIO method takes this pretrained model and trains a GP to estimate the residuals of the NN using both the output of the NN and the original input. Blue pathways are only active during the training phase. In the deployment phase, the GP provides uncertainty estimates for the predictions, while calibrating them, i.e., making point predictions more accurate. Overall, RIO transforms the standard NN regression model into a more practical probabilistic estimator.

Figure2: Capturing uncertainty of more and less accurate NNs. These figures illustrate the behavior of RIO in two cases: (left) The neural network has discovered true complex structure in the labels, so the residuals have low variance and are easy for the GP to fit with high confidence; (right) The ineffective neural network has introduced unnecessary complexity, so the residuals are modeled with high uncertainty. In both cases, RIO matches the intuition for how uncertain the NN really is.

NEURAL NETWORKS Neural Networks (NNs) learn a nonlinear transformation from input to output space based on a number of training examples. Let D ⊆ R din × R dout denote the training dataset with size n, and X = {x

Procedure of RIO Require: (X , y) = {(x i , y i )} n i=1 : training data ŷ = {ŷ i } n i=1 : NN predictions on training data x * : data to be predicted ŷ * : NN prediction on x * Ensure: ŷ * ∼ N (ŷ * + r * , var(r * )): a distribution of calibrated prediction Training Phase: 1: calculate residuals r = {r i = y i -ŷi } n i=1 2: for each optimizer step do 3:

Figure5: Quality of estimated CIs. These figures show the distribution of the percentages that testing outcomes are within the estimated 95%/90%/68% CIs over all the independent runs.

Figure 6: Distribution of CI Coverage. These figures compare the distributions of empirical coverage percentages of estimated 1%-99% CIs over all the independent runs for RIO and SVGP.

Figure7: Distribution of CI Coverage. These figures compare the distributions of empirical coverage percentages of estimated 1%-99% CIs over all the independent runs for RIO and SVGP.

Figure8: Distribution of CI Coverage. These figures compare the distributions of empirical coverage percentages of estimated 1%-99% CIs over all the independent runs for RIO and SVGP.

of ground truths and predictions

Figure 9: Distribution of Ground Truths and NN/RIO-corrected predictions. The box extends from the lower to upper quartile values of the data points, with a line at the median. The whiskers extend from the box to show the range of the data. Flier points are those past the end of the whiskers, indicating the outliers.

Figure 10: Comparison between NN and RIO-corrected outputs. Each dot represents a data point. The horizontal axis denote the original NN predictions, and the vertical axis the corresponding predictions after RIO corrections. The solid line indicates where NN predictions are same as RIOcorrected predictions (i.e., no change in output). This figure shows that RIO exhibits diverse behavior in how it calibrates output.

yacht ENB/h ENB/c airfoil CCS wine/r wine/w CCPP protein SC CT

Figure 11: Distribution of Impovement Ratio. The box extends from the lower to upper quartile values of the data (each data point represents an independent experimental run), with a line at the median. The whiskers extend from the box to show the range of the data. Flier points are those past the end of the whiskers, indicating the outliers.

Figure 12: Distribution of Ground Truths and NN/R+I-corrected predictions. The box extends from the lower to upper quartile values of the data points, with a line at the median. The whiskers extend from the box to show the of the data. Flier points are those past the end of the whiskers, indicating the outliers.

Figure 13: Comparison between NN and R+I-corrected outputs. Each dot represents a data point. The horizontal axis denote the original NN predictions, and the vertical axis the corresponding predictions after R+I corrections. The solid line indicates where NN predictions are same as R+Icorrected predictions (i.e., no change in output).

of ground truths and predictions

Figure 15: Distribution of Ground Truths and NN/R+O-corrected predictions. The box extends from the lower to upper quartile values of the data points, with a line at the median. The whiskers extend from the box to show the range of the data. Flier points are those past the end of the whiskers, indicating the outliers.

Figure 16: Comparison between NN and R+O-corrected outputs. Each dot represents a data point. The horizontal axis denote the original NN predictions, and the vertical axis the corresponding predictions after R+O corrections. The solid line indicates where NN predictions are same as R+O-corrected predictions (i.e., no change in output).

Summary of experimental results

Table content: <table><tr><td colspan="2">Dataset Method n × d</td><td>RMSE mean±std</td><td>NLPD mean±std</td><td colspan="4">Noise Variance (sec) Time Dataset Method n × d</td><td>RMSE mean±std</td><td>NLPD mean±std</td><td colspan="2">Noise Variance (sec) Time</td></tr><tr><td>yacht</td><td>NN RIO</td><td>2.30±0.93 † ‡ 1.46±0.49</td><td>-2.039±0.762</td><td>-0.82</td><td>4.02 7.16</td><td>ENB/h</td><td>NN RIO</td><td>1.03±0.51 † ‡ 0.70±0.38</td><td>-1.038±0.355</td><td>-0.46</td><td>6.65 8.18</td></tr><tr><td/><td>R+I</td><td>2.03±0.73 † ‡</td><td>2.341±0.516 † ‡</td><td>2.54</td><td>4.30</td><td/><td>R+I</td><td>0.79±0.46 † ‡</td><td>1.147±0.405 † ‡</td><td>0.63</td><td>7.52</td></tr><tr><td>308</td><td>R+O</td><td>1.88±0.66 † ‡</td><td>2.305±0.614 † ‡</td><td>1.60</td><td>6.27</td><td>768</td><td>R+O</td><td>0.80±0.43 † ‡</td><td>1.169±0.388 † ‡</td><td>0.59</td><td>7.61</td></tr><tr><td>×</td><td>Y+O</td><td>1.86±0.64 † ‡</td><td>2.305±0.639 † ‡</td><td>1.89</td><td>9.93</td><td>×</td><td>Y+O</td><td>0.88±0.48 † ‡</td><td>1.248±0.405 † ‡</td><td>0.75</td><td>11.06</td></tr><tr><td>6</td><td>Y+IO</td><td>1.58±0.52 † ‡</td><td>2.160±0.773 † ‡</td><td>1.18</td><td>9.44</td><td>8</td><td>Y+IO</td><td>0.76±0.41 † ‡</td><td>1.124±0.368 † ‡</td><td>0.56</td><td>10.64</td></tr><tr><td/><td>SVGP</td><td>4.42±0.62 † ‡</td><td>2.888±0.102 † ‡</td><td>18.56</td><td>8.96</td><td/><td>SVGP</td><td>2.13±0.18 † ‡</td><td>2.200±0.074 † ‡</td><td>4.70</td><td>10.16</td></tr><tr><td/><td>NNGP</td><td colspan="2">12.40±1.45 † ‡ 35.18±0.534 † ‡</td><td>-</td><td>7347</td><td/><td>NNGP</td><td>4.97±0.29 † ‡</td><td>32.40±0.638 † ‡</td><td>-</td><td>7374</td></tr><tr><td/><td>ANP</td><td>7.59±3.20 † ‡</td><td>1.793±0.887 † ‡</td><td>-</td><td>40.82</td><td/><td>ANP</td><td>4.08±2.27 † ‡</td><td>2.475±0.559 † ‡</td><td>-</td><td>102.3</td></tr><tr><td>ENB/c</td><td>NN RIO</td><td>1.88±0.44 † ‡ 1.48±0.33</td><td>-1.816±0.191</td><td>-1.58</td><td>6.45 8.07</td><td>airfoil</td><td>NN RIO</td><td>4.82±0.43 † ‡ 3.07±0.18</td><td>-2.554±0.053</td><td>-9.48</td><td>6.48 17.63</td></tr><tr><td/><td>R+I</td><td>1.71±0.44 † ‡</td><td>1.969±0.236 † ‡</td><td>2.22</td><td>5.02</td><td/><td>R+I</td><td>3.16±0.18 † ‡</td><td>2.583±0.051 † ‡</td><td>10.07</td><td>15.90</td></tr><tr><td>768</td><td>R+O</td><td>1.75±0.43 † ‡</td><td>2.000±0.229 † ‡</td><td>2.25</td><td>4.57</td><td>1505</td><td>R+O</td><td>4.17±0.26 † ‡</td><td>2.849±0.066 † ‡</td><td>16.64</td><td>9.97</td></tr><tr><td>×</td><td>Y+O</td><td>1.76±0.43 † ‡</td><td>2.000±0.231 † ‡</td><td>2.32</td><td>10.99</td><td>×</td><td>Y+O</td><td>4.24±0.28 † ‡</td><td>2.869±0.075 † ‡</td><td>17.81</td><td>22.72</td></tr><tr><td>8</td><td>Y+IO</td><td>1.64±0.36 † ‡</td><td>1.936±0.210 † ‡</td><td>1.96</td><td>10.56</td><td>5</td><td>Y+IO</td><td>3.64±0.53 † ‡</td><td>2.712±0.150 † ‡</td><td>14.40</td><td>24.51</td></tr><tr><td/><td>SVGP</td><td>2.63±0.23 † ‡</td><td>2.403±0.078 † ‡</td><td>6.81</td><td>10.28</td><td/><td>SVGP</td><td>3.59±0.20 † ‡</td><td>2.699±0.053 † ‡</td><td>12.67</td><td>21.74</td></tr><tr><td/><td>NNGP</td><td>4.91±0.32 † ‡</td><td>30.14±0.886 † ‡</td><td>-</td><td>7704</td><td/><td>NNGP</td><td>6.54±0.23 † ‡</td><td>33.60±0.420 † ‡</td><td>-</td><td>3355</td></tr><tr><td/><td>ANP</td><td>4.81±2.15 † ‡</td><td>2.698±0.548 † ‡</td><td>-</td><td>64.11</td><td/><td>ANP</td><td colspan="2">21.17±30.72 † ‡ 5.399±6.316 † ‡</td><td>-</td><td>231.7</td></tr><tr><td>CCS</td><td>NN RIO</td><td>6.23±0.53 † ‡ 5.97±0.48</td><td>-3.241±0.109</td><td>-24.74</td><td>9.46 13.71</td><td>wine/r</td><td>NN RIO</td><td>0.691±0.041 † ‡ 0.672±0.036</td><td>-1.094±0.100</td><td>-0.28</td><td>3.61 9.25</td></tr><tr><td>1030</td><td>R+I</td><td>6.01±0.50 † ‡</td><td>3.248±0.112 † ‡</td><td>25.40</td><td>9.52</td><td>1599</td><td>R+I</td><td colspan="2">0.669±0.036 † ‡ 1.085±0.097 † ‡</td><td>0.28</td><td>8.34</td></tr><tr><td>×</td><td>R+O</td><td>6.17±0.54 † ‡</td><td>3.283±0.120 † ‡</td><td>26.31</td><td>9.54</td><td>×</td><td>R+O</td><td colspan="2">0.676±0.035 † ‡ 1.099±0.094 ‡</td><td>0.29</td><td>5.02</td></tr><tr><td>8</td><td>Y+O</td><td>6.15±0.52 † ‡</td><td>3.279±0.117 † ‡</td><td>26.53</td><td>21.35</td><td>11</td><td>Y+O</td><td>0.676±0.034 † ‡</td><td>1.096±0.092</td><td>0.29</td><td>12.71</td></tr><tr><td/><td>Y+IO</td><td>6.06±0.49 † ‡</td><td>3.261±0.110 † ‡</td><td>25.82</td><td>23.15</td><td/><td>Y+IO</td><td>0.672±0.036 † ‡</td><td>1.094±0.098</td><td>0.28</td><td>12.48</td></tr><tr><td/><td>SVGP</td><td>6.87±0.39 † ‡</td><td>3.336±0.048 † ‡</td><td>44.55</td><td>19.85</td><td/><td colspan="3">SVGP 0.642±0.028 † ‡ 0.974±0.042 † ‡</td><td>0.40</td><td>12.17</td></tr><tr><td>wine/w</td><td>NN RIO</td><td>0.721±0.023 † ‡ 0.704±0.018</td><td>-1.090±0.038</td><td>-0.37</td><td>7.17 16.74</td><td>CCPP</td><td>NN RIO</td><td>4.96±0.53 † ‡ 4.05±0.128</td><td>-2.818±0.031</td><td>-16.30</td><td>14.52 42.65</td></tr><tr><td>4898</td><td>R+I</td><td colspan="2">0.699±0.018 † ‡ 1.081±0.037 † ‡</td><td>0.38</td><td>13.5</td><td>9568</td><td>R+I</td><td>4.06±0.13 † ‡</td><td>2.822±0.031 † ‡</td><td>16.39</td><td>39.88</td></tr><tr><td>×</td><td>R+O</td><td colspan="2">0.710±0.019 † ‡ 1.098±0.038 † ‡</td><td>0.39</td><td>6.19</td><td>×</td><td>R+O</td><td>4.32±0.15 † ‡</td><td>2.883±0.035 † ‡</td><td>18.50</td><td>18.48</td></tr><tr><td>11</td><td>Y+O</td><td colspan="2">0.710±0.019 † ‡ 1.096±0.038 † ‡</td><td>0.39</td><td>18.39</td><td>4</td><td>Y+O</td><td>4.37±0.20 † ‡</td><td>2.914±0.122 † ‡</td><td>23.98</td><td>48.27</td></tr><tr><td/><td>Y+IO</td><td>0.705±0.019 † ‡</td><td>1.090±0.038</td><td>0.38</td><td>20.06</td><td/><td>Y+IO</td><td>4.56±1.00 † ‡</td><td>2.958±0.216 † ‡</td><td>31.06</td><td>46.8</td></tr><tr><td/><td colspan="3">SVGP 0.719±0.018 † ‡ 1.081±0.022 † ‡</td><td>0.50</td><td>18.18</td><td/><td>SVGP</td><td>4.36±0.13 † ‡</td><td>2.893±0.031 † ‡</td><td>19.04</td><td>46.43</td></tr><tr><td>protein</td><td>NN RIO</td><td>4.21±0.07 † ‡ 4.08±0.06</td><td>-2.826±0.014</td><td>-15.71</td><td>151.8 149.4</td><td>SC</td><td>NN RIO</td><td>12.23±0.77 † ‡ 11.28±0.46</td><td>-3.853±0.042</td><td>-105.83</td><td>51.9 53.39</td></tr><tr><td>45730</td><td>R+I</td><td>4.11±0.06 † ‡</td><td>2.834±0.037 † ‡</td><td>15.99</td><td colspan="2">141.2 21263</td><td>R+I</td><td colspan="2">11.33±0.45 † ‡ 3.858±0.041 † ‡</td><td>107.35</td><td>47.72</td></tr><tr><td>×</td><td>R+O</td><td>4.14±0.06 † ‡</td><td>2.840±0.015 † ‡</td><td>16.18</td><td>115.1</td><td>×</td><td>R+O</td><td colspan="2">11.63±0.52 † ‡ 3.881±0.046 † ‡</td><td>112.91</td><td>30.47</td></tr><tr><td>9</td><td>Y+O</td><td>4.14±0.06 † ‡</td><td>2.840±0.015 † ‡</td><td>16.17</td><td>138.4</td><td>81</td><td>Y+O</td><td colspan="2">11.64±0.53 † ‡ 3.882±0.046 † ‡</td><td>113.61</td><td>45.35</td></tr><tr><td/><td>Y+IO</td><td>4.08±0.06</td><td>2.826±0.014</td><td>15.72</td><td>155.5</td><td/><td>Y+IO</td><td colspan="2">11.32±0.45 † ‡ 3.856±0.041 † ‡</td><td>106.93</td><td>57.74</td></tr><tr><td/><td>SVGP</td><td>4.68±0.04 † ‡</td><td>2.963±0.007 † ‡</td><td>22.54</td><td>149.5</td><td/><td>SVGP</td><td colspan="2">14.66±0.25 † ‡ 4.136±0.014 † ‡</td><td>239.28</td><td>50.89</td></tr><tr><td>CT</td><td>NN RIO</td><td>1.17±0.34 † ‡ 0.88±0.15</td><td>-1.284±0.219</td><td>-1.02</td><td>194.5 516.4</td><td>MSD</td><td>NN RIO</td><td>12.42±2.97 † ‡ 9.26±0.21</td><td>-3.639±0.022</td><td>-84.28</td><td>1136 1993</td></tr><tr><td>53500</td><td>R+I</td><td>1.17±0.34 † ‡</td><td>1.538±0.289 † ‡</td><td>1.71</td><td colspan="2">19.80 515345</td><td>R+I</td><td colspan="2">10.92±1.30 † ‡ 3.811±0.128 † ‡</td><td>135.34</td><td>282.0</td></tr><tr><td>×</td><td>R+O</td><td>0.88±0.15</td><td>1.283±0.219 † ‡</td><td>1.02</td><td>159.4</td><td>×</td><td>R+O</td><td>9.25±0.20</td><td>3.638±0.021</td><td>84.05</td><td>1518</td></tr><tr><td>384</td><td>Y+O</td><td>0.99±0.42 † ‡</td><td>1.365±0.385 † ‡</td><td>2.45</td><td>168.2</td><td>90</td><td>Y+O</td><td colspan="2">10.00±0.86 † ‡ 3.768±0.148 † ‡</td><td>169.90</td><td>1080</td></tr><tr><td/><td>Y+IO</td><td>0.91±0.16 † ‡</td><td>1.280±0.177 ‡</td><td>0.62</td><td>578.6</td><td/><td>Y+IO</td><td>9.43±0.52 ‡</td><td>3.644±0.025 † ‡</td><td>85.66</td><td>2605</td></tr><tr><td/><td>SVGP</td><td colspan="2">52.07±0.19 † ‡ 5.372±0.004 † ‡</td><td>2712</td><td>27.56</td><td/><td>SVGP</td><td>9.57±0.00 † ‡</td><td>3.677±0.000 † ‡</td><td>92.21</td><td>2276</td></tr></table>

3: Comparison among NN, RIO, and SVGP. The horizontal axis denotes the prediction RMSE of the NN, and the vertical axis the prediction RMSE of RIO (blue dots) and SVGP (yellow dots). Each dot represents an independent experimental run. Since the scales are different, the solid blue line indicates where NN and RIO/SVGP have same prediction RMSE. Thus, a dot below the line means that the method (RIO or SVGP) performs better than the NN, and vice versa. Results of SVGP on the CT dataset are not plotted because its prediction RMSE exceeded the visible scale (i.e. they were >50). RIO consistently reduces the error of the NN, and outperforms SVGP in most cases.and SVGP in terms of prediction RMSE. RIO is able to improve the NN predictions consistently regardless of how the dataset is split and how well the NN is trained. Even in the situations where original NN is much worse than SVGP, RIO still successfully improves the NN performance into a level that is comparable or better than SVGP. Appendix D.3 provides additional investigation into the behavior of RIO output corrections. RIO exhibits diverse behaviors that generally move the original NN predictions closer to the ground truth. To conclude, applying RIO to NNs not only provides additional uncertainty information, but also reliably reduces the point-prediction error.

Table content: <table/>

Spearman's correlation between RMSE and σ 2 n across RIO variants (including SVGP)

Table content: <table><tr><td/><td>yacht</td><td>ENB/h</td><td>ENB/c</td><td>airfoil</td><td>CCS</td><td>wine/r</td><td>wine/w</td><td>CCPP</td><td>protein</td><td>SC</td><td>CT</td><td>MSD</td></tr><tr><td>correlation</td><td>0.943</td><td>0.943</td><td>1.0</td><td>1.0</td><td>0.943</td><td>-0.09</td><td>0.886</td><td>1.0</td><td>0.943</td><td>1.0</td><td>0.771</td><td>0.943</td></tr><tr><td>p-value</td><td>0.005</td><td>0.005</td><td>0.0</td><td>0.0</td><td>0.005</td><td>0.872</td><td>0.02</td><td>0.0</td><td>0.005</td><td>0.0</td><td>0.072</td><td>0.005</td></tr><tr><td/><td colspan="11">Entries that are considered to indicate very strong positive monotonic correlation are marked in boldface.</td><td/></tr></table>

IMDB Age Estimation Results with DenseNet

Table content: <table><tr><td>Method</td><td>MAE</td><td>NLPD</td><td>95% CI Coverage</td><td>90% CI Coverage</td><td>68% CI Coverage</td></tr><tr><td>Pretrained DenseNet (Yang et al., 2018)</td><td>7.43</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SVGP</td><td>36.45</td><td>5.06</td><td>0.99</td><td>0.96</td><td>0.62</td></tr><tr><td>RIO</td><td>6.35</td><td>3.59</td><td>0.94</td><td>0.91</td><td>0.75</td></tr></table>

Tsun-Yi Yang, Yi-Hsuan Huang, Yen-Yu Lin, Pi-Cheng Hsiu, and Yung-Yu Chuang. Ssr-net: A compact soft stagewise regression network for age estimation. In Proc. of IJCAI, pp. 1078-1084, 2018.

Table content: <table><tr><td>Ciyou Zhu, Richard H. Byrd, Peihuang Lu, and Jorge Nocedal. Algorithm 778: L-bfgs-b: Fortran</td></tr><tr><td>subroutines for large-scale bound-constrained optimization. ACM Trans. Math. Softw., 23(4):</td></tr><tr><td>550-560, December 1997. ISSN 0098-3500. doi: 10.1145/279232.279236. URL http://doi.</td></tr><tr><td>acm.org/10.1145/279232.279236.</td></tr></table>

Summary of testing dataset The inputs to the NN are normalized to have mean 0 and standard deviation 1. The activation function is ReLU for all the hidden layers. The maximum number of epochs for training is 1000. 20% of the training data is used as validation data, and the split is random at each independent run. An early stop is triggered if the loss on validation data has not be improved for 10 epochs. The optimizer is RMSprop with learning rate 0.001, and the loss function is mean squared error (MSE).

Table content: <table><tr><td>abbreviation</td><td>full name in UCI ML repository</td><td colspan="2">dataset size dimension</td><td>note</td></tr><tr><td>yacht</td><td>Yacht Hydrodynamics Data Set</td><td>308</td><td>6</td><td>-</td></tr><tr><td>ENB/h</td><td>Energy efficiency</td><td>768</td><td>8</td><td>Heating Load as target</td></tr><tr><td>ENB/c</td><td>Energy efficiency</td><td>768</td><td>8</td><td>Cooling Load as target</td></tr><tr><td>airfoil</td><td>Airfoil Self-Noise</td><td>1505</td><td>5</td><td>-</td></tr><tr><td>CCS</td><td>Concrete Compressive Strength</td><td>1030</td><td>8</td><td>-</td></tr><tr><td>wine/r</td><td>Wine Quality</td><td>1599</td><td>11</td><td>only use winequality-red data</td></tr><tr><td>wine/w</td><td>Wine Quality</td><td>4898</td><td>11</td><td>only use winequality-white data</td></tr><tr><td>CCPP</td><td>Combined Cycle Power Plant</td><td>9568</td><td>4</td><td>-</td></tr><tr><td>CASP</td><td>Physicochemical Properties of Protein Tertiary Structure</td><td>54730</td><td>9</td><td>-</td></tr><tr><td>SC</td><td>Superconductivty Data</td><td>21263</td><td>81</td><td>-</td></tr><tr><td>CT</td><td>Relative location of CT slices on axial axis</td><td>53500</td><td>384</td><td>-</td></tr><tr><td>MSD</td><td>YearPredictionMSD</td><td>515345</td><td>90</td><td>train: first 463715, test: last 51630</td></tr><tr><td colspan="5">feed-forward NN with 4 hidden layers, each with 64 hidden neurons, is used. For all the</td></tr><tr><td colspan="5">remaining datasets, a fully connected feed-forward NN with 2 hidden layers, each with 64</td></tr><tr><td colspan="2">hidden neurons, is used.</td><td/><td/><td/></tr></table>

Summary of experimental results (with ARD)

Table content: <table><tr><td colspan="2">Dataset Method n × d</td><td>RMSE mean±std</td><td>NLPD mean±std</td><td colspan="4">Noise Variance (sec) Time Dataset Method n × d</td><td>RMSE mean±std</td><td>NLPD mean±std</td><td colspan="2">Noise Variance (sec) Time</td></tr><tr><td>yacht</td><td>NN RIO</td><td>2.35±1.16 † ‡ 1.22±0.44</td><td>-1.718±0.531</td><td>-0.46</td><td>3.41 7.73</td><td>ENB/h</td><td>NN RIO</td><td>1.01±0.38 † ‡ 0.59±0.10</td><td>-0.937±0.186</td><td>-0.22</td><td>8.47 11.98</td></tr><tr><td/><td>R+I</td><td>1.36±0.39 † ‡</td><td>1.893±0.482 † ‡</td><td>0.83</td><td>6.89</td><td/><td>R+I</td><td>0.61±0.12 † ‡</td><td>0.971±0.197 † ‡</td><td>0.25</td><td>11.07</td></tr><tr><td>308</td><td>R+O</td><td>1.88±0.65 † ‡</td><td>2.208±0.470 † ‡</td><td>1.62</td><td>6.01</td><td>768</td><td>R+O</td><td>0.77±0.32 † ‡</td><td>1.150±0.328 † ‡</td><td>0.49</td><td>10.51</td></tr><tr><td>×</td><td>Y+O</td><td>1.88±0.62 † ‡</td><td>2.227±0.439 † ‡</td><td>1.82</td><td>9.56</td><td>×</td><td>Y+O</td><td>0.83±0.37 † ‡</td><td>1.226±0.345 † ‡</td><td>0.61</td><td>15.31</td></tr><tr><td>6</td><td>Y+IO</td><td>1.09±0.30 † ‡</td><td>1.591±0.425 † ‡</td><td>0.84</td><td>9.85</td><td>8</td><td>Y+IO</td><td>0.61±0.08 † ‡</td><td>0.974±0.164 † ‡</td><td>0.25</td><td>15.96</td></tr><tr><td/><td>SVGP</td><td>0.91±0.21 † ‡</td><td>1.359±0.208 † ‡</td><td>0.77</td><td>8.93</td><td/><td>SVGP</td><td>0.77±0.05 † ‡</td><td>1.156±0.049 † ‡</td><td>0.59</td><td>14.21</td></tr><tr><td/><td>NNGP</td><td colspan="2">12.40±1.45 † ‡ 35.18±0.534 † ‡</td><td>-</td><td>7347</td><td/><td>NNGP</td><td>4.97±0.29 † ‡</td><td>32.40±0.638 † ‡</td><td>-</td><td>7374</td></tr><tr><td/><td>ANP</td><td>7.59±3.20 † ‡</td><td>1.793±0.887 † ‡</td><td>-</td><td>40.82</td><td/><td>ANP</td><td>4.08±2.27 † ‡</td><td>2.475±0.559 † ‡</td><td>-</td><td>102.3</td></tr><tr><td>ENB/c</td><td>NN RIO</td><td>1.88±0.46 † ‡ 1.42±0.18</td><td>-1.825±0.155</td><td>-1.31</td><td colspan="2">11.64 airfoil 15.27</td><td>NN RIO</td><td>4.84±0.74 † ‡ 2.49±0.15</td><td>-2.349±0.048</td><td>-6.56</td><td>8.83 17.92</td></tr><tr><td/><td>R+I</td><td>1.45±0.17 † ‡</td><td>1.837±0.134 ‡</td><td>1.45</td><td>14.02</td><td/><td>R+I</td><td>2.56±0.16 † ‡</td><td>2.376±0.054 † ‡</td><td>6.92</td><td>16.64</td></tr><tr><td>768</td><td>R+O</td><td>1.77±0.46 † ‡</td><td>2.021±0.232 † ‡</td><td>2.34</td><td>9.73</td><td>1505</td><td>R+O</td><td>4.14±0.27 † ‡</td><td>2.844±0.071 † ‡</td><td>16.31</td><td>10.97</td></tr><tr><td>×</td><td>Y+O</td><td>1.77±0.45 † ‡</td><td>2.020±0.231 † ‡</td><td>2.56</td><td>20.54</td><td>×</td><td>Y+O</td><td>4.39±1.64 † ‡</td><td>2.889±0.196 † ‡</td><td>20.08</td><td>22.90</td></tr><tr><td>8</td><td>Y+IO</td><td>1.46±0.22 † ‡</td><td>1.847±0.157 † ‡</td><td>1.42</td><td>21.65</td><td>5</td><td>Y+IO</td><td>3.73±0.35 † ‡</td><td>2.746±0.127 † ‡</td><td>14.88</td><td>25.08</td></tr><tr><td/><td>SVGP</td><td>1.73±0.25 † ‡</td><td>1.969±0.130 † ‡</td><td>3.12</td><td>18.73</td><td/><td>SVGP</td><td>3.36±0.23 † ‡</td><td>2.632±0.065 † ‡</td><td>11.11</td><td>21.80</td></tr><tr><td/><td>NNGP</td><td>4.91±0.32 † ‡</td><td>30.14±0.886 † ‡</td><td>-</td><td>7704</td><td/><td>NNGP</td><td>6.54±0.23 † ‡</td><td>33.60±0.420 † ‡</td><td>-</td><td>3355</td></tr><tr><td/><td>ANP</td><td>4.81±2.15 † ‡</td><td>2.698±0.548 † ‡</td><td>-</td><td>64.11</td><td/><td>ANP</td><td colspan="2">21.17±30.72 † ‡ 5.399±6.316 † ‡</td><td>-</td><td>231.7</td></tr><tr><td>CCS</td><td>NN RIO</td><td>6.29±0.54 † ‡ 5.81±0.50</td><td>-3.210±0.119</td><td>-23.00</td><td>6.53 9.20</td><td>wine/r</td><td>NN RIO</td><td>0.689±0.037 † ‡ 0.674±0.033</td><td>-1.091±0.082</td><td>-0.28</td><td>3.24 7.55</td></tr><tr><td>1030</td><td>R+I</td><td>5.80±0.49</td><td>3.206±0.118</td><td>23.19</td><td>8.80</td><td>1599</td><td>R+I</td><td colspan="2">0.670±0.033 † ‡ 1.085±0.084 † ‡</td><td>0.28</td><td>9.52</td></tr><tr><td>×</td><td>R+O</td><td>6.21±0.53 † ‡</td><td>3.286±0.121 † ‡</td><td>27.02</td><td>3.77</td><td>×</td><td>R+O</td><td>0.676±0.034 † ‡</td><td>1.093±0.083</td><td>0.29</td><td>5.18</td></tr><tr><td>8</td><td>Y+O</td><td>6.18±0.51 † ‡</td><td>3.278±0.114 † ‡</td><td>27.22</td><td>11.74</td><td>11</td><td>Y+O</td><td>0.675±0.033</td><td>1.088±0.081 ‡</td><td>0.29</td><td>13.56</td></tr><tr><td/><td>Y+IO</td><td>5.91±0.46 † ‡</td><td>3.228±0.108 † ‡</td><td>24.50</td><td>12.09</td><td/><td>Y+IO</td><td colspan="2">0.672±0.033 † ‡ 1.085±0.081 † ‡</td><td>0.29</td><td>14.23</td></tr><tr><td/><td>SVGP</td><td>6.20±0.40 † ‡</td><td>3.233±0.060 † ‡</td><td>34.66</td><td>11.34</td><td/><td colspan="3">SVGP 0.640±0.028 † ‡ 0.969±0.043 † ‡</td><td>0.39</td><td>13.32</td></tr><tr><td>wine/w</td><td>NN RIO</td><td>0.725±0.026 † ‡ 0.707±0.017</td><td>-1.093±0.035</td><td>-0.38</td><td>7.01 11.8</td><td>CCPP</td><td>NN RIO</td><td>4.97±0.53 † ‡ 3.99±0.13</td><td>-2.796±0.026</td><td>-15.82</td><td>10.11 26.27</td></tr><tr><td>4898</td><td>R+I</td><td colspan="2">0.702±0.018 † ‡ 1.084±0.035 † ‡</td><td>0.38</td><td>14.76</td><td>9568</td><td>R+I</td><td>3.99±0.13</td><td>2.797±0.025 † ‡</td><td>15.87</td><td>24.47</td></tr><tr><td>×</td><td>R+O</td><td colspan="2">0.711±0.019 † ‡ 1.098±0.037 † ‡</td><td>0.39</td><td>6.46</td><td>×</td><td>R+O</td><td>4.33±0.13 † ‡</td><td>2.879±0.027 † ‡</td><td>18.51</td><td>9.99</td></tr><tr><td>11</td><td>Y+O</td><td colspan="2">0.710±0.019 † ‡ 1.096±0.037 † ‡</td><td>0.39</td><td>19.77</td><td>4</td><td>Y+O</td><td colspan="2">8.94±44.78 † ‡ 2.974±0.484 † ‡</td><td>2095</td><td>27.24</td></tr><tr><td/><td>Y+IO</td><td>0.708±0.018 †</td><td>1.093±0.036</td><td>0.39</td><td>20.15</td><td/><td>Y+IO</td><td>4.57±0.97 † ‡</td><td>2.968±0.255 † ‡</td><td>36.48</td><td>27.65</td></tr><tr><td/><td colspan="3">SVGP 0.714±0.017 † ‡ 1.074±0.022 † ‡</td><td>0.50</td><td>19.42</td><td/><td>SVGP</td><td>8.80±44.83 ‡</td><td>2.917±0.468 † ‡</td><td>1455</td><td>26.61</td></tr><tr><td>protein</td><td>NN RIO</td><td>4.23±0.08 † ‡ 4.08±0.05</td><td>-2.826±0.013</td><td>-15.75</td><td>147.4 149.6</td><td>SC</td><td>NN RIO</td><td>12.41±0.84 † ‡ 11.24±0.33</td><td>-3.844±0.030</td><td>-104.6</td><td>73.11 99.91</td></tr><tr><td>45730</td><td>R+I</td><td>4.11±0.05 † ‡</td><td>2.834±0.013 † ‡</td><td>16.01</td><td colspan="2">130.9 21263</td><td>R+I</td><td colspan="2">11.27±0.32 † ‡ 3.847±0.029 † ‡</td><td>105.5</td><td>95.33</td></tr><tr><td>×</td><td>R+O</td><td>4.15±0.06 † ‡</td><td>2.843±0.015 † ‡</td><td>16.31</td><td>98.88</td><td>×</td><td>R+O</td><td colspan="2">11.68±0.42 † ‡ 3.881±0.037 † ‡</td><td>113.9</td><td>47.61</td></tr><tr><td>9</td><td>Y+O</td><td>4.15±0.06 † ‡</td><td>2.843±0.015 † ‡</td><td>16.31</td><td>120.4</td><td>81</td><td>Y+O</td><td colspan="2">11.68±0.42 † ‡ 3.882±0.036 † ‡</td><td>114.2</td><td>72.08</td></tr><tr><td/><td>Y+IO</td><td>4.08±0.05</td><td>2.826±0.013</td><td>15.75</td><td>148.6</td><td/><td>Y+IO</td><td colspan="2">11.29±0.38 † ‡ 3.848±0.034 † ‡</td><td>105.9</td><td>103.4</td></tr><tr><td/><td>SVGP</td><td>4.64±0.04 † ‡</td><td>2.955±0.008 † ‡</td><td>22.15</td><td>128.1</td><td/><td>SVGP</td><td colspan="2">14.12±0.27 † ‡ 4.090±0.015 † ‡</td><td>217.8</td><td>95.78</td></tr><tr><td>CT</td><td>NN RIO</td><td>1.12±0.29 † ‡ 0.86±0.12</td><td>-1.261±0.202</td><td>-0.95</td><td>196.6 519.6</td><td>MSD</td><td>NN RIO</td><td>12.54±1.16 † ‡ 9.79±0.22</td><td>-3.698±0.023</td><td>-94.90</td><td>777.8 3059</td></tr><tr><td>53500</td><td>R+I</td><td>1.12±0.29 † ‡</td><td>1.505±0.260 † ‡</td><td>1.55</td><td colspan="2">20.03 515345</td><td>R+I</td><td colspan="2">12.54±1.16 † ‡ 3.943±0.092 † ‡</td><td>156.6</td><td>91.01</td></tr><tr><td>×</td><td>R+O</td><td>0.86±0.12</td><td>1.261±0.202</td><td>0.95</td><td>159.4</td><td>×</td><td>R+O</td><td>9.82±0.19</td><td>3.701±0.020</td><td>96.11</td><td>2739</td></tr><tr><td>384</td><td>Y+O</td><td>0.98±0.77 ‡</td><td>1.303±0.336</td><td>2.02</td><td>180.6</td><td>90</td><td>Y+O</td><td>209.8±596.2 ‡</td><td>7.010±9.680 ‡</td><td>236.4</td><td>1491</td></tr><tr><td/><td>Y+IO</td><td>0.88±0.13 † ‡</td><td>1.256±0.145</td><td>0.55</td><td>594.5</td><td/><td>Y+IO</td><td>208.8±596.6 ‡</td><td>7.055±9.664 ‡</td><td>197.2</td><td>3347</td></tr><tr><td/><td>SVGP</td><td colspan="2">52.07±0.19 † ‡ 5.372±0.004 † ‡</td><td>2712</td><td>29.69</td><td/><td>SVGP</td><td>10.97±0.0 † ‡</td><td>3.981±0.0 † ‡</td><td>199.7</td><td>2590</td></tr></table>

Summary of experimental results (with 10,000 maximum optimizer iterations)

Table content: <table><tr><td colspan="2">Dataset Method n × d</td><td>RMSE mean±std</td><td>NLPD mean±std</td><td colspan="4">Noise Variance (sec) Time Dataset Method n × d</td><td>RMSE mean±std</td><td>NLPD mean±std</td><td>Noise Variance</td><td>Time (sec)</td></tr><tr><td>yacht</td><td>NN RIO</td><td>2.20±0.93 † ‡ 1.40±0.50</td><td>-1.883±0.568</td><td>-0.74</td><td>3.33 25.67</td><td>ENB/h</td><td>NN RIO</td><td>0.94±0.37 † ‡ 0.64±0.26</td><td>-0.968±0.273</td><td>-0.31</td><td>10.14 45.04</td></tr><tr><td/><td>R+I</td><td>1.93±0.65 † ‡</td><td>2.266±0.484 † ‡</td><td>2.19</td><td>5.59</td><td/><td>R+I</td><td>0.70±0.33 † ‡</td><td>1.043±0.317 † ‡</td><td>0.41</td><td>22.22</td></tr><tr><td>308</td><td>R+O</td><td>1.78±0.57 † ‡</td><td>2.176±0.525 † ‡</td><td>1.39</td><td>6.76</td><td>768</td><td>R+O</td><td>0.72±0.31 † ‡</td><td>1.084±0.309 † ‡</td><td>0.41</td><td>20.49</td></tr><tr><td>×</td><td>Y+O</td><td>1.78±0.56 † ‡</td><td>2.204±0.509 † ‡</td><td>1.60</td><td>23.99</td><td>×</td><td>Y+O</td><td>0.78±0.35 † ‡</td><td>1.163±0.328 † ‡</td><td>0.55</td><td>55.23</td></tr><tr><td>6</td><td>Y+IO</td><td>1.40±0.44</td><td>1.919±0.567</td><td>0.82</td><td>45.69</td><td>8</td><td>Y+IO</td><td>0.66±0.26 † ‡</td><td>1.013±0.280 † ‡</td><td>0.34</td><td>82.82</td></tr><tr><td/><td>SVGP</td><td>3.67±0.60 † ‡</td><td>2.689±0.111 † ‡</td><td>12.07</td><td>42.59</td><td/><td>SVGP</td><td>2.01±0.17 † ‡</td><td>2.145±0.071 † ‡</td><td>4.24</td><td>79.46</td></tr><tr><td/><td>NNGP</td><td colspan="2">12.40±1.45 † ‡ 35.18±0.534 † ‡</td><td>-</td><td>7347</td><td/><td>NNGP</td><td>4.97±0.29 † ‡</td><td>32.40±0.638 † ‡</td><td>-</td><td>7374</td></tr><tr><td/><td>ANP</td><td>7.59±3.20 † ‡</td><td>1.793±0.887 † ‡</td><td>-</td><td>40.82</td><td/><td>ANP</td><td>4.08±2.27 † ‡</td><td>2.475±0.559 † ‡</td><td>-</td><td>102.3</td></tr><tr><td>ENB/c</td><td>NN RIO</td><td>1.87±0.42 † ‡ 1.51±0.35</td><td>-1.852±0.198</td><td>-1.59</td><td colspan="2">11.79 airfoil 48.53</td><td>NN RIO</td><td>4.84±0.47 † ‡ 3.06±0.20</td><td>-2.551±0.058</td><td>-9.44</td><td>8.96 104.0</td></tr><tr><td/><td>R+I</td><td>1.70±0.41 † ‡</td><td>1.98±0.21 † ‡</td><td>2.17</td><td>10.96</td><td/><td>R+I</td><td>3.13±0.21 † ‡</td><td>2.573±0.059 † ‡</td><td>9.91</td><td>73.22</td></tr><tr><td>768</td><td>R+O</td><td>1.75±0.41 † ‡</td><td>2.011±0.211 † ‡</td><td>2.21</td><td>10.85</td><td>1505</td><td>R+O</td><td>4.16±0.27 † ‡</td><td>2.848±0.068 † ‡</td><td>16.58</td><td>10.92</td></tr><tr><td>×</td><td>Y+O</td><td>1.75±0.41 † ‡</td><td>2.012±0.210 † ‡</td><td>2.27</td><td>39.52</td><td>×</td><td>Y+O</td><td>4.21±0.27 † ‡</td><td>2.862±0.082 † ‡</td><td>17.89</td><td>38.69</td></tr><tr><td>8</td><td>Y+IO</td><td>1.62±0.35 † ‡</td><td>1.936±0.197 † ‡</td><td>1.86</td><td>70.94</td><td>5</td><td>Y+IO</td><td>3.19±0.30 † ‡</td><td>2.583±0.087 † ‡</td><td>10.24</td><td>119.95</td></tr><tr><td/><td>SVGP</td><td>2.52±0.21 † ‡</td><td>2.363±0.072 † ‡</td><td>6.31</td><td>84.32</td><td/><td>SVGP</td><td>3.27±0.20 † ‡</td><td>2.608±0.056 † ‡</td><td>10.56</td><td>106.0</td></tr><tr><td/><td>NNGP</td><td>4.91±0.32 † ‡</td><td>30.14±0.886 † ‡</td><td>-</td><td>7704</td><td/><td>NNGP</td><td>6.54±0.23 † ‡</td><td>33.60±0.420 † ‡</td><td>-</td><td>3355</td></tr><tr><td/><td>ANP</td><td>4.81±2.15 † ‡</td><td>2.698±0.548 † ‡</td><td>-</td><td>64.11</td><td/><td>ANP</td><td colspan="2">21.17±30.72 † ‡ 5.399±6.316 † ‡</td><td>-</td><td>231.7</td></tr><tr><td>CCS</td><td>NN RIO</td><td>6.25±0.49 † ‡ 5.96±0.47</td><td>-3.230±0.108</td><td>-25.37</td><td>6.54 14.31</td><td>wine/r</td><td>NN RIO</td><td>0.688±0.039 † ‡ 0.671±0.033</td><td>-1.088±0.087</td><td>-0.28</td><td>3.26 20.12</td></tr><tr><td>1030</td><td>R+I</td><td>5.99±0.47 † ‡</td><td>3.235±0.107</td><td>26.04</td><td>4.91</td><td>1599</td><td>R+I</td><td colspan="2">0.668±0.033 † ‡ 1.080±0.085 † ‡</td><td>0.28</td><td>12.61</td></tr><tr><td>×</td><td>R+O</td><td>6.19±0.49 † ‡</td><td>3.280±0.112 † ‡</td><td>27.43</td><td>4.04</td><td>×</td><td>R+O</td><td colspan="2">0.675±0.033 † ‡ 1.094±0.088 † ‡</td><td>0.29</td><td>4.96</td></tr><tr><td>8</td><td>Y+O</td><td>6.18±0.48 † ‡</td><td>3.276±0.109 † ‡</td><td>27.65</td><td>17.55</td><td>11</td><td>Y+O</td><td colspan="2">0.674±0.033 † ‡ 1.089±0.086 † ‡</td><td>0.29</td><td>17.01</td></tr><tr><td/><td>Y+IO</td><td>6.03±0.47 † ‡</td><td>3.246±0.107 † ‡</td><td>26.24</td><td>41.89</td><td/><td>Y+IO</td><td>0.671±0.032</td><td>1.087±0.086</td><td>0.28</td><td>34.45</td></tr><tr><td/><td>SVGP</td><td>6.62±0.37 † ‡</td><td>3.297±0.045 † ‡</td><td>41.15</td><td>73.66</td><td/><td colspan="3">SVGP 0.642±0.028 † ‡ 0.973±0.042 † ‡</td><td>0.39</td><td>70.53</td></tr><tr><td>wine/w</td><td>NN RIO</td><td>0.723±0.027 † ‡ 0.704±0.018</td><td>-1.088±0.034</td><td>-0.38</td><td>8.51 49.96</td><td>CCPP</td><td>NN RIO</td><td>4.94±0.49 † ‡ 4.03±0.13</td><td>-2.808±0.025</td><td>-16.21</td><td>17.38 151.3</td></tr><tr><td>4898</td><td>R+I</td><td colspan="2">0.700±0.017 † ‡ 1.079±0.033 † ‡</td><td>0.38</td><td>26.43</td><td>9568</td><td>R+I</td><td>4.04±0.13 † ‡</td><td>2.810±0.026 † ‡</td><td>16.28</td><td>116.0</td></tr><tr><td>×</td><td>R+O</td><td colspan="2">0.710±0.021 † ‡ 1.095±0.037 † ‡</td><td>0.39</td><td>8.87</td><td>×</td><td>R+O</td><td>4.33±0.14 † ‡</td><td>2.880±0.029 † ‡</td><td>18.56</td><td>19.02</td></tr><tr><td>11</td><td>Y+O</td><td colspan="2">0.710±0.020 † ‡ 1.093±0.037 † ‡</td><td>0.39</td><td>29.97</td><td>4</td><td>Y+O</td><td colspan="2">13.40±63.01 ‡ 3.012±0.663 † ‡</td><td>4161</td><td>100.8</td></tr><tr><td/><td>Y+IO</td><td>0.704±0.018 ‡</td><td>1.088±0.034</td><td>0.38</td><td>66.79</td><td/><td>Y+IO</td><td>4.71±1.51 † ‡</td><td>2.969±0.271 † ‡</td><td>33.58</td><td>267.1</td></tr><tr><td/><td colspan="3">SVGP 0.713±0.016 † ‡ 1.076±0.022 † ‡</td><td>0.5</td><td>158.1</td><td/><td>SVGP</td><td>4.25±0.13 † ‡</td><td>2.859±0.028 † ‡</td><td>17.94</td><td>334.6</td></tr></table>

Summary of experimental results with Random Forests (RF)

Table content: <table><tr><td colspan="2">Dataset Method n × d</td><td>RMSE mean±std</td><td>NLPD mean±std</td><td colspan="4">Noise Variance (sec) Time Dataset Method n × d</td><td>RMSE mean±std</td><td>NLPD mean±std</td><td colspan="2">Noise Variance (sec) Time</td></tr><tr><td>yacht</td><td>RF RIO</td><td>1.95±0.59 † ‡ 0.80±0.18</td><td>-1.177±0.187</td><td>-0.49</td><td>0.33 14.93</td><td>ENB/h</td><td>RF RIO</td><td>1.44±0.16 † ‡ 0.79±0.13</td><td>-1.150±0.138</td><td>-0.61</td><td>0.35 16.66</td></tr><tr><td/><td>R+I</td><td>1.62±0.67 † ‡</td><td>1.833±0.453 † ‡</td><td>1.91</td><td>12.41</td><td/><td>R+I</td><td>0.91±0.12 † ‡</td><td>1.306±0.128 † ‡</td><td>0.74</td><td>14.13</td></tr><tr><td>252</td><td>R+O</td><td>1.29±0.33 † ‡</td><td>1.658±0.267 † ‡</td><td>1.31</td><td>12.00</td><td>768</td><td>R+O</td><td>1.02±0.20 † ‡</td><td>1.401±0.187 † ‡</td><td>0.92</td><td>13.81</td></tr><tr><td>×</td><td>Y+O</td><td>1.47±0.33 † ‡</td><td>1.841±0.243 † ‡</td><td>2.29</td><td>16.43</td><td>×</td><td>Y+O</td><td>1.43±0.17 † ‡</td><td>1.792±0.146 † ‡</td><td>1.88</td><td>17.57</td></tr><tr><td>6</td><td>Y+IO</td><td>0.87±0.20 † ‡</td><td>1.268±0.193 † ‡</td><td>0.65</td><td>18.40</td><td>8</td><td>Y+IO</td><td>0.88±0.11 † ‡</td><td>1.280±0.111 † ‡</td><td>0.73</td><td>20.51</td></tr><tr><td/><td>SVGP</td><td>4.41±0.60 † ‡</td><td>2.886±0.096 † ‡</td><td>18.55</td><td>15.44</td><td/><td>SVGP</td><td>2.13±0.18 † ‡</td><td>2.199±0.074 † ‡</td><td>4.70</td><td>16.90</td></tr><tr><td/><td>NNGP</td><td colspan="2">12.40±1.45 † ‡ 35.18±0.534 † ‡</td><td>-</td><td>7347</td><td/><td>NNGP</td><td>4.97±0.29 † ‡</td><td>32.40±0.638 † ‡</td><td>-</td><td>7374</td></tr><tr><td/><td>ANP</td><td>7.59±3.20 † ‡</td><td>1.793±0.887 † ‡</td><td>-</td><td>40.82</td><td/><td>ANP</td><td>4.08±2.27 † ‡</td><td>2.475±0.559 † ‡</td><td>-</td><td>102.3</td></tr><tr><td>ENB/c</td><td>RF RIO</td><td>1.93±0.15 † ‡ 1.55±0.15</td><td>-1.854±0.087</td><td>-2.13</td><td>0.31 16.31</td><td>airfoil</td><td>RF RIO</td><td>3.68±0.17 † ‡ 2.90±0.17</td><td>-2.487±0.063</td><td>-7.37</td><td>0.36 19.80</td></tr><tr><td/><td>R+I</td><td>1.53±0.12</td><td>1.845±0.072</td><td>2.08</td><td>12.95</td><td/><td>R+I</td><td>2.95±0.17 † ‡</td><td>2.505±0.059 † ‡</td><td>7.88</td><td>17.19</td></tr><tr><td>768</td><td>R+O</td><td>1.85±0.14 † ‡</td><td>2.037±0.089 † ‡</td><td>2.93</td><td>10.62</td><td>1505</td><td>R+O</td><td>3.54±0.17 † ‡</td><td>2.693±0.059 † ‡</td><td>10.50</td><td>8.26</td></tr><tr><td>×</td><td>Y+O</td><td>1.92±0.15 † ‡</td><td>2.081±0.091 † ‡</td><td>3.20</td><td>19.82</td><td>×</td><td>Y+O</td><td>3.58±0.18 † ‡</td><td>2.714±0.085 † ‡</td><td>12.05</td><td>21.22</td></tr><tr><td>8</td><td>Y+IO</td><td>1.85±0.14 † ‡</td><td>2.039±0.087 † ‡</td><td>2.99</td><td>21.50</td><td>5</td><td>Y+IO</td><td>3.31±0.41 † ‡</td><td>2.623±0.131 † ‡</td><td>10.44</td><td>23.61</td></tr><tr><td/><td>SVGP</td><td>2.63±0.23 † ‡</td><td>2.403±0.079 † ‡</td><td>6.81</td><td>18.63</td><td/><td>SVGP</td><td>3.59±0.20 † ‡</td><td>2.698±0.054 † ‡</td><td>12.66</td><td>20.52</td></tr><tr><td/><td>NNGP</td><td>4.91±0.32 † ‡</td><td>30.14±0.886 † ‡</td><td>-</td><td>7704</td><td/><td>NNGP</td><td>6.54±0.23 † ‡</td><td>33.60±0.420 † ‡</td><td>-</td><td>3355</td></tr><tr><td/><td>ANP</td><td>4.81±2.15 † ‡</td><td>2.698±0.548 † ‡</td><td>-</td><td>64.11</td><td/><td>ANP</td><td colspan="2">21.17±30.72 † ‡ 5.399±6.316 † ‡</td><td>-</td><td>231.7</td></tr><tr><td>CCS</td><td>RF RIO</td><td>7.17±0.43 † ‡ 5.84±0.34</td><td>-3.180±0.068</td><td>-24.84</td><td>0.49 17.72</td><td>wine/r</td><td>RF RIO</td><td>0.625±0.028 † ‡ 0.642±0.028</td><td>-1.019±0.070</td><td>-0.27</td><td>0.75 21.13</td></tr><tr><td>1030</td><td>R+I</td><td>5.89±0.35 † ‡</td><td>3.190±0.069 † ‡</td><td>25.34</td><td>14.62</td><td>1599</td><td>R+I</td><td colspan="2">0.625±0.028 † ‡ 0.966±0.059 † ‡</td><td>0.31</td><td>7.57</td></tr><tr><td>×</td><td>R+O</td><td>6.90±0.41 † ‡</td><td>3.380±0.079 † ‡</td><td>34.20</td><td>9.77</td><td>×</td><td>R+O</td><td colspan="2">0.627±0.028 † ‡ 0.974±0.061 † ‡</td><td>0.30</td><td>9.26</td></tr><tr><td>8</td><td>Y+O</td><td>7.03±0.78 † ‡</td><td>3.401±0.154 † ‡</td><td>35.75</td><td>18.72</td><td>11</td><td>Y+O</td><td colspan="2">0.628±0.028 † ‡ 0.975±0.062 † ‡</td><td>0.30</td><td>20.74</td></tr><tr><td/><td>Y+IO</td><td>5.88±0.33 † ‡</td><td>3.192±0.068 † ‡</td><td>25.78</td><td>22.01</td><td/><td>Y+IO</td><td colspan="2">0.635±0.031 † ‡ 0.999±0.074 † ‡</td><td>0.28</td><td>24.88</td></tr><tr><td/><td>SVGP</td><td>6.88±0.40 † ‡</td><td>3.337±0.048 † ‡</td><td>44.62</td><td>18.09</td><td/><td>SVGP</td><td>0.642±0.028</td><td>0.974±0.042 † ‡</td><td>0.40</td><td>20.41</td></tr><tr><td>wine/w</td><td>RF RIO</td><td>0.714±0.016 † ‡ 0.702±0.015</td><td>-1.068±0.025</td><td>-0.43</td><td>1.17 36.60</td><td>CCPP</td><td>RF RIO</td><td>4.21±0.13 † ‡ 4.02±0.14</td><td>-2.801±0.029</td><td>-15.18</td><td>1.47 44.81</td></tr><tr><td>4898</td><td>R+I</td><td colspan="2">0.697±0.016 † ‡ 1.059±0.025 † ‡</td><td>0.44</td><td>33.17</td><td>9568</td><td>R+I</td><td>4.03±0.14 † ‡</td><td>2.803±0.029 † ‡</td><td>15.22</td><td>39.09</td></tr><tr><td>×</td><td>R+O</td><td colspan="2">0.715±0.016 † ‡ 1.087±0.026 † ‡</td><td>0.46</td><td>13.70</td><td>×</td><td>R+O</td><td>4.17±0.13 † ‡</td><td>2.839±0.028 † ‡</td><td>16.19</td><td>17.67</td></tr><tr><td>11</td><td>Y+O</td><td colspan="2">0.715±0.016 † ‡ 1.087±0.026 † ‡</td><td>0.46</td><td>34.16</td><td>4</td><td>Y+O</td><td colspan="2">13.26±63.05 ‡ 3.274±3.329 † ‡</td><td>2095</td><td>44.63</td></tr><tr><td/><td>Y+IO</td><td colspan="2">0.710±0.017 † ‡ 1.080±0.028 † ‡</td><td>0.45</td><td>41.52</td><td/><td>Y+IO</td><td>4.63±0.95 † ‡</td><td>2.975±0.253 † ‡</td><td>35.54</td><td>45.63</td></tr><tr><td/><td colspan="3">SVGP 0.719±0.018 † ‡ 1.081±0.022 † ‡</td><td>0.50</td><td>36.62</td><td/><td>SVGP</td><td>4.36±0.13 † ‡</td><td>2.887±0.026 † ‡</td><td>18.94</td><td>48.83</td></tr><tr><td>protein</td><td>RF RIO</td><td>5.05±0.03 † ‡ 4.55±0.03</td><td>-2.935±0.007</td><td>-20.92</td><td>19.76 298.6</td><td>SC</td><td>RF RIO</td><td>15.26±0.25 † ‡ 14.18±0.75</td><td>-4.072±0.051</td><td>-196.1</td><td>56.98 170.6</td></tr><tr><td>45730</td><td>R+I</td><td>4.57±0.04 † ‡</td><td>2.939±0.008 † ‡</td><td>21.06</td><td colspan="2">278.1 21263</td><td>R+I</td><td>14.38±0.76 ‡</td><td>4.089±0.047 † ‡</td><td>204.1</td><td>145.0</td></tr><tr><td>×</td><td>R+O</td><td>5.02±0.04 † ‡</td><td>3.033±0.007 † ‡</td><td>24.79</td><td>242.9</td><td>×</td><td>R+O</td><td colspan="2">14.98±0.40 † ‡ 4.126±0.034 † ‡</td><td>213.9</td><td>72.33</td></tr><tr><td>9</td><td>Y+O</td><td>5.04±0.04 † ‡</td><td>3.037±0.007 † ‡</td><td>25.02</td><td>206.4</td><td>81</td><td>Y+O</td><td colspan="2">14.80±0.30 † ‡ 4.121±0.059 † ‡</td><td>204.8</td><td>167.6</td></tr><tr><td/><td>Y+IO</td><td>4.56±0.03 † ‡</td><td>2.937±0.007 † ‡</td><td>21.00</td><td>300.2</td><td/><td>Y+IO</td><td>14.39±0.29</td><td>4.079±0.020</td><td>196.9</td><td>234.4</td></tr><tr><td/><td>SVGP</td><td>4.68±0.04 † ‡</td><td>2.963±0.008 † ‡</td><td>22.54</td><td>281.1</td><td/><td>SVGP</td><td colspan="2">14.66±0.25 † ‡ 4.135±0.013 † ‡</td><td>239.2</td><td>221.3</td></tr><tr><td>CT</td><td>RF RIO</td><td>9.77±0.16 † ‡ 9.32±0.49</td><td>-3.651±0.053</td><td>-86.67</td><td>302.0 1024</td><td>MSD</td><td>RF RIO</td><td>9.93±0.01 † ‡ 9.90±0.01</td><td>-3.711±0.001</td><td>-97.59</td><td>2112 2133</td></tr><tr><td>53500</td><td>R+I</td><td>9.77±0.16 † ‡</td><td>3.698±0.017 † ‡</td><td>94.04</td><td colspan="2">47.04 515345</td><td>R+I</td><td>9.93±0.01 † ‡</td><td>3.714±0.001 † ‡</td><td>98.11</td><td>136.8</td></tr><tr><td>×</td><td>R+O</td><td>9.35±0.46</td><td>3.654±0.049</td><td>86.94</td><td>228.1</td><td>×</td><td>R+O</td><td>9.92±0.01 † ‡</td><td>3.714±0.001 † ‡</td><td>97.78</td><td>1747</td></tr><tr><td>384</td><td>Y+O</td><td>9.46±0.36 † ‡</td><td>3.668±0.042 † ‡</td><td>88.00</td><td>448.7</td><td>90</td><td>Y+O</td><td>9.94±0.01 † ‡</td><td>3.715±0.001 † ‡</td><td>98.29</td><td>2198</td></tr><tr><td/><td>Y+IO</td><td>9.58±0.66 † ‡</td><td>3.676±0.094 † ‡</td><td>91.09</td><td>1439</td><td/><td>Y+IO</td><td>9.94±0.01 † ‡</td><td>3.715±0.001 † ‡</td><td>96.73</td><td>5638</td></tr><tr><td/><td>SVGP</td><td colspan="2">52.07±0.17 † ‡ 5.372±0.003 † ‡</td><td>2712</td><td>61.70</td><td/><td>SVGP</td><td>9.57±0.00 † ‡</td><td>3.677±0.000 † ‡</td><td>92.21</td><td>2276</td></tr></table>

Peer review (iclr2020):
# Summary The authors propose a method for post-hoc correction and predictive variance estimation for neural network models. The method fits a GP to the model residuals, and learns a composite kernel that combines two kernels defined on the input space and the model’s output space (called RIO, R for residual, and IO for the input-output kernel). The authors suggest that residual means and variances at test points can then be calculated explicitly using predictive distributions from the GP. The authors run a large panel of experiments across a number of datasets, and compare to a number of methods that draw connections between neural networks and GP’s. In addition, the full method is compared to a number of methods that utilize only some components of the full RIO method. In these experiments, the RIO method generally shows strong performance in both RMSE and NLPD compared to these baselines. # Feedback Overall, this is a neat method. It has the flavor of a number of other composite ML methods that have worked well in the past---e.g., boosting and platt scaling---but is different enough to stand on its own. The experimental results are quite promising. However, I am torn about the paper, because the theoretical discussion of the method is quite convoluted and seems either irrelevant or incorrect. I wish that the authors had spent more time with small demonstrations of what the procedure does in some simple settings. This would give practitioners considering the method far more intuition about when they would expect it to work and fail than the current theoretical discussion. ## Uncertainty Discussion is Lacking The motivation and discussion sell this method as an uncertainty quantification method, but almost all of the theoretical development revolves around prediction correction. The methods properties as an uncertainty quantification tool are underdeveloped. The only theoretical point made about uncertainty estimation is Theorem 2.7, which states that the scalar variance of the GP “nugget” is positively correlated with the variance of the NN’s residuals. Providing a scalar summary of noise is not particularly compelling for a method advertised as a point-prediction uncertainty quantification method. In addition, it is not clear what probability distribution the “correlation” is defined over. The argument made in the proof seems quite obvious: if a GP is used to model a noisier process (i.e., residuals with a larger variance), it will in some cases classify that variability as independent noise. If the authors wanted to focus on the properties of their method as an uncertainty quantification tool, they could discuss the assumptions underlying the GP error estimates, and when they would be likely to diverge from practical properties like predictive interval coverage. For example, because the base NN predictor is treated as fixed, it seems that this method ignore uncertainty that stems from the NN fit due to random initialization. Likewise, it seems that this method would not quantify uncertainty from resampling the data and obtaining a new NN predictor. )The coverage experiments in the appendix seem to confirm this -- generally, the predictive intervals generally under-cover the predicted values.) It’s fine if the method doesn’t quantify these types of uncertainty, but discussion of these types of issues would be far more welcome than the current convoluted theory in Section 2. This discussion might not yield theorems, but it would give practitioners useful guidelines for deciding whether the particular scheme would likely work for this application. ## Problems with the Error Correction Argument The theory section, especially 2.2, was very difficult to parse. First, as a matter of style, a sequence of Lemma and Theorem statements are given without defining most of the notation used therein, and with almost no prose providing context or intuition. In the buildup to the theorems, it is also unclear which assertions about the decompositions of y_i are assumptions about the true data generating process, and which assertions are specifications of a particular GP model. The substance also has some issues. I think the intention in this section is to get to a rather simple variance decomposition of the labels y. The question is how much variation in y or the residual is represented in the posterior predictive mean of a particular GP. It seems reasonable that in some cases, the structure in the residual may be more amenable to modeling with a stationary GP than the structure in the raw labels y. It is not clear that all of the theoretical complexity here is necessary to make this point. Instead, the authors make a convoluted argument that attempts to establish that the errors from the NN + GP approach will be smaller under very general circumstances. The argument is phrased somewhat ambiguously (it is not clear exactly what is being assumed, and what is corresponds to the specification of a working model), but depending on how one reads this section, the argument makes statements that are either too broad to be correct, or too narrow to be relevant. The argument decomposes for the raw labels and the residuals into pieces that a GP can “capture” or “represent”, and parts that it cannot. The two equations are: y_i = f(x_i) + g(x_i) + xi_i R_i = (y_i - h_NN(x_i)) = r_f(x_i) + r_g(x_i) + xi_i f(.) and r_f(.) represent the portions of the label and residual processes, respectively, that the GP *captures*. It is assumed that the GP will model this portion correctly, and leave the “epsilon-indistinguishable” portion g(.) or r_g(.) untouched. The argument then assumes that f(.) and r_f(.) will have proportional kernels, and so it is possible to show that the predictions of residuals based on r_f(.) will have smaller predictive variance than predictions based on f(.) as long as the variation represented by r_g(.) is smaller than the variation represented by g(.). On its face, this argument raises some red flags. Because h_NN(.) is allowed to be an arbitrary function, the argument here should be symmetric. Why can’t we also get a guaranteed variance reduction by adding h_NN(.) to y rather than subtracting it? Perhaps some of this is captured in the parameter delta, which quantifies the reduction in variation represented in r_g(.) vs g(.), but the argument that the kernel of r_f(.) can be no larger than the kernel f(.) in terms of trace (that is, the proportionality constant alpha is not greater than 1) does not make sense. If h_NN(x_i) is simply -f(x_i), then these arguments would not go through. At the very least, conditions need to be articulated about the properties of h_NN(.). Some of the strangeness comes from the fact that this is a poor model of most prediction problems, where the main issue with fitting a GP is not “indistinguishability”, but misspecification. Consider a process y_i that is non-stationary; say g(.) has a linear trend in some component of x. A GP with a stationary covariance kernel fit to this process (such as RBF) will attempt to explain the variation due to the linear trend with a variance kernel that encodes long-range dependence. On the other hand, if this trend were removed by a base model like an NN, the residuals would have a very different structure (perhaps they would be stationary), and in this case, the GP would fit the data with a very different covariance kernel. Unfortunately, it does not seem like the formalism here can express a notion of misspecification at all. In the theory, it is assumed that the GP will only model the portion of the labels y_i for which it is property specified (in this case, f(.)). This generally does not occur in practice, as in the example above. It might be possible for this to apply in some circumstances, but the authors give no conditions (e.g., that the process y_i be stationary). Based on this assumption, the authors assert that the fitted GP to f(.) and r_f(.) will have the same covariance kernel parameters up to some proportionality constant alpha. Much of their theoretical argument depends on this proportionality. But this proportionality cannot apply in general, and again, no conditions are given for when we might expect this to hold. It would be far more compelling if the authors proposed the very standard approach to modeling data via covariance kernels, where one first models non-stationary portions of the data with a base model, then models the correlation in the residuals with something like a GP. This is the bread-and-butter approach in, say, timeseries analysis (see, e.g., the Shumway and Stoffer textbook https://www.stat.pitt.edu/stoffer/tsa4/tsa4.htm), and the approach in this paper could be framed similarly. ## Demos I Wish I Had Seen I wish the authors had presented some demonstrations of what the GP does to the fitted values of an NN. Giving a demonstration of how the output kernel modifies predicted values, for example, would give some nice intuition the value added by this portion. I suspect that this step essentially performs something like Platt scaling, but for continuous outcomes, by shrinking predictions together so that they better match the overall distribution of observed labels. Perhaps the mechanism is different. At any rate, it would be useful to understand where the information gain is coming from, and this would be far better expressed concretely in terms of a toy data example than the theoretical arguments that are given. ## Coverage Experiments I wish the coverage experiments evaluating predictive intervals were included in the main text. As far as uncertainty quantification evaluations go, coverage is one of the few assessments that does not rely on the model itself (unlike NPLD, which uses the model’s own log-likelihood), and can be phrased as a concrete performance guarantee. Here, the goal for predictive intervals is to cover the true prediction value _at least as often_ as the nominal rate (95% intervals should cover the truth _at least_ 95% of the time), not merely that coverage be “close” to the nominal rate. This asymmetric evaluation gives you a concrete guarantee that the uncertainty estimate is conservative. The coverage experiments show that this method quite systematically under-covers compared to the end-to-end SVGP method, which generally satisfies this coverage property. I think this is important information to include about the model, and generally I think this behavior results from the fact that uncertainty is not propagated from the NN fit. This should be presented clearly in the main text.