Efficient Uncertainty Quantification and Reduction for Over-Parameterized Neural Networks
Authors: Ziyi Huang
Henry Lam
Haofeng Zhang
Venue: NeurIPS 2023 poster (2023)
Keywords: frequentist uncertainty
epistemic uncertainty
procedural variability
confidence intervals
batching
cheap bootstrap

Abstract: Uncertainty quantification (UQ) is important for reliability assessment and enhancement of machine learning models. In deep learning, uncertainties arise not only from data, but also from the training procedure that often injects substantial noises and biases. These hinder the attainment of statistical guarantees and, moreover, impose computational challenges on UQ due to the need for repeated network retraining. Building upon the recent neural tangent kernel theory, we create statistically guaranteed schemes to principally \emph{characterize}, and \emph{remove}, the uncertainty of over-parameterized neural networks with very low computation effort. In particular, our approach, based on what we call a procedural-noise-correcting (PNC) predictor, removes the procedural uncertainty by using only \emph{one} auxiliary network that is trained on a suitably labeled dataset, instead of many retrained networks employed in deep ensembles. Moreover, by combining our PNC predictor with suitable light-computation resampling methods, we build several approaches to construct asymptotically exact-coverage confidence intervals using as low as four trained networks without additional overheads.

Efficient Uncertainty Quantification and Reduction for Over-Parameterized Neural Networks
Ziyi Huang, Henry Lam, Haofeng Zhang
Columbia University
New York. NY, USA
zh2354, kh12114 humb
Abstract
Uncertainty quantification (UQ) is important for reliability assessment and enhance- ment of machine learning models. In deep learning, uncertainties arise not only from data, but also from the training procedure that often injects substantial noises and biases. These hinder the attainment of statistical guarantees and, moreover, impose computational challenges on UQ due to the need for repeated network retraining. Building upon the recent neural tangent kernel theory, we create statisti- cally guaranteed schemes to principally characterize, and remove, the uncertainty of over vreepparamttriized neural networks with very low computation effort. In par- ticular, our approach, based on what we call procedural-noise-coreee correctin (PNC) predictor, removes the procedural uncertainty by using only one auxiliary network that is trained a suitably labeled dataset, instead of many retrained networks employed in deep ensembles. Moreover, by combining our PNC predictor with suitable light- computatinn resampling methods, we build several approaches construct asymptotically exact-coverage confidence intervals using as low as four trained networks without additional overheads.
Introduction
Uncertainty quantification (UQ) concerns the dissection and estimation of various sources of errors in prediction model. It has growing importance in machine learning, as it helps assess and enhance the trustworthiness and deployment safety across many real-world tasks ranging from computer vision 67 21 and natural language processing [13
103 to autonomous driving 99 9 as well as
guiding exploration in sequential learning 9 T19 In the deep learning context, UQ encounters unique challenges on both the statistical and computational fronts. On high level, these challenges arise from the over-parametrizee and large-scale scale nature of neural networks so that, unlike classical statistical models, the prediction outcomes incur not only noises from the data, but also importantly training procedure itself 77 94 This elicits deviation from classical statistical theory that hinders the attainment of established guarantees. Moreover, because of the sizes of these models, conventional procedures such as resampling (37 30 102] demand amount of computation that could quickly become infeasible.
Our main goal of this paper propose UQ framework for over- wver-parametrized neural networks in regression that has simultaneous statistical coverage guarantee, in the sense of classical frequentist asymptotic exactness, and low computation cost, the sense requiring only few (as low as four) neural network trainings, without other extra overheads. A main driver of these strengths in our framework is a new implementable concept, which we call the Procedural-Noiss. -oiiur----ois--Crreettnggg predictor. It consists of an auxiliary network that is trained on a suitably artificially labeled dataset, with behavior mimicking the variability coming from the training procedure. To reach our goal, we
Correcting (PNC)
Authors listed alphahetically.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

synthesize and build on two recent lines of tools that appear largely segregated thus far. First is neural tangent kernel (NTK) theory 64 80 77 which provides explicit approximate formulas for well-trained infinitely wide neural networks. Importantly NTK reveals how procedural variability enters into the network prediction outcomes through, in sense, a functional shift in its corresponding kernel-based regression, which guides our PNC construction. Second is light-computation resampling methodology, including batching [47 99 [00] and the so-called cheap bootstrap method [75), which allows valid confidence interval construction using as few two model repetitions. We suitably enhance these methods to account for both data and procedural variabilities via the PNC incorporation.
We compare our framework with several major related lines of work. First, our work focuses on the quantification of epistemic uncertainty, which refers to the errors coming from the inadequacy of the model or data noises. This is different from aleatoric uncertainty, which refers to the intrinsic
stochasticity of the problem [92 96 121 62 36 or predictive uncertainty which captures the sum of epistemic and aleatoric uncertainties (but not their dissection) [94] 197 12
3 22]. Regarding epistemic uncertainty, a related line of study is deep ensemble that aggregates predictions from we will make clear later, can reduce and potentially quantify procedural variability, but a naive use would require demanding retraining effort and does not address data variability. Another line is the Bayesian UQ approach on neural networks [4 2]. This regards network weights as parameters subject to common priors such as Gaussian. Because of the computation difficulties in exact inference, an array of studies investigate efficient approximate inference approaches to estimate the posteriors [56]. While powerful, these approaches nonetheless possess inference error that could be hard to quantify, and ultimately finding rigorous guarantees on the performance of these approximate posteriors remains open our best knowledge.
multiple independent training replications [8] 73 g 94). This approach,
43 49
16
33
32
95
79
2
Statistical Framework of Uncertainty
We first describe our learning setting and define uncertainties in our framework. Suppose the input- output pair (X, ) is a random vector following an unknown probability distribution π on X' x y, where X E X C Rd is an input and Y E y C R is the response. Let the marginal distribution of X be x (x) and the conditional distribution of Y given X be Tyyx (y]æ). Given set of training data DD {(11, V1), (z2, 92), (zn: drawn independent and identically distributed (i.i.d) from T (we write æ = (Z.....,_,)) and y = (31, (y.....yy..) yx)) for short), we build a prediction model h : X' + } that best approximates Y given X. Let RD or T. denote the empirical distribution associated with the training data D, where nr is used to emphasize the sample size dependence in asymptotic results.
To this end. provided loss function y + R. we denote the population risk R2 (h) := ((x()))(((((()), Y) If h is allowed to be any possible functions, the best prediction model is the Bayes predictor [61 93]: (X) yrgmin ~xxxx [C(v, X). With finite training data of size n, classical statistical learning suggests finding hn E H. where H denotes a hypothesis class that minimizes the empirical risk, i.e., hn E argmin,ex&RRp (h). This framework fits methods such as linear or kernel ridge regression (Appendix (B). However, it is not feasible for deep learning due to non-convexity and non-identifiability Instead, in practice, gradient- gradient-based optimization methods are used, giving rise to fi.,y as a variant of fi,, where the additional variable represents the randomness in the training procedure. It is worth mentioning that this randomness generally depends on the empirical data D, and thus we use P.Ip to represent the distribution of y conditional on D.
Furthermore, we consider hr aggregate((1ngg (fimg ? PyID) where "aggregate" means an idealized aggregation approach to remove the training randomness in hn.r (known as ensemble learning 18. 17 87 45). A prime example in deep learning is deep ensemble 73) that will be detailed in the sequel. Finally, we denote h² +om hr the grand "limiting" predictor with infinite samples. The exact meaning of "lim" will be clear momentarily.
Under this framework, epistemic uncertainty, that is, errors coming from the inadequacy of the model or data noises, can be dissected into the following three sources, as illustrated in Figure Additional discussions other types of uncertainty are in the Appendix A for completeness.
Model approximation error: UQAE 11* - hig- This discrepancy between hi and h** arises from the inadequacy of the hypothesis class H. For an over-parameteriize sufficiently wide neural network
2

---

H. this error is usually negligible thanks to the universal approximation power of neural networks for any continuous functions [28 58 [33 or Lebesgue-inteee integrable functions 88
Data variability, UQnV hü ht" This measures the representativeness of the training dataset, which is the most standard epistemic uncertainty in classical statistics [TT0.
Procedural variability. UQrv = ha., - hr. This arises from the randomness in the training process for a single network hn. which is present even with deterministic infinite data. The randomness comes from the initialization of the network parameters, and also data ordering and possibly training time when running stochastic gradient descent with finite training epochs.
Y = 3o
ht
Approximation Error
Estimation Error
Data Variability
Fj
3 Quantifying Epistemic Uncertainty
har
Procedural Variability
We use frequentist framework and, for given x, we aim to construct a confidence interval for the "best" predictor h" (z).
Epistemic Uncertainty
As discussed in Section D the over-parametrized and non- convex nature of neural networks defies conventional statistical techniques and moreover introduces procedural variability that makes inference difficult.
Figure 1: Three sources of epis- temic uncertainty.
We focus on paeettrramteriiee neural networks, that is, neural networks whose width is sufficiently large, while the depth can be arbitrary such as two [120 107]. ever-parameterized neural networks give the following two theoretical advantages. First, this makes model approximation error negligible and confidence interval for hh(z)) (z) also approximately applies to hig()) Second, the NTK theory [64] implies phenomenon that the network evolves essentially as "linear model" under gradient descent, and thus the resulting predictor behaves like a shifted kernel ridge regression whose kernel is the NTK
80 54 59 118] (detailed Appendix [.. However, to our knowledge, there is no off-the-shelf
result that exactly matches our need for the epistemic uncertainty task, so we describe it below. Consider the following regression problem:
(fo, 8²) = 1
1 (fo(s)) - yi) 2- + Am 110 - 1-112-
i=1
where is the network parameters to be trained, 8 the initial network parameters, and A, > is the regularization hyper-paramede which may depend on the data size. We add regularization An, to this problem, which is slightly different from previous work (80] without the regularization An This can guarantee stable computation of the inversion of the NTK Gram matrix and can be naturally linked to kernel ridge regression, which will be introduced shortly in Proposition 3.1 We assume the network adopts the NTK parametrization with parameters randomly initialized using He initialization [55), and it is sufficiently wide to ensure the linearized neural network assumption holds; See Appendix for details. Moreover, we assume that the network is trained using the loss function in (1) via continuous-tin gradient flow by feeding the entire training data and using sufficient training time oo). In this sense, the uncertainty of data ordering and training time vanishes, making the random initialization the only uncertainty in procedural variability. The above specifications are formally summarized in Assumption C.2
With the above specifications, we have:
C holds. Then the final trained network, conditional on the initial network sgx (I)), is given by
Proposition 3.1 (Proposition C.3 Suppose that Assumption
(2) where in the subscript of h.,.or. represents training data, 8¹ represents instance of the
henee (ar) = sgn (x)(x) + K(x,x)) (K( + XnnII- (gx (x:),
initial network parameters drawn from the standard Gaussian distribution Pa= =N(0, Ip) where the dimension of 6¹ is p (He initialization); K(x, x) == (K(z,, 23)))..=1...... and K (z. x)) := (K(x, (x)x1,) (xx, 12), K(x, ITT)TTTT where K is the (population) NTK. This implies that h is the solution to the following kernel ridge regression:
a
8go + arg min
-
Σ (v
sg - g(x:))² + An Allg
9EH
n1
i=1
3

---

where H. is the reproducing kernel Hilbert space (RKHS) constructed from the NTK ((
Proposition 3.1 shows that the shifted kernel ridge regressor using NTK with a shift from an initial function 8g6 exactly the linearized neural network regressor that starts from the initial network sgs. It also reveals how procedural variability enters into the neural network prediction outcomes, highlighting the need to regard random initialization an inevitable uncertainty component in neural networks. We provide details and deviation of Proposition 3.1 in Appendix C
3.1 Challenges from the Interplay of Procedural and Data Variabilities
First, we will discuss existing challenges in quantifying and reducing uncertainty to motivate main approach based on the PNC predictor under the NTK framework. To this end, deep ensemble 8I 73 41 8 (4] is arguably the most common ensemble approach in deep learning reduce procedural variability. [8] shows that deep ensemble achieves the best performance compared with wide range of other ensemble methods, and employing more networks in deep ensemble can lead better performance. Specifically, the deep ensemble predictor 73] ke (ar) is defined as: (x) == N. Σm= (2)) where m is the number of networks in the ensemble, hn, (z) is the independently trained network with initialization 0e (with the same training data D), and 6f...., 6%. are i.i.d. samples drawn from Pør. We also introduce (x) == Er,, (hn. (x)] as the expectation and using the law of large numbers, we have limm-+so hn (z) = (z) a.s.. So h (x) behaves as an idealized deep ensemble predictor with infinitely many independent training procedures. Using Proposition 3.1 and the linearity of kernel ridge regressor with respect to data (Appendix B), we have
of ho
a0 (I) with respect to o ~ Pgs. Taking
m
(x)
= m Σ
S0x (0)(x) (x) + K (ar, x) (x) + )_nl)) y
Σ i=1
80g (2r)
m
h (x) = Ep o (x) = s(x) + K (a, x) (K (a, x) + )n (y
8 8(x))
(3)
where a(x) = g 1000 (1)) is the expectation of the initial network output 8g (z) with respect to the the distribution of the initialization parameters Pa+ = (0, Ip). It is easy to see that Ep,, [ (x)) = h; (x)) and Var pa. (f.me (x)) m Varp,, (hi,,0 (x:) where Varp,, is the variance with respect to the random initialization. As for the total variance:
Proposition 3.2. We have
Var(h: (z)) Var(fi((=)) m
[f (x))[) S (ar((( (N = Var(h, (x)).
where the variance is taken with respect to both the data 'D and the random initialization Pon-
Proposition 3.2 shows that deep ensemble improves the statistical profile of single model by reducing its procedural variability by factor 1 (but not the data variability), and achieving this reduction requires m training times. To quantify the epistemic uncertainty that contains two variabilities from data and procedure, we may employ resampling approaches such as "bootstrap on deep ensemble". This would involve two layers of sampling, the outer being the resampling of data, and the inner being the retraining of base networks with different initializations. In other words, this nested sampling amounts to a multiplicative amount of training effort in a large number of outer bootstrap resamples and the m inner retaining per resample, leading to a huge computational burden. Moreover, the data variability and procedural variability in neural networks are dependent, making the above approach delicate. More discussion about this issue is presented in Section
In the following, we introduce our PNC framework that can bypass the above issues in that:
Uncertainty reduction. We train one single network and one additional auxiliary network to completely remove the procedural variability. This is in contrast the deep ensemble approach that trains m networks and only reduces the procedural variability by an m-factor.
Uncertainty quantification. We resolve the computational challenges in "bootstrap on deep ensemble", by combining PNC predictors with low-budget inference tools that require only as low as four network trainings.

---

3.2 PNC Predictor and Procedural Variability Removal
We first develop a computationally efficient approach to obtain h; (z), the idealized deep ensemble predictor that free of procedural variability. We term our approach the procedural-noise -oieroooooii--oo (PNC) predictor, whose pseudo-code is given in Algorithm This predictor consists of a difference between essentially two neural network outcomes, h o (x)) which is the original network trained
using one initialization, and 6 (z) that is trained on a suitably artificially labeled dataset. More precisely, this dataset applies label (22) to 24, where 8 is the expected output of an untrained network with random initialization. Also note that labeling this artificial dataset does not involve any training process and, compared with standard network training, the only additional running time in the PNC predictor is to train this single artificial-labes label-trained network.
Algorithm 1 Procedural-Noits Soaieiaa-Noois---reeeeetee Correcting (PNC) Predictor
Input: Training data D = [ y1). (I22 U2),.... (I'. Un))-
Procedure: 1. Draw 0b_ Po. (O, I,) under NTK parameterization. Train standard base
network with data D and the initialization parameters 6h, which outputs ()) (-) in (3).
2. Let (1)) Epg, [sga (z)] For each zy in D, generate its "artificial" label ((x;) = Ep,, [so> (x;)]
Train an auxiliary neural network with data {(I1, 8(x1)) )), (a72, 5(xr));:-. ..> (In, {(xn))} and the ini-
tialization parameters 8 (the same one as in Step 1.), which outputs o (-). Subtracting s(-), we
obtain on,e0 (+) 2x -(().
Output: At point æ, output on (x) )n,
The development of our PNC predictor is motivated by the challenge of computing ir (a)) (detailed below). To address this challenge, we characterize the procedural noise instead, which is given by
6n (+) == o (z) - h", (x) = 3gu (2)) s(x) + ((xx)) x²" (K(x, (xx x) + Ann (((()) sgu (x)).
By Proposition 3.1 the closed-form expression of 2n.0" in 4 corresponds exactly to the artificial- label-trained network in Step 2 of [lgorithm[[] Our artificial- label-trainne network is thereby used to quantify the procedural variability directly. This observation subsequently leads to:
Theorem 3.3 (PNC). Suppose that Assumption holds. Then the oulput of the PNC predictor
(Algorithm 1 is exactly h (z) given in 3
We discuss two approaches to compute 8(x) in Step of Algorithm 1 and their implications: 1) Under He initialization, sge (z) is zero-mean Gaussian process in the infinite width limit [80), and thus we may set 5 = 0 for simplicity. Note that even we set $ = O, it does not imply that the artificil-label--rraiedd neural network in Step 2 of Algorithm T will output a zero-constant network whose parameters are all zeros. In fact, neural networks are excessively non-convex and have many nearly global minima. Starting from an initialization parameter 8¹, gradient descent this artificial- label-trained neural network will find a global minimum that is close to the H¹ but not "zero" " even if "zero" is indeed one of its nearly global minima ("nearly" in the sense of ignoring the negligible regularization term) 33 118. 23 This phenomenon can also be observed in Proposition [3. by plugging in y O. Hence, the output depends on random initialization in addition training data, and our alternative approach that does not require specific initialization to use Monte Carlo integration: s(x) = lim N Σm Sop(2) where 0% are i.i.d. from Pas. When N is finite, it introduces procedural variance that vanishes at the order N-1. Since this computation does not involve any training process and is conducted in a rapid manner practically, N > n can be applied to guarantee that the procedural variance in computing ((x) is negligible compared with the data variance at the order n -1, We practically observe that both approaches work similarly well.
network is designed to capture this procedural variability. 2) An
Finally, we provide additional remarks on why h; (I) cannot be computed easily except using our Algorithm| First, the following two candidate approaches encounter computational issues: 1) One may use deep ensemble with sufficiently many networks in the ensemble. However, this approach is tinecconu as m networks in the ensemble mean m-fold training times. m is typically as small as five in practice [73] it cannot approximate h (z) well. 2) One may use the closed-
5

---

form expression of h; in (3) which requires computing the NTK (xx, and the inversion of the NTK Gram matrix K (x.x). K(x,x') is recursively defined and does not have a simple form for computation, which might be addressed by approximating it with the empirical NTK Ko(æ, x') numerically (See Appendix [.. However, the inversion of the NTK Gram matrix gives rise to more serious computational issue: the dimension of the NTK Gram matrix is large on large datasets, making the matrix inversion very time-consuming. Another approach that seems plausible is that: 3) One may initialize one network that is equivalent to s(x) and then train it. However, this approach cannot obtain h (1:) based on Proposition 3.1 because Proposition 3.1 requires random initialization. If one starts from a deterministic network imtialization such as a zero- zero-constant network, then the NTK theory underpinning the linearized training behavior of over-parameterized neural networks breaks down and the resulting network cannot be described by Proposition 3.1
3.3 Constructing Confidence Intervals from PNC Predictors
We construct confidence intervals for h² (2: leveraging our PNC predictor in Algorithm To handle data variability, two lines of works borrowed from classical statistics may be considered. First is an analytical approach using the delta method for asymptotic normality, which involves computing the influence function [ST 38 that acts as the functional gradient of the predictor with respect to
the data distribution. was introduced in modern machine learning for understanding training
point's s effect on model's prediction [72] [3]. The second is to use resampling [77 30 [02].
such as the bootstrap or jackknife, avoid explicit variance computation. The classical resampling method requires sufficiently large resample replications and thus incurs demanding resampling effort. For instance, the jackknife approach requires the number of training times to be the same as the training data size, which is very time-consuming and barely feasible for neural networks. Standard bootstrap requires a sufficient number of resampling and retraining to produce accurate resample quantiles. Given these computational bottlenecks, consider utilizing light-computation resampling alternatives, including batching [47 99 [00] and the so-called cheap bootstrap method [75] 76).
which allows valid confidence interval construction using few two model repetitions.
Large-sample asymptotics of the PNC predictor. derive our intervals, we first gain understanding on the large- sample properties of the PNC predictor. Proposition 3.1 shows that h (a) in (3) the
solution to the following empirical risk minimization problem:
7r (.) = =(-) + min gEH n
Σ (vv s(x;) - o(xx))) Aml1
(5)
Its corresponding population risk minimization problem (i.e., removing the data variability) is: hr (-) s(-) + min E. 8(X) ((x)) Aollgl1
gEa
where An = limm An. To study the difference between the empirical and population risk mini- mization problems of kernel ridge regression, we introduce the following established result on the asymptotic normality of kernel ridge regression (See Appendix B for details):
Proposition 3.4 (Asymptotic normality of kernel ridge regression (50)). Let H be a generic RKHS.
Suppose that Assumptions B.3 and B.4 hold. Let gP,x be the solution to the following problem:
gP.A == mingen Ep[(Y g((X
g1l where P nn or π. Then
Vrllgan>n 97.20) G in 11
where G is a zero-mean Gaussian process and >> represents "converges weakly". Moreover, at point s, / M (x) - Dn.20 (x)) > (0, (xx)) where (x) = Seexy I IF (2: T) (x))dx((2) and IF is the influence function of statistical functional gr,
Next, apply the above proposition to problems about h; Let Ti (P) be the solution of the following problem min [Y - sX) { + >ollall? that is evaluated at a point æ. Then we have the following surge-sample asymptotic of the PNC predictor, providing the thcoretical foundation for our subsequent interval construction approaches.
Theorem 3.5 (Large-sample asymptotics of the PNC predictor). Suppose that Assumption C.2 holds.
Suppose that Assumptions B.3 and B.4 hold when Y is replaced by Y - 3(X). Input the training

---

data D into Algorithm
1 to obtain (x) - @n We have
Vn ( -0, - h((x)) + (x),
(7)
where
d²()) JEEXX
(2; (
(8)
Thus, asymptotically (in the sense of + 00) exact (1 - al-eevel confidence interval of h" (z) is vn 91 --- where Qa is the a-quantile of
2(w) [hmn.o on (x) &n, vn 91 g (x) on + [(z)
the standard Gaussian distribution 10.
Theorem 3.5 does not indicate the value of o² (z). In general, (x() is unknown and needs to be estimated. It is common to approximate (2; with (2; T1, mn) and set (a:) E TF2 (2; T,, Tn (x)). This method known as the infinitesimal jackknife variance estimator
1102.. In Appendix B we derive the exact closed-form expression of the infinitesimal jackknife variance estimation o² (a), and further prove its consistency as well as the statistical coverage guarantee of confidence intervals built upon it. These results could be of theoretical interest. Yet in practice, the computation of (z(x) requires the evaluation of the NTK Gram matrix and its inversion, which is computationally demanding for large and thus not recommended for practical implementation on large datasets.
In the following, we provide two efficient approaches that avoid explicit estimation of the asymptotic variance as in the infinitesimal jackknife approach, and the computation of the NTK Gram matrix inversion.
PNC-enhanced batching. We propose an approach for constructing confidence interval that is particularly useful for large datasets, termed PNC-enhanced batching. The pseudo-code is given in Algorithm 2 Originating from simulation analysis [47 99 [OO), the key idea of batching is to construct a se self-normalizing (-statistic that "cancels out" the unknown variance, leading to a valid confidence interval without explicitly needing to compute this variance. can be used to conduct inference on serially dependent simulation outputs where the standard error is difficult to compute analytically. Previous studies have demonstrated the effectiveness of batching on the use of inference for Markov chain Monte Carlo [46] 40 66] and also the so-called input uncertainty problem (48). Its application deep learning uncertainty quantification was not revealed in previous work, potentially due to the additional procedural variability. Integrating it with the PNC predictor, PNC-eenanced batching very efficient and meanwhile possesses asymptotically exact coverage of its confidence interval, as stated below.
Algorithm 2 PNC-Enhanced Batching
Input:
Training dataset DD of size n. The number of batches m' 2 2.
Procedure: Split the training data D into m' batches and input each batch in Algorithm to output i. Re (x) 22.0n (ar) for [mm], where n'
ml
2. Compute (z())) ni Ejz (Kar Ha -0, (x)).
and Sy(x)² = m/ En'i (2)) - 22 o0 (x) (E(t))?
Output: At point x, output (x) and Sz(x))
Theorem 3.6 (Exact coverage of PNC PNC-enhanced batching confidence interval). Suppose that As- sumption C.2 holds. Suppose that Assumptions B.3 and B.4l when Y is replaced by Y 7 -(X).. Choose any m' 2 2 in Algorithm 2 Then an asymptotically exact (1 - al-eevel confidence interval of htt (x)))) is (Z)()) Sg(z)
Vmm 41-1 + 1 (B(r) + Sala)
vmn 21 -1- $], where go is the a-quantile of the distribution
bm with degree of freedom ml 1.
PNC-enhanced cheap bootstrap. We propose an alternative approach for constructing a confidence interval that works for large datasets but is also suitable for smaller ones, termed PNC-enhunced cheap bootstrap. The pseudo-code is given in Algorithm 3 Cheap bootstrap 75
76
74
is
modified bootstrap procedure with substantially less retraining effort than conventional bootstrap methods, via leveraging the asymptotic independence between the original and resample estimators
7

---

and asymptotic normality. Note that our proposal is fundamentally different from the naive use of bootstrap or bagging when additional randomness appears [73 8T 45], which mixes the procedural and data variabilities and does not directly provide confidence intervals. Like PNC-enhanced batching, PNC-enhanced cheap bootstrap also avoids the explicit estimation of the asymptotic variance. The difference between the above two approaches is that PNC-enhanced batching divides data into a small number of batches, and thus is suggested for large datasets while PNC-enhanced cheap bootstrap re-selects samples from the entire dataset, hence suited also for smaller datasets. On the other hand, in terms running time, when Rt = m' - 1, PNC-enhanced cheap bootstrap and PNC-enhanced batching share the same number of network training, but since batching is trained subsets of the data, the individual network training in PNC-enhanced batching is faster than PNC-enhanced cheap bootstrap. PNC-enhanced cheap bootstrap also enjoys asymptotically exact coverage of its confidence interval, as stated below.
Algorithm 3 PNC-EEhanced Cheap Bootstrap
Input: Training dataset ID of size n. The number of replications 1.
Procedure: 1. Input 'D in Algorithm to output h (a) 00.00 (2x).
2. For each replication E [R], resample D, i.e. independently and uniformly sample with replacement from (1r1 v1), ..... (Iin vn)) n times obtain D {(=]" "i"),,.,, )- (17)) ux2) Input
D in Algorithm to output h (0,0(()) (xx) - 02 (x).
Compute (x((x) = hn, 10 (xr) (x))
and Sc(x)) = 1/2 Eこ (hit (x) 3nd no (x) vt(()))
Output: At point z, output (1( and (c))?
Theorem 3.7 (Exact coverage of PNC-enhanced cheap bootstrap confidence interval). Suppose that Assumption C.2 holds. Suppose that Assumptions [B.3 and hold when Y is replaced by Y s(X). Choose any
in Algorithm 3 Then an asymptotically exact (1 - a-v-eeel confidence interval of he (2) is [c( - Sc(r) 4, (.Vc())) ((c(z))S(()))))))) Sc(r) $] where do is the a-quantile of the distribution tR with degree of freedom R.
Experiments
We conduct numerical experiments demonstrate the effectiveness of our approaches Our proposed approaches are evaluated on the following two tasks: 1) construct confidence intervals and 2) reduce procedural variability improve prediction. With known ground-truth regression function, training data are regenerated from the underlying synthetic data generative process. According to the NTK parameterization in Section C our base network is formed with two fully connected layers with n x 32 neurons in each hidden layer to ensure the network is sufficiently wide and over pre-pparamttrrieed. Detailed optimization specifications are described in Proposition C.3 Our synthetic datasets #1 are generated with the following distributions: X ~ Unif([0, 0.2l*) and
Y "> ()) + (0, (O0,0000011) 0.0013.. The training set D {(z,, y4) : i 1,..., n) is formed by drawing i.i.d. samples of (X, Y) from the above distribution with sample size n. We consider multiple dimension settings d = 2, 4, 8, 16 and data size settings
128, 256, 512, 1024 to study the
effects on different dimensionalities and data sizes. Additional experimental results on more datasets are presented in Appendix F The implementation details of our experiments are also provided in Appendix F
Constructing confidence intervals. We use xo (0.1...., 0.1) as the fixed test point for confidence intervals construction. Let yo Σi= (0.1) be the ground-truth label for To without aleatoric noise. Our goal is to construct confidence interval at zo for yo- To evaluate the performance of confidence intervals, we set the number of experiments J = 100. In each repetition j E [J], we generate a new training dataset from the same synthetic distribution and construct a new confidence interval [Lg (2o), 1 (xo)] with 95% or 90% confidence level, and then check the coverage rate (CR):
CR = 3 プド j 1 (o) Sm Sm (wo)- The primary evaluation of the confidence interval is based on whether its coverage rate is equal to or larger than the desired confidence level. In addition to CR,
"The source code for experiments is available at httpp:::///gthab.b.oo///////////
8

---

Table I: Confidence interval construction on synthetic datasets #1 with different data sizes = 128, 256, 512, 1024 and different data dimensions d 2, 4,8, 8, 16. The CR results that attain the desired confidence level 55% are in bold.
PNC-eahanced batching
PNC-enhanced cheap bootstrap
DropoutUO
95%CRCCR/W) SP%CIICR/W)
MP
95ECKCR/WW 90CCRC))W
MP
95%CCICC/W WSCCCCRWW)
MP
[d = 2)
0.98/0.443
0,95/0.0323
0.1998
0.98A0.0571
0.95/U.D433
0.1983
0.93/0.0564
0.87/0.440
0.2119
256
0.99/0.0301
0.96/0.0222
0.200M4
1.000.03766
0,95/0,0289
0.1997
0.99/0.0390
0.96/0.0127
0.2045
n = 512
0,97/0.9211
0.95/0.0156
0.2004
0.9800.024
0,96/0.0226
0.1992
0.95/0.0313
0,93/0.0267
0.2049
1024
0.96/0.0152
0,92/0.0112
0.2003
0,9700.2005
0,94/0.057
0,201
0,96/0.0287
0,89/0.0244
0.2052
[d = 4)
0.92/
0,94/0.0460
0.4012
1.000,08000
0,95/0.0629
0.4013
0,96/0.0868
0,90N.07422
0.4157
256
0.98/0.0411
0,95/0.0304
0.7991
0.9900.0569
0.96/1.0437
0.3988
0,97/0.0566
0,92/0.0481
0.4071
512
0.97/0.0296
0295
0.94/0.0218
0.3988
0.98/0.0396
0.96/0.0304
0.3989
0.96/0.0375
0,900000316
0.4045
1024
0,97/0.0213
0,95/0.0158
0.3983
0,98/0.0296
0,98/0.0227
0,3988
0,93/0.03-40
0,8410.0290
0.4055
[d 8)
128
0.98/0.0865
0,95/0.0639
0.7999
0.98/0.104
0,98/0.0832
0.7980
0.91/0.1345
0.87
0.8245
256
0.99/0.0594
0.96/0.0439
0.8014
1.000.07766
0.99/1.0611
0.7957
0.88/0.0791
0.81/0.0656
0.8152
= 512
0.98/0.0393
0.94/0.0290
0.7983
0.99/0.0588
0,99/0,0452
0,8012
0.99/0 0632
0,97/0.0538
0.7998
1024
0,95/0.0270
0,92/0.0200
0.7981
0,99000093
0,96/U.0302
0.7997
0.88/0.0421
0.82/9.0356
0,8040
[d I6)
128
0.95/0.1068
0,94/0.0790
1,5946
0.99
0,47/0.1 1044
1.6057
0.990.1555
0.960.1311
16093
256
0,98/0.0730
0,94/0.0540
1,5966
1,000. 1137
0,98/0.0833
1,5954
0. 0.90/O. 1072
0.86/0.0909
161777
512
0,98/0.0543
0.93/0.0401
1.5966
0.99/0.0788
0,98/0,0605
1.5972
0.86/0.099 0.00
0.800007990
1.6132
1024
0,95/0.0388
0.92/0.0287
1.5976
0,97X00050
0.95/U.0422
1,5980
0.87/0.0760
0,83/0.0647
1.6079
we also report the median point of the interval (MP): MP = E ((()(())+L and the interval width (IW): IW = (U,)()) (zo)) MP reflects the most likely point prediction for Zo. while IW reflects the conservativeness of the confidence interval, as a wide interval can easily achieve 95% coverage rate but is not practically preferable. In the confidence interval tasks, we use the tropouu-based Bayesian approximate inference (DropoutUQ) 43 as the baseline for comparison since it is one of the most widely-used representatives in Bayesian uncertainty quantification.
+ (zoo)
Table reports the CR, IW, and MP of 95% and 90% confidence intervals from our proposals and baselines. We have the following observations:
1) The CR values: For our two proposed approaches, in the majority of experiments, CR 2 95% for 95% confidence intervals, and CR 2 90% for 90% confidence intervals, thus satisfying the coverage requirements, while in very few experiments (in Appendix D the performances are degraded. Nonetheless, this occasional degeneration appears insignificant and can be potentially explained by the fact that the statistical guarantee of confidence intervals generated by PNC-enhanced batching and cheap bootstrap in the asymptotic sense. In contrast, DQ does not have such statistical guarantees. The intervals from DropoutUQ cannot normally maintain satisfactory coverage rate when they are narrow, although they have a similar or larger interval width as PNC-enhanced batching or cheap bootstrap. Moreover, we notice that the dropout rate has a significant impact on the interval width of DropoutUQ, and thus, additional tuning efforts are demanded for this baseline, while our approaches do not need this level of tuning procedure. These observations demonstrate the robustness and effectiveness of our proposals.
2) The IW values: Overall, the IW values are relatively small for both of our approaches, indicating that our framework is able to successfully generate narrow confidence intervals with high coverage rates. Moreover, as shown, the IW values from both of our approaches decrease roughly at the rate n along with an increased training data size which corroborates with our theoretical results well (Section[.3) In comparison with PNC-cnhanced cheap bootstrap, we observe that PNC-cnhanced batching tends to generate narrower intervals and has less training time per network, only subset of training data is involved in one training trial. In general, PNC-enhanced batching requires data splits, making its accuracy vulnerable to small sample sizes, and thus leading to less stable performance compared with PNC-enhanced cheap bootstrap. Therefore, we recommend employing PNC-enhanced cheap bootstrap small datasets and applying PNC-enhanced batching on datasets with sufficient training samples.
3) The MP values: MP appears consistent with the ground-truth label yo in all experiments, showing that our confidence intervals can accurately capture the ground-truth label on average under multiple problem settings.

---

Table 2: Reducing procedural variability to improve prediction on synthetic datasets #1 with different data sizes n 128, 256, 512, 1024 and different data dimensions d 2,4, 8, 16. The best MSE results are in bold.
MSE
One base network
PNC predictor
Deep ensemble (5 networks)
Deep ensemble (2 networks)
2)
(6.68 + 2.74) 10-
(3.68 +: 1.84) x 10-4
(4.22 1 1.74)
10-4
(6.28 +: 2.74) 10-4
256
(2.38 + 1.06)
(6.22+2.00 x 10
(8.86
4.94)
(1:25 + 0.67)
612
(1.03 # 0.98)
10
(2.06
0.72)) x x0
(3.32 1.10)
10
(5.11 3.33) 10-5
1024
(6.98 1 4.77)
10-5
(8.92 5.69) x 10
(1.72 0.77)
10
(5.75 1 2.28)
10-5
(d 4)
128
(2.11 # 1.49) x 10
(1.18
x 10
(1.53 0.60)
x 10 8
(1.83 # 0.80) x 10
a 256
(8.82 1 3.26) x 10"4
H 0.88)) x 10
(4.22 1 2.01)
10
(5.47 1 1,91) x 10"4
n 512
(5.35 1 1.91)
x 10"4
(1.92 0.53) x 10
(2.88 1 1.80)
10
(3.99 1 1.87) x 10"4
1 1024
(2:23
+ 0.83)
x 10
(4.22
(4.22
x 10
(8.50 2.39)
10
(1.65 + 0.37) x 10
8)
128
(4.07 1 1.04) 10-3
(2.54
10
(2.73 0.86)
(3.27 1 0.96) 10-3
256
0.703
10-3
(1.05 0.188 x 10-3
(1.31 0.33)
10
(1.48 + 0.38) 10-3
612
(1.34
10
0.700)
x 10
(7.40 1.35)
10
(1.08 # 0.40)
10
1024
(8.54 10"*
(2.02
x 10
(3.79 1.01) 10
(5.91 + 1.87) 10
(d 16)
n 128
(9.03 x 10-3 (6.00 0.83) x 10-3
(6.37 : 0.67) 10
(7.47 + 1.29) x 10-3
256
(6.76
x 10 (4.19 : 0.94) X 10
(4.83 1.14)
(5.46 + 1.49) x
(4.19
10
(L.60 (.35)) x 10-
(2.05
0.30)
(2.87 1 0.34) 10-3
1021
(3.16 0.35)
10 (7.41 0.99) x 10
(1.20 0.13) 10
(1.91 0.20) 10-3
Reduce procedural variability to improve prediction. In this part, we illustrate that our proposed PNC preditor (Algorithm [ can achieve better prediction by reducing procedural variability. We use empirical mean square error (MSE) to evaluate the prediction performance. Specifically, the empirical MSE is computed as MSE(h) == N. Ee (h(x²) - v1)² where (z/, v(), i= 1..., No. are i.i.d. test data independent of the training. We set Nte- = 2048 for test performance evaluation in all experiments. We compare the prediction performance of our proposed PNC predictor with the following approaches: base network, i.e., a standard network training with a single random initialization, and the deep ensemble approach with different numbers (m) of networks in the ensemble 77.. In deep ensemble, we consider m = 2 as it has a similar running time as our PNC (with one additional network training), and m = 5 as it is widely used in previous work 73 94.. All networks share the same hyperparameters and specifications.
Table 2 reports the sample mean and sample standard deviation of MSE from our proposed PNC predictor and other baseline approaches with 10 experimental repetition times. As shown, the proposed PNC achieves notably smaller MSE than baselines with similar computational costs (the base network and deep ensemble with two networks) in all experiments. Moreover, it achieves better or analogous results compared with the deep ensemble with 5 networks, but the running time of the latter is 2.5 times as much as that of PNC. It also has relatively small sample standard deviations, showing its capacity to reduce variability in networks effectively. These results verify that our proposed PNC can successfully produce better prediction and simultaneously bears significantly reduced computational costs compared to the classical deep ensemble.
5 Concluding Remarks
In this study, we propose a systematic epistemic uncertainty assessment framework with simultancous statistical coverage guarantees and low computation costs for over- werrpaaameetrized neural networks in regression. Benefiting from our proposed PNC approach, our study shows promise to characterize and eliminate procedural uncertainty by introducing only one artificial-label-teeee trained network. Integrated with suitable tight-computation resampling methods, we provide two effective approaches to construct asymptotically exact-coverage confidence intervals using few model retrainings. Our evaluation results different settings corroborate our theory and show that our approach can generate confidence intervals that are narrow and have satisfactory coverages simultaneously. In the future, we will extend our current framework to general loss functions that can potentially handle classification tasks to open up opportunities for broader applications.
10

---

Acknowledgments and Disclosure of Funding
This work has been supported in part by the National Science Foundation under grants CAREER CMMI-1834710 and IIS- 1849280, and the Cheung-Kon Kong Innovation Doctoral Fellowship. The authors thank the anonymous reviewers for their constructive comments which have helped greatly improve the quality of our paper.
References
[1] Y. Abbasi-Yadkori, D. Pál, and C. Szepesvári. Improved algorithms for linear stochastic bandits. Advances in Neural Information Processing Systems, 24:2312-2320, 2011.
[2] M. Abdar, F. Pourpanah, S. Hussain, D. Rezazadegan, L. Liu, M. Ghavamzadeh, P. Fieguth,
X. Cao, A. Khosravi, U. R. Acharya, et al. A review of uncertainty quantification in deep learning: Techniques, applications and challenges. Information Fusion, 76:243-299 297, 2021.
[3] Alaa and M. Van Der Schaar, Frequentist uncertainty in recurrent neural networks via blockwise influence functions. In International Conference on Machine Learning, pages 175-190. PMLR, 2020.
[4] Z. Allen-Zhu, Y. Li, and Y. Liang. Learning and generalization in overparameterized neural networks, going beyond two layers. Advances in Neural Information Processing Systems, 2019.
[5] z. Allen-Zhu, Y. Li, and Z. Song. A convergence theory for deep learning via over- parameterization. In International Conference Machine Learning, pages 242-252. PMIR, 2019.
[6] Z. Allen Zhu, Li, and Z. Song. On the convergence rate of training recurrent neural networks. Advances in Neural Information Processing Systems, 32:6676 668666666888666 2019.
[7] Arora, S. Du, W. Hu, Z. Li, R. R. Salakhutdinov, and Wang. On exact computation with infinitely wide neural net. Advances in Neural Information Processing Systems, 32:8141-8150, 2019.
[8] A. Ashukha, A. Lyzhov, D. Molchanov, and D. Vetrov. Pitfalls of in-domain uncertainty estimation and ensembling in deep learning. arXiv preprint arXiv:2002.064770 2020.
[9] P. Auer, N. Cesa-Bianchi, and P. Fischer. Finite-time analysis of the multiarmed bandit problem. Machine Learning, (7(2):235-255) 2002.
[10] Y. Bai, S. Mei, H. Wang, and C. Xiong. Understanding the under-coverage bias in uncertainty estimation. arXiv preprint arXiv:2106.05513 2021.
[11] F. Barber, E. Candes, A. Ramdas, and R. Tibshirani. The limits of distribution-fret 2019.
conditional predictive inference. arXiv preprint
[12] R. F. Barber, E. J. Candes, A. Ramdas, and R. J. Tibshirani. Predictive inference with the jackknife+. arXiv preprint arXiv:1905.022 02928, 2019.
[13] S. Basu, Pope, and S. Feizi. Influence functions in deep learning are fragile. arXiv preprint arXiv:2006. 2et.20006.4651.. 2020.
[14] A. Berlinet and C. Thomas-Agnan. Reproducing kernel Hilbert spaces in probability and statistics. Springer Science & Business Media, 2011.
[15] A. Bietti and J. Mairal. On the inductive bias of neural tangent kernels. Advances in Neural Information Processing Systems, 32, 2019.
[16] Blundell, J. Cornebise, Kavukcuoglu, and D. Wierstra. Weight uncertainty in neural network. In International Conference on Machine Learning, pages 1613-1622, PMLR, 2015.
[17] L. Breiman. Bagging predictors. Machine Learning, 24:123-140, 1996.
11

---

[18] L. Breiman. Random forests. Machine Learning. 45:5-34 -2,,,,2220 2001.
[19] Y. Cao and Q Gu. Generalization bounds of stochastic gradient descent for wide and deep neural networks. Advances in Neural Information Processing Systems,
2019.
[20] Y. Cao and Gu. Generalization enror bounds of gradient descent for learning over- parameterized deep relu networks. Proceedings of the AAAI Conference on Artificial In- telligence, (4(04))3349-3333 3356, 2020.
[21] D. Carvalho, R. Clark, A. Nicastro, and H. Kelly. Scalable uncertainty for computer vision with functional variational inference. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 12003-12013, 2020.
[22] H. Chen, Huang, Lam, H. Qian, and H. Zhang. Learning prediction intervals for regres- sion: Generalization and calibration. In International Conference on Artificial Intelligence and Statistics, 2021.
[23] L. Chizat, E. Oyallon, and F. Bach, On lazy training in differentiable programming. Advances in neural information processing systems, 32, 2019.
[24] A. Christmann and R. Hable. On the consistency of the bootstrap approach for support vector machines and related kernel-based methods. Empirical inference: Festschrift in honor of vladimir n. vapnik, pages 231-244, 2013.
[25] A. Christmann and I. Steinwart. Consistency and robustness of kernel-basee regression in convex risk minimization. Bernoulli, 13(3):799-811 2007.
[26] Christmann and I. Steinwart. How svms can estimate quantiles and the median. Advances in Neural Information Processing Systems, 90:305-312, 2007.
[27] F. Cucker and S. Smale. On the mathematical foundations of learning. Bulletin of the American Mathematical Society, 39(1):1 -9, 2002.
[28] G. Cybenko. Approximation by superpositions of a sigmoidal function. Mathematics of Control, Signals and Systems, 2(4):303-31 -14, 1989.
[29] N. Dalmasso, T. Pospisil, A. B. Lee, R. Izbicki, P. E. Freeman, and A. I. Malz. Condi- tional density estimation tools in python and with applications to photometric redshifts and likelihood-free cosmological inference. Astronomy and Computing, 30: 100362, 2020.
[30] A. C. Davison and D. V. Hinkley. Bootstrap methods and their application. Cambridge university press, 1997.
[31] M. Debruyne, M. Hubert, J. A. Suykens. Model selection in kernel based regression using the influence function. Journal of Machine Learning Research, 9(10), 2008
[32] S. Depeweg, J. M. Hemááneez-Lobtto F. Doshi-Velez, and S. Udluft. Learning and pol- icy search stochastic dynamical systems with bayesian neural networks. arXiv preprint 1605.:6005.77
2016.
[33] S. Depeweg, J.-M. Hernandez-Lobato, F. Doshi- and S. Udluft. Decomposition of uncertainty in bayesian deep learning for efficient and disk-sensitive learning. International Conference on Machine Learning, pages 1184-1193. PMLR, 2018.
[34] S. Du, J. Lee, H. Li, L. Wang, and X. Zhai. Gradient descent finds global minima of deep neural networks. In International Conference on Machine Learning, pages 1675-1685. PMLR, 2019.
[35] S. S. Du, X. Zhai, B. Poczos, and A. Singh. Gradient descent provably optimizes over- parameterized neural networks. In International Conference on Learning Representations, 2018.
[36] V. Dutordoir, H. Salimbeni, J. Hensman, and M. Deisenroth. Gaussian process conditional density estimation. Advances in Neural Information Processing Systems, pages 2385-2395, 2018.
12

---

[37] B. Efron and R. J. Tibshirani. An introduction to the bootstrap. CRC press, 1994.
[38] L. T. Fernholz. Von Mises calculus for statistical functionals, volume 19. Springer Science & Business Media, 2012.
[39] R. A. Fisher. Inverse probability. Mathematical Proceedings of the Cambridge Philosophical Society, 26()):228-35, 1930.
[40] J. M. Flegal, G. L.. Jones, et al. Batch means and spectral variance estimators in Markov chain Monte Carlo. The Annals of Statistics, 38(2):1034-107), 2010.
[41] S. Fort, H. Hu, and B. Lakshminarayanan. Deep ensembles: A loss landscape perspective. arXiv preprint arXiv:1912.02235 2019.
[42] P. Freeman, R. Izbicki, and Lee. A unified framework for constructing, tuning and assessing photometric redshift density estimates in selection bias setting. Monthly Notices of the Royal Astronomical Society, 468(4):4556 4565, 2017.
[43] Y. Gal and Z. Ghahramani. Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In International Conference on Machine Learning, pages 1050- 1059, 2016.
[44] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin. Bayesian data analysis. CRC Press, 2013.
[45] P. Geurts, D. Ernst, and L. Wehenkel. Extremely randomized trees. Machine Learning, 63:3-42, 2006.
[46] C. J. Geyer. Practical Markov chain Monte Carlo. Statistical Science, 7(4):473-483, 1992.
[47] P. W. Glynn and D. L. Iglehart. Simulation output analysis using standardized time series. Mathematics of Operations Research, 15(1):1-16, 1990.
[48] P. W. Glynn and H. Lam. Constructing simulation output intervals under input uncertainty via data sectioning. In 2018 Winter Simulation Conference (WSC), pages 1551-1562. IEEE, 2018.
[49] A. Graves. Practical variational inference for neural networks. Advances in Neural Information Processing Systems, 24, 2011.
[50] R. Hable. Asymptotic normality of support vector machine variants and other regularized kernel methods. Journal of Multivariate Analysis, 106:92-117, 2012.
[51] F. R. Hampel. The influence curve and its role in robust estimation. Journal of the American Statistical Association, 69(346):383-393, 1974.
[52] F. R. Hampel. Potential surprises. International Symposium on Imprecise Probability: Theories and Applications (ISIPTA), page 209, 2011.
[53] B. Hanin and M. Sellke, Approximating continuous functions by relu nets of minimal width. arXiv preprint arXiv: 1278, 2017.
[54] B. He, B. Lakshminarayanan, and Y. W. Teh. Bayesian deep ensembles via the neural tangent kernel. arXiv preprint arXiv:2007.05864, 2020.
[55] K. He, X. Zhang, S. Ren, and J. Sun. Delving deep into rectifiers: Surpassing human-level performance imagenet classification. In Proceedings of the IEEE International Conference on Computer Vision, pages 1026-1034, 2015.
[56] J. M. Hernández-Lobaut and R. Adams. Probabilistic backpropagation for scalable learning of bayesian neural networks. In International Conference on Machine Learning, pages 1861- 1869, 2015.
[57] M. P. Holmes, A. G. Gray, and C. L. Isbell Jr. Fast nonparametric conditional density estimation. In Proceedings of the Twenty-Thin Third Conference on Uncertainty Artificial Intelligence, pages 175-182, 2007.
13

---

[58] K. Hornik. Approximation capabilities of multilayer feedforward networks. Neural Networks, 4(2):251-257, 1991.
[59] W. Hu, Z. Li, and D. Yu. Simple and effective regularization methods for training on nois- ily labeled data with generalization guarantee. In International Conference on Learning Representations, 2020.
[60] Z. Huang, H. Lam, and H. Zhang. Quantifying epistemic uncertainty in deep learning. arXiv preprint a 10. /2/22, 2021.
[61] E. Hüllermeier and W. Waegeman. Aleatoric and epistemic uncertainty in machine learning: An introduction to concepts and methods. Machine Learning, 110(3):457-566, 2021.
[62] R. Izbicki and A. B. Lee. Nonparametric conditional density estimation in high-dimensional regression setting. Journal of Computational and Graphical Statistics, 25(4):1297-1316, 2016.
[63] R. Izbicki, A. B. Lee, and P. E. Freeman. Photo-z estimation: An example of nonparametric conditional density estimation under selection bias. Ann. Appl. Stat., 11(2):698-724, 06 2017.
[64] A. Jacot, F. Gabriel, and C. Hongler. Neural tangent kernel: Convergence and generalization in neural networks. arXiv preprint arXiv: 1806. 07572, 2018
[65] L. Jaeckel. The infinitesimal jackknife. memorandum. Technical report, MM 72-1215-11, Bell Lab. Murray Hill, NJ, 1972.
[66] G. L. Jones, M. Haran, B. S. Caffo, and R. Neath. Fixed-width output analysis for Markoy chain Monte Carlo. Journal of the American Statistical Association, 101(47))1537-1147, 2006.
[67] A. Kendall and Y. Gal. What uncertainties do need in bayesian deep learning for computer vision? arXiv preprint arXiv: 1703.04977 2017.
[68] A. Khosravi, S. Nahavandi, D. Creighton, and A. F. Atiya. Lower upper bound estimation method for construction of neural network-based prediction intervals. IEEE Transactions on Neural Networks, (2(3):337-346 2010.
[69] A. Khosravi, S. Nahavandi, D. Creighton, and A. F. Atiya. Comprehensive review of neural network-basee prediction intervals and new advances. IEEE Transactions Neural Networks, 22(9):1341-133 1356, 2011.
[70] R. Koenker and K. F. Hallock. Quantile regression. Journal of Economic Perspectives, 15(4):143-156, 2001.
[71] P. W. Koh, K.-S. Ang, H. H. Teo, and P. Liang. On the accuracy of influence functions for measuring group effects. arXiv preprint 1rXiv:1905.13289 2019.
[72] P. W. Koh and P. Liang. Understanding black-box predictions via influence functions. In International Conference on Machine Learning, pages 1885-1894. PMLR, 2017.
[73] B. Lakshminarayanan. A. Pritzel, and C. Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. In Advances in Neural Information Processing Systems, pages 6402 6413, 2017.
[74] H. Lam. Cheap bootstrap for input uncertainty quantification. In 2022 Winter Simulation Conference (WSC), pages 2318-2329. IEEE, 2022.
[75] H. Lam. A cheap bootstrap method for fast inference. arXiv preprint arXiv:2202.00090) 2022.
[76] H. Lam and Z. Liu. Bootstrap in high dimension with low computation. International Conference on Machine Learning (ICML), 2023.
[77] Lam and H. Qian. Subsampling to enhance efficiency in input uncertainty quantification. Operations Research, 7003):189111133 2022.
14

---

[78] H. Lam and H. Zhang. Doubly robust stein-kerellized monte carlo estimator: Simultane- ous birs-variance reduction and supercanonical convergence. Journal of Machine Learning Research, (1885)1-558. 2023.
[79] J. Lee, Y. Bahri, R. Novak, S. S. Schoenholz, J. Pennington, and J. Sohl-Dickstein. Deep neural networks as gaussian processes. arXiv preprint arXiv:1711.00165, 2017.
[80] J. Lee, L. Xiao, S. Schoenholz, Y. Bahri, R. Novak, J. Sohl-Dickstein, and J. Pennington. Wide neural networks of any depth evolve as linear models under gradient descent. Advances in Neural Information Processing Systems, 32:8572-8558
2019.
[81] S. Lee, S. Purushwalkam, M. Cogswell, D. Crandall, and D. Batra. Why m heads are better than one: Training a diverse ensemble of deep networks. arXiv preprint arXiv:151 151.06314 2015.
[82] J. Lei, M. G'Sell, A. Rinaldo, R. J. Tibshirani, and L. Wasserman. Distribution-fer predictive inference for regression. Journal of the American Statistical Association, 1136523);1094-1111 2018.
[83] J. Lei, A. Rinaldo, and L.. Wasserman. A conformal prediction approach to explore functional data. Annals of Mathematics and Artificial Intelligence, 74(1-2):29 43, 2015.
[84] J. Lei and L. Wasserman. Distribution-free prediction bands for non- non- porammercitticcccc regression. Journal of the Royal Statistical Society: Series (Statistical Methodology), 76(1):71-96, 2014.
[85] Y. Li and Y. Liang, Learning overparameterized neural networks via stochastic gradient descent on structured data. Advances in Neural Information Processing Systems, 2018.
[86] Z. Li, R. Wang, D. Yu, S. S. Du, W. Hu, R. Salakhutdinov, and S. Arora. Enhanced convolu- tional neural tangent kernels. arXiv preprint arXiv: 1r11vv1911.00009, 2019.
[87] N. Littlestone and M. K. Warmuth. The weighted majority algorithm. Information and computation, 108(2):112-611, 1994.
[88] Z. Lu, H. Pu, F. Wang, Z. Hu, and L. Wang. The expressive power of neural networks: A view from the width. In Proceedings of the 31st International Conference on Neural Information Processing Systems, pages
2017.
[89] N. Meinshausen. Quantile regression forests. Journal of Machine Learning Research, 7(Jun):833999, 2006.
[90] R. Michelmore, M. Kwiatkowska, and Y. Gal. Evaluating uncertainty quantification in end-to- end autonomous driving control. arXiv preprint 1rXiv:1811.0881 2018.
[91] Michelmore, M. Wicker, L. Laurenti, Cardelli, Y. Gal, and M. Kwiatkowska. Uncertainty quantification with statistical guarantees in end-to-end autonomous driving control. In 2020 IEEE International Conference on Robotics and Automation, pages 7344-7350. IEEE, 2020.
[92] M. Mirza and S. Osindero. Conditional generative adversarial nets. arXiv preprint arXiv: 1784,4, 2014.
[93] M. Mohri, Rostamizadeh, and Talwalkar. Foundations of machine learning. MIT Press, 2018.
[94] T. Pearce, M. Zaki, A. Brintrup, and A. Neely. High-quality prediction intervals for deep learning: A distribution-fres free, ensembled approach. In International Conference on Machine Learning, PMLR: Volume 80, 2018.
[95] K. Posch and J. Pilz. Correlated parameters to accurately measure uncertainty in deep neural networks. IEEE Transactions on Neural Networks and Learning Systems, 1037-10077- 2020.
[96] Y. Ren, Zhu, J. Li, and Luo. Conditional generative moment-matching networks, Advances in Neural Information Processing Systems, 29, 2016.
15

---

[97] Y. Romano, E. Patterson, and E. Candes. Conformalized quantile regression. In Advances in Neural Information Processing Systems, pages 3543-3553, 2019.
[98] N. Rosenfeld, Y. Mansour, and E. Yom-Tov. Discriminative learning of prediction intervals. In International Conference on Artificial Intelligence and Statistics, pages 347-355, 2018.
[99] Schmeiser. Batch size effects the analysis of simulation output. Operations Research, 30(3):556-568, 1982.
[100] L. Schruben. Confidence interval estimation using standardized time series. Operations Research, 3166):1000-1008, 1983.
[101] R. Senge, S. Bösner, K. Dembezyński, J. Haasenritter, O. Hirsch, N. DonnerBanzhoff, and E. Hüllermeier. Reliable classification: Learning classifiers that distinguish aleatoric and epistemic uncertainty. Information Sciences, 255:16-29, 2014.
[102] J. Shao and D. Tu. The jackknife and bootstrap. Springer Science & Business Media, 2012.
[103] A. Siddhant and Z. C. Lipton. Deep bayesian active learning for natural language processing: Results of large-scale empirical study. arXiv preprint 18080 2018.
[104] S. Smale and D.-X. Zhou. Shannon sampling and function reconstruction from point values. Bulletin of the American Mathematical Society, 41(3):279-306 2004.
[105] S. Smale and D.-X. Zhou. Shannon sampling II: Connections to learning theory. Applied and Computational Harmonic Analysis, 19(3):285-302 2005.
[106] S. Smale and D.-X. Zhou. Learning theory estimates via integral operators and their approxi- mations. Constructive Approximation, 26(2):153-177 2007.
[107] S. Sooodd,, I. Ishikawa, and M. Ikeda. Ridge regression with over- wer-parametrized two-layer networks converge to ridgelet spectrum. In International Conference on Artificial Intelligence and Statistics, pages 2674-2666 PMLR, 2021.
[108] I. Steinwart, A. Christmann, et al. Estimating conditional quantiles with the help of the pinball loss. Bernoulli, 17(1):211-225, 2011.
[109] H. Sun and Q. Wu. Application of integral operator for regularized least-square regression. Mathematical and Computer Modelling, (0(12)):76--885 2009.
[110] W. Van der Vaart. Asymptotic statistics, volume Cambridge University Press, 2000.
[111] V. Vovk, A. Gammerman, and G. Shafer. Algorithmic learning in a random world. Springer Science & Business Media, 2005.
[112] K. Wang. Pseudo-labeling for kernel ridge regression under covariate shift, arXiv preprint 1rXiv:2302.10160, 2023.
[113] Y. Xiao and W. Y. Wang. Quantifying uncertainties in natural language processing tasks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 7322-7329, 2019.
[114] G. Yang. Scaling limits of wide neural networks with weight sharing: Gaussian process behavior, gradient independence, and neural tangent kernel derivation. arXiv preprint arXiv:1902.04760. 2019.
[115) G. Yang. Tensor programs ii: Neural tangent kernel for any architecture. arXiv preprint arXiv:2006. /4548, 2020.
[116] H. Zhang, J. Zimmerman, D. Nettleton, and D. J. Nordman. Random forest prediction intervals. The American Statistician, pages 1-15, 2019.
[117] W. Zhang, D. Zhou, I.. Li, and Q. Gu. Neural thompson sampling. International Conference Learning Representations, 2021.
16

---

[118] Y. Zhang, Z. Xu, T. Luo, and Z. Ma. A type of generalization enror induced by initial- ization in deep neural networks. In Mathematical and Scientific Machine Learning, pages 144-164, 164. PMLR, R, 2020.
[119] D. Zhou, L. Li, and Q. Gu. Neural contextual bandits with ucb-based exploration. In International Conference on Machine Learning, pages 11492-11502. PMLR, 2020.
[120] M. Zhou, R. Ge, and C. Jin. A local convergence theory for mildly over-parameterized two-layer neural network. In Conference Learning Theory, pages
PMLR, 2021.
[121] Zhou, Y. Jiao, J. Liu, and Huang. A deep generative approach conditional sampling. Journal of the American Statistical Association, pages 1-12, 2022.
[122] D. Zou, Y. Cao, D. Zhou, and Q. Gu. Gradient descent optimizes over- parameterrzzed deep relu networks. Machine Learning, 109(3):467 2020.
17

---

Appendices
We provide further results and discussions in this supplementary material. Section A discusses aleatoric uncertainty and predictive uncertainty. Section B discusses statistical inference for kernel- based regression and, in particular, the asymptotic normality of kernel-based regression. Section C discusses the training of the linearized neural networks based on the NTK theory. In particular, we show in Proposition (...ssttoonnnnthhh.. that the shifted kemel ridge regressor using NTK with a shift from an initial function Sgb is exactly the linearized neural network regressor that starts from the initial network Sgh. Section D provides additional remarks explain some of the details in the main paper. Section E| presents the proofs for results in the paper. Section presents experimental details and more experimental results.
Appendix A
Other Types of Uncertainty
Back in 1930, [39] first introduced formal distinction between aleatory and epistemic uncertainty in statisties [52]. while in modern machine learning, their distinction and connection were investigated in [I01] and further extended to deep learning models 67 33]. In 194], the differences between procedural and data variabilities in epistemic uncertainty were intuitively described; however, no rigorous definition or estimation method was provided. Our previous draft [60] discussed intuitive approaches for quantifying procedural and data variabilities, which some the ideas in this work originate from. In the following, we present an additional discussion other types of uncertainty for the sake of completeness. Section A.T presents aleatoric uncertainty and Section A.2 presents predictive uncertainty.
A.1 Aleatoric Uncertainty
We note that if could remove all the epistemic uncertainty, i.e., letting UQ O. the best predictor we could get is the Bayes predictor hig- However, h as a point estimator, is not able to capture the randomness in x if it is not point mass distribution.
The easiest way to think of aleatoric uncertainty is that is captured by Zy|x. At the level of the realized response or label value, this uncertainty is represented by
UQAU = (g(
If the connection between X and Y is non-deterministic, the description of new prediction problem involves a conditional probability distribution Ty Standard neural network predictors can only provide a single output y. Thus, even given full information of the distribution the uncertainty in the prediction single output y remains. This uncertainty canmot be removed by better modeling more data.
There are multiple work aiming estimate the aleatoric uncertainty, i.e., learn Ty|x. For instance, conditional quantile regression aims to learn each quantile of the distribution Tylx 770
89
108]
Conditional density estimation aims to approximately describe the density of Ty
57
36
42
63
However, we remark that these approaches also face their own epistemic uncertainty in the for recent studies on the epistemic bias in quantile regression.
estimation. See 226 [08 10
A.2
Predictive Uncertainty
Existing work on uncertainty measurement in deep learning models mainly focuses on prediction sets (predictive uncertainty), which captures the of epistemic and aleatoric uncertainties 94 97
12
(22).
In certain scenarios, is not necessary estimate each uncertainty separately. A distinction between aleatoric and epistemic uncertainties might appear less significant. The user may only concern about the overall uncertainty related to the prediction, which is called predictive uncertainty and can be thought as the summation of the epistemic and aleatoric uncertainties.
The most common way to quantify predictive uncertainty is the prediction set: We aim to find map C. X + 2 which maps input to subset of the output space that for each test point
18

---

Zo E A' the prediction set ((z)) is likely to cover the true outcome yo [ prediction set also called prediction interval in the case of regression 68 69 The prediction set communicates predictive uncertainty via a statistical guarantee on the marginal coverage, i.e.,
83 82)
This
P(yu E ((xu)) 2 8
for a small threshold 8 > O where the probability is taken with respect to both training data Drr (epistemic uncertainty) for learning O and the test data (Zo,yo) (aleatoric uncertainty). It is more tempting obtain a statistical guarantee with only aleatoric uncertainty by considering the probablity conditional on the prediction set, i.e.,
P(yo E ((x))()) 2 5
However, this guarantee is in general very hard to achieve in the finite-sample case. Even asymp- totically. (9) easy achieve unless we have a very simple structure on the data distribution [98] [16]]. A recent study show that (9) could hold in the finite-sample sense if we could leverage a set of validation data [22
Appendix B Statistical Inference for Kernel-Based Regression
B.1 Classical Statistical Learning Framework and Kernel Ridge Regression
Following Section 2 assuming that the input-output pair (X, (X,Y) a random vector following an unknown probability distribution on X' x y where X X Rd is input and E E the corresponding output. Let the marginal distribution of X be x (x) and the conditional distribution of Y given be (y)=).
Suppose learner has access to set of training data D {(a11, (x2, v2), (Ii, Un)). which is assumed to be independent and identically distributed (i.i.d)) according to the data distribution æ. Let (I1, .... )².)" and (y1, )² for short. Let TD or TT denote the empirical distribution associated with the training data D, where Ty is used to emphasize the sample size dependence in asymptotic results:
n
TD nT
((),))
The goal of the learner is to obtain function h X' + y such that (=x) is "good" predictor for the response if X = x is observed for any E X. This is typically done by minimizing the following (ground-truth) population risk
RA (h)
(Chx))
where L y R is the loss function. For instance, L can be the square error in regression or cross-entropy loss in classification. If h is allowed to be any possible functions, the best predictions in the sense of minimizing the risk are described by the Bayes predictor hi [61] 93
h (X) arg min Ey xx [C(j, Y)
jEY
Note that hi cannot be obtained in practice, since the conditional distribution Ty xx (w/x) is unknown. In particular, the least squares loss [(h(X), = (h(X) Y)² yields the gground-rruth regression function defined by
gf (x) = B(x.Y)~> Y = x).
As the ground-truth distribution æ unknown, we can neither compute nor minimize the population risk R (h) directly. We define the risk under general distribution P.
(h == [Exx)>>(((((()),
In particular, for P nn (the empirical distribution associated with the training data D), an empirical
risk is derived:
Rr, (h) =
c(h(x:).9))
19

---

In practice, minimizing the risk over h is restricted to a certain function class. Let H. be hypothesis class, which is set of functions {h + V/h H). For instance, 1) H. could be nonparametric class such as reproducing kernel Hilbert space (RKHS). The resulting method is known as kernel- based regression; See below. 2) could be a parametric class (he X + y/0 E )) where 0 the parameter, and θ is the set of all possible parameters, e.g., the linear coefficients in a linear regression model, or the set of all network parameters in neural network model.
In classical statistical learning, one is interested in finding a hypothesis In E H that minimizes the population risk
9x arg min R (h). LEH
(10)
which is called the true risk minimizer. We remark that , is the best choice in the sense of minimizing the risk within the hypothesis set H. and different choices of H may lead to different hg. As π is unknown, the learner may consider minimizing the empirical risk:
Ix.
1= arg min Rx, (h)
(11)
hEH
which is called the empirical risk minimizer. More generally, to avoid overfitting in the finite sample regime, the learner may consider a regularized empirical risk minimization problem:
I arg min R.,,, (h) + >m Aml1
(12)
and its corresponding population problem
An arg min Ro (h) + Ao |/
(13)
NEH
which are called the true regularized risk minimizer and the empirical regularized risk minimizer, respectively, Here, An 2 0 is the regularization parameter, which may depend on the data size n, and we assume it has a limit Åy lim, An In general, the target 9x is not equal to Or..A, or gn.. We omit 0 in the subscript if X = and write 9r., == 9a..0 and gx = 9x.0 for simplicity.
The above framework fits classical machine learning approaches such as linear regression or kerne-- based convex regression (such as kernel ridge regression), since (12) as well as (13) has a unique solution in this setting. However, this framework is not suitable for deep learning: the empirical
(regularized) risk minimizer 9r.,,,, not unique and cannot be precisely obtained due to the non- convex nature of neural networks. Therefore, the gradient- based methods in deep learning can only find an approximate solution fr..,., subjected to procedural variability.
Kernel ridge regression. Next, we apply the above framework to kernel ridge regression and review some basic existing results about kernel ridge regression, which can be found, e.g., in 04 78 [04]
106 27 109].
The based regression means that the hypothesis class H. in statistical learning is chosen to be an RKHS. Formally, a Hilbert space H. consisting of functions h : X + R with an inner product ( H H + R is a RKHS if there exists a symmetric positive definite function k X' x X - R called a (reproducing) kernel, such that for all x E X, we have ((;; E H and for all x E X' and h H. we have h(x) (n(-), k(; 2)). We use the above l to denote the kernel associated with H. Note that for any symmetric positive definite function k. we can naturally construct, from k; an RKHS H whose kernel is exactly k [14).
When H is a RKHS, Ir..A., in 12 can be computed for a number of convex loss functions L. In particular, for the least squares loss L(h(X), Y) (h(X) - Y) the resulting problem (12 is well known as the kernel ridge regression, and there exists a closed-form expression for both fr,.,,, and Ir.Ag= as we discuss below.
Formally, the kernel ridge regression problem is given by
Im, (a) == arg min
15 g(0,))2 + Am
A1 1=1
where >,, > 0 a regularization hyperparameter that may depend on the cardinality of the training data set. There is a closed-form formula for its solution, stated below.
20

---

Proposition B.1. Let
k(x, x) (k(I,, 2,))nn R
be the kernel Gram matrix and
k(x, x) = (k( k(x, xT
Then the (unique) kernel ridge regression solution is given as Am (4)) k(x, ))" (k( Aanl)
is worth mentioning the linearity of kernel ridge regression: When we have two outputs = (yn,... yn and y (9/1 "..1"191)) yn 2/ we have
k (k<(x, ar) + Ann (y +y)
= k x) (k(x, x) + A,nd) + k(x, x)² (k(x, + Xn
which means the solution to
arg min sEH
Uy + v/S g(0,))² + + Anllelie
the summation of the solution to and the solution to
n
arg min gEH
-
Σ
(1) g(x;))² Anl19l12
11
j=1
n
arg min gEH
Σ
wu/ - g(x;))) + An
n
Define O₂ I + (Txx the kernel integral operator
(O) = / k(x, xt g')((()))x)xx( E. DE (mxx).
where [22(Tx)) = [s: ( co).
109 shows that O2 is compact and positive self-adjoint linear operator on L²( Note that since
LE is positive: O 2 O. we have that O₂ + A strictly positive for any A > 0: Og + A > O, and thus its inverse operator (Ox >n exists and is unique. Note that (Ox, AnT) -1 is also linear operator on L (xx))
Next, consider the population risk minimization problem corresponding to Ir..A, as follows:
ON,Ao = arg min
[E= - g( ((x)) {o1lg112
9E
It is easy to see that this problem is equivalent to
Ag = arg min {
- 9(x)) +
9EH
where (xr) )x V = x is the (ground-truth) regression function. This problem has the following explicit closed-fie -losed-form expression of the solution (a proof can be found in (27):
Proposition B.2. The solution of gr.>0 is given as G=.Ac = (O2 + -
The linearity of kernel ridge regression also holds for the population-level solution:
(g + g d') =(O, + An -0
=(O, + >o - +
=(Os + -1 + (Ox + -o
-gn + 9'r
for any two functions g. (Tx) since both (Ox + 1) and O₁ are linear operators on L (1Tx).
21

---

B.2 Asymptotic Normality of Kernel-Based Regression
In this section, we review existing results on the asymptotic normality of kernel-based regression, which is established using the influence function concept [25] 24| 50].
Let H. be generic RKHS with the associated kernel k( x>). Let || || denote the norm H. Note that the feature map k₂ = k(z;)) :) is a function in H. We consider the connection between the true regularized risk minimizer gs Ao and the empirical regularized risk minimizer 9xn.An introduced in (11) and (12) when H. is a RKHS.
The asymptotic normality of kemel--aaee regression roughly states that under some mild conditions, yn(gr..A. (x.Ag) is asymptotically normal as n + oc. This result provides the theoretical foundations to conduct asymptotic statistical inference for kernel-based regression, in particular, building asymptotic confidence interval for the ground-trume truttttuthhhhhh 9n.Ao
To be rigorous, we first list the following assumptions introduced by 500
Assumption B.3. Let (Ω, .A, Q) be a probability space. Suppose that A' C F is closed and bounded with Morel-o--aleebra (x)). and y I is closed with horel-o-algebre B()). The worel-o-algeebe of X' x y is denoted by ( x }). Let H. be an RKHS with kernel & and let be probability measure on (X' x y, B(X x V)). Assume that the kernel of H. ks X X' an m-times continuously differentiable kernel i Rd Rd - IR such that m > d/2 and + 0. Let
the restriction of
Å E (0, +x)) be any positive number. Suppose that the sequence An satisfy that An λo ol )
Assumption B.4. Let C y x R + [0, +oo) be loss function satisfying the following conditions:
C is a convex loss, i.c., z + L(y, 2) is convex in a for every y E y.
The partial derivatives
8 Cly 2)
C'(y,2) = Oz ((y,2,, ("(y, = 02
exist for every (v,z) E/ IR.
The maps
(u,z) + L'(,,)), (v,z) + f'
are continuous.
There a E (,y),, and for every a E (0, +oo), there is a by E I (mTy and a constant 8 E (0, +00) such that, for every V,
[C(y, 2) < b(y) + |2| Vz,
sup =E|-a,a]
(',, < (,(2),
sup =E[-o,a]
(C" 6%
where 2 is a constant.
Note that the conditions on the loss function Assumption B.4 are satisfied, c.g., for the logistic loss for classification or least-square loss for regression with π such that Exy [Y+]] < co. Therefore, Assumption B.4 is reasonable for the kernel ridge regression problem we consider in this work.
Theorem B.5 (Theorem 1.1 in (50). Suppose that Assumptions B.3 and B.4 Then there is tight, meeeabbe Gaussian process G Ω + H, w + ((w) such that
VnlgT,,>, gn. + G in H
where + represents "converges weakly". The Gaussian process G is zero-mean; i.e., Bigg G)xl = for every g H.
Note that Theorem B.5 implies that for any H. the random variable
Ω R. - (g. (Gw))x
has zero-mean normal distribution. In particular, letting = kig the reproducing property of k implies that,
Vrlggn. (2) Ao G(z)
22

---

To obtain the variance of this limiting distribution, denoted as 5 (z), we review the tool of influence functions follows.
Let P he a general distribution on the domain X' x V. Let 2 = (Zx, E X x y be a general point. Let T be general statistical functional that maps distribution P to real value, that is, T; P + T(P) E IR. The Gâteaux derivative of ((),, T-(-), is defined as
T(P + =Q) T(P)
(Q(Q)
lim 50+0
Let
Ps. (1- + = =(6, -
where 8 denotes the point mass distribution in 2. Then the influence function of ((P) at the point a special Gâteaux derivative defined as
T(P -
[F = lim c-0
Tp(6 P))
Under some mild conditions (e.g.. is Hadamard differentiable 38 [TO]), for the statistical functional T holds:
the central limit theorem
(n(T(m))-T(T) + ( (x, (1:
as n + oo where nn is the empirical distribution associated with n samples i.i.d drawn from x. Applying the above concepts to the kernel-based regression, we let T(P) be the solution to the following problem:
min Er g(x)) Aolggl
gEH
Note that (50 has shown that Vn(gr,, (x) Gr,Ao (x)) has the same limiting behavior as
Vnlgg,,.oo (x:) Ax (x)) so only the limiting regularization hyper-parameter Ao is used T(P).
The proof of Theorem 3.1 in [50) provides the doseddffrred expression of the Gâteaux derivative Ty, (Q) as follows:
= (xx
for Q lin(Bs) where lin By) corresponds to a subset of finite measures on x x y, x ))) (See Proposition A.8. in [50) and Sp H + H is defined by
Sp(f) 2hof Ep (Y (X)))
where is the feature map of H. (e.g., we can simply take (=x((-) k(x,,) :)) Note that Sp is continuous linear operator that is invertible (See Proposition A.5. in (50). In particular, letting Q P. we obtain that
- = (
- L' Sp. (zx) (zx)) (zz). (z.))
(14)
Note that the above special Gâtcaux derivative (influence function) of T(P), (p P). had been derived initially in (25
The proof of Theorem 3.1 in (50) shows that T(P) Hadamard differentiable, and the asymptotic normality holds. Using the above expression of the influence function, conclude that:
Vn(gTa.Aa (z) (z)) = (000g(())))
where
(x(x) = /
SSXXX (2)
and IF(2; is given by (14
IF(2:T, (ō m) -
(E, (T, (X)
Summarizing the above discussion, Proposition 3.4 follows.
23

---

B.3 Infinitesimal Jackknife for Kernel-Based Regression
In this section, we follow up on our discussion in Section 3.3 on the infinitesimal jackknife variance estimation. We derive the closed-formed expression of the infinitesimal jackknife variance estimation for (a)) in Section [.2 in the kernel ridge regression. We also show the consistency infinitesimal jackknife variance estimation in general kemel-based regression. Our consistency result appears new in the literature.
First, we estimate (2o) as follows:
1 Σ (2; 2,ED
(z0) =
where D {I11 91). (22. 12)..... (In, 8n)) and Tr is the empirical distribution associated with the data D. This estimator is known the infinitesimal jackknife jackkinff variance estimator [65]. In this section, we use To instead æ in the variance or influence function at test point Zn to avoid confusion between To and 2z. The latter Zn is referred to the z-componenet of 2.
Based on 25 follows.
31
50
we derive the closed-form expression of the variance estimation (2( as
Theorem B.6 (Expression of infinitesimal jackknife for kernel ridge regression). Suppose that Assumptions B.3 and B.4 hold. For the least squares loss [(h(X), Y) (h(X) -Y)* in kernel ridge regression, IF F 2; Tn) ) is given by
IF T, (zo) k(zo. x) (k(x, x)) -ond)- Mc(x) Me(Fo)
(15)
where Mz(x) = (M, (z) M2(x,)) and
M. Mx, = Irn ( 1 (zy Irn (2x)))
2o
for H, = Fo, F1,
Zn-
Note that although Hadamard differentiability is able to guarantee the central limit theorem of T. the consistency of variance estimation generally requires more than Hadamard differentiability. In general, we need some additional continuity condition (such as continuously Gâteaux differentiability) to guarantee that F T, nT) is indeed "close" IF2 (2; T, T) and EE.eD IF2 T, mn (250 indeed "close" SΛxAxY IF2 (2; T. (x) (T (z). We show that this is achievable for the infinites- imal jackknife variance estimator for kernel ridge regression by only imposing a weak assumption.
Theorem B.7 (Consistency of infinitesimal jackknife for kernel-based regression). Suppose that Assumptions and B. hold. Moreover, assume that y C IR is bounded, and b and by in Assumption B.4 ar4are& bounded on V. Then we have
1 Σ
(x0)
(2) = (20)) a.s.
=EXXY
2;ED
as n + oc. Hence, an asymptotically exact (1 - lev-evel confidence interval of 9n,2o (x)) is
((x0) 91 2. Ir (xv) {( 91 2
Inn.An (oo)
where qa is the a-quantile of the standard Gaussian distribution N (0, 1).
The proof of Theorems B.6 and B.7 is given in Appendix E
B.4 PNC-Enhanced Infinitesimal Jackknife
In this section, we return to our problem setting in Section B.3 to our PNC predictor and develop the PNC-enhanced infinitesimal jackknife confidence interval for par-parameterizee neural networks.
3.3 We apply the general results in Section
24

---

Recall that the statistical functional Ti is associated with the following problem:
1 Ev, s(z,) - g(zx)) An.llll
on (-) (n) s(-) min 9 72
Consider the PNC-enhanced infinitesimal jackknife variance estimator:
0² (x0) πE
1
(F; 2:ED
Σ
We obtain that
Theorem B.8 (Expression of PNC-enhanced infinitesimal jackknife). Suppose that Assumption
C.2
holds. Suppose that Assumptions B.3 and B.4 hold when Y is replaced by Y 8(X). Then
2; T1, Ta) given by
IF(2; (20)) = K (Fo, x)² (K (K(x, x) + Aond) (Ma((x) -M (20)
(16)
where M,(x) = (M₂ (x1) M2 (xn))) and
1 M2(1;) == hn (I;; - 2n - s(z;)
- (2v - (hn,or(z)) - (2 (2n)
do
for v, $1, In. Hence
(x00) Σ
1F Ti Ti,tTi (F))
2ED
K (No, x)² x)) + ) -M,, (x) - M2,
LED
Theorem 8.8
immediately follows from Theorem B.6
Theorem B.9 (Exact coverage PNC-enhanced infinitesimal jackknife confidence interval). Suppose that Assumption C.2 holds. Suppose that Assumptions B.3 and B.4 hold when Y is replaced by - 3(X). Moreover, assume that y C bounded. Then we have
8² (x(0)
Σ (2; TI 2ED
1
J (2: ()))dd = a
a.s.
as n + oc. Hence, an asymptotically exact (1 - ae)-leel confidence interval of h" (av) is (xx)
( 2n (xo)
(xr) a1 2 on (xT) - 2n,,s (2o) +
91
n
where qa the a-quaantie of the standard Gaussian distribution (O, 1) and the computation of (2x) is given by Theorem B.&
Note that in our problem setting, we can set b(y) = 2 and 1/ = a + |y| for kernel ridge regression. Since y is bounded, b and b2 are naturally bounded on y. Therefore, Theorem B.9 follows from Theorem B.7
Appendix C Theory of Neural Tangent Kernel
In this section, we provide brief review of the theory of neural tangent kernel (NTK). Then discuss particularly one result employed by this paper, i.e.,
The shifted kernel ridge regressor using NTK with shift from initial function Sgs is exactly the linearized neural network regressor that starts from the initial network sge.
has attracted lot of attention since it provides a new perspective on training dynamics, generalization, and expressibility of over-parameterized neural networks. Recent papers show that when training an over-parametrized neural network, the weight matrix at each layer is close to its
NTK 64
25

---

An over wwrr-parameteiized neural network can rapidly reduce training enror to zero via gradient descent, and thus finds global minimum despite the objective function being non-convex Moreover, the trained network also exhibits good generalization property [19 20) These observations are implicitly described by a notion, NTK, suggested by (4) This kernel is able to characterize the training behavior of sufficiently wide fully-connected neural networks and build
initialization 85 35|
34
4 122
a new connection between neural networks and kernel methods. Another line of work is to extend NTK for fully-connected neural networks to CNTK for convolutional neural networks 7 [4 86
NTK is a fruitful and rapidly growing area.
We focus on the result that is used as the foundation of our epistemic uncertainty quantification: Since the gradient of over- parameterizee neural networks is nearly constant and close its initialization, training networks is very similar to shifted kernel ridge regression where the feature map of the
kernel is given by the gradient of networks [D]. Multiple previous works 7 80 54
59 118 L5]
have established some kind of equívalences between the neural network regressor and the kernel ridge regressor based on wide variety of assumptions or settings, such as the full-rank assumption of NTK, or zero-valued initialization, or introducing small multiplier K, normalized training inputs ||æ:||2 = 1. These results are also of different forms. Since a uniform statement is not available, we do not intend to dig into those detailed assumptions or theorems as they are not the focus of this paper (Interested readers may refer to the above papers).
In the following, we provide some discussions on how to obtain the equivalence in Proposition 3.1 adopted in the main paper, based on a less (technically) rigorous assumption. This assumption is borrowed from the linearized neural network property introduced in [80], where they replace the outputs of the neural network with their first-order Taylor expansion, called the linearized neural network. They show that in the infinite width limit, the outputs of the neural network are the same as the linearized neural network. [80] does not add regularization in the loss. Instead, we introduce a regularization An > 0 in the loss since regularization is common and useful in practice. Moreover, it can guarantee the stable computation of the inversion of the NTK Gram matrix and can be naturally linked to kernel ridge regression. In the following, we will review the linearized neural network assumption (Assumption C.I) as well as other network specifications. Based on them, we show that Proposition C.3 holds, which provides the starting point for uncertainty quantification framework. We also provide some additional remarks in Section D
NTK parameterization. We consider fully- collecccneeeted neural network with any depth defined formally as follows. Suppose the network consists of L hidden layers. Denote g(d) = x and dg = d. Let
(z) = ww g0) = =V V
where ((()) pd is the weight matrix in the (-th layer, is coordinate-ve wise activation function, co = 1o. and / E [4]. The output of the neural network is
fo(x) == (((+) (x) = w(+ (2)
where ((++1) E gixdr, and (W()).., ..., (L+1) represents all the parameters in the network. Note that here we use NTK parametrization with a width-dependent scaling factor, which is thus slightly different from the standard parametrization. Unlike the standard parameterization which only normalizes the forward dynamics of the network, the NTK parameterization also normalizes its backward dynamics.
He Initialization. The NTK theory depends on the random initialization of the network. Suppose the dimension of network parameters 0 is p. We randomly initialize all the weights to be i.i.d. A(0,, 1) random variables. In other words, we set Po, = (0, In) and let 0³ be an instantiation drawn from Po.. This initialization method is essentially known as the He initialization [55].
NTK and Linearized neural network. With the above NTK parameterization and He initialization, the population NTK expression is defined recursively as follows: (64 71 For / [L].
(00 = TTa ER,
(e,xx') = (20-)) (x-1)
E 22
(x,)) {d-1)(d'
(Z,z)) Co CcB(unn)-
26

---

We also define derivative covariance: y "0 x') co (uBaaa))wrrrannnoo/od(/
The final population NTK is defined as
L+1
K(x, x')
I1-1) x')
II E
where ()(++))(rr x²) = for convenience. Let (. be the standard inner product in RP.
Let x(0; Vafe(z) E IR denote the gradient of the network. The empirical (since the initialization random) NTK matrix is defined as
Ko(z, x') J(O;
In practice, we may use the empirical NTK matrix numerically compute the population NTK matrix.
One of the most important results in the NTK theory is that the empirical NTK matrix converges to the population NTK matrix as the width of the network increases (64 114 715.7:
Ko = (()
K(x, x')
where + represents "converge in probability" 7 or "almost surely" "
The NTK theory shows that =(;;)) = Vofe(r) is approximately a constant independent of the network parameter 0 (but still depends on x) when the network is sufficiently wide. Therefore, when we treat ((;; x) as constant independent of the network parameter 9, we can obtain a model obtained from the first-order Taylor expansion of the network around its initial parameters. which called linearized neural network [80]. Our study on epistemic uncertainty is based on the linearized neural network.
To be rigorous, we adopt the following linearized neural network assumption based on the NTK theory:
Assumption C.1. [Linearized neural network assumption] Suppose that J(0; .r) = Vofo(x) E 1 independent of 0 during network training, which thus denoted as J(x). The population NTK is then given by K (xx,,' x²) J())) [xx)
We introduce the NTK vector/matrix as follows. For = (x1,,..; we let
K( x) := (K (zi, (2),,,
[x))((x) g=x=
be the NTK Gram matrix evaluated on data For an arbitrary point x E Rd, we let R" be the kernel value evaluated between the point i and data æ, i.c.,
<(x,x) (K(x, (2), K(x, x,,)) Jx)J[(x) Elx
Training dynamics of the linearized network. Based on Assumption C.I the training dynamics of the linearized network can be derived as follows.
We consider the regression problem with the following regularized empirical loss
R(fe)
((fo(z)). A..1 8
(17)
n
where 6h ((0) is the initialization of the network. In the following, we use fe(o) to represent the network with parameter ((L) that evolves with the time l. Let y be the learning rate. Suppose the network parameter is trained via continuous-tim gradient flow. Then the evolution of the parameters (t) and fo(e) can be written as
d0(t) n
(1 ( Core C+ (0() 8(0))
(18)
(x) dt
d0(t) dt
V Vofuce)
(19)
27

---

where (t) and (2) are the network parameters and network output at time t.
Under Assumption C.1
have that
(x())(x) = fa(o [()](0) 8(0)),
which gives
Velecn J(x). Suppose the loss function is an MSE loss, i.e. ((U.9) = (ÿ -y) Then (n) (z).g) 2(f((z) (a) y)
Therefore, under Assumption
(18 becomes
d0(t)
[xx) [ 2A, (0(t) e(0))
-- 2 [xx) +())) 0(0)) y) + 22()
-2n (z)) (fo()(x) J(z)0(0)) y) 00) + (=) J(x) + Anx 0(t)
7
Note that / [ (Joro) (ar) - )()) y) X, (0) and (a)( T [(x) + An are both inde-
pendent of 1. Solving this ordinary differential equation, obtain
0(t) =0()) [x=)
(=) ((x) + An1)
(z) - yy (1 (1-exp exp (-2m (-)(())T( ((x)) )(())+>n1))) )nd)))
Hence the network output at time i becomes
(x) (x)
-J(x)(0(1) - 0(0))
=- --=3(=))(a))
[=) J(xx) + An1)
(m
I (I-exp exp (-2ml (_)(Z)) [(x) )n1) xnn))) (z)(0)(2)) - y)
K(x,x) (xx(xx,))+)) + n
(227p (1 (I I (x-xxp/
K(x,z) + Xn ))) (f)r0)()) - y))
n2
where the last inequality used the notation K (xx, æ) = J (x) and K K x) (x)) Therefore, with sufficient time of training (t oo), the final trained network is
(z) lim (z)
=e(oo (x:) 1 K(r,z)) x) K(x, Xn11 (z) - v)
(z) K(x, x) (K xz) + I) -1 (y (x))
soo (r) + K (z, x (xx [x, x) + I) (g Sas (x))
where in the last inequality, in the main paper. Summarizing the above discussion, we conclude that
Sge represent the network initialized with the parameter 0 as
Assumption C.2. Suppose that the network training is specified as follows:
The network adopts the NTK parametrization and parameters are randomly initialized using He initialization.
28

---

2. The network is sufficiently (infinitely) wide so that the linearized neural network assumption
(Assumption C.I holds.
3. The network is trained using the loss function in ① via continuous-time gradient flow by feeding the entire training data and using sufficient training time (t + oo).
Proposition C.3. Suppose that Assumption C.2 holds. Then the final trained network given by
hn (2) Sgb (z)
(xx, x) (K (x, x)) + Xn) (y
Sgr (x))
(20)
For the rest part of Proposition 3.1 please refer to Lemma
C.4 below.
Shifted kernel ridge regression. On the other hand, we build shifted kernel ridge regression based on the NTK as follows.
Let H. be the reproducing kernel Hilbert space constructed from the kernel function K (æx x²) (the population NTK); See Appendix B Consider the following kernel ridge regression problem on the RKHS H:
1 Σ i=1
8gb arg min
sgi (2.) g(x;))² 2 + An 11g 1 (y
(21)
gEH
where An is a regularization hyper- parameere and Sgr is a known given function. Lemma C.4. 20 is the solution to the following kernel ridge regression problem
Sgh + arg min
(ye (x;) g(x:))² Alll
(22)
gE
1=1
Proof of Lemma C.4
We first let dJ
8p0 (1x)) be the shifted label. Now consider the kernel
ridge regression problem
1 n i=1
arg min
Σ (ū: g(xx))²= An l9
gER
Standard theory in kernel ridge regression (Proposition B.I shows that the closed-form solution to the above problem is given by:
K (1r, (K a)) An K(x x>) (K (a, x)) + X,, sgo (ax)
as K is the kernel function of H. This equation corresponds to the definition of in 800 (x) in 20
From Proposition
C.3 and Lemma C.4 we immediately conclude that
The shifted kernel ridge regressor using NTK with shift from initial function 8gs exactly the linearized neural network regressor that starts from the initial network sga.
Appendix D Additional Remarks
We make the following additional remarks explain some of the details in the main paper:
1. Regarding Proposition 3.1 Proposition| C.3 (.3) in Section Proposition 3.1 is the theoretical foundation of this paper develop our framework for efficient uncertainty quantification and reduction for neural networks. We recognize the limitation of this proposition: The lincarized neural network assumption (Assumption C.I must hold guarantee that the exact equality (2) holds in Proposition 3.1 In reality, the network can never be infinitely wide, and thus an additional error will appear in (2). Unfortunately, with finite network width, this error might be roughly estimated up to certain order but still involves some unknown constants that hide beneath the data and network [4 80), which extremely difficult to quantify precisely in practice. This work does not deal with such an eтror from finite network width. Instead, we assume that the linearized neural network assumption readily holds as the starting point for developing our uncertainty quantification framework. Equivalently, one may view our work as uncertainty quantification for linearized neural networks.
29

---

2. Regarding the regularization parameters in Section We introduce the regularization parameters A, with the regularization term centering at the randomly initialized network
parameter in T throughout this paper. This form of regularized loss is also adopted in
[119 [17.. has some advantages: 1) common assumption in NTK theory is to assume that the smallest eigenvalue of the NTK Gram matrix is bounded away from zero [4 7 80] With An > 0 introduced, there is an additional term nd in the NTK Gram matrix so that this assumption is always valid naturally. 2) It helps to stabilize the inversion of the NTK Gram matrix any computation is needed, 3) It can naturally linked to the well-known kernel ridge regression, which we leverage in this work develop uncertainty quantification approaches.
On the other hand, we note that with Ao = lim An > 0, the "best" predictor h² in (6)
given in our framework is not exactly the same as the ground-truth regression function g
E (xx,Y)> (Y| = (OK + Aqd) x (o # gf. However, previous work [109] 106 1 shows that
x]. In fact, Proposition B.2 shows that h 3 3+ (Oxxxx + to (0)
with some conditions, h - 9 + 0 in some metrics when A - O. Therefore, we may use very small limit value An in practice to make the error h - negligible,
did in
experiments.
3. Regarding naive resampling approaches such as "bootstrap on a deep ensemble" mentioned in Section 3.1 We show that bootstrap on the deep ensemble estimator hW (x) face computational challenges in constructing confidence intervals. Since hn (x) involves two randomnesses, data variability and procedural variability, we need to bring up two asymptotic normalities to handle them. One approach is first to consider procedural vari- ability given data and then to consider data variability, as described as follows. Another approach is first consider data variability given procedure and then consider procedural variability. This will encounter similar challenges as the first approach, so we only discuss the first approach here.
We intuitively present the technical discussions to explain the computational challenges. First, conditional on the training data of size n, by the central limit theorem, we have
vm(hh (x) h (x:) & (10 (h) as m - oc where & represents "approximately equal in distribution" " and ch Varr,, (hm, or (I:) which depends on the training data. Second, we assume that the central limit theorem holds for hr (a)) (which is rigorously justified in Section 3.3): vT(h; h- (x)) = _N (²) as n - co. Moreover, we assume for simplicity that we can add up the above two convergences in distribution to obtain
Vn(h? (x) (x)))) & NN - -.) m
Note that (. (x) and hr h (x) are clearly dependent, both depending on the training data of size n. Therefore, the statistical correctness of the above summation of two convergences in distribution needs to be justified, e.g., using techniques in 448 [7] However, even we assume that all the technical parts can be addressed and the asymptotic normality of 1. (x)) is established, computation issues still arise the "bootstrap on a deep ensemble" approach: 1)" c2" is meaningful in the limit if we have the limit 2 This typically requires the ensemble size m to be large, at least in the order (1), which increases the computational cost of the deep ensemble predictor. However, this is not the only computationally demanding part. 2) The estimation of the asymptotic variance "Ç² + e²? requires additional computational effort. When using resampling techniques such as the bootstrap, we need to resample R times to obtain (hr (x))*² where E [R] This implies that need Rm network training times in total, and according to Point 1), it is equivalent to at least R x w(1) network training times, which is extremely computationally demanding. Therefore, naive use of the "bootstrap on deep ensemble" approach is barely feasible in neural networks.
Sa
Appendix E Proofs
In this section, we provide technical proofs of the results in the paper.
30

---

Proof of Theorem 3.2 In terms of estimation bias, we note that
m E (z) (z) = m Σ Blhn, (x)
Elhn (x) h (x)).
shown in This implics that the deep ensemble predictor and the single network predic- tor have the same conditional procedural mean given the data. terms of variance, note that (x)) hn,or, (z),, derived in (2, are conditionally independent given D. Therefore the deep ensemble predictor gives
Var (x)
Var(hn (2)())) _ Var (x))
while the single model prediction gives Var(hn, (x) Now using conditioning, we have
Var(h" (a:)
-ar((B(h" (x) + E(Var(A"
-Var(h( B[Var(h (c) m
(z)) + B[Vrr (x)
-ar(hhn,000 (x))
as desired.
Proof Theorem 3.3 Recall that ,= (.) is obtained by training an auxiliary neural network with
data ((11,6(1)) (oz. (x2)), (222).... (In, (xx))) and the initialization parameters 6¹. Then Proposition
3.1 immediately implies that
(K(,, a)) + (ind) so>(x)). 3 e (z) son (x) K(xx)) (K),))
Note that this auxiliary network starts from the same initialization 0 as in h (x). Then we have
on (x) = H'n.e on (I) - -((x)
Sob (x)) ={xx) K(x,, (K(K() AnmIt-((8(0) so (a:)
hn (x) (x)
by and Hence,
h;((x) = Nn (x) - (6' (x) 3(x)) = R Ao (x)) - En
as desired.
Proof of Theorem 3.5
Theorem 3.3
implies that
h (I) (x) - hi, (x) (x)
which is the solution to the following problem
( = x(-) + min 9ER m
> (ws 8(xx)) - g(xx))) + Ao A011
Therefore, its corresponding population risk minimization problem is:
h (-) = 5(-) + min [Y s(X) g(x)) + Ao 11
sEH
Applying Proposition 3.4 to ig and h" we have that
(())()))) s(x)) (li 8(x))) + ~0(000(x))
31

---

In other words,
vn (on,o0())) (x) ~((0,
This shows that asymptotically as n - 0o we have
91 h (x) )x (x)
S 9- + +1
where ya is the a-quantile standard Gaussian distribution (V, Hence
Tinoor (x) - Φn,en (x))
o(x) vn
q Mi, (x) Sn (x) +
( vn. 11
is an asymptotically exact (1 a--leeel confidence interval of (x))
Proof of Theorem 3.6 Note that the statistics
Re R (x)
Jú ne, (a:)
from the j-th batch i.i.d. for j E [m']. Moreover, by Theorem 3.5
we have the asymptotic
normality:
Vnl (o Sn o0 (I) hz(z) NO 2 (z)),
meaning n' n/ml + oc. Therefore by the property of Gaussian distribution and the principle of batching, we have
Vn'
m
(bB(z) h² hx))) <(x)
Zi
m
and
vn(x) (x)
m
Zi
771
SB(z)/Vml
2
172
(Zj
7,)2
(m/-111
(1
j==
for i...d. (0, 1) random variables Zi...., ..., Zm', where we use the continuous mapping theorem to deduce the weak convergence. Note that
197 1 23
tm
V m/ (m/ j=1 (Z
5m
2,)
ml 1==1
Hence, asymptotically as > Do have
4B(() (x)
IP --
SB(x)/VmL 5Q
+ +1
where is the a-quantile of the distribution 11 with degree of freedom m' - 1. Hence
Sp(z)
SB(x) 41
()(())-
Vm 91- (b(())
Vm
asymptotically exact (1 levlevel confidence interval of h" hx:)..
Proof of Theorem 3.7 Note that for j E [R... the statistics
i (x) 3z (Ir)
from the j-th bootstrap replication are i.i.d. conditional on the dataset D. By Theorem 3.5 we have
the asymptotic normality:
/ (2)) - (x) - -1(())) = ~N(0,0²(x))
32

---

as 7 + oc. By Theorem 2 in [24]. we have the following asymptotic normality: For j E [R]
hhd Vn p, (a)
3 n, (x) (hin.e - on or (2)) (z)),
conditional on the training data D. Therefore Assumption 1 in 75 is satisfies and thus Theorem [7] holds, showing that asymptotically as n - co, we have
vc(x)) - h (z)
S q--
+ +1
Sc(r)
where is the a-quantile the distribution In with degree of freedom R. Hence
( (c(a) (0( + (c())q
an asymptotically exact (1 levee confidence interval of h" (z))
Proof of Theorem B.6 Recall that & is the feature map associated with the kernel & of the RKHS H. We apply the influence function formula in [25] (see also [31 50) to obtain the infinitesimal jackknife:
(F(2:T,m,) = (2 + (' Dr (2x))) ((zx)
where Sp : 7 - 76 is defined by
Sp(f) = 2f + (y (X))
This can be seen from [4) by setting P n. and using the closed-form solution of the kernel ridge regression.
To compute the exact formula, we need to obtain Sx, and S Since the loss function L(0.)) v)
(j y) we have
2 n
Sm. (f) = 220l + (
(x)))))_)2x))))((X)))))))))))))) ( 2No
n j=1 (x,))((x))
Σ
Suppose ST.1 (2>08t.: 91. Then zo, we have
220 2o (50) Sr. (ō) (zo)) 2dom1 (5o) + Σ
(1,))k(
==1
Hence,
1 n
(1())) Io (Fo)
Σ
őn (2y (k(o,, 2.).
don
This implies that we need to evaluate ğ (Iz) on training data first, which is straightforward by letting To Fi Fn.
1 Agn
I1 (x) = ST..do (x))
k(a, an)r (m)
so
(xx) = (k(x, x) + Agnl) -1 Dm
and
ğn (xa) = I5nnAa (21))
k(oo, ar) T (x))
don
= In (0 (50) k( xr)" (k(a x.) + 2on1) -1 "annoo (x)
Next we compute
S-10(2,) = 92. At Zn, we have
k(Zz, (zx) -r. (ā2(T0)) 2A0ā2(Ho)
Ta (z) )k(20 xz))
33

---

Hence
(2(20)
k(zz, 20)
Σ (1x0,
2Ao
Xon
This implies that need to evaluate (2(Zz) on training data first, which is straightforward by letting To = 21, Zu.
1
- k( xt)
k(æ, x))2(()
2No
Ann
so
(2(x)) = x) + Apnd) (5/1) k k(zrs x))
and
1
1
(2(0))
k(2z, 20)
k( (2(x)
2A0
Apn
1
k(Zz; 20)
k x) (k(xxx) + Aorl) "k(2, x)
2Ao
2Ao
Combing previous results, we obtain
[F(;; (2o)
=- Imn (oo) + k(o0, (k((,, æ) + >oom))
In do do (x)
+ 2(22 (2x))
k(22,
k(10, x)² (k(x,, + Aonl) klzz, x))
2Ao
2λo
=k(zo, )²² (k(xx xx) + SonI)
Iota (x)
(Zy Irrnndo x))
Ao
1
- -an (10) + - (2n - Drndoo ( (2z)) 20)
-k(sp, x) (k((x, \pml) ( (z
as desired.
Proof of Theorem B.7 Recall from Section B.2 we use the equivalence notations:
TF(2;T, Tp PP)
First, we prove the following Claim:
sup
| -T) T", (6= Tn) + 0, a.s.
(23)
=EX'xy'
From Lemma A.6. and Lemma A.7. in 50) we have
(T(() = --S7'(WP(Q))
In this equation,
Wp: Di:(inss)) - H, Q (PPA)
a continuous linear operator lin Bg) where the space lin (B)) is defined in [50] Moreover, the operator norm of Wp satisfies ||Wr| < 1. In addition,
H. - H. f + 2hof + Ep (X))x
an invertible continuous linear operator on where the feature map of H. This implies that Sp' 71 - H is also continuous linear operator H.
Next, we apply Lemma A.7 [50] by considering the following two sequences: The first sequence is 6_ n, E lin(Bs) which obviously satisfies
( π.,) - (6_ - T)ll = || kno + O, a.s.
The second sequence nn E Bg which obviously satisfies
π knlloo 0,
34

---

Following the proof of Lemma A.7 in [50], we have that
[T - T- (6. - Tn)l
(W,(2 - =)) - ST. (W nn (6g 1 aT)))
<||S; (0, m)) S (W nn (dg - m + || (6g T)) - '5'(War. (62 - x
'
+ +115,-(W.,, (62 mn)) 5T. SW.,,, (6= xn))ll
<||s="|||W,,0_ x) WT. (6= = h + 1-1s= [[= π) (6= TTa)) loo
+ +| -- W |11|, nn lloe
S|IS; || WT, (T)llæ + ||s- [|| T. ||||| ||| 2o (since WT, is linear operator)
+
( - =|/0o + || - -| 1lo)
Therefore, only need study the three terms in the last equation. Note that the first two terms independent of 2, and thus it follows from Step and Step 4 in the proof of Lemma A.7 in (50) that
sup =EXxy
11| W., (T) + ||S -Il I|l|
Tullo
= -|SS [ W.., (T) + || II[ an |||| - xnlloo - 0, a.s.
To show that the third term also satisfies that
||S_ - -- - [., (115, a110o + || - xnoloo +0 0, a.s., =EX'X''
sup
we only need to note the following fact:
1) ||S - + O, a.s., by Step 2 in the proof of Lemma A.7 in [50]. Moreover, this equation is
independent of 2.
2) ||W,, n, by Step in the proof of Lemma A.7 in [50]. Moreover, this equation is independent of z.
3) || oo +00, a.s., since || Talloo a.s.. Moreover, this equation is independent of 2.
4) sup:Exxy ||0g 11100 +00. To see this, note that by definition
sup || 1 sup sup =EX'Xy'
/ gd(6. m)
sup sup g(z) / gdv S sup l|gll xEX geg
=EXXX gEg
Deg
where the last term is independent of and the space supgeg, o/ < by the definition of Gi. supge(t) llal < +oo by our additional assumption on
91 U S2 (b} defined in (50]
Since both X' and y are bounded and closed, suppex / /k(r, x) is bounded above by, say x < co.
Thus, for every h E with ||
Ci, we have | 60 S C k. By definition of G2. for any g G2. we can write g(z, y) = L' (y, fo(x)), with |1/o æ S co and [ S 1. Note that [ S 1
implies that || I k and [| S Co implies that [| S cps. Thus Assumption B.4 shows that ||C' (g. fu(x))]l S b/enn (y) which is bounded on y by our additional assumption (uniformly for fo). Hence we conclude that sup lI S k.su, Von (v) < +00. Summarizing the above discussion, obtain that
sup ||0 - -o S 2sup Ilol < 400 =EX'X
SEG
Combining the above points 1)-4), we conclude that
sup xEX'x
||
- - IT (116- T/1o0 + || Tulloo) + O, a.s.
Hence, we obtain our Claim 23
Applying Kno to the reproducing property of l implies that
sup
[ T)(z)) T"d (8. n. (2o)] + O, a.s.
xEXxy
35

---

Note that since X' and y are bounded and closed by our assumptions, we have that
Vexw (T= - (xo)) dm (2) - 2xEXxy (T* (62 mn)()))) dr(=))
<| x x x sup
[6= - (2)))) -T²,, (6g - xx)(
=EXXY
x 1 Vrexx (T= - (xro)) + TT,, (6= -n)((o)) dr(())]
-0 0.8.
where use the fact that Sxe SEe (T (xo)) dr(z)= =22 +00. On the other hand, it follows from the strong law of large numbers that
1 Σ (6, (6,, - (xn)) - JEEKYY (Ta, (62 - -m)((2o))))) dr(=)) + 0, 0.8. 71 2EDD
Hence, we conclude that
T Nn (6z, (T) 1 JSEXXY (T((5=-))<zo) dx(z) 0, a.s.
Σ
7
ADD
In other words,
((x0) -- 1 Σ (F
f J
IF - = (=)
a.s.
ni 2EE
For confidence intervals, the proof straightforward. Theorem B.5 shows that n + oo we have
9 A. (zo) In Ao (x0) = (0, {(xa)/ /n
By Slutsky's theorem and {(zo) E(ro)
1, a.8., have
yn., A., (20) Iv (200 = 10,
5(xo)\ vn
This implies that asymptotically
00, we have
IT (a (50) (0
P
S 91
+ a
{(00)/ v
where da is the a-quantile of standard Gaussian distribution N'0,,) 1). Hence
Innidn (zo) É(zo) Vn
n1 -1 IT,,A (2o) + É(xO)
/n n1 +
an asymptotically exact (1 al-level confidence interval of Ar (I)).
Appendix F
Experiments: Details and More Results
F.1 Experimental Details
We provide more details about our experimental implementation in Section
Throughout our experiments, we use a two-layer fully-connected neural network as the base predictor based on the NTK specifications in Section C (Proposition C.3) However, we need resolve the conflicts between the theoretical assumptions therein (e.g., continuous- gradient flow) and practical implementation (e.g., discrete-time gradient descent), and at the same time, guarantee that the training procedure indeed operates in the NTK regime so that the first-order Taylor expansion (linearized neural network assumption) works well following specifications in all experiments.
Therefore, we will use the
36

---

1. The network adopts the NTK parametrization, and its parameters are randomly initialized using He initialization. The ReLU activation function is used in the network.
2. The network has 32 x n hidden neurons in its hidden layer where n is the size of the entire training data. The network should be sufficiently wide so that the NTK theory holds.
3. The network is trained using the regularized square loss ( with regularization hyper- parameter ', = 0.1
4. The network is trained using the (full) batch gradient descent (by feeding the whole dataset).
5. The learning rate and training epochs are properly tuned based on the specific dataset. The epochs should not be too small since the training needs to converge to good solution, but the learning rate should also not be too large because we need to stipulate that the training procedure indeed operates in the NTK regime (an area around the initialization). Note that in practice. we cannot use the continuous- time gradient flow, and the network can never be infinitely wide. Therefore, with the fixed learning rate in gradient descent, we do not greatly increase the number of epochs so that the training procedure will likely find a solution that is not too far from the initialization.
6. We set m 4 in the PNC-enhanced batching approach and R = 4 in the PNC-enhanced cheap bootstrap approach.
7. In DropoutUQ, the dropout rate is a crucial hyper-parameter. We find that the dropout rate has significant impact on the interval width of DropoutUQ: large dropout rate produces a wide interval since the dropout rate is linked to the variance of the prior Gaussian distribution 443. Therefore, to make fair comparisons between different approaches, we adjust the dropout rate in DropoutU so that they have a similar or larger interval width as PNC-enhanced batching or PNC-enhanced cheap bootstrap.
8. All experiments are conducted on a single GeForce RTX 2080 Ti GPU.
F22 Additional Experiments
In this section, we present additional experimental results on more datasets. These experimental results further support our observations and claims in Section 4 demonstrating the robustness and effectiveness of our proposals.
First, we consider additional synthetic datasets.
Synthetic Datasets #2: X ~ Unif( 0.2)d) and Y ~ ΣL xin (i)((()))))))) + WO0 0.001 The training set D = {x 3:) : I,....n] n is formed by i.i.d. samples of (X,Y) with sample size n from the above data generating process. We use 2i = 00.1,0.1,... 0.1) and yo = Σi= 0.1 s0n(0.1) the fixed test point in the confidence interval task. The rest of the setup is the same Section 4 The implementation specifications are the same as in Section [.I] The results for constructing confidence intervals are displayed in Table 3 and the results for reducing procedural variability are displayed in Table
Next, we consider real-world datasets.
is worth mentioning the main reason for adopting synthetic datasets in our experiments. That is, precise evaluation of confidence interval requires the following two critical components, which nevertheless are not available real-world data: 1) One needs to know the ground-truth regression function to check the coverage of confidence intervals, while the label in real-world data typically contains some aleatoric noise. 2) To estimate the coverage rate of the confidence interval, we need to repeat the experiments multiple times by regenerating new independent datasets from the same data distribution. In practice, we cannot regenerate new real-world datasets.
However, to conduct simulative experiments on realistic real-world data, we provide possible way to mimic an evaluation of the confidence interval on real-world datasets as follows.
Step 1. For a given real-world dataset, select a fixed test data point (zo. yo) and a subset of the real-world data (ri, y:) (i E Z) that does not include (zo, yo)-
Step For each experimental repetition j E [J], add certain artificial independent noise on the label to obtain "simulated" real-world dataset (z4, Yt + es,z) (i I) where ei are all independent
37

---

Table 3: Confidence interval construction on synthetic datasets #2 with different data sizes = 128, 256, 512, 1024 and different data dimensions d 2, 4,8, 16. The CR results that attain the desired confidence level 55% are in bold.
PNC-eahanced batching
PNC-enhanced cheap bootstrap
DropoutUO
95%CRCCR/W) SP%CIICR/W)
MP
95SCKCR/WW) 90RCRCR/W)
MP
95%CCICC/W WSCCCCRWW)
MP
[d = 2)
0098/0.0075
0.91/0.0056
0.0201
1,000,0333
0.99/0.024
0.0184
0,990.03392
0.97/0.0330
0.0216
256
0.97/0.0058
0.930.00443
0.0197
0.97/0.0215
0,93/.0165
0.0185
0,96/0.0213
0.94/0.0169
0.0214
n = 512
0,96/0.0050
0.95/0.0037
0.0199
0.9400.0148
0,90/0.0114
0.0186
0.96/0.0164
0,900000135
0.0216
1024
0.96/0.0041
0,94/0.0030
0.0199
0,950.00155
0,93000008
0,0203
0,91/0.0124
0,86/0,0102
0.0210
[d = 4)
0
0,97/0.0022
0.0395
1,000.0443
0,98/0.0340
0.0370
0,9600.062
0.93/U.3375
0.0441
256
0.98/0.0088
0.92/0.0065
0.0401
0.96/0.0366
0.94/.02666
0.0374
0,98/0.0344
0,94/0.0282
0.0426
512
0.94/0.0073 0073
0.91/0.0054
0.0403
0.99/0.0236
0.94/0.0181
0.0373
0,94/0.0310
0,88/0.0260
0.0427
1024
10000.0055
0,91/0.00-0
0.0399
0,94/0.0167
0.87/0.0128 28
0.0373
0,95/0.0241
0.90/900203
0.0417
[d 8)
128
0.97/0.0145
0,930.0107
0.0798
1,000.0683
L00/0.0524
0.0849
0,94/0.0661
0.86/0.0548
0.0876
256
0.96/0.01 1O
0,94/0.0082
0.0801
0.99/0.0400
0.97/1.041
0.0846
0.92/0.0538
0.91/0.0446
0.0853
= 512
0,97/0.0092
0.92/0.0006 000.00668
0.0793
0.99/0.0393
0,96/0.0302
0.0846
0.91/0.0427
0.85/0.0359
0.0862
1024
0,95/0.0084
0,90/0.0022
0,0792
0,9770.084
0,95/0.0288
0,0820
0,85/0.0295
0.79/0.0248
0.0809
[d I6)
128
0.95/0.0214
0,94/0.0158
0.1608
1,000.
0.99/1.771
0.1568
0.92/0.1124
0.87/0.0958
0.1814
256
0,97/0.0174
0.91/0.0129
0,1595
1.000.07255
100/0.0556
0,1578
0,95/0.0800
0.88/0.0677
0.12
512
0,96/0.0140 40
0.92/0.0104
0.1604
0.99/0.0516
0,98/0.0396
0.1569
0.88/0.0565
0,80N0.0475
0.1689
1024
0,95/0.0114
0.88/0.0084
0,1599
0,9800.0350
0.95/U.0269
0,1568
0.86/0.0584
0.83/0.0498
0.1685
Table 4: Reducing procedural variability to improve prediction on synthetic datasets #2 with different data sizes 128, 256, 512, 1024 and different data dimensions d = 2, 8, 16. The best MSE results are in bold.
MSE
One network
PNC predictor
Deep ensembll networks)
Deep ensemble (2 networks)
2)
128
(1.02 + 1,16) 10 (1.38 + 0.18) x 10
(2.90 1.72)
10
(5.39 + 3.01) 10
256
(7.36
+ 4.84) 10-5
(1.28 0.100 x 10-
(2.39
1.17)
10-3
(3.67 # 1.61) 10-5
512
(4.41
2.34)
x 10
i 0.10) x 10- (1.11
(1.79 1 5.19)
10
(2.52 1 1.55) x 10-
n 1024
(3.46
1.91)
X 10
(9.78 + 0.94) x 10-
(1.03 + 0.39)
x 10-
(1.82 + 1.27) X 10-5
[d 4)
(1.89 1 1.23) 10
(1.98 + 0.25) x 10
(4.64 1 1.49) 10
(1.02 1 0.31) 10"4
256
(1.79
10
0.08) x 10- + (1.49
(455 1.58)
(1.02 0.37) 10-4
512
(1.72
(1.29
+ 0.05) x 10"
(3.73 0.88
(7.78 1.72)
1 1024
(9.75
10-3
(1.05 + 0.03) x 10-
(3.70 0.79)
10
(8.92 # 2.09) 10-3
(d 8)
128
(8.00 + 1.36) x 10 (5.07 + 0.65) x 10-3
(2.18 0.59) x 10-4
(4.77 : 1.32) X 10-4
n 256
(4.50
+ 0.30) x 10 (2.30
+ 0.24) x 10
(9.44 1.44)
x 10 6
(2.08 0.19) x 10
a 512
(3.80
1 0.35)
X 10 (1.47
: 0.07) x 10
(8.44 1 1.11)
10
(2.02 1 0.23) X 10-4
n 1024
(3.46
1 0.50) x 10
(.12:0004)
x 10
(8.11 1 1.27)
10
(1.71 1 0.23) x 10
I6)
128
(2.32 + 0.24) 10
(1.39
10
(5.71 0.74)
(1.36 # 0.30) 10
256
(1.86
1 0.21)
10-3
(6.55 0.75) x 10
(3.96 0.34)
10
(9.49 1 0.67) 10
512
(1.42
: 0.08)
10-3
x 10-3
(3.00 0.16)
10
(7.14 # 0.42)
10-4
1024
(1.18
+ 0.06)
10
(2.58 1 0.24)
x 10
(2.56
181 10
(6.30 + 0.39)
10
Gaussian random variables. Construct a confidence interval based on this "regenerating" training data (2i, yr + ei.g) (i E Z) and evaluate it on (zo; 30)-
In the above setting, Yo is treated as the true mean response of zo without aleatoric noise, and e represents the only data variability. As discussed above, precise evaluation is impossible for real- world data; the above procedure, although not precise, is the best we can do to provide a resembling evaluation. Using this setting, we conduct experiments on real-world benchmark regression datasets Boston, Concrete, and Energy. These results are shown in Table
from UCI datasets:
From the results, we observe that our approaches also work well for these "simulated" real-world datasets. Under various problem settings, our proposed approach can robustly provide accurate confidence intervals that satisfy the coverage requirement of the confidence level. In contrast, DropoutUQ does not have such statistical guarantees. Overall, our approaches not only work for synthetic datasets but also is scalable to be applied potentially in benchmark real-world datasets.
38

---

Finally, we conduct experiments of reducing procedural variability to show the performance improve- ment capability of proposed models and report MSE results on the "simulated" real-world datasets: Boston, Concrete, and Energy. In each experiment, we use an 80%/20% random split for training data and test data in the original dataset, and then report MSE on its corresponding "simulated" dataset from our proposed PNC predictor and other baseline approaches with 10 experimental repetition times, These results are shown in G
Table 5: Confidence interval construction on "simulated" real-world benchmark datasets: Boston, Concrete, and Energy. J 40. The CR results that attain the desired confidence level 95% are in bold.
PNC-onhanced batching
PNC-onhanced cheap bootstrap
DropoutUQ
95%CKCR/WW)
90RCI(CR/W)
MP
95%CCI(CR/W) 90%CIICR/W)
MP
95%CICCR/WW) 90%CIICR/WW)
MP
Boston
Test
L0/1172
0,975/0.866
0.486
0.925/1.024
0.875/0.786
0,85/09997
0.7/0.832
0.109
Test 2
1080.873
0.95/3.645
0.100
0.975X0.869
0.925/0.667
0.042
0.75X1.981
0.625/0.42
0.089
Test
1,00.896
1.0006622
-0.047
0.975/1.473
0.85/1.131
0.0875
0.725/1.00 .000
0.575/0.866
-0.119
Test
L0N0.915
1.0/6666
0.219
0,925A0.829
0,875/U.637
0,372
0.90/1.164
0.825/09981
0.131
Text 5
LO00.587
1.0/0.344
0.145
1.00.531
0,975/0.007
0.272
0,75/1.093
0.725/0.920
0.454
Concrete
Test 1
L0/1.778
100/10119
0.416
0.975/1.013
0.975/0.778
0.379
0.725/1. 176
0.675/1.016
0.088
Test
1.01.1171
LON0667
0.362
0.675/1.1 112
0.65/0.952
Test 3
1001,673
1.0/1.237
0,0455
1,/0,970
0.97500.745
0,516
0.875/1 .061
0,75A18887
0.597
Test 4
1501
1.0/0855
0.484
L.00,960
0,95/0.737
0.466
0.875/1. 101
0.825/0.949
0.552
Tost
10/1.208
1.00.8933
0.216
L.00.601
LOV0.462
0.309
0.825/1 000
0,8/0,800
0.684
Energy
Test
0,975/0.566
0.975/0
LOV0.458
0.152
0.75/0.854
0.725/0.733
0.037
Test 2
1,0/0.422
1.0/0.112
0.294
0,95/0,570
0.875/0.437
0.191
0,8010.611
0.725/0.504
0.231
Tea
1,000.435
1.00.3211
.0.734
L.00.724
LO00.556
-0.766
0.75X0. 749
0.675/0.630
-0.908
Tost 4
1.000.572
1,0/0,433
-0.723
1,006613
10/0,470
-0.792
0.82540.651
0.775/0.549
-0.847
Test
0,975/0,512
0.975/0.379
-0.679
0,975/0.550
0.975/0.422
-0.702
0.825/0.644
0.80/0.550
-0.807
Table Reducing procedural variability to improve prediction on "simulated" real-world benchmark datasets: Boston, Concrete, and Energy. The best MSE results are in bold.
MSE
One
PNC predictor
Deep ensemble networks)
Deep ensemble networks)
Boston
Test
(3.52 +: 0.78) x 10
(1.04 1 0.03) x 10
(1.56 + 0.28)
10
(2.24 + 0.51) x 10-1
Text
(3.15 +: 0.59) x 10
(1.01 + 0.04) 0.044 x 10-1
(1.38 0.26)
10-1
(2.03 : 0.47) x 10-1
Test
(4.18 +: 0.69) x 10-1
(1.93 + 0.05) x 10-1
12.88 0.33) 10
(2.68 # 0.37) x 10
Teat
(3.33
0.81) x 10
(1.17 147700.00 0.06) x 10-1
(1.18 x 10
(1.84 0.24) x 10
Test
(3.46 1.46:: x 10
x 10-1
(1.99 (1.90:0..1) x 10"1
(2.45 1 0.32) x 10
Concrete
Test
(2.47 10
: x 10
(2.02
10-1
(2.21 10
Test 2
(3.00 0.29) 10
(2.18 x 10-1
(2.2811.15)
10
(2.43 10
3est
(2.39 0.16)
1011
(1.85 x 10-1
(2.01 10-1
(2.18 10
Text
(2.85 :: 0.53) 10
(1.44
x 10
(1.92 10
(2.17 0.38) 10
5
(2.51 0
10
(1.75 10
(1.86 0.06) 10-1
(2.16 10
Energy
Test
(1.09 :.14) 10
(7.33:0.08) x 10-2
(8.56
10-2
(9.91 1 0.86) X 10
T
(1.08 +0100) 10
(7.46
x 10
(8.29
(9.57 = 0.84 10
Test
(1.23 (1.00)
10
X (0.12)
(8.45
(1.01 1 0.12) 10
Test
(1.25 (1.25:0.12)) 10 (7.24
x 10
(..31
10
(L07:0..1) 10
5
(1.12 0.09)
(6.81
1 0.23) x 10
07.47
(9.69 0.37)
F.3 Confidence Intervals for Coverage Results
When evaluating the performance of confidence intervals, the CR value is computed based finite number of repetitions, which will incur a binomial error on the (population) coverage estimate. Therefore, in addition to reporting CR point estimate, we also report confidence interval for the (population) coverage. Table repporss the Clopper-Pearson exact binomial confidence interval. Note that this interval can be computed straightforwardly by the point estimate (CR) and the number of repetitions (J). Therefore, we can refer to Table 7 to additionally obtain a confidence interval of coverage for results, e.g., in Tables aadd]
39

---

Table 7: Clopper-Pearson exact binomial confidence interval for coverage estimates at confidence level 95%. J 100.
Confidence interval coverage
coverage
coverage
1.00
[0.964, 1.000
0.99
[0.946, 1.000)
0.98
00.330, 0.998
0.97
00.915, 0.984
0.96
0.901. 0.989
0.95
0.887, 0.984
00.874,..33 0.978
0.93
0.861, 0.971
0.92
0.848, 0.960
0.91
[0.836, 0.958
0.90
008211, 0.951
0.89
00.812 0.911
[0.800,
0.7766, 0.921
0.85
(0.765, 0.914
0.84
0.755, 0.906
0.83
0.742, 0.898
0.890
0..719, 0.882]
0.708, 0.873
40


Reviews:
[{"review_id": "mWmqEUl4NC", "reviewer": "Reviewer_VEbh", "summary": "In the paper, authors present a new approach to quantify and mitigate a specific aspect of epistemic uncertainty in model predictions. They identify and quantify \"procedural variability,\" a type of epistemic uncertainty that arises from noise in the training process. Based on the Neural Tangent Kernel (NTK) theory, the authors introduce a \"procedure-noise-correction\" method. They offer theoretical support and experimental evidence to justify their approach.", "questions": "Can the method be somehow adapted for classification?\n\n\n\nEdit: I would like to thank the authors for their answers during rebuttal period.", "limitations": "As I can see, this paper has no potential negative societal impact.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Generally, the majority of research papers propose methods that quantify the overall epistemic uncertainty, but a few discuss the estimation and evaluation of its different components. The process of splitting epistemic uncertainty into different parts, while not a novel concept, remains significantly under-studied. Therefore, the problem addressed by this paper is promising and brings potential value for further studies.\n\nThe strengths of this paper are:\n\n1) The authors built a mathematical framework for their proposed method.\n2) The authors' approach to quantify and remove procedural variability is indeed novel and opens new perspectives for future works.\n3) Combined with the PNC-Enhanced batching, the presented approach seems practical, and practitioners in the field may find it beneficial.\n4) The supplementary material is verbose and resembles a textbook, providing a comprehensive guide to the details of derivations, proofs, and certain theoretical aspects like the NTK.", "weaknesses": "There are some aspects of this paper I consider as weaknesses:\n\n1. Narrative of the Paper: Despite its overall good organization and clarity, Section 3 in the main document is challenging to understand. Reorganizing this section could enhance the paper's readability. Some propositions could be moved to the Appendix, freeing up space for the experimental section (see Weakness 2). Adding illustrative elements, such as a schematic diagram of the dataset construction and overall architecture, could further improve the readers' understanding.\n2. Experiments: The set of experiments presented appear overly simplified and unrepresentative. Even though the multi-layer perceptron (MLP) models with two layers and varying unit numbers technically qualify as over-parameterized for d=1,2,4,8, the example is overly simplistic. The only experiment is with a sinusoidal signal form, and a y-label variance of 0.001 isn't explained or justified. The paper could be enriched by the inclusion of generally accepted benchmarks, such as regression benchmarks containing UCI datasets[1]. Moreover, it's clearly seen from your results that an ensemble of 5 models substantially outperforms an ensemble of 2 models, and training these ensembles isn't overly computationally demanding considering the complexity of data and models. Hence, comparing with larger ensembles, say 50 or 100 members, could provide more insightful results.\n3. The initial part of the paper seems to promote the enhancement of **uncertainty quantification** by disentangling epistemic uncertainty into different components, but the actual evaluation of this uncertainty is lacking. The experiment appears to improve the **optimization process** rather than the quantification of uncertainty. An experimental suggestion would be to select a measure for epistemic uncertainty and apply it to both in-distribution (validation split) and out-of-distribution data. Computing the uncertainty measure for each data point across independently trained models several times and comparing the resulting ROC-AUC with the one obtained through your method could demonstrate whether the elimination of procedural variability truly improves the estimates of epistemic uncertainty.\n4. There is no code available, so I can not check the reproduction of the experiments.\n\n[1] Gal Y., Ghahramani Z. Dropout as a bayesian approximation: Representing model uncertainty in deep learning //international conference on machine learning. – PMLR, 2016. – С. 1050-1059.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We greatly appreciate the reviewer for the careful reading and valuable feedback. We address the reviewer's concerns below.\n\nQ1. (Experiments) For additional experiments on real-world datasets, please refer to our global response to all reviewers and our discussions there. We hope this would alleviate your concern that our data set is overly simplistic (and not motivated as well). For larger ensembles, the emphasis of our paper is that our PNC-predictor only requires 2 networks training. Therefore, a fair comparison with the same running time should be 2 networks in the deep ensemble approach, and our approach clearly outperforms it. Deep ensemble with any networks of more than 3 (5 or 50) does not change our conclusion even if they perform similarly to ours, since they require much more running times than ours.\n\nQ2. (Evaluation) Our work focuses on the quantification of $\\textit{epistemic uncertainty}$, which refers to the errors coming from the inadequacy of the model or training data noises. This is different from $\\textit{aleatoric uncertainty}$ (the randomness in the conditional data distribution). We provided a discussion on the difference between them in Appendix A. In particular, epistemic uncertainty is typically used on in-distribution data, while detecting out-of-distribution data generally requires aleatoric uncertainty or predictive uncertainty that is different from our study scope, even though this is certainly an interesting future direction for us to explore.\n\nWhen speaking of the “improve the optimization process” experiment, our understanding is that the reviewer probably means Table 2. However, we also have Table 1 for evaluation of our epistemic uncertainty quantifier. We conduct the evaluation of epistemic uncertainty statistically via the task of constructing confidence intervals, which is our task at the beginning of Section 4 (Table 1). More precisely, we aim to construct an interval that can cover the ground-truth value $h^*(x)$ with high probability with respect to the training randomness. To evaluate the confidence interval performance, we repeat our experiments 40 times to obtain 40 confidence intervals and check their coverage rate. For 95% confidence level, the coverage rate of 40 confidence intervals should be around or larger than 95%. Our performance appears competitive in this task as well.\n\nQ3. (Code.) We will definitely make our code publicly available. Also, and at the time being, Section E.1 Experimental Details provides our experimental details and configurations.\n\nQ4. (Classification) The theoretical developments needed for regression tasks and classification tasks appear different and not easily adaptable. NTK is the essential component of our theory, providing the theoretical foundation to establish the statistical properties of our algorithms. In previous literature, NTK is generally employed for mean-square-error (MSE) loss, such as \n\nDu, Simon, et al. \"Gradient descent finds global minima of deep neural networks.\" International conference on machine learning. PMLR, 2019. \n\nLee, Jaehoon, et al. \"Wide neural networks of any depth evolve as linear models under gradient descent.\" Advances in neural information processing systems 32 (2019).\n\nHe, Bobby, Balaji Lakshminarayanan, and Yee Whye Teh. \"Bayesian deep ensembles via the neural tangent kernel.\" Advances in neural information processing systems 33 (2020): 1010-1022.\n\nas well as other references for NTK in our paper. The major technical discrepancy related to NTK between the regression tasks and the classification tasks is that the softmax operation in the network final layer and the cross-entropy-type loss used in the classification does not have an explicit statistical characterization of neural network training. Therefore, most NTK-based theories consider MSE loss. It would be an interesting but different direction to study the statistical guarantee for classification neural networks, which we recognize as our future research direction. \n\n\n\n\n"}, {"review_id": "zHIJuTRF2o", "reviewer": "Reviewer_cgMr", "summary": "The authors focus on the task of uncertainty quantification for neural networks. \nThey contribute (i) a procedure to remove procedural uncertainty, an uncertainty that arises due to\nrandomness in the training procedure, and (ii) an approach to cheaply construct confidence intervals with asymptotic coverage guarantees. The proposal is evaluated on two small synthetic experiments.\n", "questions": "- The evaluation is limited to a single function in a very low-dimensional setting. Can the authors provide a further experimental evaluation with respect to the scalability of the proposed approach?\n- Wrt the paragraph around l174-l193. I am not sure I understand this. $\\bar s \\equiv 0.2$, or $\\bar s \\equiv 0$, gives us a constant target, that is independent of the input, i.e., $\\hat \\phi_{n, \\theta^b}'(\\cdot)$ should learn to be constant, giving us $\\bar \\phi_{n,\\theta^b} \\equiv 0.2$, i.e., step 2 of Algorithm 1 becomes irrelevant. Which part am I misunderstanding?\n", "limitations": "Discussion on limitations and broader impact are missing.\n", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The proposed approach offers a cheap (compared to ensemble models) and principled approach to remove procedural uncertainty and construct confidence intervals. \n- The paper is well written.\n", "weaknesses": "The main weakness of the paper is its empirical evaluation, which consists of a single simple data set, a sum of sine functions, evaluated with small data sets (up to 1024) in very low dimensional spaces (up to eight dimensions). Any kind of scalability guarantee is not given/evaluated, which is concerning, especially since the theory assumes access to the full training data (l118). \nA second synthetic function (sum of exponentials and squared covariates) is evaluated in the appendix, again without too many details. \n\n\n### Minor\n- Throughout the paper, the term _extremely_ is used repeatedly (extremely fast, extremely time-consuming, extremely few) instead of providing proper quantitative numbers (if only in the order of magnitudes, or relative to other approaches).\n- Prop 3.4 $\\pi_n$ is only defined in the appendix\n- When submitting to NeurIPS, please follow the NeurIPS style guide, i.e., place your captions above tables, no vertical lines in tables, etc.\n- The regularizing $\\lambda$ is fixed to $\\lambda \\equiv 0.1^{10}$, i.e., essentially zero? Why keep it at all?\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We greatly appreciate the reviewer for the careful reading and valuable feedback. We address the reviewer's concerns below.\n\nQ1. (“extremely”) Thanks for pointing this out. We agree that “extremely” should be clarified. In the revised version, we will avoid using “extremely” but provide more concrete numbers to elaborate our point. For instance, we will change,\n\n“Given these computational bottleneck, we will utilize light-computation resampling alternatives, including batching and the so-called cheap bootstrap method, which allows valid confidence interval construction using extremely few model retrainings.”\n\nto\n\n“Given these computational bottleneck, we will utilize light-computation resampling alternatives, including batching and the so-called cheap bootstrap method, which allows valid confidence interval construction using as few as one additional model retraining.”\n\nQ2. ($\\pi_n$) We will add this definition to the main body.\n\nQ3. (Format) We will revise the format of our tables accordingly.\n\nQ4. (Regularization) Our theory applies to any regularization greater than zero, regardless of how small it is. We introduce it in our theory so that the eigenvalue of the NTK Gram matrix will be bounded away from zero and thus the inversion of the regularized NTK Gram matrix can be computed stably. In practice, the regularization is typically very small or simply zero, e.g., in \n\nLee, Jaehoon, et al. \"Wide neural networks of any depth evolve as linear models under gradient descent.\" Advances in neural information processing systems 32 (2019).\n\nWe maintain a small but nonzero regularization in the experiment to match our theoretical results. \n\nQ5. (Experiments) We understand this concern. We have carried out additional experiments on more sophisticated problems using “simulated” real-world data. Please refer to our global response to all reviewers.\n\nQ6. (Constant label) Even if every training data has label 0, the output network is not a zero-constant network: The major difference between neural networks and classical models such as linear regression, is that neural networks are non-convex and have many local/global minima. Starting from an initialization parameter $\\theta^b$, gradient descent will find a global minimum that is very close to $\\theta^b$ (but not $0$ in general although $0$ is one of the global minima). This phenomenon has been characterized by the NTK theory or lazy training, e.g.,\n\nDu, Simon, et al. \"Gradient descent finds global minima of deep neural networks.\" International conference on machine learning. PMLR, 2019. \n\nChizat, Lenaic, Edouard Oyallon, and Francis Bach. \"On lazy training in differentiable programming.\" Advances in neural information processing systems 32 (2019).\n\nAnother way to see this is that in Proposition 3.1, if we plug in $\\textbf{y}=0$ (zero-label on all training data), we will not get a zero-constant network. Instead, the output network depends on the initialization network. This phenomenon also highlights the need to consider random initialization as an important uncertainty component in the UQ task. \n\n\n\n\n"}, {"review_id": "zQZL2wKgI0", "reviewer": "Reviewer_bPtq", "summary": "The paper focuses on uncertainty quantification, specifically on estimating the statistical range of predictions made by a deep neural network (DNN). This is achieved through the use of a novel DNN called Procedural Noise Correcting (PNC), which is capable of estimating the variability inherent in the training process, particularly the variability arising from different initializations of the DNN. In order to obtain a more accurate estimate of the idealized DNN, denoted as h^*, the authors employ the Neural Tangent Kernel (NTK) and train the PNC DNN. To construct the confidence interval for the PNC, they utilize either batching or bootstrap methods.", "questions": "Please see Weaknesses", "limitations": "The authors' paper lacks emphasis on the societal impact of their work. Typically, I wouldn't expect authors in the field of Uncertainty Quantification to delve into this aspect. However, I do feel that the claims made in the paper are overly assertive. For instance, on line 29, the paper mentions general uncertainty quantification, yet it only focuses on regression tasks. Additionally, the claim of working with Overparametrized Neural Networks seems exaggerated, considering that the DNN employed in the study consists of just two fully connected layers. It would be advisable for the authors to exercise caution when making such claims.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "In this paper, the authors introduce a remarkable deep neural network (DNN) named PNC, which effectively addresses uncertainty quantification. The presentation of the paper is well written, although there are a few missing notations. However, the supplementary materials provided by the authors greatly enhance the overall completeness of the paper. The problem tackled by the authors is intriguing, and to the best of my knowledge, I have not come across any other papers that address this particular issue using a similar approach.\n", "weaknesses": "There are several instances of missing notations in the paper, which make it challenging to comprehend. I'm curious about the authors' decision to only consider the random initialization of network parameters when discussing procedural uncertainty. Why did they not take into account random batches or the non-deterministic backpropagation? Algorithm 1 is also difficult for me to understand. If only \\bar(s) is predicted, how can h_n*(x) be obtained? Is it necessary to apply the NTK? I did not see it explicitly mentioned in algorithm 1.\n\nRegarding the experimental protocol, I believe it could benefit from a more robust approach. It might be valuable to compare the proposed method to other approaches, such as:\n\n[1] Blundell, Charles, et al. \"Weight uncertainty in neural network.\" International conference on machine learning. PMLR, 2015.\n[2] Maddox, W. J., Izmailov, P., Garipov, T., Vetrov, D. P., & Wilson, A. G. (2019). A simple baseline for bayesian uncertainty in deep learning. Advances in neural information processing systems, 32.\n\n[3] He, Bobby, Balaji Lakshminarayanan, and Yee Whye Teh. \"Bayesian deep ensembles via the neural tangent kernel.\" Advances in neural information processing systems 33 (2020): 1010-1022.\n\nAdditionally, I find it intriguing to include experiments on additional datasets such as MNIST, CIFAR-10, and CIFAR-100, not limited to regression tasks but also for classification purposes. It would provide valuable insights into how the proposed DNN performs in the presence of epistemic uncertainty. Why not evaluate its behavior by applying out-of-distribution (OOD) data?\n\nThe authors highlight that PNC involves minimal computation; however, they only apply their technique to relatively simple datasets. This raises questions about whether it truly qualifies as a low-budget technique.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We greatly appreciate the reviewer for the careful reading and valuable feedback. We address the reviewer's concerns below.\n\nQ1. (NTK and algorithms) Most of the essential notations are introduced in the statement of Proposition 3.1 and at the beginning of Section 3.1, while we defer the details related to NTK in Appendix C due to the page limit. Regarding the missing notations, we would appreciate it if the reviewer could give more guidance.\n\nNTK is the essential component of our theory, which provides the theoretical foundation to establish the statistical properties of all our algorithms, including justifying Algorithm 1. Therefore, we need the NTK theory to derive the statistical guarantee of Algorithm 1.\n\nThe major difference between our work and previous literature is that our UQ framework has simultaneous $\\textit{statistical coverage guarantee}$ and $\\textit{low computation cost}$. Some Bayesian approaches are fast-to-implement empirically but do not (or have not been shown to) enjoy the statistical coverage guarantee that ours does. The major theoretical difficulty in neural networks, compared with classical models such as linear regression, is that neural networks are non-convex and have many global minima, making standard statistical theory impossible to be applied. The recent NTK theory gives a statistical characterization of neural network training, which provides a possible foundation for theoretical UQ. With this, our theoretical framework can characterize the uncertainty from the random initialization, build connections between UQ for neural networks and UQ for kernel-based regression, and develop large-sample asymptotic theory, which takes the very first step in the UQ framework with statistical coverage guarantee, even only considering random initialization. It is more appealing to incorporate random batches or SGD in the framework, but then it becomes unclear how to obtain statistical coverage guarantees, which is left as our future research direction.\n\n$\\bar{s}$ is not a “prediction”; it is the average of the initial network (without using data or training). $\\bar{s}$ is used to generate the procedural-shifted label in Step 2, which is then used to train the second network. $\\hat{h}^*_{n}(x)$ is the output of Algorithm 1, which is the subtraction of two neural networks from Steps 1 and 2, respectively. \n\nQ2. (Compare with Bayesian approaches) We cited these works and discussed them in Section 1: “While powerful, these approaches nonetheless possess inference error that could be hard to quantify, and ultimately finding rigorous guarantees on the performance of these approximate posteriors remains open to our best knowledge.” The major difference between our work and these works is that our UQ framework has a statistical coverage guarantee. Some Bayesian approaches are fast-to-implement empirically but do not (or have not been shown to) enjoy the statistical coverage guarantee that ours does. For the illustration in our experiments, DropoutUQ is probably the most cited (well-known) work in deep-learning-based UQ, and could be a representative of these Bayesian methods without statistical guarantees.\n\nQ3. (Classification) NTK is the essential component of our theory, providing the theoretical foundation to establish the statistical properties of our algorithms. In previous literature, NTK is generally employed for mean-square-error (MSE) loss, e.g.,\n\nDu, Simon, et al. \"Gradient descent finds global minima of deep neural networks.\" International conference on machine learning. PMLR, 2019.\n\nLee, Jaehoon, et al. \"Wide neural networks of any depth evolve as linear models under gradient descent.\" Advances in neural information processing systems 32 (2019).\n\nas well as other references for NTK in our paper. The major technical discrepancy related to NTK between the regression tasks and the classification tasks is that the softmax operation in the network final layer and the cross-entropy-type loss used in the classification does not have an explicit statistical characterization of neural network training. Therefore, most NTK-based theories consider MSE loss. It would be an interesting but different direction to study the statistical guarantee for classification neural networks, which we recognize as our future research direction. \n\nQ4. (Computation) First of all, we have added more experiments on “simulated” real-world datasets to validate our performances; please refer to our global response to all reviewers and our discussions there. We also want to emphasize that our rigorous theory holds for any $m’\\ge 2$ in PNC-enhanced batching and $R \\ge 1$ in PNC-enhanced cheap bootstrap, meaning that our approaches can construct asymptotically exact-coverage confidence intervals using as low as two PNC procedures (four network training in total) without additional overheads. In terms of running times, it is as low as 4 times the standard network (point prediction) training time.\n\nQ5. (Terminology) We understand the reviewer’s concern and would definitely remove or clarify the overclaims in the paper (e.g., “general” uncertainty quantification, etc). Nonetheless, we should point out that “Overparametrized” used in deep learning theory, by convention, means the width is sufficiently large, while the depth can be as small as 2, e.g., \n\nLi, Yuanzhi, Tengyu Ma, and Hongyang R. Zhang. \"Learning over-parametrized two-layer neural networks beyond NTK.\" Conference on learning theory. PMLR, 2020.\n\nLi, Yuanzhi, and Yingyu Liang. \"Learning overparameterized neural networks via stochastic gradient descent on structured data.\" Advances in neural information processing systems 31 (2018).\n\nNTK theory and our theory are applied to neural networks with any depth. In experiments, we use a two-layer neural network for illustration. So, this terminology appears fine to us, but we would clarify it in the paper as suggested by the reviewer.\n"}, {"review_id": "k3m3gqpp0L", "reviewer": "Reviewer_nEzW", "summary": "This paper proposes an epistemic uncertainty assessment framework which comes with statistical coverage guarantees and low computation costs for over-parametrized neural networks.  This approach seeks to remove procedural uncertainty by using one auxiliary network. In addition, approaches are provided to construct confidence intervals using a small number of retrained networks on different problem settings.\n", "questions": "How wide is wide enough that the NTK theory used is approximately valid?\n\nHave you looked into studying more realistic problems (even if still only simulated)?\n", "limitations": "Some discussion should be provided on limitations. For example, how large does the NN need to be? ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors introduce a Procedural-Noise-Correcting (PNC) predictor, which uses one auxiliary network trained on the data but where all labels set to zero (instead of +-1).  This mimicks the variability coming from the training procedure leveraging asymptotic results from NTK theory. The idea is simple and intuitive and appears to work on some toy experiments.", "weaknesses": "Experiments are only provided on low-dimensional toy problems.\n\nIn addition, for the coverage studies, only 40 samples are used which results in relatively large binomial errors on the reported coverage properties. This issue isn't discussed at all. It would be better to both increase the number of samples and to report the results in a way that respects the uncertainty due to finite samples.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We greatly appreciate the reviewer for the careful reading and valuable feedback. We address the reviewer's concerns below.\n\nQ1. (Toy problems) To fully evaluate our approaches, we conducted additional experiments on the real-world dataset. For additional experiments, please refer to our global response to all reviewers and our discussions there.\n\nQ2. (Selection on repeat times) The 40 here represents 40 experimental repetitions: Essentially, we regenerate the entire datasets 40 times and re-run our approaches on new datasets 40 times to generate 40 confidence intervals. We tried more repetition times, such as 80 times in the experiment: the plot of the coverage rate and interval width against the number of experimental repetitions is shown in the attached PDF file (in our global response to all reviewers). As shown, 40 times is sufficiently large to guarantee a robust evaluation of the coverage rate and interval width.\n\nQ3. (Width in NTK) There are some rough lower bounds of the width in the literature so that the NTK holds approximately; See, e.g.,  [*], [**].\n\nEssentially, they show that the width should be at least a polynomial function of the sample size as well as other constants that are intrinsic in the data. In the experiment, we choose a 32$\\times$sample size (main paper) or 128$\\times$sample size (Appendix) to reflect this dependence. Note that the bounds derived from theory are rough and typically too large to be implemented in practice, so we need to balance it with our actual computation units. \n\nQ4. (Experiments) For additional experiments on real-world datasets, please refer to our global response to all reviewers and our discussions there.\n\n[*]: Arora, Sanjeev, et al. \"On exact computation with an infinitely wide neural net.\" Advances in neural information processing systems 32 (2019). \n\n[**]: Du, Simon, et al. \"Gradient descent finds global minima of deep neural networks.\" International conference on machine learning. PMLR, 2019. \n\n\n\n"}]

Author rebuttal:
We greatly appreciate the reviewers for their careful reading and valuable feedback. 

In this global response, we will show some additional experiments to showcase our approach to more challenging problems than in the paper via “simulated” real-world datasets. In the following, we will first discuss why synthetic datasets are employed in our experiments and then provide our experimental results and rationale.

The major reason for using synthetic datasets in our experiments is that precise evaluation of a confidence interval requires the following two critical components, which nevertheless are not available in real-world data: 1) We need to know the ground-truth regression function, while the label in real-world data typically contains some aleatoric noise. 2) To estimate the coverage rate of the CI, we need to repeat the experiments multiple times (40 times in our paper) by keep regenerating new independent datasets (from the same data distribution). In practice, we cannot regenerate new real-world datasets.

Thus, to address the reviewer’s concern about experiments on larger, more realistic and preferably real-world data, we provide a way to mimic an evaluation of the CI on real-world datasets as follows.

Step 1. Select a subset of the real-world data $(x_i, y_i)$ ($i \in \mathcal{I}$) and a test data point $(x_{k}, y_{k})$ ($k \notin \mathcal{I}$).

Step 2. For each experimental repetition j, add some “artificial” independent noise on the label to obtain a “simulated” real-world dataset $(x_i, y_i+\epsilon_{i,j})$ ($i \in \mathcal{I}$) where $\epsilon_{i,j}$ are all independent variables. Construct CI based on this “regenerating” training data $(x_i, y_i+\epsilon_{i,j})$ ($i \in \mathcal{I}$) and evaluate the coverage on $(x_{k}, y_{k})$.

In the above setting, $y_{k}$ becomes the “true” mean response of $x_{k}$ without aleatoric noise, and $\epsilon_{i,j}$ represents the only data variability. Again, as we discussed above, precise evaluation is impossible for real-world data; the above procedure, although not precise, is the best we can do to provide a resembling evaluation.

Using this setting, we conduct experiments on real-world benchmark regression datasets from UCI datasets: Boston, Concrete, and Energy. These results are shown in the Table in our attached PDF file. 

From the results, we can see that our approaches also work well for these “simulated” real-world datasets. They can provide accurate confidence intervals that satisfy the coverage requirement of the confidence level. In contrast, DropoutUQ does not have such statistical guarantees. Therefore, our approaches not only work for synthetic datasets but also is scalable well to be applied in benchmark real-world datasets.


Meta-review: Before rebuttal the main critic of the reviewers was the limmited experimental analysis on synthetic data. During rebuttal authors added additional experimental results which convinced all reviewers to vote for acceptance. The novel results should definitely be added to a revised version of the paper.
Decision: Accept (poster)