Efficient Uncertainty Quantification and Reduction for
  Over-Parameterized Neural Networks
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.
Ziyi Huang, Henry Lam, Haofeng Zhang
# Efficient Uncertainty Quantification and Reduction for Over-Parameterized Neural Networks

###### 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 _quantify_, and _remove_, the procedural 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 _one_ auxiliary network that is trained on a suitably labeled data set, 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.

## 1 Introduction

Uncertainty quantification (UQ) concerns the dissection and estimation of various sources of errors in a 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 [64; 21] and natural language processing [106; 98] to autonomous driving [85; 86], as well as guiding exploration in sequential learning [9; 1; 111]. In the deep learning context, UQ encounters unique challenges on both the statistical and computational fronts. On a high level, these challenges arise from the over-parametrized and large-scale nature of neural networks so that, unlike classical statistical models, the prediction outcomes incur not only noises from the data, but also importantly the training procedure itself [70; 89]. This elicits a deviation from the classical statistical theory that hinders the attainment of established guarantees. Moreover, because of the sizes of these models, conventional procedures such as resampling [97] demand an amount of computation that could quickly become infeasible.

Our main goal of this paper is to propose a UQ framework for epistemic uncertainty for over-parametrized neural networks that has simultaneous _statistical coverage guarantee_, in the sense of classical frequentist asymptotic exactness, and _low computation cost_, in the sense of 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-Noise-Correcting (PNC)_ predictor. It consists of an auxiliary network that is trained on a suitably artificially labeled data set, with behavior mimicking the variability coming from the training procedure. To reach our goal, we synthesize and build on two recent lines of tools that appear largely segregated thus far. First is neural tangent kernel (NTK) theory [61; 75; 7], which provides explicit approximate
formulas for well-trained infinitely wide neural networks. Importantly, NTK reveals how procedural variability enters into the prediction outcomes through, in a sense, an approximately shifted kernel regression that guides our PNC construction. Second is light-computation resampling methodology, including batching [44; 94; 95] and the so-called cheap bootstrap method [71], which allows valid confidence interval construction using extremely few model retrainings. 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 [87; 91; 11; 59; 34], or predictive uncertainty which captures the sum of epistemic and aleatoric uncertainties (but not their dissection) [89; 92; 12; 3; 22]. Regarding epistemic uncertainty, a related line of study is deep ensemble that aggregates predictions from multiple independent training replications [76; 70; 38; 8; 89]. This approach, as 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 [41; 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 [40; 46; 16; 31; 30; 90; 74; 53]. 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.

## 2 Statistical Framework of Uncertainty

We first describe our learning setting and define uncertainties in our framework. Suppose the input-output pair \((X,Y)\) is a random vector following an unknown probability distribution \(\pi\) on \(\mathcal{X}\times\mathcal{Y}\), where \(X\in\mathcal{X}\subset\mathbb{R}^{d}\) is an input and \(Y\in\mathcal{Y}\subset\mathbb{R}\) is the response. Let the marginal distribution of \(X\) be \(\pi_{X}(x)\) and the conditional distribution of \(Y\) given \(X\) be \(\pi_{Y|X}(y|x)\). Given a set of training data \(\mathcal{D}=\{(x_{1},y_{1}),(x_{2},y_{2}),...,(x_{n},y_{n})\}\) drawn i.i.d. from \(\pi\) (we write \(\mathbf{x}=(x_{1},...,x_{n})^{T}\) and \(\mathbf{y}=(y_{1},...,y_{n})^{T}\) for short), we build a prediction model \(h:\mathcal{X}\rightarrow\mathcal{Y}\) that best approximates \(Y\) given \(X\). Let \(\pi_{\mathcal{D}}\) be the empirical distribution associated with the training data \(\mathcal{D}\).

To this end, provided a loss function \(\mathcal{L}:\mathcal{Y}\times\mathbb{R}\rightarrow\mathbb{R}\), we denote the population risk \(R_{\pi}(h):=\mathbb{E}_{(X,Y)\sim\pi}[\mathcal{L}(h(X),Y)]\). If \(h\) is allowed to be any possible functions, the best prediction model is the _Bayes predictor_[58; 88]: \(h^{*}_{B}(X)\in\text{argmin}_{y\in\mathcal{Y}}\mathbb{E}_{Y\sim\pi_{Y|X}}[ \mathcal{L}(y,Y)|X]\). With finite training data of size \(n\), classical statistical learning suggests finding \(\hat{h}_{n}\in\mathcal{H}\), where \(\mathcal{H}\) denotes a hypothesis class that minimizes the empirical risk, i.e., \(\hat{h}_{n}\in\text{argmin}_{h\in\mathcal{H}}R_{\pi_{\mathcal{D}}}(h)\). This framework fits methods such as linear or kernel ridge regression (Appendix B), but because of non-convexity and other issues, it is not feasible for deep learning. Instead, in practice, gradient-based optimization methods are used, giving rise to \(\hat{h}_{n,\gamma}\) as a variant of \(\hat{h}_{n}\), where the additional variable \(\gamma\) represents the randomness in the training procedure. It is worth mentioning that this randomness generally depends on the empirical data \(\mathcal{D}\), and thus we use \(P_{\gamma|\mathcal{D}}\) to represent the distribution of \(\gamma\) conditional on \(\mathcal{D}\).

Furthermore, we consider \(\hat{h}_{n}^{*}=\text{aggregate}(\{\hat{h}_{n,\gamma}:\gamma\sim P_{\gamma| \mathcal{D}}\})\) where "aggregate" means an idealized aggregation approach to remove the training randomness in \(\hat{h}_{n,\gamma}\) (known as ensemble learning [70; 18; 17; 82; 42]). A prime example in deep learning is deep ensemble [70] that we will detail in the sequel. Finally, we denote \(h^{*}=\lim_{n\rightarrow\infty}\hat{h}_{n}^{*}\) as the grand "limiting" predictor when we have infinite samples. The exact meaning of "lim" will be clear momentarily. Under this framework, we can dissect epistemic uncertainty into three sources (illustrated in Figure 1):

**Model approximation error**\(\text{UQ}_{AE}=h^{*}-h_{B}^{*}\). This discrepancy between \(h_{B}^{*}\) and \(h^{*}\) arises from the inadequacy of the hypothesis class \(\mathcal{H}\). For an over-parameterized sufficiently wide neural network \(\mathcal{H}\), this error is usually negligible because of the universal approximation power of neural networks for any continuous functions [27; 55; 50] or Lebesgue-integrable functions [83].

Figure 1: Three sources of epistemic uncertainty.
**Data variability** UQ\({}_{DV}=\hat{h}_{n}^{*}-h^{*}\). This measures the representativeness of the training dataset, which is the most standard epistemic uncertainty in classical statistics [104].

**Procedural variability** UQ\({}_{PV}=\hat{h}_{n,\gamma}-\hat{h}_{n}^{*}\). This arises from the randomness in the training process for a single network \(\hat{h}_{n,\gamma}\), which is present even with deterministic or infinite data. The randomness comes from the initialization of the network parameters, and also data ordering and possibly training time when we run stochastic gradient descent with finite training epochs.

## 3 Quantifying Epistemic Uncertainty

We use a frequentist framework and, for a given \(x\), we aim to construct a confidence interval for the Bayes predictor \(h_{B}^{*}(x)\). As discussed in the introduction, the over-parametrized nature of neural networks defies conventional statistical techniques and moreover introduces procedural variability that makes inference difficult. We focus on sufficiently wide neural networks, which gives two theoretical advantages. First, this makes model approximation error negligible and a confidence interval for \(h^{*}(x)\) also applies to \(h_{B}^{*}(x)\). Second, the NTK theory [61] implies a phenomenon that the network evolves essentially as a "linear model" under gradient descent, and thus the resulting predictor behaves like a (shifted) kernel ridge regression whose kernel is the NTK [7; 75; 51; 56; 110] (more details in Appendix C). However, we do not find off-the-shelf results that exactly match our need for the epistemic uncertainty task, so we describe our result below. Consider the following regression problem:

\[\hat{R}_{n}(f_{\theta},\theta^{b})=\frac{1}{n}\sum_{i=1}^{n}(f_{ \theta}(x_{i})-y_{i})^{2}+\lambda_{n}\|\theta-\theta^{b}\|_{2}^{2}. \tag{1}\]

where \(\theta^{b}\) is an initialization network parameter and \(\lambda_{n}>0\) is the regularization hyper-parameter which may depend on the data size. We add regularization \(\lambda_{n}\) to this problem, which is slightly different from previous work [75] without the regularization \(\lambda_{n}\), since it can guarantee the stable computation of the inversion of the NTK Gram matrix and can be naturally linked to kernel ridge regression, as we will show shortly in Proposition 3.1. We assume the network adopts the NTK parametrization, and its parameters are randomly initialized using the He initialization, and it is sufficiently (infinitely) wide so that the linearized neural network assumption holds; See Appendix C for details. Moreover, we assume that the network is trained using the loss function in (1) via continuous-time gradient flow by feeding the entire training data and using sufficient training time (\(t\rightarrow\infty\)). In this sense, the uncertainty of data ordering and training time vanishes, and thus the random initialization is the only uncertainty in procedural variability. The above specifications are formally summarized in Assumption C.2. With them, we have

**Proposition 3.1** (Proposition C.3).: _Suppose that Assumption C.2 holds. Then the final trained network, conditional on the initial network \(s_{\theta^{b}}(x)\), is given by_

\[\hat{h}_{n,\theta^{b}}(x)=s_{\theta^{b}}(x)+\mathbf{K}(x,\mathbf{x})^{T} (\mathbf{K}(\mathbf{x},\mathbf{x})+\lambda_{n}n\mathbf{I})^{-1}(\mathbf{y}-s_{\theta^{b}}(\mathbf{x})), \tag{2}\]

_where in the subscript of \(\hat{h}_{n,\theta^{b}}\), \(n\) represents \(n\) training data, \(\theta^{b}\) represents an instance of the initial network parameters drawn from the standard Gaussian distribution \(P_{\theta^{b}}=\mathcal{N}(0,\mathbf{I}_{p})\) where the dimension of \(\theta^{b}\) is \(p\) (He initialization); \(\mathbf{K}(\mathbf{x},\mathbf{x}):=(K(x_{i},x_{j}))_{i,j=1,...,n}\) and \(\mathbf{K}(x,\mathbf{x}):=(K(x,x_{1}),K(x,x_{2}),...,K(x,x_{n}))^{T}\) where \(K\) is the (population) NTK. This implies that \(\hat{h}_{n,\theta^{b}}\) is the solution to the following kernel ridge regression: \(s_{\theta^{b}}+\underset{g\in\bar{\mathcal{H}}}{\arg\min}\ \frac{1}{n}\sum_{i=1}^{n}(y_{i}-s_{\theta^{b}}(x_{i})-g(x_{i}))^{2}+ \lambda_{n}\|g\|_{\bar{\mathcal{H}}}^{2}\) where \(\bar{\mathcal{H}}\) is the reproducing kernel Hilbert space constructed from the NTK \(K(x,x^{\prime})\)._

Proposition 3.1 shows that the shifted kernel ridge regressor using NTK with a shift from an initial function \(s_{\theta_{b}}\) is exactly the linearized neural network regressor that starts from the initial network \(s_{\theta_{b}}\). We provide details and deviation of Proposition 3.1 in Appendix C.

### Challenges from the Interplay of Procedural and Data Variabilities

Next, to motivate our main approach based on the PNC predictor, let us discuss the existing challenges in quantifying and reducing uncertainty under the NTK framework. To this end, deep ensemble [76; 70; 38; 8; 89] is arguably the most common ensemble approach in deep learning to reduce procedural variability. [8] shows that deep ensemble achieves the best performance compared with a wide
range of other ensemble methods, and more networks in deep ensemble lead to better performance. Specifically, the _deep ensemble predictor_[70]\(\hat{h}_{n}^{m}(x)\) is defined as: \(\hat{h}_{n}^{m}(x):=\frac{1}{m}\sum_{i=1}^{m}\hat{h}_{n,\theta_{i}^{k}}(x)\) where \(m\) is the number of networks in the ensemble, \(\hat{h}_{n,\theta_{i}^{b}}(x)\) is the independently trained network with initialization \(\theta_{j}^{b}\) (with the same training data \(\mathcal{D}\)) and \(\theta_{1}^{b},...,\theta_{m}^{b}\) are i.i.d. samples drawn from \(P_{\theta^{b}}\). We also introduce \(\hat{h}_{n}^{*}(x):=\mathbb{E}_{P_{\theta^{b}}}[\hat{h}_{n,\theta^{b}}(x)]\) as the expectation of \(\hat{h}_{n,\theta^{b}}(x)\) with respect to \(\theta^{b}\sim P_{\theta^{b}}\). Taking \(m\rightarrow\infty\) and using the law of large numbers, we have \(\lim_{m\rightarrow\infty}\hat{h}_{n}^{m}(x)=\hat{h}_{n}^{*}(x)\;a.s..\). So \(\hat{h}_{n}^{*}(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

\[\hat{h}_{n}^{m}(x)=\frac{1}{m}\sum_{i=1}^{m}s_{\theta_{i}^{b}}(x)+\mathbf{K}(x, \mathbf{x})^{T}(\mathbf{K}(\mathbf{x},\mathbf{x})+\lambda_{n}n\mathbf{I})^{-1}\left(\mathbf{y}-\frac{1 }{m}\sum_{i=1}^{m}s_{\theta_{i}^{b}}(\mathbf{x})\right)\]

\[\hat{h}_{n}^{*}(x)=\mathbb{E}_{P_{\theta^{b}}}[\hat{h}_{n,\theta^{b}}(x)]= \bar{s}(x)+\mathbf{K}(x,\mathbf{x})^{T}(\mathbf{K}(\mathbf{x},\mathbf{x})+\lambda_{n}n\mathbf{I})^{-1} (\mathbf{y}-\bar{s}(\mathbf{x})) \tag{3}\]

where \(\bar{s}(x)=\mathbb{E}_{P_{\theta^{b}}}[s_{\theta^{b}}(x)]\) is the expectation of \(s_{\theta^{b}}(x_{0})\) with respect to the the distribution of the initialization parameters \(P_{\theta^{b}}=\mathcal{N}(0,\mathbf{I}_{p})\). It is easy to see that \(\mathbb{E}_{P_{\theta^{b}}}[\hat{h}_{n}^{m}(x)]=\hat{h}_{n}^{*}(x)\) and \(\text{Var}_{P_{\theta^{b}}}(\hat{h}_{n}^{m}(x))=\frac{1}{m}\text{Var}_{P_{ \theta^{b}}}(\hat{h}_{n,\theta^{b}}(x))\) where \(\text{Var}_{P_{\theta^{b}}}\) is the variance with respect to the random initialization. As for the total variance,

**Proposition 3.2**.: _We have_

\[\text{Var}(\hat{h}_{n}^{m}(x))=\text{Var}(\hat{h}_{n}^{*}(x))+\frac{1}{m} \mathbb{E}[\text{Var}(\hat{h}_{n,\theta^{b}}(x)|\mathcal{D})]\leq\text{Var}( \hat{h}_{n}^{*}(x))+\mathbb{E}[\text{Var}(\hat{h}_{n,\theta^{b}}(x)|\mathcal{D })]=\text{Var}(\hat{h}_{n,\theta^{b}}(x)),\]

_where the variance is taken with respect to both the data \(\mathcal{D}\) and the random initialization \(P_{\theta^{b}}\)._

Proposition 3.2 shows that deep ensemble improves the statistical profile of a single model by reducing its procedural variability by a factor \(\frac{1}{m}\) (but not the data variability), and achieving this reduction requires \(m\) trainings. To quantify the epistemic uncertainty in a deep ensemble, we may employ resampling approaches such as "bootstrap on a 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 the large number of outer bootstrap resamples and the \(m\) inner retaining per resample. This leads to a huge computational burden.

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 to deep ensemble that trains \(m\) networks and only reduces the procedural variability by an \(m\)-factor.

**Uncertainty quantification.** We resolve the computational challenges in "bootstrap on a deep ensemble", by combining PNC predictors with low-budget inference tools that require only as low as _four_ network trainings.

### PNC Predictor and Procedural Variability Removal

We first develop a computationally efficient approach to obtain \(\hat{h}_{n}^{*}(x)\), the idealized deep ensemble predictor that is free of procedural variability. We term our approach the procedural-noise-correcting (PNC) predictor, whose pseudo-code is given in Algorithm 1. This predictor consists of a difference between essentially two neural network outcomes, \(\hat{h}_{n,\theta^{b}}(x)\) which is the original network trained using one initialization, and \(\hat{\phi}_{n,\theta^{b}}(x)\) that is trained on a suitably artificially labeled dataset. More precisely, this dataset applies label \(\bar{s}(x_{i})\) to \(x_{i}\), the expected outcome of an _untrained_ network with random initialization. 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-label-trained network.

To motivate this predictor, we first note that, instead of directly computing \(\hat{h}_{n}^{*}(x)\), we can compute the procedural noise, which is given by

\[\hat{\phi}_{n,\theta^{b}}(\cdot):=\hat{h}_{n,\theta^{b}}(x)-\hat{h}_{n}^{*}(x)=s _{\theta^{b}}(x)-\bar{s}(x)+\mathbf{K}(x,\mathbf{x})^{T}(\mathbf{K}(\mathbf{x},\mathbf{x})+\lambda _{n}n\mathbf{I})^{-1}(\bar{s}(\mathbf{x})-s_{\theta^{b}}(\mathbf{x})). \tag{4}\]
By Proposition 3.1, this closed-form expression of \(\hat{\phi}_{n,\theta^{b}}\) corresponds exactly to the auxiliary network described above. This subsequently leads to:

**Theorem 3.3** (Pnc).: _Suppose that Assumption C.2 holds. Then the output of the PNC predictor (Algorithm 1) is exactly \(\hat{h}_{n}^{*}(x)\) given in (3)._

There are two approaches to computing \(\bar{s}(x)\) in Algorithm 1: 1) Under He initialization, \(s_{\theta^{b}}(x)\) is a zero-mean Gaussian process in the infinite width limit [75], and thus we can numerically set \(\bar{s}\equiv 0\). 2) An alternative approach that does not require specific initialization is to use Monte Carlo integration: \(\bar{s}(x)=\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}s_{\theta^{b}_{i}}(x)\) where \(\theta^{b}_{i}\) are i.i.d. from \(P_{\theta^{b}}\). When \(N\) is finite, it introduces procedural variance that vanishes at the order \(N^{-1}\). However, since this computation does not require any training and is extremely fast in practice, we may use \(N\gg n\) to guarantee that compared with data variance at the order \(n^{-1}\), procedural variance in computing \(\bar{s}(x)\) is negligible. In experiments, both approaches work well, and we will use the first approach for efficiency. Note that although \(\bar{s}\equiv 0\) in the infinite width limit [75], it does not imply that the trained auxiliary neural network in Algorithm 1 will be a zero network (all parameters are zero). In fact, over-parameterized neural networks have many global minima, and the auxiliary neural network will find one of the global minima that is very close to the \(\theta^{b}\)[33; 10].

Compared with Algorithm 1, the following candidate approaches encounter computational issues: 1) Using deep ensemble with sufficiently many networks in the ensemble. However, this approach is time-consuming as \(m\) networks in the ensemble mean \(m\)-fold training times. \(m\) is typically as small as 5 in practice [70] so it cannot approximate \(\hat{h}_{n}^{*}(x)\) well. 2) Using the closed-form expression of \(\hat{h}_{n}^{*}\) in (3), which requires computing the NTK \(K(x,x^{\prime})\) and the inversion of the NTK Gram matrix. \(K(x,x^{\prime})\) is recursively defined and does not have a simple form for computation, which might be addressed by approximating it with the empirical NTK \(K_{\theta}(x,x^{\prime})\) numerically (See Appendix C). However, the inversion of the NTK Gram matrix is another issue: When the training data size is large, the NTK Gram matrix is large, and the matrix inversion becomes time-consuming.

### Constructing Confidence Intervals from PNC Predictors

We construct confidence intervals for \(h^{*}(x)\) leveraging our PNC predictor in Algorithm 1. 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 [48; 35] that acts as the functional gradient of the predictor with respect to the data distribution. It was initially introduced in modern machine learning for understanding a training point's effect on a model's prediction [69; 68; 13]. The second is to use resamplings [97], such as bootstrap or jackknife, to avoid explicit variance computation. The classical resampling method requires enough resample replications and thus incurs demanding resampling effort. For instance, the jackknife requires the same number of training times as the training data size, which is extremely time-consuming. Standard bootstrap requires a sufficient number of resampling and retraining to produce accurate resample quantiles. Given these computational bottleneck, we will utilize light-computation resampling alternatives, including batching [44; 94; 95] and the so-called cheap bootstrap method [71; 72], which allows valid confidence interval construction using extremely few model retrainings.

**Large-sample asymptotics of the PNC predictor.** To derive our intervals, we first gain understanding on the large-sample properties of the PNC predictor. Proposition 3.1 shows that \(\hat{h}_{n}^{*}(x)\) in (3) is the
solution to the following empirical risk minimization problem:

\[\hat{h}_{n}^{*}(\cdot)=\bar{s}(\cdot)+\min_{g\in\mathcal{H}}\frac{1}{n}\sum_{i=1 }^{n}[(y_{i}-\bar{s}(x_{i})-g(x_{i}))^{2}]+\lambda_{n}\|g\|_{\mathcal{H}}^{2} \tag{5}\]

Therefore, its corresponding population risk minimization problem (i.e., removing the data variability) is:

\[h^{*}(\cdot)=\bar{s}(\cdot)+\min_{g\in\mathcal{H}}\mathbb{E}_{\pi}[(Y-\bar{s}( X)-g(X))^{2}]+\lambda_{0}\|g\|_{\mathcal{H}}^{2} \tag{6}\]

where \(\lambda_{0}=\lim_{n\to\infty}\lambda_{n}\). To study the difference between the empirical and population risk minimization 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 [47]).: _Suppose that Assumptions B.3 and B.4 hold. Let \(f_{P}\) be the solution to the following problem: \(\min_{g\in\mathcal{H}}\mathbb{E}_{P}[(Y-g(X))^{2}]+\lambda_{P}\|g\|_{\mathcal{ H}}^{2}\) where \(P=\pi_{n}\) or \(\pi\). Then \(\sqrt{n}(f_{\pi_{n}}-f_{\pi})\Rightarrow\mathbb{G}\) in \(\mathcal{\bar{H}}\) where \(\mathbb{G}\) is a zero-mean Gaussian process and \(\Rightarrow\) represents "converges weakly". Moreover, at point \(x\), \(\sqrt{n}(f_{\pi_{n}}(x)-f_{\pi}(x))\Rightarrow\mathcal{N}(0,\sigma_{\pi}^{2}(x))\) where \(\sigma_{\pi}^{2}(x)=\int_{z\in\mathcal{X}\times\mathcal{Y}}IF^{2}(z;f_{P},\pi) (x)d\pi(z)\) and \(IF\) is the influence function of statistical functional \(f_{P}\)._

Next, we apply the above proposition to our problems about \(\hat{h}_{n}^{*}\). Let \(T_{1}(P)(x)\) be the solution of the following problem \(\min_{g\in\mathcal{H}}\mathbb{E}_{P}[(Y-\bar{s}(X)-g(X))^{2}]+\lambda_{P}\|g \|_{\mathcal{H}}^{2}\) that is evaluated at a point \(x\). Here \(\lambda_{\pi_{n}}=\lambda_{n}\) and \(\lambda_{\pi}=\lambda_{0}\). Then we have the following large-sample asymptotics of the PNC predictor, which provides the theoretical foundation for constructing an asymptotic confidence interval.

**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-\bar{s}(X)\). Input the training data \(\mathcal{D}\) into Algorithm 1 to obtain \(\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)\). We have_

\[\sqrt{n}\left(\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)-h^{*}(x) \right)\Rightarrow\mathcal{N}(0,\sigma_{\theta^{b}}^{2}(x)), \tag{7}\]

_where_

\[\sigma_{\theta^{b}}^{2}(x)=\int_{z\in\mathcal{X}\times\mathcal{Y}}IF^{2}(z;T_ {1},\pi)(x)d\pi(z). \tag{8}\]

_Thus, an asymptotically (in the sense of \(n\to\infty\)) exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x)\) is \([\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)-q_{1-\frac{\alpha}{2}},\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)+q_{1-\frac{\alpha}{2}}]\) where \(q_{\alpha}\) is the \(\alpha\)-quantile of the distribution \(\mathcal{N}(0,\frac{1}{n}\sigma_{\theta^{b}}^{2}(x))\)._

In general, \(\sigma_{\theta^{b}}^{2}\) is unknown and must be estimated. It is common to approximate \(IF^{2}(z;T_{1},\pi)\) with \(IF^{2}(z;T_{1},\pi_{n})\) and \(\hat{\sigma}_{\theta^{b}}^{2}(x)=\frac{1}{n}\sum_{z_{i}\in\mathcal{D}}IF^{2}(z _{i};T_{1},\pi_{n})(x)\). This method is known as the infinitesimal jackknife variance estimator [97]. We provide the closed-form expression of the variance estimation \(\hat{\sigma}_{\theta^{b}}^{2}(x)\) in Appendix B that could be of theoretical interest. Yet in practice, the computation of \(\hat{\sigma}_{\theta^{b}}^{2}(x)\) requires the evaluation of the NTK Gram matrix and its inversion, which is computationally demanding for large \(n\) and thus not recommended for practical implementation on a large dataset. In the following, we provide two approaches that avoid the explicit estimation of the asymptotic variance in the infinitesimal jackknife approach, especially the Gram matrix inversion.

**PNC-enhanced batching.** We propose an approach for constructing a confidence interval that is especially useful for large datasets, termed _PNC-enhanced batching_. The pseudo-code is given in Algorithm 2. Originating from simulation analysis [94, 95, 44], the key idea of batching is to construct a self-normalizing \(t\)-statistic that "cancels out" the unknown variance, leading to a valid confidence interval without explicitly needing to compute the variance. It 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 [43, 37, 63] and also the so-called input uncertainty problem [45]. Its application in deep learning uncertainty quantification was not revealed in previous work, potentially due to the additional procedural variability. Integrating it with the PNC predictor, PNC-enhanced batching is very efficient and enjoys the asymptotically exact coverage of its confidence interval, as stated below.
**Theorem 3.6** (Exact coverage of PNC-enhanced batching confidence interval).: _Suppose that Assumption C.2 holds. Suppose that Assumptions B.3 and B.4 hold when \(Y\) is replaced by \(Y-\bar{s}(X)\). Then an asymptotically exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x)\) is \([\psi_{B}(x)-\frac{S_{B}(x)}{\sqrt{m^{\prime}}}q_{1-\frac{\alpha}{2}},\psi_{B}( x)+\frac{S_{B}(x)}{\sqrt{m^{\prime}}}q_{1-\frac{\alpha}{2}}]\) where \(q_{\alpha}\) is the \(\alpha\)-quantile of the t distribution \(t_{m^{\prime}-1}\) with degree of freedom \(m^{\prime}-1\)._

**PNC-enhanced cheap bootstrap.** We propose an alternative approach for constructing a confidence interval that is also suitable for small datasets, termed _PNC-enhanced cheap bootstrap_. The pseudo-code is given in Algorithm 3. Cheap bootstrap [71, 72] is a modified bootstrap procedure with substantially less retraining effort than conventional bootstrap methods, via leveraging the independence between resample and original estimators. Note that our proposal is fundamentally different from the naive use of bootstrap or bagging when the additional randomness appears [70, 76, 42], which mixes the procedural and data variability, 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 them is that the 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 the samples from the entire dataset, which is also suitable for small datasets. In terms of running time, when \(R=m^{\prime}-1\), PNC-enhanced cheap bootstrap and PNC-enhanced batching share the same number of network training, but since batching is trained on the 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 the asymptotically exact coverage of its confidence interval, as stated below.

**Input:** Training dataset \(\mathcal{D}\) of size \(n\). The number of replications \(R\geq 1\).

**Procedure: 1.** Input \(\mathcal{D}_{1}\cup\mathcal{D}_{2}\) in Algorithm 1 to output \(\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)\).

**2.** For each replication \(b\in[R]\), resample \(\mathcal{D}\), i.e., independently and uniformly sample with replacement from \(\{(x_{1},y_{1}),...,(x_{n},y_{n})\}\)\(n\) times to obtain \(\mathcal{D}^{*b}=\{(x_{1}^{*b},y_{1}^{*b}),...,(x_{n}^{*b},y_{n}^{*b})\}\). Input \(\mathcal{D}^{*b}\) in Algorithm 1 to output \(\hat{h}_{n,\theta^{b}}^{*b}(x)-\hat{\phi}_{n,\theta^{b}}^{*b}(x)\).

**3.** Compute \(\psi_{C}(x)=\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)\),

and \(S_{C}(x)^{2}=\frac{1}{R}\sum_{b=1}^{R}\left(\hat{h}_{n,\theta^{b}}^{*b}(x)- \hat{\phi}_{n,\theta^{b}}^{*b}(x)-\psi_{C}(x)\right)^{2}.\)

**Output:** At point \(x\), output \(\psi_{C}(x)\) and \(S_{C}(x)^{2}\).

**Algorithm 3** PNC-Enhanced Cheap Bootstrap

**Theorem 3.7** (Exact coverage of PNC-enhanced cheap bootstrap confidence interval).: _Suppose that Assumption C.2 holds. Suppose that Assumptions B.3 and B.4 hold when \(Y\) is replaced by \(Y-\bar{s}(X)\). Then an asymptotically exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x)\) is \([\psi_{C}(x)-S_{C}(x)q_{1-\frac{\alpha}{2}},\psi_{C}(x)+S_{C}(x)q_{1-\frac{ \alpha}{2}}]\) where \(q_{\alpha}\) is the \(\alpha\)-quantile of the t distribution \(t_{R}\) with degree of freedom \(R\)._

## 4 Experiments

We conduct numerical experiments to demonstrate the effectiveness of our approaches. Our code will be publicly available. Our proposed approaches are evaluated on the following two tasks: 1) construct confidence intervals and 2) remove procedural variability. With a 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\times 32\) neurons in each hidden layer to ensure the network is sufficiently wide and over-parameterized. Detailed optimization specifications are described in Proposition C.3. Our synthetic dataset is generated with the following distribution: \(X\sim\text{Unif}([0,0.2]^{d})\) and \(Y\sim\sum_{i=1}^{d}\sin(X^{(i)})+\mathcal{N}(0,0.001^{2})\). The training set \(\mathcal{D}=\{(x_{i},y_{i}):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=1,2,4,8\) and data size settings \(n=128,256,512,1024\) to study the effects on different dimensionalities and data sizes. Additional experimental results are presented in Appendix E.

**Constructing confidence intervals.** We use \(x_{0}=(0.1,...,0.1)\) as the fixed test point for confidence intervals construction. Let \(y_{0}=\sum_{i=1}^{d}\sin(0.1)\) be the ground-truth label for \(x_{0}\) without aleatoric noise. To evaluate the construction performance on confidence intervals, we set the repeat time \(J=40\). In each repetition \(j\in[J]\), we generate a new training dataset from the same synthetic distribution and construct a new confidence interval \([L_{j}(x_{0}),U_{j}(x_{0})]\) with \(95\%\) or \(90\%\) confidence level, and then check the coverage rate (CR): \(\text{CR}=\frac{1}{J}\sum_{j=1}^{J}\mathbf{1}_{L_{j}(x_{0})\leq y_{0}\leq U_{j }(x_{0})}\). 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, we also report the median point of the interval (MP): \(\text{MP}=\frac{1}{J}\sum_{j=1}^{J}\frac{1}{2}(U_{j}(x_{0})+L_{j}(x_{0}))\) and the interval width (IW): \(\text{IW}=\frac{1}{J}\sum_{j=1}^{J}(U_{j}(x_{0})-L_{j}(x_{0}))\). MP reflects the most likely point prediction for \(x_{0}\), 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 construction tasks, we use the dropout-based Bayesian approximate inference (DropoutUQ) [40] as the baseline for comparison since it is the most widely-used baseline in Bayesian uncertainty quantification. The implementation details are provided in Appendix E. Table 1 reports the CR, IW, and Mean of \(95\%\) and \(90\%\) confidence intervals from our proposals and baselines. We have the following observations:

1) The CR values: In the majority of experiments, \(\text{CR}\geq 95\%\) for \(95\%\) confidence intervals and \(\text{CR}\geq 90\%\) for \(90\%\) confidence intervals satisfy the requirement of confidence intervals, while in very few experiments, the performances are degraded. This occasional degeneration is insignificant and can be potentially explained by the fact that the statistical guarantee of confidence intervals generated by PNC-enhanced batching and cheap bootstrap is in the asymptotic sense. Overall, the CR for \(95\%\)/\(90\%\) confidence intervals are always above \(92.5\%\)/\(87.5\%\) from both proposed methods. In contrast, DropoutUQ does not have such statistical guarantees. The intervals from DropoutUQ cannot maintain a satisfactory coverage rate when the interval is narrow, although they have a similar or larger interval width as PNC-enhanced batching or cheap bootstrap. Moreover, we find that the

\begin{table}
\begin{tabular}{l|c c c|c c c|c c c c} \hline \hline  & \multicolumn{3}{c|}{PNC-enhanced batching} & \multicolumn{3}{c|}{PNC-enhanced cheap bootstrap} & \multicolumn{3}{c}{DropoutUQ} \\  & 95\%(CIR/W) & 90\%(CIR/W) & MP & 95\%(CIR/W) & 90\%(CIR/CIR/W) & MP & 95\%(CIR/W) & 90\%(CIR/W) & MP \\ \hline \((d=1)\) & \(\mathbf{0.95}\)/0.00598 & \(\mathbf{0.95}\)/0.00042 & 0.10031 & \(\mathbf{1.0}\)/0.00557 & \(\mathbf{0.95}\)/0.00428 & \(\mathbf{0.9987}\)/**0.00565** & \(\mathbf{0.925}\)/0.00341 & \(\mathbf{0.09972}\) \\ \(n=256\) & \(\mathbf{1.0}\)/0.00437 & \(\mathbf{0.95}\)/0.00323 & \(\mathbf{0.09990}\) & \(\mathbf{0.95}\)/0.00303 & \(\mathbf{0.87}\)/0.00233 & \(\mathbf{0.09986}\) & \(\mathbf{1.0}\)/0.00336 & \(\mathbf{0.925}\)/0.00211 & \(\mathbf{0.09973}\) \\ \(n=512\) & \(\mathbf{0.975}\)/0.00241 & \(\mathbf{0.95}\)/0.00178 & \(\mathbf{0.9994}\) & \(\mathbf{1.0}\)/0.00262 & \(\mathbf{1.0}\)/0.00173 & \(\mathbf{0.09986}\) & \(\mathbf{0.975}\)/0.00243 & \(\mathbf{0.875}\)/0.00153 & \(\mathbf{0.09983}\) \\ \(n=1024\) & \(\mathbf{0.975}\)/0.00163 & \(\mathbf{0.90}\)/0.00120 & \(\mathbf{0.09983}\) & \(\mathbf{0.975}\)/0.00145 & \(\mathbf{0.975}\)/0.00107 & \(\mathbf{0.99968}\) & \(\mathbf{0.925}\)/0.00145 & \(\mathbf{0.600}\)/0.000815 & \(\mathbf{0.09978}\) \\ \hline \((d=2)\) & \(\mathbf{0.975}\)/0.01810 & \(\mathbf{0.99}\)/0.01344 & \(\mathbf{0.20205}\) & \(\mathbf{0.95}\)/0.01252 & \(\mathbf{0.925}\)/0.00408 & \(\mathbf{0.19984}\) & \(\mathbf{0.95}\)/0.01853 & \(\mathbf{0.95}\)/0.01466 & \(\mathbf{0.20131}\) \\ \(n=256\) & \(\mathbf{1.0}\)/0.01260 & \(\mathbf{0.90}\)/0.009317 & \(\mathbf{0.2010}\)/1.000473 & \(\mathbf{1.00}\)/0.00656 & \(\mathbf{1.1991}\)/0.07950 & \(\mathbf{0.990}\)/0.007128 & \(\mathbf{0.19903}\) \\ \(n=512\) & \(\mathbf{1.0}\)/0.005555 & \(\mathbf{0.925}\)/0.004847 & \(\mathbf{1.9999}\)/0.005612 & \(\mathbf{0.925}\)/0.004347 & \(\mathbf{0.9197}\)/0.0250/0.01644 & \(\mathbf{0.850}\)/0.00185 & \(\mathbf{0.19929}\) \\ \(n=1024\) & \(\mathbf{0.95}\)/0.003825 & \(\mathbf{0.875}\)/0.002829 & \(\mathbf{1.0974}\) & \(\mathbf{1.0}\)/0.004187 & \(\mathbf{1.0}\)/0.003214 & \(\mathbf{1.0973}\)/0.85005160 & \(\mathbf{0.750}\)/0.003905 & \(\mathbf{0.19893}\) \\ \hline \((d=4)\) & \(\mathbf{0.4}\) & \(\mathbf{1.0}\)/0.05706 & \(\mathbf{0.975}\)/0.04219 & \(\mathbf{0.409}\)/0.00929 & \(\mathbf{0.975}\)/0.03487 & \(\mathbf{0.925}\)/0.02677 & \(\mathbf{0.39605}\)/0.04288 & \(\mathbf{0.875}\)/0.03489 & \(\mathbf{0.40271}\) \\ \(n=256\) & \(\mathbf{0.95}\)/0.02495 & \(\mathbf{0.925}\)/0.01845 & \(\mathbf{0.9640}\)/0.02553 & \(\mathbf{0.925}\)/0.01730 & \(\mathbf{0.3961}\)/0.03691 & \(\mathbf{0.95}\)/0.02524 & \(\mathbf{0.925}\)/0.02067 & \(\mathbf{0.40111}\) \\ \(n=512\) & \(\mathbf{0.925}\)/0.01805 & \(\mathbf{0.875}\)/0.01334 & \(\mathbf{0.400}\)/0.0950 & \(\mathbf{0.95}\)/0.019147 & \(\mathbf{0.990}\)/0.01471 & \(\mathbf{0.97367}\)/0.900191 & \(\mathbf{0.825}\)/0.01556 & \(\mathbf{0.39820}\) \\ \(n=1024\) & \(\mathbf{0.975}\)/0.01412 & \
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 evidently demonstrate the robustness and effectiveness of our proposals. 2) The IW values: As shown, the IW values reduce roughly at the rate \(n^{-\frac{1}{2}}\) along with an increased training data size, which corroborates with our theory (Section 3.3) well. In comparison with PNC-enhanced cheap bootstrap, the PNC-enhanced batching tends to generate wider intervals when the sample size is small, as it requires data splits, and thus its accuracy might be hindered on small datasets. Herein, we recommend using PNC-enhanced batching for larger datasets while applying PNC-enhanced cheap bootstrap for smaller datasets. Overall, our IW values are small for both approaches, indicating that our framework is able to successfully generate narrow confidence intervals with high coverage rates. 3) The MP values: MP appears highly consistent to the ground-truth label \(y_{0}\) in all experiments, showing that our confidence intervals can accurately capture the ground-truth label on average under multiple problem settings.

**Removing procedural variability.** We compare our proposed PNC predictor (Algorithm 1) with the following approaches: a base network, that is, a standard network training with a single random initialization, and the deep ensemble approach with different numbers (\(m\)) of networks in the ensemble [70]. In deep ensemble, we consider \(m=2\) as it has a similar running time as our PNC (with one auxiliary network), and \(m=5\) as it is widely used in previous work [70; 89]. All networks share the same hyperparameters and specifications. Table 2 reports the mean and standard deviation of mean square error (MSE) from our proposed PNC predictor and other baseline approaches with \(10\) repetition times. The MSE is computed as \(\text{MSE}(h):=\frac{1}{N_{te}}\sum_{i=1}^{N_{te}}(h(x^{\prime}_{i})-y^{\prime }_{i})^{2}\) where \((x^{\prime}_{i},y^{\prime}_{i}),i=1,...,N_{te}\) are i.i.d. testing data independent of the training. We set \(N_{te}=2048\) for test performance evaluation in all experiments. From Table 2, we observe that our 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/analogous results as the deep ensemble with 5 networks, but the latter has 2.5 times the running time as much as PNC. These results indicate that our proposed PNC can achieve compatible results with deep ensemble but with significantly reduced computational costs.

## 5 Concluding Remarks

In this study, we propose a novel epistemic uncertainty assessment framework with simultaneous statistical coverage guarantees and low computation costs for over-parametrized neural networks. Benefiting from our proposed PNC approach, our study shows promise to remove procedural uncertainty by using only one auxiliary network. Integrated with suitable light-computation resampling methods, we provide two effective approaches to construct asymptotically exact-coverage confidence intervals using as low as two retrained networks on different problem settings. Our evaluation results fully corroborate our theory and show that our approach can generate confidence intervals that are narrow and has satisfactory confidence levels simultaneously.

\begin{table}
\begin{tabular}{c|c c c} \hline \hline MSE & One base network & PNC predictor & Deep ensemble (5 networks) & Deep ensemble (2 networks) \\ \hline \(d=1\), \(n=128\) & \((6.77\pm 4.17)\times 10^{-5}\) & \((3.11\pm 1.79)\times 10^{-5}\) & \((3.58\pm 2.93)\times 10^{-5}\) & \((4.55\pm 2.87)\times 10^{-5}\) \\ \(d=1\), \(n=256\) & \((1.53\pm 1.58)\times 10^{-5}\) & \((5.74\pm 3.66)\times 10^{-6}\) & \((6.34\pm 4.35)\times 10^{-6}\) & \((1.05\pm 0.83)\times 10^{-5}\) \\ \(d=1\), \(n=512\) & \((3.97\pm 3.37)\times 10^{-6}\) & \((1.85\pm 0.67)\times 10^{-6}\) & \((2.19\pm 1.86)\times 10^{-6}\) & \((2.50\pm 3.12)\times 10^{-6}\) \\ \(d=1\), \(n=1024\) & \((1.58\pm 0.57)\times 10^{-6}\) & \((5.59\pm 2.30)\times 10^{-7}\) & \((4.50\pm 3.01)\times 10^{-7}\) & \((0.74\pm 1.01)\times 10^{-6}\) \\ \hline \(d=2\), \(n=128\) & \((6.68\pm 2.74)\times 10^{-4}\) & \((3.68\pm 1.84)\times 10^{-4}\) & \((4.22\pm 1.74)\times 10^{-4}\) & \((6.28\pm 2.74)\times 10^{-4}\) \\ \(d=2\), \(n=256\) & \((2.38\pm 1.06)\times 10^{-4}\) & \((6.22\pm 2.01)\times 10^{-5}\) & \((8.86\pm 4.94)\times 10^{-5}\) & \((1.25\pm 0.67)\times 10^{-4}\) \\ \(d=2\), \(n=512\) & \((1.03\pm 0.98)\times 10^{-4}\) & \((2.06\pm 0.72)\times 10^{-5}\) & \((3.32\pm 1.10)\times 10^{-5}\) & \((5.11\pm 3.33)\times 10^{-5}\) \\ \(d=2\), \(n=1024\) & \((6.98\pm 4.77)\times 10^{-5}\) & \((8.92\pm 5.69)\times 10^{-6}\) & \((1.72\pm 0.77)\times 10^{-5}\) & \((5.75\pm 2.28)\times 10^{-5}\) \\ \hline \(d=4\), \(n=128\) & \((2.11\pm 1.49)\times 10^{-3}\) & \((1.18\pm 0.25)\times 10^{-3}\) & \((1.53\pm 0.60)\times 10^{-3}\) & \((1.83\pm 0.80)\times 10^{-3}\) \\ \(d=4\), \(n=256\) & \((8.82\pm 3.26)\times 10^{-4}\) & \((4.09\pm 0.85)\times 10^{-4}\) & \((4.22\pm 2.01)\times 10^{-4}\) & \((5.47\pm 1.91)\times 10^{-4}\) \\ \(d=4\), \(n=512\) & \((5.35\pm 1.91)\times 10^{-4}\) & \((1.92\pm 0.53)\times 10^{-4}\) & \((2.88\pm 1.80)\times 10^{-4}\) & \((3.99\pm 1.87)\times 10^{-4}\) \\ \(d=4\), \(n=1024\) & \((2.23\pm 0.83)\times 10^{-4}\) & \((4.22\pm 0.43)\times 10^{-5}\) & \((8.50\pm 2.39)\times 10^{-5}\) & \((1.65\pm 0.57)\times 10^{-4}\) \\ \hline \(d=8\), \(n=128\) & \((4.07\pm 1.04)\times 10^{-3}\) & \((2.54\pm 0.29)\times 10^{-3}\) & \((2.73\pm 0.86)\times 10^{-3}\) & \((3.27\pm 0.96)\times 10^{-3}\) \\ \(d=8\), \(n=256\) & \((2.13\pm 0.70)\times 10^{-3}\) & \((1.05\pm 0.18)\times 10^{-3}\) & \((1.34\pm 0.33)\times 10^{-3}\) & \((1.48\pm 0.38)\times 10^{-3}\) \\ \(d=8\), \(n=512\) & \((1.36\pm 0.34)\times 10^{-3}\) & \((**5.04\pm 0.70)\times 10^{-4}\) & \((7.40\pm 1.35)\times 10^{-4}\) & \((1.08\pm 0.40)\times 10^{-3}\) \\ \(d=8\), \(n=1024\) & \((8.54\pm 2.23)\times 10^{-4}\) & \((**2.02\pm 0.24)\times 10^{-4}\) & \((3.79\pm 1.01)\times 10^{-4}\) & \((5.91\pm 1.87)\times 10^{-4}\) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Reducing procedural variability on synthetic datasets with different data sizes \(n=128,256,512\), \(1024\) and different data dimensions \(d=1,2,4,8\). The best results are in **bold**.
## References

* [1] Y. Abbasi-Yadkori, D. Pal, and C. Szepesvari. 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-297, 2021.
* [3] A. 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 overparameterization. In _International Conference on Machine Learning_, pages 242-252. PMLR, 2019.
* [6] Z. Allen-Zhu, Y. Li, and Z. Song. On the convergence rate of training recurrent neural networks. _Advances in Neural Information Processing Systems_, 32:6676-6688, 2019.
* [7] S. Arora, S. S. Du, W. Hu, Z. Li, R. R. Salakhutdinov, and R. Wang. On exact computation with an 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.06470_, 2020.
* [9] P. Auer, N. Cesa-Bianchi, and P. Fischer. Finite-time analysis of the multiarmed bandit problem. _Machine Learning_, 47(2):235-256, 2002.
* [10] Y. Bai, S. Mei, H. Wang, and C. Xiong. Understanding the under-coverage bias in uncertainty estimation. _arXiv preprint arXiv:2106.05515_, 2021.
* [11] R. F. Barber, E. J. Candes, A. Ramdas, and R. J. Tibshirani. The limits of distribution-free conditional predictive inference. _arXiv preprint arXiv:1903.04684_, 2019.
* [12] R. F. Barber, E. J. Candes, A. Ramdas, and R. J. Tibshirani. Predictive inference with the jackknife+. _arXiv preprint arXiv:1905.02928_, 2019.
* [13] S. Basu, P. Pope, and S. Feizi. Influence functions in deep learning are fragile. _arXiv preprint arXiv:2006.14651_, 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] C. Blundell, J. Cornebise, K. 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.
* [18] L. Breiman. Random forests. _Machine learning_, 45:5-32, 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_, 32:10836-10846, 2019.
* [20] Y. Cao and Q. Gu. Generalization error bounds of gradient descent for learning over-parameterized deep relu networks. _Proceedings of the AAAI Conference on Artificial Intelligence_, 34(04):3349-3356, 2020.
* [21] E. D. Carvalho, R. Clark, A. Nicastro, and P. 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, Z. Huang, H. Lam, H. Qian, and H. Zhang. Learning prediction intervals for regression: Generalization and calibration. In _International Conference on Artificial Intelligence and Statistics_, 2021.
* [23] 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.
* [24] A. Christmann and I. Steinwart. Consistency and robustness of kernel-based regression in convex risk minimization. _Bernoulli_, 13(3):799-819, 2007.
* [25] A. Christmann and I. Steinwart. How svms can estimate quantiles and the median. _Advances in Neural Information Processing Systems_, 20:305-312, 2007.
* [26] F. Cucker and S. Smale. On the mathematical foundations of learning. _Bulletin of the American Mathematical Society_, 39(1):1-49, 2002.
* [27] G. Cybenko. Approximation by superpositions of a sigmoidal function. _Mathematics of Control, Signals and Systems_, 2(4):303-314, 1989.
* [28] N. Dalmasso, T. Pospisil, A. B. Lee, R. Izbicki, P. E. Freeman, and A. I. Malz. Conditional density estimation tools in python and r with applications to photometric redshifts and likelihood-free cosmological inference. _Astronomy and Computing_, 30:100362, 2020.
* [29] M. Debruyne, M. Hubert, and J. A. Suykens. Model selection in kernel based regression using the influence function. _Journal of Machine Learning Research_, 9(10), 2008.
* [30] S. Depeweg, J. M. Hernandez-Lobato, F. Doshi-Velez, and S. Udluft. Learning and policy search in stochastic dynamical systems with bayesian neural networks. _arXiv preprint arXiv:1605.07127_, 2016.
* [31] S. Depeweg, J.-M. Hernandez-Lobato, F. Doshi-Velez, and S. Udluft. Decomposition of uncertainty in bayesian deep learning for efficient and risk-sensitive learning. In _International Conference on Machine Learning_, pages 1184-1193. PMLR, 2018.
* [32] 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.
* [33] 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.
* [34] V. Dutordoir, H. Salimbeni, J. Hensman, and M. Deisenroth. Gaussian process conditional density estimation. In _Advances in Neural Information Processing Systems_, pages 2385-2395, 2018.
* [35] L. T. Fernholz. _Von Mises calculus for statistical functionals_, volume 19. Springer Science & Business Media, 2012.
* [36] R. A. Fisher. Inverse probability. _Mathematical Proceedings of the Cambridge Philosophical Society_, 26(4):528-535, 1930.
* [37] 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-1070, 2010.
* [38] S. Fort, H. Hu, and B. Lakshminarayanan. Deep ensembles: A loss landscape perspective. _arXiv preprint arXiv:1912.02757_, 2019.
* [39] P. E. Freeman, R. Izbicki, and A. B. Lee. A unified framework for constructing, tuning and assessing photometric redshift density estimates in a selection bias setting. _Monthly Notices of the Royal Astronomical Society_, 468(4):4556-4565, 2017.
* [40] 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.
* [41] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin. _Bayesian data analysis_. CRC press, 2013.
* [42] P. Geurts, D. Ernst, and L. Wehenkel. Extremely randomized trees. _Machine learning_, 63:3-42, 2006.
* [43] C. J. Geyer. Practical Markov chain Monte Carlo. _Statistical Science_, 7(4):473-483, 1992.
* [44] P. W. Glynn and D. L. Iglehart. Simulation output analysis using standardized time series. _Mathematics of Operations Research_, 15(1):1-16, 1990.
* [45] 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.
* [46] A. Graves. Practical variational inference for neural networks. _Advances in neural information processing systems_, 24, 2011.
* [47] R. Hable. Asymptotic normality of support vector machine variants and other regularized kernel methods. _Journal of Multivariate Analysis_, 106:92-117, 2012.
* [48] F. R. Hampel. The influence curve and its role in robust estimation. _Journal of the American Statistical Association_, 69(346):383-393, 1974.
* [49] F. R. Hampel. Potential surprises. _International Symposium on Imprecise Probability: Theories and Applications (ISIPTA)_, page 209, 2011.
* [50] B. Hanin and M. Sellke. Approximating continuous functions by relu nets of minimal width. _arXiv preprint arXiv:1710.11278_, 2017.
* [51] B. He, B. Lakshminarayanan, and Y. W. Teh. Bayesian deep ensembles via the neural tangent kernel. _arXiv preprint arXiv:2007.05864_, 2020.
* [52] K. He, X. Zhang, S. Ren, and J. Sun. Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In _Proceedings of the IEEE International Conference on Computer Vision_, pages 1026-1034, 2015.
* [53] J. M. Hernandez-Lobato and R. Adams. Probabilistic backpropagation for scalable learning of bayesian neural networks. In _International Conference on Machine Learning_, pages 1861-1869, 2015.
* [54] M. P. Holmes, A. G. Gray, and C. L. Isbell Jr. Fast nonparametric conditional density estimation. In _Proceedings of the Twenty-Third Conference on Uncertainty in Artificial Intelligence_, pages 175-182, 2007.
* [55] K. Hornik. Approximation capabilities of multilayer feedforward networks. _Neural Networks_, 4(2):251-257, 1991.
* [56] W. Hu, Z. Li, and D. Yu. Simple and effective regularization methods for training on noisily labeled data with generalization guarantee. In _International Conference on Learning Representations_, 2020.
* [57] Z. Huang, H. Lam, and H. Zhang. Quantifying epistemic uncertainty in deep learning. _arXiv preprint arXiv:2110.12122_, 2021.
* [58] E. Hullermeier and W. Waegeman. Aleatoric and epistemic uncertainty in machine learning: An introduction to concepts and methods. _Machine Learning_, 110(3):457-506, 2021.
* [59] R. Izbicki and A. B. Lee. Nonparametric conditional density estimation in a high-dimensional regression setting. _Journal of Computational and Graphical Statistics_, 25(4):1297-1316, 2016.
* [60] 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.
* [61] A. Jacot, F. Gabriel, and C. Hongler. Neural tangent kernel: Convergence and generalization in neural networks. _arXiv preprint arXiv:1806.07572_, 2018.
* [62] L. Jaeckel. The infinitesimal jackknife. memorandum. Technical report, MM 72-1215-11, Bell Lab. Murray Hill, NJ, 1972.
* [63] G. L. Jones, M. Haran, B. S. Caffo, and R. Neath. Fixed-width output analysis for Markov chain Monte Carlo. _Journal of the American Statistical Association_, 101(476):1537-1547, 2006.
* [64] A. Kendall and Y. Gal. What uncertainties do we need in bayesian deep learning for computer vision? _arXiv preprint arXiv:1703.04977_, 2017.
* [65] 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_, 22(3):337-346, 2010.
* [66] A. Khosravi, S. Nahavandi, D. Creighton, and A. F. Atiya. Comprehensive review of neural network-based prediction intervals and new advances. _IEEE Transactions on Neural Networks_, 22(9):1341-1356, 2011.
* [67] R. Koenker and K. F. Hallock. Quantile regression. _Journal of Economic Perspectives_, 15(4):143-156, 2001.
* [68] P. W. Koh, K.-S. Ang, H. H. Teo, and P. Liang. On the accuracy of influence functions for measuring group effects. _arXiv preprint arXiv:1905.13289_, 2019.
* [69] P. W. Koh and P. Liang. Understanding black-box predictions via influence functions. In _International Conference on Machine Learning_, pages 1885-1894. PMLR, 2017.
* [70] 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.
* [71] H. Lam. A cheap bootstrap method for fast inference. _arXiv preprint arXiv:2202.00090_, 2022.
* [72] H. Lam and Z. Liu. Bootstrap in high dimension with low computation. _International Conference on Machine Learning (ICML)_, 2023.
* [73] H. Lam and H. Zhang. Doubly robust stein-kernelized monte carlo estimator: Simultaneous bias-variance reduction and supercanonical convergence. _Journal of Machine Learning Research_, 24(85):1-58, 2023.
* [74] 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.
* [75] 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-8583, 2019.
* [76] 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:1511.06314_, 2015.
* [77] J. Lei, M. G'Sell, A. Rinaldo, R. J. Tibshirani, and L. Wasserman. Distribution-free predictive inference for regression. _Journal of the American Statistical Association_, 113(523):1094-1111, 2018.
* [78] 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.
* [79] J. Lei and L. Wasserman. Distribution-free prediction bands for non-parametric regression. _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_, 76(1):71-96, 2014.
* [80] Y. Li and Y. Liang. Learning overparameterized neural networks via stochastic gradient descent on structured data. _Advances in Neural Information Processing Systems_, 2018.
* [81] Z. Li, R. Wang, D. Yu, S. S. Du, W. Hu, R. Salakhutdinov, and S. Arora. Enhanced convolutional neural tangent kernels. _arXiv preprint arXiv:1911.00809_, 2019.
* [82] N. Littlestone and M. K. Warmuth. The weighted majority algorithm. _Information and computation_, 108(2):212-261, 1994.
* [83] 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 6232-6240, 2017.
* [84] N. Meinshausen. Quantile regression forests. _Journal of Machine Learning Research_, 7(Jun):983-999, 2006.
* [85] R. Michelmore, M. Kwiatkowska, and Y. Gal. Evaluating uncertainty quantification in end-to-end autonomous driving control. _arXiv preprint arXiv:1811.06817_, 2018.
* [86] R. Michelmore, M. Wicker, L. Laurenti, L. 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 (ICRA)_, pages 7344-7350. IEEE, 2020.
* [87] M. Mirza and S. Osindero. Conditional generative adversarial nets. _arXiv preprint arXiv:1411.1784_, 2014.
* [88] M. Mohri, A. Rostamizadeh, and A. Talwalkar. _Foundations of machine learning_. MIT press, 2018.
* [89] T. Pearce, M. Zaki, A. Brintrup, and A. Neely. High-quality prediction intervals for deep learning: A distribution-free, ensembled approach. In _International Conference on Machine Learning, PMLR: Volume 80_, 2018.
* [90] K. Posch and J. Pilz. Correlated parameters to accurately measure uncertainty in deep neural networks. _IEEE Transactions on Neural Networks and Learning Systems_, 32(3):1037-1051, 2020.
* [91] Y. Ren, J. Zhu, J. Li, and Y. Luo. Conditional generative moment-matching networks. _Advances in Neural Information Processing Systems_, 29, 2016.
* [92] Y. Romano, E. Patterson, and E. Candes. Conformalized quantile regression. In _Advances in Neural Information Processing Systems_, pages 3543-3553, 2019.
* [93] 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.
* [94] B. Schmeiser. Batch size effects in the analysis of simulation output. _Operations Research_, 30(3):556-568, 1982.
* [95] L. Schruben. Confidence interval estimation using standardized time series. _Operations Research_, 31(6):1090-1108, 1983.
* [96] R. Senge, S. Bosner, K. Dembczynski, J. Haasenritter, O. Hirsch, N. Donner-Banzhoff, and E. Hullermeier. Reliable classification: Learning classifiers that distinguish aleatoric and epistemic uncertainty. _Information Sciences_, 255:16-29, 2014.
* [97] J. Shao and D. Tu. _The jackknife and bootstrap_. Springer Science & Business Media, 2012.
* [98] A. Siddhant and Z. C. Lipton. Deep bayesian active learning for natural language processing: Results of a large-scale empirical study. _arXiv preprint arXiv:1808.05697_, 2018.
* [99] 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.
* [100] S. Smale and D.-X. Zhou. Shannon sampling II: Connections to learning theory. _Applied and Computational Harmonic Analysis_, 19(3):285-302, 2005.
* [101] S. Smale and D.-X. Zhou. Learning theory estimates via integral operators and their approximations. _Constructive Approximation_, 26(2):153-172, 2007.
* [102] I. Steinwart, A. Christmann, et al. Estimating conditional quantiles with the help of the pinball loss. _Bernoulli_, 17(1):211-225, 2011.
* [103] H. Sun and Q. Wu. Application of integral operator for regularized least-square regression. _Mathematical and Computer Modelling_, 49(1-2):276-285, 2009.
* [104] A. W. Van der Vaart. _Asymptotic statistics_, volume 3. Cambridge university press, 2000.
* [105] V. Vovk, A. Gammerman, and G. Shafer. _Algorithmic learning in a random world_. Springer Science & Business Media, 2005.
* [106] 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.
* [107] 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.
* [108] G. Yang. Tensor programs ii: Neural tangent kernel for any architecture. _arXiv preprint arXiv:2006.14548_, 2020.
* [109] H. Zhang, J. Zimmerman, D. Nettleton, and D. J. Nordman. Random forest prediction intervals. _The American Statistician_, pages 1-15, 2019.
* [110] Y. Zhang, Z.-Q. J. Xu, T. Luo, and Z. Ma. A type of generalization error induced by initialization in deep neural networks. In _Mathematical and Scientific Machine Learning_, pages 144-164. PMLR, 2020.
* [111] 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.
* [112] X. Zhou, Y. Jiao, J. Liu, and J. Huang. A deep generative approach to conditional sampling. _Journal of the American Statistical Association_, pages 1-12, 2022.
* [113] D. Zou, Y. Cao, D. Zhou, and Q. Gu. Gradient descent optimizes over-parameterized deep relu networks. _Machine Learning_, 109(3):467-492, 2020.
## 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 Proposition 3.1 that the shifted kernel ridge regressor using NTK with a shift from an initial function \(s_{\theta_{b}}\) is exactly the linearized neural network regressor that starts from the initial network \(s_{\theta_{b}}\). Section D presents the proofs for results in the paper. Section E presents experimental details and more experimental results. Section F provides additional methodology and results, extending our framework to semi-supervised learning by introducing the independent PNC predictor.

## Appendix A Other Types of Uncertainty

Back in 1930, [36] proposed a formal distinction between aleatory and epistemic uncertainty [49], while in modern machine learning, their distinction and connection were further investigated in [96] and then extended to deep learning models [64, 31]. In [89], the differences between procedural and data variabilities in epistemic uncertainty were intuitively described; however, no rigorous definition or estimation method was provided. We present an additional discussion on other types of uncertainty for the sake of completeness. Section A.1 presents aleatoric uncertainty and Section A.2 presents predictive uncertainty.

### Aleatoric Uncertainty

We note that if we could remove all the epistemic uncertainty, i.e., letting \(\text{UQ}_{EU}=0\), the best predictor we could get is the Bayes predictor \(h_{B}^{*}\). However, \(h_{B}^{*}\), as a point estimator, is not able to capture the randomness in \(\pi_{Y|X}\) if it is not a point mass distribution.

The easiest way to think of aleatoric uncertainty is that it is captured by \(\pi_{Y|X}\). At the level of the realized response or label value, this uncertainty is represented by

\[\text{UQ}_{AU}=y-h_{B}^{*}(x)\]

If the connection between \(X\) and \(Y\) is non-deterministic, the description of a new prediction problem involves a conditional probability distribution \(\pi_{Y|X}\). Standard neural network predictors can only provide a single output \(y\). Thus, even given full information of the distribution \(\pi\), the uncertainty in the prediction of a single output \(y\) remains. This uncertainty cannot be removed by better modeling or more data.

There are multiple work aiming to estimate the aleatoric uncertainty, i.e., to learn \(\pi_{Y|X}\). For instance, conditional quantile regression aims to learn each quantile of the distribution \(\pi_{Y|X}\)[67, 84, 102]. Conditional density estimation aims to approximately describe the density of \(\pi_{Y|X}\)[54, 34, 59, 28, 39, 60]. However, we remark that these approaches also face their own epistemic uncertainty in the estimation. See [25, 102, 10] for recent studies on the epistemic bias in quantile regression.

### Predictive Uncertainty

Existing work on uncertainty measurement in deep learning models mainly focuses on _prediction sets_ (_predictive uncertainty_), which captures the sum of epistemic and aleatoric uncertainties [89, 92, 12, 3, 22].

In certain scenarios, it is not necessary to 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 a map \(\hat{C}:\mathcal{X}\to 2^{\mathcal{Y}}\) which maps an input to a subset of the output space so that for each test point \(x_{0}\in\mathcal{X}\), the prediction set \(\hat{C}(x_{0})\) is likely to cover the true outcome \(y_{0}\)[105, 11, 12, 79, 78, 77]. This
prediction set is also called a prediction interval in the case of regression [65; 66]. The prediction set communicates predictive uncertainty via a statistical guarantee on the marginal coverage, i.e.,

\[\mathbb{P}(y_{0}\in\hat{C}(x_{0}))\geq 1-\delta\]

for a small threshold \(\delta>0\) where the probability is taken with respect to both training data \(\mathcal{D}_{tr}\) (epistemic uncertainty) for learning \(\hat{C}\) and the test data \((x_{0},y_{0})\) (aleatoric uncertainty). It is more tempting to obtain a statistical guarantee with only aleatoric uncertainty by considering the probablity conditional on the prediction set, i.e.,

\[\mathbb{P}(y_{0}\in\hat{C}(x_{0})|\hat{C})\geq 1-\delta. \tag{9}\]

However, this guarantee is in general very hard to achieve in the finite-sample case. Even asymptotically, (9) is not easy to achieve unless we have a very simple structure on the data distribution [93; 109]. 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

### Classical Statistical Learning Framework and Kernel Ridge Regression

Following Section 2, assuming that the input-output pair \((X,Y)\) is a random vector following an unknown probability distribution \(\pi\) on \(\mathcal{X}\times\mathcal{Y}\) where \(X\in\mathcal{X}\subset\mathbb{R}^{d}\) is an input and \(Y\in\mathcal{Y}\subset\mathbb{R}\) is the corresponding output. Let the marginal distribution of \(X\) be \(\pi_{X}(x)\) and the conditional distribution of \(Y\) given \(X\) be \(\pi_{Y|X}(y|x)\).

Suppose a learner has access to a set of training data \(\mathcal{D}=\{(x_{1},y_{1}),(x_{2},y_{2}),...,(x_{n},y_{n})\}\), which is assumed to be independent and identically distributed (i.i.d.) according to the data distribution \(\pi\). Let \(\mathbf{x}=(x_{1},...,x_{n})^{T}\) and \(\mathbf{y}=(y_{1},...,y_{n})^{T}\) for short. Let \(\pi_{n}\) or \(\pi_{\mathcal{D}}\) be the empirical distribution associated with the training data \(\mathcal{D}\):

\[\pi_{\mathcal{D}}=\pi_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{(x_{i},y_{i})}.\]

The goal of the learner is to obtain a function \(h:\mathcal{X}\rightarrow\mathcal{Y}\) such that \(h(x)\) is a "good" predictor for the response \(y\) if \(X=x\) is observed for any \(x\in\mathcal{X}\). This is typically done by minimizing the following (ground-truth) _population risk_

\[R_{\pi}(h):=\mathbb{E}_{(X,Y)\sim\pi}[\mathcal{L}(h(X),Y)]\]

where \(\mathcal{L}:\mathcal{Y}\times\mathbb{R}\rightarrow\mathbb{R}\) is the _loss function_. For instance, \(\mathcal{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_\(h_{B}^{*}\)[88; 58]:

\[h_{B}^{*}(X):=\operatorname*{arg\,min}_{\bar{y}\in\mathcal{Y}}\,\mathbb{E}_{Y \sim\pi_{Y|X}}[\mathcal{L}(\hat{y},Y)|X].\]

Note that \(h_{B}^{*}\) cannot be obtained in practice, since the conditional distribution \(\pi_{Y|X}(y|x)\) is unknown. In particular, the least squares loss \(\mathcal{L}(h(X),Y)=(h(X)-Y)^{2}\) yields the _(ground-truth) regression function_ defined by

\[g_{\pi}^{*}(x):=\mathbb{E}_{(X,Y)\sim\pi}[Y|X=x].\]

As the ground-truth distribution \(\pi\) is unknown, we can neither compute nor minimize the population risk \(R_{\pi}(h)\) directly. We define the risk under a general distribution \(P\),

\[R_{P}(h):=\mathbb{E}_{(X,Y)\sim P}[\mathcal{L}(h(X),Y)].\]

In particular, for \(P=\pi_{n}\) (the empirical distribution associated with the training data \(\mathcal{D}\)), an _empirical risk_ is derived:

\[R_{\pi_{n}}(h):=\mathbb{E}_{(X,Y)\sim\pi_{n}}[\mathcal{L}(h(X),Y)]=\frac{1}{n }\sum_{i=1}^{n}\mathcal{L}(h(x_{i}),y_{i})\]
In practice, minimizing the risk over \(h\) is restricted to a certain function class. Let \(\mathcal{H}\) be a hypothesis class, which is a set of functions \(\{h:\mathcal{X}\to\mathcal{Y}|h\in\mathcal{H}\}\). For instance, 1) \(\mathcal{H}\) could be a nonparametric class such as a _reproducing kernel Hilbert space (RKHS)_. The resulting method is known as _kernel-based regression_; See below. 2) \(\mathcal{H}\) could be a parametric class \(\{h_{\theta}:\mathcal{X}\to\mathcal{Y}|\theta\in\Theta\}\) where \(\theta\) is the parameter, and \(\Theta\) 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 a neural network model.

In classical statistical learning, one is interested in finding a hypothesis \(g_{\pi}\in\mathcal{H}\) that minimizes the population risk

\[g_{\pi}:=\operatorname*{arg\,min}_{h\in\mathcal{H}}\,R_{\pi}(h). \tag{10}\]

which is called the true risk minimizer. We remark that \(h_{\pi}\) is the best choice in the sense of minimizing the risk within the hypothesis set \(\mathcal{H}\) and different choices of \(\mathcal{H}\) may lead to different \(h_{\pi}\). As \(\pi\) is unknown, the learner may consider minimizing the empirical risk:

\[g_{\pi_{n}}:=\operatorname*{arg\,min}_{h\in\mathcal{H}}\,R_{\pi_{n}}(h) \tag{11}\]

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:

\[g_{\pi_{n},\lambda_{n}}:=\operatorname*{arg\,min}_{h\in\mathcal{H}}\,R_{\pi_{ n}}(h)+\lambda_{n}\|h\|_{\mathcal{H}}^{2} \tag{12}\]

and its corresponding population problem

\[g_{\pi,\lambda_{0}}:=\operatorname*{arg\,min}_{h\in\mathcal{H}}\,R_{\pi}(h)+ \lambda_{0}\|h\|_{\mathcal{H}}^{2} \tag{13}\]

which are called the true regularized risk minimizer and the empirical regularized risk minimizer, respectively. Here, \(\lambda_{n}\geq 0\) is the regularization parameter, which may depend on the data size \(n\), and we assume it has a limit \(\lambda_{0}=\lim_{n\to\infty}\lambda_{n}\). In general, the target \(g_{\pi}\) is not equal to \(g_{\pi_{n},\lambda_{n}}\) or \(g_{\pi_{n}}\) (We omit \(0\) in the subscript if \(\lambda=0\)).

The above framework fits classical machine learning approaches such as linear regression or kernel-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 \(g_{\pi_{n},\lambda_{n}}\) or \(h_{\pi_{n}}\) is 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 \(g_{\pi_{n},\lambda_{n}}\) or \(h_{\pi_{n}}\), 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 [14, 101, 73, 99, 103, 14].

The kernel-based regression means that the hypothesis class \(\mathcal{H}\) in statistical learning is chosen to be an RKHS. Formally, a Hilbert space \(\mathcal{H}\) consisting of functions \(h:\mathcal{X}\to\mathbb{R}\) with an inner product \(\langle\cdot,\cdot\rangle:\mathcal{H}\times\mathcal{H}\to\mathbb{R}\) is a RKHS if there exists a symmetric positive definite function \(k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}\), called a (reproducing) kernel, such that for all \(x\in\mathcal{X}\), we have \(k(\cdot,x)\in\mathcal{H}\) and for all \(x\in\mathcal{X}\) and \(h\in\mathcal{H}\), we have \(h(x)=\langle h(\cdot),k(\cdot,x)\rangle\). We use the above \(k\) to denote the kernel associated with \(\mathcal{H}\). Note that for any symmetric positive definite function \(k\), we can naturally construct, from \(k\), an RKHS \(\mathcal{H}\) whose kernel is exactly \(k\)[14].

When \(\mathcal{H}\) is a RKHS, \(g_{\pi_{n},\lambda_{n}}\) in (12) can be efficiently computed for a number of convex loss functions \(\mathcal{L}\). In particular, for the least squares loss \(\mathcal{L}(h(X),Y)=(h(X)-Y)^{2}\), the resulting problem (12) is well known as the kernel ridge regression, and there exists a closed-form expression for both \(g_{\pi_{n},\lambda_{n}}\) and \(g_{\pi,\lambda_{0}}\), as we discuss below.

Formally, the kernel ridge regression problem is given by

\[g_{\pi_{n},\lambda_{n}}(x):=\operatorname*{arg\,min}_{g\in\mathcal{H}}\left\{ \frac{1}{n}\sum_{j=1}^{n}(y_{j}-g(x_{j}))^{2}+\lambda_{n}\|g\|_{\mathcal{H}}^ {2}\right\}\]

where \(\lambda_{n}>0\) is a regularization hyperparameter that may depend on the cardinality of the training data set. There is a closed-form formula for its solution, as stated below.
**Proposition B.1**.: _Let_

\[\mathbf{k}(\mathbf{x},\mathbf{x})=(k(x_{i},x_{j}))_{n\times n}\in\mathbb{R}^{n\times n}\]

_be the kernel Gram matrix and_

\[\mathbf{k}(x,\mathbf{x})=(k(x,x_{1}),\cdots,k(x,x_{n}))^{T}.\]

_Then the (unique) kernel ridge regression solution is given as_

\[g_{\pi_{n},\lambda_{n}}(x)=\mathbf{k}(x,\mathbf{x})^{T}(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_ {n}n\mathbf{I})^{-1}\mathbf{y}\]

It is worth mentioning the linearity of kernel ridge regression: When we have two outputs \(\mathbf{y}=(y_{1},...,y_{n})^{T}\) and \(\mathbf{y}^{\prime}=(y^{\prime}_{1},...,y^{\prime}_{n})^{T}\), we have

\[\mathbf{k}(x,\mathbf{x})^{T}(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_{n}n\mathbf{I})^{- 1}(\mathbf{y}+\mathbf{y}^{\prime})\] \[= \mathbf{k}(x,\mathbf{x})^{T}(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_{n}n\mathbf{I})^{- 1}\mathbf{y}+\mathbf{k}(x,\mathbf{x})^{T}(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_{n}n\mathbf{I})^{-1} \mathbf{y}\]

which means the solution to

\[\operatorname*{arg\,min}_{g\in\mathcal{H}}\left\{\frac{1}{n}\sum_{j=1}^{n}(y_ {j}+y^{\prime}_{j}-g(x_{j}))^{2}+\lambda_{n}\|g\|_{\mathcal{H}}^{2}\right\}\]

is the summation of the solution to

\[\operatorname*{arg\,min}_{g\in\mathcal{H}}\left\{\frac{1}{n}\sum_{j=1}^{n}(y_ {j}-g(x_{j}))^{2}+\lambda_{n}\|g\|_{\mathcal{H}}^{2}\right\}\]

and the solution to

\[\operatorname*{arg\,min}_{g\in\mathcal{H}}\left\{\frac{1}{n}\sum_{j=1}^{n}(y^{ \prime}_{j}-g(x_{j}))^{2}+\lambda_{n}\|g\|_{\mathcal{H}}^{2}\right\}.\]

Define \(\mathcal{O}_{k}:L^{2}(\pi_{X})\to L^{2}(\pi_{X})\) as the kernel integral operator

\[(\mathcal{O}_{k}g)(x):=\int_{\mathcal{X}}k(x,x^{\prime})g(x^{\prime})\pi_{X}( x^{\prime})dx^{\prime},\;x\in\mathcal{X},\;g\in L^{2}(\pi_{X}).\]

where \(L^{2}(\pi_{X}):=\{f:\int_{\mathcal{X}}f(x)^{2}\pi_{X}(x)dx<\infty\}\).

[103] shows that \(\mathcal{O}_{k}\) is a compact and positive self-adjoint linear operator on \(L^{2}(\pi_{X})\). Note that since \(L_{k}\) is positive: \(\mathcal{O}_{k}\geq 0\), we have that \(\mathcal{O}_{k}+\lambda_{0}I\) is strictly positive for any \(\lambda_{0}>0\): \(\mathcal{O}_{k}+\lambda_{0}I>0\), and thus its inverse operator \((\mathcal{O}_{k}+\lambda_{0}I)^{-1}\) exists and is unique. Note that \((\mathcal{O}_{k}+\lambda_{0}I)^{-1}\) is also a linear operator on \(L^{2}(\pi_{X})\).

Next, consider the population risk minimization problem corresponding to \(g_{\pi_{n},\lambda_{n}}\) as follows:

\[g_{\lambda_{0}}:=\operatorname*{arg\,min}_{g\in\mathcal{H}}\left\{\mathbb{E}_ {\pi}[(Y-g(X))^{2}]+\lambda_{0}\|g\|_{\mathcal{H}}^{2}\right\}.\]

It is easy to see that this problem is equivalent to

\[g_{\lambda_{0}}:=\operatorname*{arg\,min}_{g\in\mathcal{H}}\left\{\mathbb{E}_ {\pi}[(g_{\pi}^{*}(X)-g(X))^{2}]+\lambda_{0}\|g\|_{\mathcal{H}}^{2}\right\}.\]

where \(g_{\pi}^{*}(x)=\mathbb{E}_{(X,Y)\sim\pi}[Y|X=x]\) is the (ground-truth) regression function. This problem has the following explicit closed-form expression of the solution (a proof can be found in [26]):

**Proposition B.2**.: _The solution of \(g_{\lambda_{0}}\) is given as \(g_{\lambda_{0}}=(\mathcal{O}_{k}+\lambda_{0}I)^{-1}\mathcal{O}_{k}g_{\pi}^{*}\)._

The linearity of kernel ridge regression also holds for the population-level solution:

\[(g+g^{\prime})_{\lambda_{0}}= (\mathcal{O}_{k}+\lambda_{0}I)^{-1}\mathcal{O}_{k}(g+g^{\prime})_{\pi}^ {*}\] \[= (\mathcal{O}_{k}+\lambda_{0}I)^{-1}\mathcal{O}_{k}(g_{\pi}^{*}+g_ {\pi}^{*})\] \[= (\mathcal{O}_{k}+\lambda_{0}I)^{-1}\mathcal{O}_{k}g_{\pi}^{*}+( \mathcal{O}_{k}+\lambda_{0}I)^{-1}\mathcal{O}_{k}g_{\pi}^{\prime*}\] \[= g_{\lambda_{0}}+g_{\lambda_{0}}^{\prime}\]

for any two functions \(g,g^{\prime}\in L^{2}(\pi_{X})\) since both \((\mathcal{O}_{k}+\lambda_{0}I)^{-1}\) and \(\mathcal{O}_{k}\) are linear operators on \(L^{2}(\pi_{X})\).
### 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 [24, 23, 47].

Let \(\mathcal{H}\) be a generic RKHS with the associated kernel \(k(x,x^{\prime})\). Let \(\|\cdot\|_{\mathcal{H}}\) denote the norm on \(\mathcal{H}\). Note that the feature map \(k_{x}=k(x,\cdot)\) is a function in \(\mathcal{H}\). We consider the connection between the true regularized risk minimizer \(g_{\pi,\lambda_{0}}\) and the empirical regularized risk minimizer \(g_{\pi_{n},\lambda_{n}}\) introduced in (11) and (12) when \(\mathcal{H}\) is a RKHS.

The asymptotic normality of kernel-based regression roughly states that under some mild conditions, \(\sqrt{n}(g_{\pi_{n},\lambda_{n}}-g_{\pi,\lambda_{0}})\) is asymptotically normal as \(n\to\infty\). This result provides the theoretical foundations to conduct asymptotic statistical inference for kernel-based regression, in particular, building asymptotic confidence interval for the ground-truth \(g_{\pi,\lambda_{0}}\).

To be rigorous, we first introduce the following assumptions:

**Assumption B.3**.: Let \((\Omega,\mathcal{A},\mathcal{Q})\) be a probability space. Suppose that \(\mathcal{X}\subset\mathbb{R}^{d}\) is closed and bounded with Borel-\(\sigma\)-algebra \(\mathcal{B}(\mathcal{X})\), and \(\mathcal{Y}\subset\mathbb{R}\) is closed with Borel-\(\sigma\)-algebra \(\mathcal{B}(\mathcal{Y})\). The Borel-\(\sigma\)-algebra of \(\mathcal{X}\times\mathcal{Y}\) is denoted by \(\mathcal{B}(\mathcal{X}\times\mathcal{Y})\). Let \(\mathcal{H}\) be an RKHS with kernel \(k\) and let \(\pi\) be a probability measure on \((\mathcal{X}\times\mathcal{Y},\mathcal{B}(\mathcal{X}\times\mathcal{Y}))\). Assume that the kernel of \(\mathcal{H}\), \(k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}\) is the restriction of an \(m\)-times continuously differentiable kernel \(\tilde{k}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}\) such that \(m>d/2\) and \(k\neq 0\). Let \(\lambda_{0}\in(0,+\infty)\) be any positive number. Suppose that the sequence \(\lambda_{n}\) satisfy that \(\lambda_{n}-\lambda_{0}=o(\frac{1}{\sqrt{n}})\).

**Assumption B.4**.: Let \(\mathcal{L}:\mathcal{Y}\times\mathbb{R}\to[0,+\infty)\) be a loss function satisfying the following conditions:

* \(\mathcal{L}\) is a convex loss, i.e., \(z\mapsto L(y,z)\) is convex in \(z\) for every \(y\in\mathcal{Y}\).
* The partial derivatives \[\mathcal{L}^{\prime}(y,z)=\frac{\partial}{\partial z}\mathcal{L}(y,z),\quad \mathcal{L}^{\prime\prime}(y,z)=\frac{\partial^{2}}{\partial^{2}z}\mathcal{L }(y,z)\] exist for every \((y,z)\in\mathcal{Y}\times\mathbb{R}\).
* The maps \[(y,z)\mapsto\mathcal{L}^{\prime}(y,z),\quad(y,z)\mapsto\mathcal{L}^{\prime \prime}(y,z)\] are continuous.
* There is a \(b\in L^{2}(\pi_{Y})\), and for every \(a\in(0,+\infty)\), there is a \(b^{\prime}_{a}\in L^{2}(\pi_{Y})\) and a constant \(b^{\prime\prime}_{a}\in(0,+\infty)\) such that, for every \(y\in\mathcal{Y}\), \[|\mathcal{L}(y,z)|\leq b(y)+|z|^{p}\ \forall z,\quad\sup_{z\in[-a,a]}| \mathcal{L}^{\prime}(y,z)|\leq b^{\prime}_{a}(y),\quad\sup_{z\in[-a,a]}| \mathcal{L}^{\prime\prime}(y,z)|\leq b^{\prime\prime}_{a}\] where \(p\geq 1\) is a constant.

Note that the conditions on the loss function \(\mathcal{L}\) in Assumption B.4 are satisfied, e.g., for the logistic loss for classification or least-square loss for regression with \(\pi\) such that \(\mathbb{E}_{\pi_{Y}}[Y^{4}]<\infty\). Therefore, Assumption B.4 is reasonable for the kernel ridge regression problem we consider in this work.

**Theorem B.5** (Theorem 3.1 in [47]).: _Suppose that Assumptions B.3 and B.4 hold. Then there is a tight, Borel-measurable Gaussian process \(\mathbb{G}:\Omega\to\mathcal{H},\omega\mapsto\mathbb{G}(\omega)\) such that_

\[\sqrt{n}(g_{\pi_{n},\lambda_{n}}-g_{\pi,\lambda_{0}})\Rightarrow\mathbb{G}\text { in }\mathcal{H}\]

_where \(\Rightarrow\) represents "converges weakly". The Gaussian process \(\mathbb{G}\) is zero-mean; i.e., \(\mathbb{E}[\langle g,\mathbb{G}\rangle_{\mathcal{H}}]=0\) for every \(g\in\mathcal{H}\)._

Note that Theorem B.5 implies that for any \(g\in\mathcal{H}\), the random variable

\[\Omega\to\mathbb{R},\quad\omega\mapsto\langle g,\mathbb{G}(\omega)\rangle_{ \mathcal{H}}\]

has a zero-mean normal distribution. In particular, letting \(g=k_{x}\), the reproducing property of \(k\) implies that,

\[\sqrt{n}(g_{\pi_{n},\lambda_{n}}(x)-g_{\pi,\lambda_{0}}(x))\Rightarrow\mathbb{ G}(x)\]
To obtain the variance of this limiting distribution, denoted as \(\zeta^{2}(x)\), we review the tool of influence functions as follows.

Let \(P\) be a general distribution on the domain \(\mathcal{X}\times\mathcal{Y}\). Let \(z=(z_{x},z_{y})\in\mathcal{X}\times\mathcal{Y}\) be a general point. Suppose \(P_{\varepsilon,z}=(1-\varepsilon)P+\varepsilon\delta_{z}\) where \(\delta_{z}\) denotes the point mass distribution in \(z\). Let \(T\) be a general statistical functional that maps a distribution \(P\) to a real value: \(T:P\to T(P)\). Then the influence function of \(T\) at \(P\) at the point \(z\) is defined as

\[IF(z;T,P)=\lim_{\varepsilon\to 0}\frac{T(P_{\varepsilon,z})-T(P)}{ \varepsilon}=T^{\prime}_{P}(\delta_{z}-P).\]

Under some mild conditions (e.g., \(T\) is \(\rho_{\infty}\)-Hadamard differentiable [35]), the central limit theorem holds: \(\sqrt{n}(T(\pi_{n})-T(\pi))\) is asymptotically normal with mean 0 and variance \(\int_{z}IF^{2}(z;T,\pi)d\pi(z)\) where \(\pi_{n}\) is the empirical distribution associated with \(n\) samples i.i.d. drawn from \(\pi\).

In particular, let \(T(P)\) be the solution to target problem:

\[\min_{g\in\mathcal{H}}\mathbb{E}_{P}[(Y-g(X))^{2}]+\lambda_{P}\|g\|_{\mathcal{ H}}^{2}\]

where \(P=\pi_{n}\) or \(\pi\), \(\lambda_{\pi_{n}}=\lambda_{n}\) and \(\lambda_{\pi}=\lambda_{0}\). Then the proof of Theorem 3.1 in [47] provides the closed-formed expression of the Gateaux derivative \(T^{\prime}_{P}(Q)\) as follows:

\[T^{\prime}_{P}(Q)=-S_{P}^{-1}(\mathbb{E}_{Q}[\mathcal{L}^{\prime}(Y,g_{P, \lambda_{0}}(X))\Phi(X)])\]

for \(Q\in\text{lin}(B_{S})\) where \(\text{lin}(B_{S})\) corresponds to a subset of finite measures on \((\mathcal{X}\times\mathcal{Y},\mathcal{B}(\mathcal{X}\times\mathcal{Y}))\) (See Proposition A.8. in [47]) and \(S_{P}:\mathcal{H}\rightarrow\mathcal{H}\) is defined by

\[S_{P}(f)=2\lambda_{0}f+\mathbb{E}_{P}[\mathcal{L}^{\prime\prime}(Y,g_{P, \lambda_{0}}(X))f(X)\Phi(X)]\]

where \(\Phi\) is the feature map of \(\mathcal{H}\) (e.g., we can simply take \(\Phi(x)(\cdot)=k(x,\cdot)\)). Note that \(S_{P}\) is a continuous linear operator which is invertible (See Proposition A.5. in [47]). In particular, letting \(Q=\delta_{z}-P\), we obtain that

\[T^{\prime}_{P}(\delta_{z}-P)=S_{P}^{-1}(\mathbb{E}_{P}[\mathcal{L}^{\prime}(Y,g_{P,\lambda_{0}}(X))\Phi(X)])-\mathcal{L}^{\prime}(z_{y},g_{P,\lambda_{0}}( z_{x}))S_{P}^{-1}\Phi(X) \tag{14}\]

Note that the above special Gateaux derivative (influence function) of \(T(P)\), \(T^{\prime}_{P}(\delta_{z}-P)\), was initially derived by [24].

The proof of Theorem 3.1 in [47] shows that \(T(P)\) is Hadamard differentiable, and thus the asymptotic normality follows. Using the above expression of the influence function, we conclude that:

\[\sqrt{n}(g_{\pi_{n},\lambda_{n}}(x)-g_{\pi,\lambda_{0}}(x))\Rightarrow\mathcal{ N}(0,\xi^{2}(x))\]

where

\[\xi^{2}(x)=\int_{z\in\mathcal{X}\times\mathcal{Y}}IF^{2}(z;T,\pi)(x)d\pi(z)\]

and \(IF(z;T,\pi)\) is given by (14):

\[IF(z;T,\pi)=T^{\prime}_{\pi}(\delta_{z}-\pi)=S_{\pi}^{-1}(\mathbb{E}_{\pi}[ \mathcal{L}^{\prime}(Y,g_{\pi,\lambda_{0}}(X))\Phi(X)])-\mathcal{L}^{\prime}( z_{y},g_{\pi,\lambda_{0}}(z_{x}))S_{\pi}^{-1}\Phi(X).\]

Summarizing the above discussion, we have Proposition 3.4.

### 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 \(\xi^{2}(x)\) in Section B.2 in the kernel ridge regression. We also show the consistency of infinitesimal jackknife variance estimation in general kernel-based regression. Our consistency result appears new in the literature.

First, we estimate \(\xi^{2}(x_{0})\) as follows:

\[\hat{\xi}^{2}(x_{0})=\frac{1}{n}\sum_{z_{i}\in\mathcal{D}}IF^{2}(z_{i};T,\pi_{ n})(x_{0})\]

where \(\mathcal{D}=\{(x_{1},y_{1}),(x_{2},y_{2}),...,(x_{n},y_{n})\}\) and \(\pi_{n}\) is the empirical distribution associated with the data \(\mathcal{D}\). This estimator is known as the infinitesimal jackknife variance estimator [62]. In this section,
we use \(x_{0}\) instead of \(x\) in the variance or influence function at a test point \(x_{0}\) to avoid confusion between \(x_{0}\) and \(z_{x}\). The latter \(z_{x}\) is referred to the \(x\)-componenet of \(z\).

Based on [24, 29, 47], we derive the closed-form expression of the variance estimation \(\hat{\xi}^{2}(x_{0})\) as follows.

**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 \(\mathcal{L}(h(X),Y)=(h(X)-Y)^{2}\) in kernel ridge regression, \(IF(z;T,\pi_{n})\) is given by_

\[IF(z;T,\pi_{n})(x_{0})=\mathbf{k}(x_{0},\mathbf{x})^{T}(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_ {0}n\mathbf{I})^{-1}M_{z}(\mathbf{x})-M_{z}(x_{0}) \tag{15}\]

_where \(M_{z}(\mathbf{x}):=(M_{z}(x_{1}),...,M_{z}(x_{n}))^{T}\) and_

\[M_{z}(x_{i}):=g_{\pi_{n},\lambda_{0}}(x_{i})-\frac{1}{\lambda_{0}}(z_{y}-g_{ \pi_{n},\lambda_{0}}(z_{x}))k(z_{x},x_{i})\]

_for \(x_{i}=x_{0},x_{1},\cdots,x_{n}\)._

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 Gateaux differentiability) to guarantee that \(IF^{2}(z;T,\pi_{n})\) is indeed "close" to \(IF^{2}(z;T,\pi)\) and \(\frac{1}{n}\sum_{z_{i}\in\mathcal{D}}IF^{2}(z_{i};T,\pi_{n})(x_{0})\) is indeed "close" to \(\int_{z\in\mathcal{X}\times\mathcal{Y}}IF^{2}(z;T,\pi)(x_{0})d\pi(z)\). We show that this is achievable for the infinitesimal 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 B.3 and B.4 hold. Moreover, assume that \(\mathcal{Y}\subset\mathbb{R}\) is bounded, and \(b\) and \(b^{\prime}_{a}\) in Assumption B.4 are bounded on \(\mathcal{Y}\). Then we have_

\[\hat{\xi}^{2}(x_{0})=\frac{1}{n}\sum_{z_{i}\in\mathcal{D}}IF^{2}(z_{i};T,\pi_ {n})(x_{0})\rightarrow\int_{z\in\mathcal{X}\times\mathcal{Y}}IF^{2}(z;T,\pi)( x)d\pi(z)=\xi^{2}(x_{0}),\quad a.s.\]

_as \(n\rightarrow\infty\). Hence, an asymptotically exact \((1-\alpha)\)-level confidence interval of \(g_{\pi,\lambda_{0}}(x_{0})\) is_

\[\left[g_{\pi_{n},\lambda_{n}}(x_{0})-\frac{\hat{\xi}(x_{0})}{n}q_{1-\frac{ \alpha}{2}},g_{\pi_{n},\lambda_{n}}(x_{0})+\frac{\hat{\xi}(x_{0})}{n}q_{1-\frac {\alpha}{2}}\right]\]

_where \(q_{\alpha}\) is the \(\alpha\)-quantile of the standard Gaussian distribution \(\mathcal{N}(0,1)\)._

The proof of Theorems B.6 and B.7 is given in Appendix D.

### PNC-Enhanced Infinitesimal Jackknife

In this section, we return to our problem setting in Section 3.3. We apply the general results in Section B.3 to our PNC predictor and develop the PNC-enhanced infinitesimal jackknife confidence interval for over-parameterized neural networks.

Recall that the statistical functional \(T_{1}\) is associated with the following problem:

\[\hat{h}_{n,\theta^{\text{b}}}(\cdot)-\hat{\phi}_{n,\theta^{\text{b}}}(\cdot)- \bar{s}(\cdot)=\min_{g\in\mathcal{H}}\frac{1}{n}\sum_{i=1}^{n}[(y_{i}-\bar{s} (x_{i})-g(x_{i}))^{2}]+\lambda_{n}\|g\|_{\mathcal{H}}^{2}\]

Consider the PNC-enhanced infinitesimal jackknife variance estimator:

\[\hat{\sigma}_{\theta^{\text{b}}}^{2}(x_{0})=\frac{1}{n}\sum_{z_{i}\in\mathcal{ D}}IF^{2}(z_{i};T_{1},\pi_{n})(x_{0}).\]

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-\bar{s}(X)\). Then \(IF(z;T_{1},\pi_{n})\) is given by_

\[IF(z;T_{1},\pi_{n})(x_{0})=\mathbf{K}(x_{0},\mathbf{x})^{T}(\mathbf{K}(\mathbf{x},\mathbf{x})+ \lambda_{0}n\mathbf{I})^{-1}M_{z}(\mathbf{x})-M_{z}(x_{0}) \tag{16}\]
_where \(M_{z}(\mathbf{x}):=(M_{z}(x_{1}),...,M_{z}(x_{n}))^{T}\) and_

\[M_{z}(x_{i}):=\hat{h}_{n,\theta^{b}}(x_{i})-\hat{\phi}_{n,\theta^{b}}(x_{i})- \bar{s}(x_{i})-\frac{1}{\lambda_{0}}(z_{y}-(\hat{h}_{n,\theta^{b}}(z_{x})-\hat{ \phi}_{n,\theta^{b}}(z_{x})))K(z_{x},x_{i})\]

_for \(x_{i}=x_{0},x_{1},\cdots,x_{n}\)._

Theorem B.8 immediately follows from Theorem B.6.

**Theorem B.9** (Exact coverage of 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 \(Y-\bar{s}(X)\). Moreover, assume that \(\mathcal{Y}\subset\mathbb{R}\) is bounded. Then we have_

\[\hat{\sigma}_{\theta^{b}}^{2}(x_{0})=\frac{1}{n}\sum_{z_{i}\in\mathcal{D}}IF^{ 2}(z_{i};T_{1},\pi_{n})(x_{0})\to\int_{z\in\mathcal{X}\times\mathcal{Y}}IF^{2 }(z;T_{1},\pi)(x)d\pi(z)=\sigma_{\theta^{b}}(x_{0}),\quad a.s.\]

_as \(n\to\infty\). Hence, an asymptotically exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x_{0})\) is_

\[\left[\hat{h}_{n,\theta^{b}}(x_{0})-\hat{\phi}_{n,\theta^{b}}(x_{0})-\frac{ \hat{\sigma}_{\theta^{b}}(x_{0})}{n}q_{1-\frac{a}{2}},\hat{h}_{n,\theta^{b}}( x_{0})-\hat{\phi}_{n,\theta^{b}}(x_{0})+\frac{\hat{\sigma}_{\theta^{b}}(x_{0})}{n}q _{1-\frac{a}{2}}\right]\]

_where \(q_{\alpha}\) is the \(\alpha\)-quantile of the standard Gaussian distribution \(\mathcal{N}(0,1)\) and the computation of \(\hat{\sigma}_{\theta^{b}}(x_{0})\) is given by Theorem B.8._

Note that in our problem setting, we can set \(b(y)=y^{2}\) and \(b^{\prime}_{a}=a+|y|\) for kernel ridge regression. Since \(\mathcal{Y}\) is bounded, \(b\) and \(b^{\prime}_{a}\) are naturally bounded on \(\mathcal{Y}\). Therefore, Theorem B.9 follows from Theorem B.7.

## Appendix C Theory of Neural Tangent Kernel

In this section, we provide a brief review of the theory of neural tangent kernel (NTK). We will then discuss one of the results employed by this paper, i.e.,

_The shifted kernel ridge regressor using NTK with a shift from an initial function \(s_{\theta_{b}}\) is exactly the linearized neural network regressor that starts from the initial network \(s_{\theta_{b}}\)._

NTK [61] has attracted a 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 initialization [80, 33]. An over-parameterized neural network can rapidly reduce training error to zero via gradient descent, and thus finds a global minimum despite the objective function being non-convex [32, 4, 5, 6, 113]. Moreover, the trained network also exhibits good generalization property [19, 20]. These observations are implicitly described by a notion, NTK, suggested by [61]. This kernel is able to characterize the training behavior of sufficiently wide fully-connected neural networks and build 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, 107, 81]. 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-parameterized neural networks is nearly constant and close to its initialization, training networks is very similar to a shifted kernel ridge regression where the feature map of the kernel is given by the gradient of networks [75]. Multiple previous works [7, 75, 51, 56, 110, 15] have established some kind of equivalences between the neural network regressor and the kernel ridge regressor based on a wide variety of assumptions or settings, such as the full-rank assumption of NTK, or zero-valued initialization, or introducing a small multiplier \(\kappa\), or normalized training inputs \(\|x_{i}\|_{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 [75], 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. [75] does not add regularization in the loss. Instead, we introduce a regularization \(\lambda_{n}>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.1) as well as other network specifications. Based on them, we show that Proposition C.3 holds, which provides the starting point for our uncertainty quantification framework. We also provide some additional remarks in Section C.1.

**NTK parameterization.** We consider a fully-connected neural network with any depth defined formally as follows. Suppose the network consists of \(L\) hidden layers. Denote \(g^{(0)}(x)=x\) and \(d_{0}=d\). Let

\[f^{(l)}(x)=W^{(l)}g^{(l-1)}(x),\;g^{(l)}(x)=\sqrt{\frac{c_{\sigma}}{d_{l}}} \sigma(f^{(l)}(x))\]

where \(W^{(l)}\in\mathbb{R}^{d_{l}\times d_{l-1}}\) is the weight matrix in the \(l\)-th layer, \(\sigma\) is a coordinate-wise activation function, \(c_{\sigma}=\mathbb{E}_{z\sim\mathcal{N}(0,1)}[\sigma^{2}(z)]^{-1}\), and \(l\in[L]\). The output of the neural network is

\[f_{\theta}(x):=f^{(L+1)}(x)=W^{(L+1)}g^{(L)}(x)\]

where \(W^{(L+1)}\in\mathbb{R}^{1\times d_{L}}\), and \(\theta=(W^{(1)},...,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 \(\theta\) is \(p\). We randomly initialize all the weights to be i.i.d. \(\mathcal{N}(0,1)\) random variables. In other words, we set \(P_{\theta_{b}}=\mathcal{N}(0,\mathbf{I}_{p})\) and let \(\theta^{b}\) be an instantiation drawn from \(P_{\theta_{b}}\). This initialization method is essentially known as the He initialization [52].

**NTK and Linearized neural network.** With the above NTK parameterization and He initialization, the population NTK expression is defined recursively as follows: [61, 7]: For \(l\in[L]\),

\[\Sigma^{(0)}(x,x^{\prime})=x^{T}x^{\prime}\in\mathbb{R},\]

\[\Lambda^{(l)}(x,x^{\prime})=\begin{pmatrix}\Sigma^{(l-1)}(x,x)&\Sigma^{(l-1)} (x,x^{\prime})\\ \Sigma^{(l-1)}(x,x^{\prime})&\Sigma^{(l-1)}(x^{\prime},x^{\prime})\end{pmatrix} \in\mathbb{R}^{2\times 2},\]

\[\Sigma^{(l)}(x,x^{\prime})=c_{\sigma}\mathbb{E}_{(u,v)\sim\mathcal{N}(0, \Lambda(l))}[\sigma(u)\sigma(v)]\in\mathbb{R}.\]

We also define a derivative covariance:

\[\Sigma^{(l)\prime}(x,x^{\prime})=c_{\sigma}\mathbb{E}_{(u,v)\sim\mathcal{N}(0, \Lambda(l))}[\sigma^{\prime}(u)\sigma^{\prime}(v)]\in\mathbb{R},\]

The final population NTK is defined as

\[K(x,x^{\prime})=\sum_{l=1}^{L+1}\left(\Sigma^{(l-1)}(x,x^{\prime})\prod_{s=l}^ {L+1}\Sigma^{(s)\prime}(x,x^{\prime})\right)\]

where \(\Sigma^{(L+1)\prime}(x,x^{\prime})=1\) for convenience. Let \(\langle\cdot,\cdot\rangle\) be the standard inner product in \(\mathbb{R}^{p}\).

Let \(J(\theta;x):=\nabla_{\theta}f_{\theta}(x)\in\mathbb{R}^{1\times p}\) denote the gradient of the network. The empirical (since the initialization is random) NTK matrix is defined as

\[K_{\theta}(x,x^{\prime})=\langle J(\theta;x);J(\theta;x^{\prime})\rangle.\]

In practice, we may use the empirical NTK matrix to 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 [61, 107, 108, 7]:

\[K_{\theta}(x,x^{\prime})=\langle J(\theta;x),J(\theta;x)\rangle\to K(x,x^{ \prime})\]

where \(\rightarrow\) represents "converge in probability" [7] or "almost surely" [107].
The NTK theory shows that \(J(\theta;x)=\nabla_{\theta}f_{\theta}(x)\) is approximately a constant independent of the network parameter \(\theta\) (but still depends on \(x\)) when the network is sufficiently wide. Therefore, when we treat \(J(\theta;x)\) as a constant independent of the network parameter \(\theta\), we can obtain a model obtained from the first-order Taylor expansion of the network around its initial parameters, which is called a linearized neural network [75]. 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(\theta;x)=\nabla_{\theta}f_{\theta}(x)\in\mathbb{R}^{1\times p}\) is independent of \(\theta\) during network training, which is thus denoted as \(J(x)\). The population NTK is then given by \(K(x,x^{\prime})=\langle J(x),J(x)\rangle=J(x)J(x)^{\top}\).

We introduce the NTK vector/matrix as follows. For \(\mathbf{x}=(x_{1},...,x_{j})^{\top}\), we let

\[\mathbf{K}(\mathbf{x},\mathbf{x}):=(K(x_{i},x_{j}))_{i,j=1,...,n}=J(\mathbf{x})J(\mathbf{x})^{\top }\in\mathbb{R}^{n\times n}\]

be the NTK Gram matrix evaluated on data \(\mathbf{x}\). For an arbitrary point \(x\in\mathbb{R}^{d}\), we let \(\mathbf{K}(x,\mathbf{x})\in\mathbb{R}^{n}\) be the kernel value evaluated between the point \(x\) and data \(\mathbf{x}\), i.e.,

\[\mathbf{K}(x,\mathbf{x}):=(K(x,x_{1}),K(x,x_{2}),...,K(x,x_{n}))^{T}=J(x)J(\mathbf{x})^{ \top}\in\mathbb{R}^{1\times n}.\]

**Training dynamics of the linearized network.** Based on Assumption C.1, the training dynamics of the linearized network can be derived as follows.

We consider the regression problem with the following regularized empirical loss

\[\hat{R}(f_{\theta})=\frac{1}{n}\sum_{i=1}^{n}\mathcal{L}(f_{\theta}(x_{i}),y_{ i})+\lambda_{n}\|\theta-\theta^{b}\|_{2}^{2}. \tag{17}\]

where \(\theta^{b}=\theta(0)\) is the initialization of the network. In the following, we use \(f_{\theta(t)}\) to represent the network with parameter \(\theta(t)\) that evolves with the time \(t\). Let \(\eta\) be the learning rate. Suppose the network parameter is trained via continuous-time gradient flow. Then the evolution of the parameters \(\theta(t)\) and \(f_{\theta(t)}\) can be written as

\[\frac{d\theta(t)}{dt}=-\eta\left(\frac{1}{n}\nabla_{\theta}f_{\theta(t)}(\mathbf{ x})^{\top}\nabla_{f_{\theta(t)}(\mathbf{x})}\mathcal{L}+2\lambda_{n}(\theta(t)- \theta(0))\right) \tag{18}\]

\[\frac{df_{\theta(t)}(\mathbf{x})}{dt}=\nabla_{\theta}f_{\theta(t)}(\mathbf{x})\frac{d \theta(t)}{dt} \tag{19}\]

where \(\theta(t)\) and \(f_{\theta(t)}(x)\) are the network parameters and network output at time \(t\).

Under Assumption C.1, we have that

\[f_{\theta(t)}(x)=f_{\theta(0)}(x)+J(x)(\theta(t)-\theta(0)),\]

which gives

\[\nabla_{\theta}f_{\theta(t)}(x)=J(x).\]

Suppose the loss function is an MSE loss, i.e., \(\mathcal{L}(\tilde{y},y)=(\tilde{y}-y)^{2}\). Then

\[\nabla_{f_{\theta(t)}(\mathbf{x})}\mathcal{L}(f_{\theta(t)}(\mathbf{x}),\mathbf{y})=2(f_{ \theta(t)}(\mathbf{x})-\mathbf{y})\]

Therefore, under Assumption C.1, (18) becomes

\[\frac{d\theta(t)}{dt} =-\eta\left(\frac{1}{n}J(\mathbf{x})^{\top}\nabla_{f_{\theta(t)}(\bm {x})}\mathcal{L}+2\lambda_{n}(\theta(t)-\theta(0))\right)\] \[=-\eta\left(\frac{2}{n}J(\mathbf{x})^{\top}\left(f_{\theta(0)}(\mathbf{x} )+J(\mathbf{x})(\theta(t)-\theta(0))-\mathbf{y}\right)+2\lambda_{n}(\theta(t)-\theta( 0))\right)\] \[=-2\eta\left(\frac{1}{n}J(\mathbf{x})^{\top}\left(f_{\theta(0)}(\mathbf{x })-J(\mathbf{x})\theta(0)-\mathbf{y}\right)-\lambda_{n}\theta(0)+\left(\frac{1}{n}J( \mathbf{x})^{\top}J(\mathbf{x})+\lambda_{n}\mathbf{I}\right)\theta(t)\right)\]
Note that \(\frac{1}{n}J(\mathbf{x})^{\top}\left(f_{\theta(0)}(\mathbf{x})-J(\mathbf{x})\theta(0)-\mathbf{y} \right)-\lambda_{n}\theta(0)\) and \(\left(\frac{1}{n}J(\mathbf{x})^{\top}J(\mathbf{x})+\lambda_{n}\mathbf{I}\right)\) are both independent of \(t\). Solving this ordinary differential equation, we obtain

\[\theta(t)= \theta(0)-\frac{1}{n}J(\mathbf{x})^{\top}\left(\frac{1}{n}J(\mathbf{x})^{ \top}J(\mathbf{x})+\lambda_{n}\mathbf{I}\right)^{-1}\] \[\left(\mathbf{I}-\exp\left(-2\eta t\left(\frac{1}{n}J(\mathbf{x})^{\top}J (\mathbf{x})+\lambda_{n}\mathbf{I}\right)\right)\right)\left(f_{\theta(0)}(\mathbf{x})- \mathbf{y}\right)\]

Hence the network output at time \(t\) becomes

\[f_{\theta(t)}(x)-f_{\theta(0)}(x)\] \[= J(x)(\theta(t)-\theta(0))\] \[= -\frac{1}{n}J(x)J(\mathbf{x})^{\top}\left(\frac{1}{n}J(\mathbf{x})^{\top }J(\mathbf{x})+\lambda_{n}\mathbf{I}\right)^{-1}\] \[= -\frac{1}{n}\mathbf{K}(x,\mathbf{x})\left(\frac{1}{n}\mathbf{K}(\mathbf{x},\mathbf{x })+\lambda_{n}\mathbf{I}\right)^{-1}\left(\mathbf{I}-\exp\left(-2\eta t\left(\frac{1} {n}\mathbf{K}(\mathbf{x},\mathbf{x})+\lambda_{n}\mathbf{I}\right)\right)\right)\left(f_{\theta (0)}(\mathbf{x})-\mathbf{y}\right)\]

where the last inequality used the notation \(\mathbf{K}(x,\mathbf{x})=J(x)J(\mathbf{x})^{\top}\) and \(\mathbf{K}(x_{0},\mathbf{x})=J(\mathbf{x})J(\mathbf{x})^{\top}\).

Therefore, with sufficient time of training (\(t\to\infty\)), the final trained network is

\[f_{\theta(\infty)}(x) =\lim_{t\to\infty}f_{\theta(t)}(x)\] \[=f_{\theta(0)}(x)-\frac{1}{n}\mathbf{K}(x,\mathbf{x})\left(\frac{1}{n}\bm {K}(\mathbf{x},\mathbf{x})+\lambda_{n}\mathbf{I}\right)^{-1}\left(f_{\theta(0)}(\mathbf{x})- \mathbf{y}\right)\] \[=f_{\theta(0)}(x)+\mathbf{K}(x,\mathbf{x})\left(\mathbf{K}(\mathbf{x},\mathbf{x})+n \lambda_{n}\mathbf{I}\right)^{-1}\left(\mathbf{y}-f_{\theta(0)}(\mathbf{x})\right)\] \[=s_{\theta^{b}}(x)+\mathbf{K}(x,\mathbf{x})\left(\mathbf{K}(\mathbf{x},\mathbf{x})+n \lambda_{n}\mathbf{I}\right)^{-1}\left(\mathbf{y}-s_{\theta^{b}}(\mathbf{x})\right)\]

where in the last inequality, we use \(s_{\theta^{b}}\) to represent the network initialized with the parameter \(\theta^{b}\) as in the main paper. Summarizing the above discussion, we conclude that

**Assumption C.2**.: Suppose that the network training is specified as follows:

1. The network adopts the NTK parametrization and its parameters are randomly initialized using He initialization.
2. The network is sufficiently (infinitely) wide so that the linearized neural network assumption (Assumption C.1) holds.
3. The network is trained using the loss function in (1) via continuous-time gradient flow by feeding the entire training data and using sufficient training time (\(t\to\infty\)).

**Proposition C.3**.: _Suppose that Assumption C.2 holds. Then the final trained network is given by_ \[\hat{h}_{n,\theta^{b}}(x)=s_{\theta^{b}}(x)+\mathbf{K}(x,\mathbf{x})\left(\mathbf{K}(\bm {x},\mathbf{x})+\lambda_{n}n\mathbf{I}\right)^{-1}\left(\mathbf{y}-s_{\theta^{b}}(\mathbf{x}) \right).\] (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 a shifted kernel ridge regression based on the NTK as follows.

Let \(\bar{\mathcal{H}}\) be the reproducing kernel Hilbert space constructed from the kernel function \(K(x,x^{\prime})\) (the population NTK); See Appendix B. Consider the following kernel ridge regression problem on the RKHS \(\bar{\mathcal{H}}\):

\[s_{\theta^{b}}+\operatorname*{arg\,min}_{g\in\bar{\mathcal{H}}}\frac{1}{n}\sum _{i=1}^{n}(y_{i}-s_{\theta^{b}}(x_{i})-g(x_{i}))^{2}+\lambda_{n}\|g\|_{\bar{ \mathcal{H}}}^{2} \tag{21}\]

where \(\lambda_{n}\) is a regularization hyper-parameter and \(s_{\theta^{b}}\) is a known given function.
**Lemma C.4**.: (20) _is the solution to the following kernel ridge regression problem_

\[s_{\theta^{b}}+\operatorname*{arg\,min}_{g\in\bar{\mathcal{H}}}\,\frac{1}{n}\sum _{i=1}^{n}(y_{i}-s_{\theta^{b}}(x_{i})-g(x_{i}))^{2}+\lambda_{n}\|g\|_{\bar{ \mathcal{H}}}^{2}. \tag{22}\]

Proof of Lemma c.4.: We first let \(\tilde{y}_{i}=y_{i}-s_{\theta^{b}}(x_{i})\) be the shifted label. Now consider the kernel ridge regression problem

\[\operatorname*{arg\,min}_{g\in\bar{\mathcal{H}}}\,\frac{1}{n}\sum_{i=1}^{n}( \tilde{y}_{i}-g(x_{i}))^{2}+\lambda_{n}\|g\|_{\bar{\mathcal{H}}}^{2}.\]

Standard theory in kernel ridge regression (Proposition B.1) shows that the closed-form solution to the above problem is given by:

\[\mathbf{K}(x,\mathbf{x})^{T}(\mathbf{K}(\mathbf{x},\mathbf{x})+\lambda_{n}n\mathbf{I})^{-1}\tilde{\bm {y}}=\mathbf{K}(x,\mathbf{x})^{T}(\mathbf{K}(\mathbf{x},\mathbf{x})+\lambda_{n}n\mathbf{I})^{-1}(\mathbf{ y}-s_{\theta^{b}}(\mathbf{x}))\]

as \(K\) is the kernel function of \(\bar{\mathcal{H}}\). This equation corresponds to the definition of \(\hat{h}_{n,\theta^{b}}(x)-s_{\theta^{b}}(x)\) in (20). 

From Proposition C.3 and Lemma C.4, we immediately conclude that

_The shifted kernel ridge regressor using NTK with a shift from an initial function \(s_{\theta_{b}}\) is exactly the linearized neural network regressor that starts from the initial network \(s_{\theta_{b}}\)._

### Additional Remarks.

We make the following additional remarks:

1. We introduce the regularization parameters \(\lambda_{n}>0\) in this paper, which has some advantages: 1) A common assumption in NTK theory is to assume that the smallest eigenvalue of the NTK Gram matrix is bounded away from zero [32, 7, 75]. With \(\lambda_{n}>0\) is introduced, there is an additional term \(\lambda_{n}n\mathbf{I}\) in the NTK Gram matrix so that this assumption is valid naturally. 2) It helps to stabilize the inversion of the NTK Gram matrix if any computation is needed. 3) It can be naturally linked to the well-known kernel ridge regression, which we leverage in this work to develop uncertainty quantification approaches. On the other hand, we note that with \(\lambda_{0}=\lim_{n\to\infty}\lambda_{n}>0\), the "best" predictor \(h^{*}\) in the neural network is not exactly the same as the ground-truth regression function \(g^{*}_{\pi}:=\mathbb{E}_{(X,Y)\sim\pi}[Y|X=x]\). In fact, Proposition B.2 shows that \(h^{*}=(\mathcal{O}_{K}+\lambda_{0}I)^{-1}\mathcal{O}_{K}g^{*}_{\pi}\neq g^{*} _{\pi}\). However, previous work [103, 101] shows that with some mild conditions, \(h^{*}-g^{*}_{\pi}\to 0\) in some sense when \(\lambda_{0}\to 0\). Therefore, we may use a very small \(\lambda_{0}\) in practice to make the error \(h^{*}-g^{*}_{\pi}\) negligible, as we did in our experiments.
2. Proposition 3.1 (Proposition C.3) is the theoretical foundation of this paper to develop our framework for efficient uncertainty quantification and reduction for neural networks. We recognize the limitation of this proposition: The linearized neural network assumption (Assumption C.1) must hold to 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 a certain order but still involves some unknown constants that hide beneath the data and network [32, 7, 75], which is thus extremely difficult to quantify precisely in practice. Hence, this work does not deal with such an error 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.

## Appendix D Proofs

In this section, we provide technical proof of the results in the paper.
Proof of Theorem 3.2.: In terms of estimation bias, we note that

\[\mathbb{E}[\hat{h}_{n}^{m}(x)|\mathcal{D}]=\frac{1}{m}\sum_{i=1}^{m}\mathbb{E}[ \hat{h}_{n,\theta_{i}^{b}}(x)|\mathcal{D}]=\mathbb{E}[\hat{h}_{n,\theta^{b}}(x )|\mathcal{D}]=\hat{h}_{n}^{*}(x).\]

as shown in (3). This implies that the deep ensemble predictor and the single network predictor have the same conditional procedural mean given the data. In terms of variance, note that \(\hat{h}_{n,\theta_{i}^{b}}(x),...,\hat{h}_{n,\theta_{m}^{b}}(x)\), derived in (2), are conditionally independent given \(\mathcal{D}\). Therefore the deep ensemble predictor gives

\[\text{Var}(\hat{h}_{n}^{m}(x)|\mathcal{D})=\frac{1}{m^{2}}\sum_{i=1}^{m}\text{ Var}(\hat{h}_{n,\theta_{i}^{b}}(x)|\mathcal{D})=\frac{1}{m}\text{Var}(\hat{h}_{n, \theta^{b}}(x)|\mathcal{D})\]

while the single model prediction gives \(\text{Var}(\hat{h}_{n,\theta^{b}}(x)|\mathcal{D})\). Now using conditioning, we have

\[\text{Var}(\hat{h}_{n}^{m}(x))\] \[= \text{Var}(\mathbb{E}[\hat{h}_{n}^{m}(x)|\mathcal{D}])+\mathbb{E }[\text{Var}(\hat{h}_{n}^{m}(x)|\mathcal{D})]\] \[= \text{Var}(\hat{h}_{n}^{*}(x))+\frac{1}{m}\mathbb{E}[\text{Var}( \hat{h}_{n,\theta^{b}}(x)|\mathcal{D})]\] \[\leq \text{Var}(\hat{h}_{n}^{*}(x))+\mathbb{E}[\text{Var}(\hat{h}_{n, \theta^{b}}(x)|\mathcal{D})]\] \[= \text{Var}(\hat{h}_{n,\theta^{b}}(x))\]

as desired. 

Proof of Theorem 3.3.: Recall that \(\hat{\phi}^{\prime}_{n,\theta^{b}}(\cdot)\) is obtained by training an auxiliary neural network with data \(\{(x_{1},\bar{s}(x_{1})),(x_{2},\bar{s}(x_{2})),...,(x_{n},\bar{s}(x_{n}))\}\) and the initialization parameters \(\theta^{b}\). Then Proposition 3.1 immediately shows that

\[\hat{\phi}^{\prime}_{n,\theta^{b}}(x)=s_{\theta^{b}}(x)+\mathbf{K}(x,\mathbf{x})^{T}( \mathbf{K}(\mathbf{x},\mathbf{x})+\lambda_{n}n\mathbf{I})^{-1}(\bar{s}(\mathbf{x})-s_{\theta^{b}} (\mathbf{x})).\]

Note that this auxiliary network starts from the same initialization \(\theta^{b}\) as in \(\hat{h}_{n,\theta^{b}}(x)\). Then we have

\[\hat{\phi}_{n,\theta^{b}}(x) =\hat{\phi}^{\prime}_{n,\theta^{b}}(x)-\bar{s}(x)\] \[=s_{\theta^{b}}(x)-\bar{s}(x)+\mathbf{K}(x,\mathbf{x})^{T}(\mathbf{K}(\mathbf{x}, \mathbf{x})+\lambda_{n}n\mathbf{I})^{-1}(\bar{s}(\mathbf{x})-s_{\theta^{b}}(\mathbf{x}))\] \[=\hat{h}_{n,\theta^{b}}(x)-\hat{h}_{n}^{*}(x)\]

This shows that

\[\hat{h}_{n}^{*}(x)=\hat{h}_{n,\theta^{b}}(x)-(\hat{\phi}^{\prime}_{n,\theta^{b }}(x)-\bar{s}(x))=\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)\]

as desired. 

Proof of Theorem 3.5.: Theorem 3.3 implies that

\[\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)=\hat{h}_{n}^{*}(x)\]

which is the solution to the following problem

\[\hat{h}_{n}^{*}(\cdot)=\bar{s}(\cdot)+\min_{g\in\mathcal{H}}\frac{1}{n}\sum_{i =1}^{n}[(y_{i}-\bar{s}(x_{i})-g(x_{i}))^{2}]+\lambda_{n}\|g\|_{\mathcal{H}}^{2}\]

Therefore, its corresponding population risk minimization problem is:

\[h^{*}(\cdot)=\bar{s}(\cdot)+\min_{g\in\mathcal{H}}\mathbb{E}_{\pi}[(Y-\bar{s}( X)-g(X))^{2}]+\lambda_{0}\|g\|_{\mathcal{H}}^{2}\]

Applying Proposition 3.4 to \(\hat{h}_{n}^{*}\) and \(h^{*}\), we have that

\[\sqrt{n}\left((\hat{h}_{n}^{*}(x)-\bar{s}(x))-(h^{*}(x)-\bar{s}(x))\right) \Rightarrow\mathcal{N}(0,\sigma_{\theta^{b}}^{2}(x))\]
In other words,

\[\sqrt{n}\left(\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)-h^{*}(x) \right)\Rightarrow\mathcal{N}(0,\sigma_{\theta^{b}}^{2}(x))\]

This shows that asymptotically as \(n\rightarrow\infty\) we have

\[\mathbb{P}\left(-q_{1-\frac{\alpha}{2}}\leq\frac{\hat{h}_{n,\theta^{b}}(x)- \hat{\phi}_{n,\theta^{b}}(x)-h^{*}(x)}{\sigma_{\theta^{b}}(x)/\sqrt{n}}\leq q_ {1-\frac{\alpha}{2}}\right)\to 1-\alpha\]

where \(q_{\alpha}\) is the \(\alpha\)-quantile of standard Gaussian distribution \(\mathcal{N}(0,1)\). Hence

\[\left[\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)-\frac{\sigma_{ \theta^{b}}(x)}{\sqrt{n}}q_{1-\frac{\alpha}{2}},\hat{h}_{n,\theta^{b}}(x)- \hat{\phi}_{n,\theta^{b}}(x)+\frac{\sigma_{\theta^{b}}(x)}{\sqrt{n}}q_{1- \frac{\alpha}{2}}\right]\]

is an asymptotically exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x)\). 

Proof of Theorem 3.6.: Note that the statistics

\[\hat{h}_{n^{\prime},\theta^{b}}^{j}(x)-\hat{\phi}_{n^{\prime},\theta^{b}}^{j} (x)\]

from the \(j\)-th batch is i.i.d. for \(j\in[m^{\prime}]\). Moreover, by Theorem 3.5, we have the asymptotic normality:

\[\sqrt{n^{\prime}}\left(\hat{h}_{n^{\prime},\theta^{b}}(x)-\hat{\phi}_{n^{ \prime},\theta^{b}}(x)-h^{*}(x)\right)\Rightarrow\mathcal{N}(0,\sigma_{\theta ^{b}}^{2}(x)),\]

as \(n\rightarrow\infty\) meaning \(n^{\prime}=n/m^{\prime}\rightarrow\infty\). Therefore by the property of Gaussian distribution and the principle of batching, we have

\[\frac{\sqrt{n^{\prime}}}{\sigma_{\theta^{b}}(x)}(\psi_{B}(x)-h^{*}(x)) \Rightarrow\frac{1}{m^{\prime}}\sum_{i=1}^{m^{\prime}}Z_{i}\]

and

\[\frac{\psi_{B}(x)-h^{*}(x)}{S_{B}(x)/\sqrt{m^{\prime}}}\Rightarrow\frac{ \frac{1}{m^{\prime}}\sum_{i=1}^{m^{\prime}}Z_{i}}{\sqrt{\frac{1}{m^{\prime}( m^{\prime}-1)}\sum_{j=1}^{m^{\prime}}(Z_{j}-\frac{1}{m^{\prime}}\sum_{i=1}^{m^{ \prime}}Z_{i})^{2}}}\]

for i.i.d. \(\mathcal{N}(0,1)\) random variables \(Z_{1},...,Z_{m^{\prime}}\), where we use the continuous mapping theorem to deduce the weak convergence. Note that

\[\frac{\frac{1}{m^{\prime}}\sum_{i=1}^{m^{\prime}}Z_{i}}{\sqrt{\frac{1}{m^{ \prime}(m^{\prime}-1)}\sum_{j=1}^{m^{\prime}}(Z_{j}-\frac{1}{m^{\prime}}\sum_ {i=1}^{m^{\prime}}Z_{i})^{2}}}\stackrel{{\text{\small d}}}{{=}}t_ {m^{\prime}-1}\]

Hence, asymptotically as \(n\rightarrow\infty\) we have

\[\mathbb{P}\left(-q_{1-\frac{\alpha}{2}}\leq\frac{\psi_{B}(x)-h^{*}(x)}{S_{B}( x)/\sqrt{m^{\prime}}}\leq q_{1-\frac{\alpha}{2}}\right)\to 1-\alpha\]

where \(q_{\alpha}\) is the \(\alpha\)-quantile of the t distribution \(t_{m^{\prime}-1}\) with degree of freedom \(m^{\prime}-1\). Hence

\[\left[\psi_{B}(x)-\frac{S_{B}(x)}{\sqrt{m^{\prime}}}q_{1-\frac{\alpha}{2}}, \psi_{B}(x)+\frac{S_{B}(x)}{\sqrt{m^{\prime}}}q_{1-\frac{\alpha}{2}}\right]\]

is an asymptotically exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x)\). 

Proof of Theorem 3.7.: Note that the statistics

\[\hat{h}_{n^{\prime},\theta^{b}}^{*b}(x)-\hat{\phi}_{n^{\prime},\theta^{b}}^{*b} (x)\]

from the \(b\)-th replication is i.i.d. for \(b\in[R]\) conditional on the dataset \(\mathcal{D}\). By Theorem 3.5, we have the asymptotic normality:

\[\sqrt{n}\left(\hat{h}_{n,\theta^{b}}(x)-\hat{\phi}_{n,\theta^{b}}(x)-h^{*}(x) \right)\Rightarrow\mathcal{N}(0,\sigma_{\theta^{b}}^{2}(x)),\]
as \(n\to\infty\). By Theorem 2 in [23], we have the following asymptotic normality:

\[\sqrt{n}\left(\hat{h}_{n,\theta^{8}}^{sb}(x)-\hat{\phi}_{n,\theta^{8}}^{sb}(x)- \left(\hat{h}_{n,\theta^{8}}(x)-\hat{\phi}_{n,\theta^{8}}(x)\right)\right) \Rightarrow\mathcal{N}(0,\sigma_{\theta^{8}}^{2}(x)),\]

as \(n\to\infty\) conditional on the data \(\mathcal{D}\). Therefore Assumption 1 in [71] is satisfies and thus Theorem 1 in [71] holds, showing that asymptotically as \(n\to\infty\), we have

\[\mathbb{P}\left(-q_{1-\frac{n}{2}}\leq\frac{\psi_{C}(x)-h^{*}(x)}{S_{C}(x)} \leq q_{1-\frac{n}{2}}\right)\to 1-\alpha\]

where \(q_{\alpha}\) is the \(\alpha\)-quantile of the t distribution \(t_{R}\) with degree of freedom \(R\). Hence

\[\left[\psi_{C}(x)-S_{C}(x)q_{1-\frac{n}{2}},\psi_{C}(x)+S_{C}(x)q_{1-\frac{n}{ 2}}\right]\]

is an asymptotically exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x)\). 

Proof of Theorem b.6.: Recall that \(\Phi\) is the feature map associated with the kernel \(k\) of the RKHS \(\mathcal{H}\). We apply the influence function formula in [24] (see also [29, 47]) to calculate the infinitesimal jackknife:

\[IF(z;T,\pi_{n})=-S_{\pi_{n}}^{-1}(2\lambda_{0}g_{\pi_{n},\lambda_{0}})+ \mathcal{L}^{\prime}(z_{y}-g_{\pi_{n},\lambda_{0}}(z_{x}))S_{\pi_{n}}^{-1} \Phi(z_{x})\]

where \(S_{P}:\mathcal{H}\to\mathcal{H}\) is defined by

\[S_{P}(f)=2\lambda_{0}f+\mathbb{E}_{P}[\mathcal{L}^{\prime\prime}(Y,g_{P, \lambda_{0}}(X))f(X)\Phi(X)]\]

To compute this, we need to obtain \(S_{\pi_{n}}\) and \(S_{\pi_{n}}^{-1}\). Since the loss function is \(\mathcal{L}(\hat{y},y)=(\hat{y}-y)^{2}\), we have

\[S_{\pi_{n}}(f)=2\lambda_{0}f+\mathbb{E}_{\pi_{n}}[\mathcal{L}^{\prime\prime}( Y,g_{\pi_{n},\lambda_{0}}(X))f(X)\Phi(X)]=2\lambda_{0}f+\frac{2}{n}\sum_{j=1}^{n}f (x_{j})\Phi(x_{j}).\]

Suppose \(S_{\pi_{n}}^{-1}(2\lambda_{0}g_{\pi_{n},\lambda_{0}})=\tilde{g}_{1}\). Then at \(x_{0}\), we have

\[2\lambda_{0}g_{\pi_{n},\lambda_{0}}(x_{0})=S_{\pi_{n}}(\tilde{g}_{1}(x_{0}))= 2\lambda_{0}\tilde{g}_{1}(x_{0})+\frac{2}{n}\sum_{j=1}^{n}\tilde{g}_{1}(x_{j} )k(x_{0},x_{j}).\]

Hence,

\[\tilde{g}_{1}(x_{0})=g_{\pi_{n},\lambda_{0}}(x_{0})-\frac{1}{\lambda_{0}n} \sum_{j=1}^{n}\tilde{g}_{1}(x_{j})k(x_{0},x_{j}).\]

This implies that we need to evaluate \(\tilde{g}_{1}(x_{j})\) on training data first, which is straightforward by letting \(x_{0}=x_{1},...,x_{n}\):

\[\tilde{g}_{1}(\mathbf{x})=g_{\pi_{n},\lambda_{0}}(\mathbf{x})-\frac{1}{\lambda_{0}n} \mathbf{k}(\mathbf{x},\mathbf{x})\tilde{g}_{1}(\mathbf{x})\]

so

\[\tilde{g}_{1}(\mathbf{x})=(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_{0}n\mathbf{I})^{-1}(\lambda _{0}n\mathbf{I})g_{\pi_{n},\lambda_{0}}(\mathbf{x})\]

and

\[\tilde{g}_{1}(x_{0}) =g_{\pi_{n},\lambda_{0}}(x_{0})-\frac{1}{\lambda_{0}n}\mathbf{k}(x_{0 },\mathbf{x})^{T}\tilde{g}_{1}(\mathbf{x})\] \[=g_{\pi_{n},\lambda_{0}}(x_{0})-\mathbf{k}(x_{0},\mathbf{x})^{T}(\mathbf{k}( \mathbf{x},\mathbf{x})+\lambda_{0}n\mathbf{I})^{-1}g_{\pi_{n},\lambda_{0}}(\mathbf{x})\]

Next we compute \(S_{\pi_{n}}^{-1}\Phi(z_{x})=\tilde{g}_{2}\). At \(x_{0}\), we have

\[k(z_{x},x_{0})=\Phi(z_{x})(x_{0})=S_{\pi_{n}}(\tilde{g}_{2}(x_{0}))=2\lambda_{0 }\tilde{g}_{2}(x_{0})+\frac{2}{n}\sum_{j=1}^{n}\tilde{g}_{2}(x_{j})k(x_{0},x_{ j})\]

Hence

\[\tilde{g}_{2}(x_{0})=\frac{1}{2\lambda_{0}}k(z_{x},x_{0})-\frac{1}{\lambda_{0} n}\sum_{j=1}^{n}\tilde{g}_{2}(x_{j})k(x_{0},x_{j})\]
This implies that we need to evaluate \(\tilde{g}_{2}(x_{j})\) on training data first, which is straightforward by letting \(x_{0}=x_{1},...,x_{n}\):

\[\tilde{g}_{2}(\mathbf{x})=\frac{1}{2\lambda_{0}}\mathbf{k}(z_{x},\mathbf{x})-\frac{1}{ \lambda_{0}n}\mathbf{k}(\mathbf{x},\mathbf{x})\tilde{g}_{2}(\mathbf{x})\]

so

\[\tilde{g}_{2}(\mathbf{x})=(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_{0}n\mathbf{I})^{-1}\left( \frac{n}{2}\mathbf{I}\right)\mathbf{k}(z_{x},\mathbf{x})\]

and

\[\tilde{g}_{2}(x_{0}) =\frac{1}{2\lambda_{0}}k(z_{x},x_{0})-\frac{1}{\lambda_{0}n}\mathbf{ k}(x_{0},\mathbf{x})^{T}\tilde{g}_{2}(\mathbf{x})\] \[=\frac{1}{2\lambda_{0}}k(z_{x},x_{0})-\frac{1}{2\lambda_{0}}\mathbf{ k}(x_{0},\mathbf{x})^{T}(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_{0}n\mathbf{I})^{-1}\mathbf{k}(z_{x}, \mathbf{x})\]

Combing previous results, we obtain

\[IF(z;T,\pi_{n})(x_{0})\] \[= -g_{\pi_{n},\lambda_{0}}(x_{0})+\mathbf{k}(x_{0},\mathbf{x})^{T}(\mathbf{k}( \mathbf{x},\mathbf{x})+\lambda_{0}n\mathbf{I})^{-1}g_{\pi_{n},\lambda_{0}}(\mathbf{x})\] \[+2(z_{y}-g_{\pi_{n},\lambda_{0}}(z_{x}))\left(\frac{1}{2\lambda_ {0}}k(z_{x},x_{0})-\frac{1}{2\lambda_{0}}\mathbf{k}(x_{0},\mathbf{x})^{T}(\mathbf{k}(\mathbf{ x},\mathbf{x})+\lambda_{0}n\mathbf{I})^{-1}\mathbf{k}(z_{x},\mathbf{x})\right)\] \[= \mathbf{k}(x_{0},\mathbf{x})^{T}(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_{0}n\mathbf{I })^{-1}\left(g_{\pi_{n},\lambda_{0}}(\mathbf{x})-\frac{1}{\lambda_{0}}(z_{y}-g_{ \pi_{n},\lambda_{0}}(z_{x}))\mathbf{k}(z_{x},\mathbf{x})\right)\] \[-g_{\pi_{n},\lambda_{0}}(x_{0})+\frac{1}{\lambda_{0}}(z_{y}-g_{ \pi_{n},\lambda_{0}}(z_{x}))k(z_{x},x_{0})\] \[= \mathbf{k}(x_{0},\mathbf{x})^{T}(\mathbf{k}(\mathbf{x},\mathbf{x})+\lambda_{0}n\mathbf{I })^{-1}M_{z}(\mathbf{x})-M_{z}(x_{0})\]

as desired. 

Proof of Theorem b.7.: Recall from Section B.2, we use the equivalence notations:

\[IF(z;T,P)=T^{\prime}_{P}(\delta_{z}-P).\]

First, we prove the following Claim:

\[\sup_{z\in X\times\mathcal{Y}}\|T^{\prime}_{\pi}(\delta_{z}-\pi)-T^{\prime}_{ \pi_{n}}(\delta_{z}-\pi_{n})\|_{\mathcal{H}}\to 0,\quad a.s. \tag{23}\]

From Lemma A.6. and Lemma A.7. in [47], we have

\[T^{\prime}_{P}(Q)=-S_{P}^{-1}(W_{P}(Q)).\]

In this equation,

\[W_{P}:\text{lin}(B_{S})\to\mathcal{H},\quad Q\mapsto\mathbb{E}_{Q}[\mathcal{ L}^{\prime}(Y,g_{P,\lambda_{0}}(X))\Phi(X)]\]

is a continuous linear operator on \(\text{lin}(B_{S})\) where the space \(\text{lin}(B_{S})\) is defined in [47]. Moreover, the operator norm of \(W_{P}\) satisfies \(\|W_{P}\|\leq 1\). In addition,

\[S_{P}:\mathcal{H}\to\mathcal{H},\quad f\mapsto 2\lambda_{0}f+\mathbb{E}_{P}[ \mathcal{L}^{\prime\prime}(Y,g_{P,\lambda_{0}}(X))f(X)\Phi(X)]\]

is an invertible continuous linear operator on \(\mathcal{H}\) where \(\Phi\) is the feature map of \(\mathcal{H}\). This implies that \(S_{P}^{-1}:\mathcal{H}\to\mathcal{H}\) is also a continuous linear operator on \(\mathcal{H}\).

Next, we apply Lemma A.7 in [47] by considering the following two sequences: The first sequence is \(\delta_{z}-\pi_{n}\in\text{lin}(B_{S})\) which obviously satisfies

\[\|(\delta_{z}-\pi_{n})-(\delta_{z}-\pi)\|_{\infty}=\|\pi-\pi_{n}\|_{\infty}\to 0,\quad a.s.\]

The second sequence is \(\pi_{n}\in B_{S}\) which obviously satisfies

\[\|\pi-\pi_{n}\|_{\infty}\to 0,\quad a.s.\]
Following the proof of Lemma A.7 in [47], we have that

\[\|T_{\pi}^{\prime}(\delta_{z}-\pi)-T_{\pi_{n}}^{\prime}(\delta_{z}- \pi_{n})\|_{\mathcal{H}}\] \[= \|S_{\pi}^{-1}(W_{\pi}(\delta_{z}-\pi))-S_{\pi_{n}}^{-1}(W_{\pi_{n }}(\delta_{z}-\pi_{n}))\|_{\mathcal{H}}\] \[\leq \|S_{\pi}^{-1}(W_{\pi}(\delta_{z}-\pi))-S_{\pi}^{-1}(W_{\pi_{n}}( \delta_{z}-\pi))\|_{\mathcal{H}}+\|S_{\pi}^{-1}(W_{\pi_{n}}(\delta_{z}-\pi))-S_ {\pi}^{-1}(W_{\pi_{n}}(\delta_{z}-\pi_{n}))\|_{\mathcal{H}}\] \[+\|S_{\pi}^{-1}(W_{\pi_{n}}(\delta_{z}-\pi_{n}))-S_{\pi_{n}}^{-1}( W_{\pi_{n}}(\delta_{z}-\pi_{n}))\|_{\mathcal{H}}\] \[\leq \|S_{\pi}^{-1}\|\|W_{\pi}(\delta_{z}-\pi)-W_{\pi_{n}}(\delta_{z} -\pi)\|_{\mathcal{H}}+\|S_{\pi}^{-1}\|\|W_{\pi_{n}}\|\|(\delta_{z}-\pi)-( \delta_{z}-\pi_{n})\|_{\infty}\] \[+\|S_{\pi}^{-1}-S_{\pi_{n}}^{-1}\|\|W_{\pi_{n}}\|\|\delta_{z}-\pi _{n}\|_{\infty}\] \[\leq \|S_{\pi}^{-1}\|W_{\pi}(\pi)-W_{\pi_{n}}(\pi)\|_{\mathcal{H}}+\|S _{\pi}^{-1}\|\|W_{\pi_{n}}\|\|\pi-\pi_{n}\|_{\infty}\quad\text{(since $W_{\pi_{n}}$ is a linear operator)}\] \[+\|S_{\pi}^{-1}-S_{\pi_{n}}^{-1}\|\|W_{\pi_{n}}\|\left(\|\delta_{z }-\pi\|_{\infty}+\|\pi-\pi_{n}\|_{\infty}\right)\]

Therefore, we only need to study the three terms in the last equation. Note that the first two terms are independent of \(z\), and thus it follows from Step 3 and Step 4 in the proof of Lemma A.7 in [47] that

\[\sup_{z\in\mathcal{X}\times\mathcal{Y}}\|S_{\pi}^{-1}\|\|W_{\pi} (\pi)-W_{\pi_{n}}(\pi)\|_{\mathcal{H}}+\|S_{\pi}^{-1}\|\|W_{\pi_{n}}\|\|\pi- \pi_{n}\|_{\infty}\] \[= \|S_{\pi}^{-1}\|W_{\pi}(\pi)-W_{\pi_{n}}(\pi)\|_{\mathcal{H}}+\| S_{\pi}^{-1}\|\|W_{\pi_{n}}\|\|\pi-\pi_{n}\|_{\infty}\to 0,\quad a.s.\]

To show that the third term also satisfies that

\[\sup_{z\in\mathcal{X}\times\mathcal{Y}}\|S_{\pi}^{-1}-S_{\pi_{n}}^{-1}\|\|W_{ \pi_{n}}\|\left(\|\delta_{z}-\pi\|_{\infty}+\|\pi-\pi_{n}\|_{\infty}\right) \to 0,\quad a.s.,\]

we only need to note the following fact:

1) \(\|S_{\pi}^{-1}-S_{\pi_{n}}^{-1}\|\to 0,\ a.s.\), by Step 2 in the proof of Lemma A.7 in [47]. Moreover, this equation is independent of \(z\).

2) \(\|W_{\pi_{n}}\|\leq 1\) by Step 1 in the proof of Lemma A.7 in [47]. Moreover, this equation is independent of \(z\).

3) \(\|\pi-\pi_{n}\|_{\infty}<+\infty,\ a.s.\), since \(\|\pi-\pi_{n}\|_{\infty}\to 0,\ a.s.\). Moreover, this equation is independent of \(z\).

4) \(\sup_{z\in\mathcal{X}\times\mathcal{Y}}\|\delta_{z}-\pi\|_{\infty}<+\infty\). To see this, we note that by definition

\[\sup_{z\in\mathcal{X}\times\mathcal{Y}}\|\delta_{z}-\pi\|_{\infty}=\sup_{z\in \mathcal{X}\times\mathcal{Y}}\sup_{g\in\mathcal{G}}\left|\int gd(\delta_{z} -\pi)\right|=\sup_{z\in\mathcal{X}\times\mathcal{Y}}\sup_{g\in\mathcal{G}} \left|g(z)-\int gd\pi\right|\leq 2\sup_{g\in\mathcal{G}}\|g\|_{\infty}\]

where the last term is independent of \(z\) and the space \(\mathcal{G}=\mathcal{G}_{1}\cup\mathcal{G}_{2}\cup\{b\}\) is defined in [47]. \(\sup_{g\in\mathcal{G}_{1}}\|g\|_{\infty}\leq 1\) by the definition of \(\mathcal{G}_{1}\). \(\sup_{g\in\{b\}}\|g\|_{\infty}<+\infty\) by our additional assumption on \(b\). Since both \(\mathcal{X}\) and \(\mathcal{Y}\) are bounded and closed, \(\sup_{x\in\mathcal{X}}\sqrt{k(x,x)}\) is bounded above by, say \(\kappa<\infty\). Thus, for every \(h\in\mathcal{H}\) with \(\|h\|_{\mathcal{H}}\leq C_{1}\), we have \(\|h\|_{\infty}\leq C_{1}\kappa\). By definition of \(\mathcal{G}_{2}\), for any \(g\in\mathcal{G}_{2}\), we can write \(g(x,y)=\mathcal{L}^{\prime}(y,f_{0}(x))f_{1}(x)\) with \(\|f_{0}\|_{\mathcal{H}}\leq c_{0}\) and \(\|f_{1}\|_{\mathcal{H}}\leq 1\). Note that \(\|f_{1}\|_{\mathcal{H}}\leq 1\) implies that \(\|f_{1}\|_{\infty}\leq\kappa\) and \(\|f_{0}\|_{\mathcal{H}}\leq c_{0}\) implies that \(\|f_{0}\|_{\infty}\leq c_{0}\kappa\). Thus Assumption B.4 shows that \(\|\mathcal{L}^{\prime}(y,f_{0}(x))\|_{\infty}\leq b^{\prime}_{\rm co\kappa}(y)\) which is bounded on \(\mathcal{Y}\) by our additional assumption (uniformly for \(f_{0}\)). Hence we conclude that \(\sup_{g\in\mathcal{G}_{2}}\|g\|_{\infty}\leq\kappa\sup_{y\in\mathcal{Y}}b^{ \prime}_{\rm co\kappa}(y)<+\infty\). Summarizing the above discussion, we obtain that

\[\sup_{z\in\mathcal{X}\times\mathcal{Y}}\|\delta_{z}-\pi\|_{\infty}\leq 2\sup_{g\in \mathcal{G}}\|g\|_{\infty}<+\infty\]

Combining the above points 1)-4), we conclude that

\[\sup_{z\in\mathcal{X}\times\mathcal{Y}}\|S_{\pi}^{-1}-S_{\pi_{n}}^{-1}\|\|W_{ \pi_{n}}\|\left(\|\delta_{z}-\pi\|_{\infty}+\|\pi-\pi_{n}\|_{\infty}\right)\to 0, \quad a.s.\]

Hence, we obtain our Claim (23).

Applying \(k_{x_{0}}\) to (23), the reproducing property of \(k\) implies that

\[\sup_{z\in\mathcal{X}\times\mathcal{Y}}\left|T_{\pi}^{\prime}(\delta_{z}-\pi)(x _{0})-T_{\pi_{n}}^{\prime}(\delta_{z}-\pi_{n})(x_{0})\right|\to 0,\quad a.s.\]
Note that since \(\mathcal{X}\) and \(\mathcal{Y}\) are bounded and closed by our assumptions, we have that

\[\left|\int_{z\in\mathcal{X}\times\mathcal{Y}}\left(T^{\prime}_{\pi} (\delta_{z}-\pi)(x_{0})\right)^{2}d\pi(z)-\int_{z\in\mathcal{X}\times\mathcal{Y }}\left(T^{\prime}_{\pi_{n}}(\delta_{z}-\pi_{n})(x_{0})\right)^{2}d\pi(z) \right|\] \[\leq |\pi(\mathcal{X}\times\mathcal{Y})|\times\sup_{z\in\mathcal{X} \times\mathcal{Y}}\left|T^{\prime}_{\pi}(\delta_{z}-\pi)(x_{0})-T^{\prime}_{ \pi_{n}}(\delta_{z}-\pi_{n})(x_{0})\right|\] \[\times\left|\int_{z\in\mathcal{X}\times\mathcal{Y}}\left(T^{ \prime}_{\pi}(\delta_{z}-\pi)(x_{0})+T^{\prime}_{\pi_{n}}(\delta_{z}-\pi_{n}) (x_{0})\right)d\pi(z)\right|\] \[\to 0,\quad a.s.\]

where we use the fact that \(\int_{z\in\mathcal{X}\times\mathcal{Y}}\left(T^{\prime}_{\pi}(\delta_{z}-\pi) (x_{0})\right)^{2}d\pi(z)=\xi^{2}(x_{0})<+\infty\). On the other hand, it follows from the strong law of large numbers that

\[\frac{1}{n}\sum_{z_{i}\in\mathcal{D}}\left(T^{\prime}_{\pi_{n}}(\delta_{z_{i}} -\pi_{n})(x_{0})\right)^{2}-\int_{z\in\mathcal{X}\times\mathcal{Y}}\left(T^{ \prime}_{\pi_{n}}(\delta_{z}-\pi_{n})(x_{0})\right)^{2}d\pi(z)\to 0,\quad a.s.\]

Hence, we conclude that

\[\frac{1}{n}\sum_{z_{i}\in\mathcal{D}}\left(T^{\prime}_{\pi_{n}}(\delta_{z_{i}} -\pi_{n})(x_{0})\right)^{2}-\int_{z\in\mathcal{X}\times\mathcal{Y}}\left(T^{ \prime}_{\pi}(\delta_{z}-\pi)(x_{0})\right)^{2}d\pi(z)\to 0,\quad a.s.\]

In other words,

\[\hat{\xi}^{2}(x_{0})=\frac{1}{n}\sum_{z_{i}\in\mathcal{D}}IF^{2}(z_{i};T,\pi_{ n})(x_{0})\to\int_{z\in\mathcal{X}\times\mathcal{Y}}IF^{2}(z;T,\pi)(x)d\pi(z)= \xi^{2}(x_{0}),\quad a.s.\]

For confidence intervals, the proof is straightforward. Theorem B.5 shows that as \(n\to\infty\) we have

\[\frac{g_{\pi_{n},\lambda_{n}}(x_{0})-g_{\pi,\lambda_{0}}(x_{0})}{\xi(x_{0})/ \sqrt{n}}\Rightarrow\mathcal{N}(0,1)\]

By Slutsky's theorem and \(\frac{\xi(x_{0})}{\xi(x_{0})}\to 1,\ a.s.,\) we have

\[\frac{g_{\pi_{n},\lambda_{n}}(x_{0})-g_{\pi,\lambda_{0}}(x_{0})}{\hat{\xi}(x_{ 0})/\sqrt{n}}\Rightarrow\mathcal{N}(0,1)\]

This implies that asymptotically as \(n\to\infty\), we have

\[\mathbb{P}\left(-q_{1-\frac{\alpha}{2}}\leq\frac{g_{\pi_{n},\lambda_{n}}(x_{0 })-g_{\pi,\lambda_{0}}(x_{0})}{\hat{\xi}(x_{0})/\sqrt{n}}\leq q_{1-\frac{ \alpha}{2}}\right)\to 1-\alpha\]

where \(q_{\alpha}\) is the \(\alpha\)-quantile of standard Gaussian distribution \(\mathcal{N}(0,1)\). Hence

\[\left[g_{\pi_{n},\lambda_{n}}(x_{0})-\frac{\hat{\xi}(x_{0})}{\sqrt{n}}q_{1- \frac{\alpha}{2}},g_{\pi_{n},\lambda_{n}}(x_{0})+\frac{\hat{\xi}(x_{0})}{ \sqrt{n}}q_{1-\frac{\alpha}{2}}\right]\]

is an asymptotically exact \((1-\alpha)\)-level confidence interval of \(g_{\pi,\lambda_{0}}(x)\). 

## Appendix E Experiments: Details and More Results

### Experimental Details

We provide more details about our experimental implementation in Section 4.

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 to resolve the conflicts between the theoretical assumptions therein (e.g., continuous-time 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 [61, 32, 75, 110]. Therefore, we will use the following specifications in all experiments.
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\times 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 (1) with regularization hyper-parameter \(\lambda_{n}\equiv 0.1^{10}\).
4. The network is trained using the (full) batch gradient descent (by feeding the whole dataset).
5. The learning rate is properly tuned based on the specific dataset, typically around \(0.5\) and between \(0.1\) and \(0.7\).
6. The training epochs are properly tuned based on the specific dataset, typically around 300-2000. The epochs should not be too small since the training needs to converge to a good solution, but it should also not be too large because we need to stipulate that the training procedure indeed operates in the NTK regime (which is 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.
7. We set \(m^{\prime}=4\) in the PNC-enhanced batching approach and \(R=4\) in the PNC-enhanced cheap bootstrap approach.
8. In DropoutUQ, the dropout rate is a crucial hyper-parameter. We find that the dropout rate has a significant impact on the interval width of DropoutUQ: A large dropout rate produces a wide interval since the dropout rate is linked to the variance of the prior Gaussian distribution [40]. Therefore, to make fair comparisons between different approaches, we adjust the dropout rate in DropoutUQ so that they have a similar or larger interval width as PNC-enhanced batching or PNC-enhanced cheap bootstrap.
9. All experiments are conducted on a single GeForce RTX 2080 Ti GPU.

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

Synthetic Datasets 2: \(X\sim\text{Unif}([0,0.2]^{d})\) and \(Y\sim\sum_{i=1}^{d}\exp(X^{(i)})+(X^{(i)})^{2}+\mathcal{N}(0,0.01^{2})\). The training set \(\mathcal{D}=\{(x_{i},y_{i}):i=1,...,n\}\) is formed by i.i.d. samples of \((X,Y)\) with sample size \(n\). We use \(x_{0}=(0.1,0.1,...,0.1)\) as the fixed test point in the confidence interval task. The rest setup is the same as Section 4. For this dataset, the base network has \(128\times n\) hidden neurons in its hidden layer. More hidden neurons are used in these datasets to guarantee that the network is sufficiently wide and the training procedure operates in the NTK regime. The rest implementation specifications are the same as in Section E.1. The results for constructing confidence intervals are displayed in Table 3, and the results for removing procedural variability are displayed in Table 4.

## Appendix F Additional Methodology and Results

In semi-supervised learning, the unlabeled data appear in the training dataset and typically constitute the majority of the training dataset. Our original PNC algorithm only involves labeled data without

\begin{table}
\begin{tabular}{c|c c c|c c c|c c c} \hline \hline  & \multicolumn{2}{c|}{PNC-enhanced batching} & \multicolumn{2}{c|}{PNC-enhanced cheap bootstrap} & \multicolumn{2}{c}{DropoutUQ} \\  & 95\%CI(CR/TW) & 90\%CI(CR/TW) & MP & 95\%CI(CR/TW) & 90\%CI(CR/TW) & MP & 95\%CI(CR/TW) & 90\%CI(CR/TW) & MP \\ \hline \((d=16)\) & & & & & & & & & \\ \(n=128\) & **0.95**/0.2201 & **0.90**/1.627 & 17.858 & **0.95**/0.1905 & **0.90**/1.462 & 17.830 & **0.95**/0.2227 & **0.95**/0.1919 & 17.834 \\ \(n=256\) & 0.9250.1604 & **0.90**/1.186 & 17.857 & **0.975**/0.1293 & **0.90**/0.0993 & 17.849 & **0.95**/0.1602 & **0.925**/0.1368 & 17.831 \\ \(n=512\) & **0.975**/0.1080 & **0.90**/0.798 & 17.846 & **0.95**/0.1077 & **0.90**/0.0827 & 17.843 & 0.80,10.088 & 0.775,0.0929 & 17.836 \\ \(n=1024\) & **0.95**/0.0855 & 0.8750.0632 & 17.845 & **0.975**/0.0790 & **0.95**/0.606 & 17.836 & 0.8250,0.0950 & 0.800,0.0794 & 17.832 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Confidence interval construction on synthetic datasets 2 with different data sizes \(n=128,256,512,1024\) and dimension \(d=16\). The CR results that attain the desired confidence level \(95\%/90\%\) are in **bold**.
utilizing unlabeled data. To handle semi-supervised learning, we propose a variant of PNC algorithm, which is especially useful for leveraging unlabeled data. We term it an Independent Procedural-Noise-Correcting (IPNC) algorithm whose pseudo-code is given in Algorithm 4. Similarly to our work in Section 3.3, we also develop the large-sample asymptotic results for the IPNC predictor in Theorem F.1. This theorem is then used to guide the development of the IPNC-enhanced batching or IPNC-enhanced cheap bootstrap.

In IPNC, we use two independent datasets, one with labeled data to train the base network and one with possibly unlabeled data to train the artificial-label-trained network. This independent training idea stipulated the independence between the two networks in IPNC.

**Input:** Training dataset \(\mathcal{D}_{1,2}=\mathcal{D}_{1}\cup\mathcal{D}_{2}\) with size \(n_{1}+n_{2}\), which is divided into two disjoint (independent) datasets \(\mathcal{D}_{1}\) and \(\mathcal{D}_{2}\) where \(|\mathcal{D}_{1}|=n_{1}\) and \(|\mathcal{D}_{2}|=n_{2}\). \(\mathcal{D}_{1}\) is labeled while \(\mathcal{D}_{1}\) may not be labeled.

**Procedure:**

**1.** Draw \(\theta^{b}\sim P_{\theta^{b}}=\mathcal{N}(0,\mathbf{I}_{p})\) under NTK parameterization. Train a standard base network with data \(\mathcal{D}_{1}\) and the initialization parameters \(\theta^{b}\), which outputs \(\hat{h}_{n_{1},\theta^{b}}(\cdot)\).

**2.** For each \(x_{i}\in\mathcal{D}_{2}\), generate its procedural-shifted label \(\bar{s}(x_{i})=\mathbb{E}_{P_{\theta^{b}}}[s_{\theta^{b}}(x_{i})]\). Train an auxiliary neural network with data \(\mathcal{D}^{\prime}_{2}=\{(x_{i},\bar{s}(x_{i})):x_{i}\in\mathcal{D}_{2}\}\) and the initialization parameters \(\theta^{b}\) (the same one as in Step 1.), which outputs \(\hat{\phi}^{\prime}_{n_{2},\theta^{b}}(\cdot)\). Subtracting \(\bar{s}\), we obtain \(\hat{\phi}_{n_{2},\theta^{b}}=\hat{\phi}^{\prime}_{n_{2},\theta^{b}}-\bar{s}\) in (4).

**Output:** At point \(x\), output \(\hat{h}_{n_{1},\theta^{b}}(x)-\hat{\phi}_{n_{2},\theta^{b}}(x)\).

**Algorithm 4** Independent Procedural-Noise-Correcting (IPNC) Predictor

To derive the results for IPNC algorithm, we start with the following decomposition:

\[\hat{h}_{n,\theta^{b}}(x)-h^{*}(x)=\hat{h}_{n,\theta^{b}}(x)-\hat{h}^{*}_{n}( x)+\hat{h}^{*}_{n}(x)-h^{*}(x).\]

For the first part \(\hat{h}_{n,\theta^{b}}(x)-\hat{h}^{*}_{n}(x)=\hat{\phi}_{n,\theta^{b}}(x)\), as we discussed in Section 3.1, it is the solution to the following empirical risk minimization problem

\[\hat{\phi}_{n,\theta^{b}}(\cdot)=s_{\theta^{b}}(\cdot)-\bar{s}(\cdot)+\min_{g \in\mathcal{H}}\frac{1}{n}\sum_{i=1}^{n}(\bar{s}(x_{i})-s_{\theta^{b}}(x_{i})- g(x_{i}))^{2}+\lambda_{n}\|g\|_{\tilde{\mathcal{H}}}^{2}\]

Its corresponding population risk minimization problem is as follows:

\[\phi_{\theta^{b}}(\cdot):=s_{\theta^{b}}(\cdot)-\bar{s}(\cdot)+\min_{g\in \mathcal{H}}\mathbb{E}_{\pi}[(\bar{s}(X)-s_{\theta^{b}}(X)-g(X))^{2}]+\lambda_ {0}\|g\|_{\tilde{\mathcal{H}}}^{2}.\]

Note that \(\phi_{\theta^{b}}(\cdot)\) is data-free but depends on the random initialization \(\theta^{b}\). Recall that

\[\hat{h}^{*}_{n}(\cdot)=\bar{s}(\cdot)+\min_{g\in\mathcal{H}}\frac{1}{n}\sum_{ i=1}^{n}[(y_{i}-\bar{s}(x_{i})-g(x_{i}))^{2}]+\lambda_{n}\|g\|_{\tilde{\mathcal{H}}}^ {2}\]

and

\[h^{*}(\cdot)=\bar{s}(\cdot)+\min_{g\in\mathcal{H}}\mathbb{E}_{\pi}[(Y-\bar{s}( X)-g(X))^{2}]+\lambda_{0}\|g\|_{\tilde{\mathcal{H}}}^{2}.\]

For \(\hat{\phi}_{n,\theta^{b}}(x)=\hat{h}_{n,\theta^{b}}(x)-\hat{h}^{*}_{n}(x)\) in (4), it is the solution to the following empirical risk minimization problem

\[\hat{\phi}_{n,\theta^{b}}(\cdot)=s_{\theta^{b}}(\cdot)-\bar{s}(\cdot)+\min_{g \in\mathcal{H}}\frac{1}{n}\sum_{i=1}^{n}(\bar{s}(x_{i})-s_{\theta^{b}}(x_{i}) -g(x_{i}))^{2}+\lambda_{n}\|g\|_{\tilde{\mathcal{H}}}^{2} \tag{24}\]

\begin{table}
\begin{tabular}{c|c c c} \hline MSE & One base network & PNC predictor & Deep ensemble (5 networks) Deep ensemble (2 networks) \\ \hline \(d=16,n=128\) & \((1.70\pm 0.37)\times 10^{-2}\) & \(\mathbf{(1.33\pm 0.19)\times 10^{-2}}\) & \((1.41\pm 0.29)\times 10^{-2}\) & \((1.57\pm 0.32)\times 10^{-2}\) \\ \(d=16,n=256\) & \((1.07\pm 0.16)\times 10^{-2}\) & \(\mathbf{(8.51\pm 1.67)\times 10^{-3}}\) & \((8.54\pm 1.71)\times 10^{-3}\) & \((9.75\pm 1.70)\times 10^{-3}\) \\ \(d=16,n=512\) & \((9.83\pm 1.74)\times 10^{-3}\) & \(\mathbf{(6.17\pm 0.93)\times 10^{-3}}\) & \((6.72\pm 0.91)\times 10^{-3}\) & \((7.51\pm 1.03)\times 10^{-3}\) \\ \(d=16,n=1024\) & \((8.49\pm 1.20)\times 10^{-3}\) & \(\mathbf{(4.90\pm 0.36)\times 10^{-3}}\) & \((5.47\pm 0.53)\times 10^{-3}\) & \((6.42\pm 0.74)\times 10^{-3}\) \\ \hline \end{tabular}
\end{table}
Table 4: Reducing procedural variability on synthetic datasets 2 with different data sizes \(n=128,256,512,1024\) and dimension \(d=16\). The best results are in **bold**.
We consider its corresponding population risk minimization problem as follows:

\[\phi_{\theta^{b}}(\cdot):=s_{\theta^{b}}(\cdot)-\bar{s}(\cdot)+\min_{g\in\mathcal{ H}}\mathbb{E}_{\pi}[(\bar{s}(X)-s_{\theta^{b}}(X)-g(X))^{2}]+\lambda_{0}\|g\|_{\mathcal{H}}^ {2}. \tag{25}\]

Proposition B.2 shows that \(\phi_{\theta^{b}}\) is given by

\[\phi_{\theta^{b}}(\cdot)=s_{\theta^{b}}(\cdot)-\bar{s}(\cdot)+(\mathcal{O}_{K }+\lambda_{0}I)^{-1}\mathcal{L}(\bar{s}-s_{\theta^{b}}).\]

For the second part \(\hat{h}_{n}^{*}(x)-h^{*}(x)\), we have (4) and (6). Therefore we have

\[\hat{h}_{n,\theta^{b}}(\cdot)-h^{*}(\cdot)-\phi_{\theta^{b}}(\cdot)=\text{( \ref{eq:h_
**Theorem F.1** (Large-sample asymptotics of the IPNC predictor).: _Suppose that Assumption C.2 holds. Suppose that Assumptions B.3 and B.4 hold when \(Y\) is replaced by \(Y-s_{\theta^{b}}(X)\) or \(\bar{s}(X)-s_{\theta^{b}}(X)\). Input the training data \(\mathcal{D}\) into Algorithm 4. We have_

\[\sqrt{n_{1}}\left(\hat{h}_{n_{1},\theta^{b}}(x)-h^{*}(x)-\phi_{ \theta^{b}}(x)\right)\Rightarrow\mathcal{N}(0,\beta_{2,\theta^{b}}^{2}(x)), \tag{32}\] \[\sqrt{n_{2}}\left(\hat{\phi}_{n_{2},\theta^{b}}(x)-\phi_{\theta^{ b}}(x)\right)\Rightarrow\mathcal{N}(0,\beta_{3,\theta^{b}}^{2}(x)), \tag{33}\]

_where_

\[\beta_{2,\theta^{b}}^{2}(x)=\int_{z}IF^{2}(z;T_{2},\pi)(x)d\pi(z), \;\beta_{3,\theta^{b}}^{2}(x)=\int_{z}IF^{2}(z;T_{3},\pi)(x)d\pi(z). \tag{34}\]

_Let \(n_{2}=\zeta n_{1}\) where \(\zeta>0\) is a fixed constant. An asymptotically (in the sense of \(n_{1}\rightarrow\infty\)) exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x)\) is_

\[\left[\hat{h}_{n_{1},\theta^{b}}(x)-\hat{\phi}_{n_{2},\theta^{b}}( x)-\left(\frac{\beta_{2,\theta^{b}}(x)}{\sqrt{n_{1}}}+\frac{\beta_{3,\theta^{b}}(x) }{\sqrt{n_{2}}}\right)q_{1-\frac{\alpha}{2}},\right.\] \[\left.\hat{h}_{n_{1},\theta^{b}}(x)-\hat{\phi}_{n_{2},\theta^{b}}( x)+\left(\frac{\beta_{2,\theta^{b}}(x)}{\sqrt{n_{1}}}+\frac{\beta_{3,\theta^{b}}(x) }{\sqrt{n_{2}}}\right)q_{1-\frac{\alpha}{2}}\right]\]

_where \(q_{\alpha}\) is the \(\alpha\)-quantile of the standard Gaussian distribution._

Compared with Theorem 3.5, the major difference is that \(\mathcal{D}_{1}\) and \(\mathcal{D}_{2}\) are independent, which results in the independence between (32) and (33), and thus we can directly sum these two equations. Based on the above theorem, the IPNC-enhanced batching or IPNC-enhanced cheap bootstrap can be similarly derived, as shown below.

**IPNC-enhanced batching.** See Algorithm 5.

**Input:** Training dataset \(\mathcal{D}_{1,2}=\mathcal{D}_{1}\cup\mathcal{D}_{2}\) with size \(n_{1}+n_{2}\), which is divided into two disjoint (independent) datasets \(\mathcal{D}_{1}\) and \(\mathcal{D}_{2}\) where \(|\mathcal{D}_{1}|=n_{1}\) and \(|\mathcal{D}_{2}|=n_{2}\). \(\mathcal{D}_{1}\) is labeled while \(\mathcal{D}_{1}\) may not be labeled. The number of batches \(m^{\prime}\geq 2\).

**Procedure: 1.** Split \(\mathcal{D}_{1}\) into \(m^{\prime}\) batches \(\mathcal{D}_{1,1},...,\mathcal{D}_{1,m^{\prime}}\). Split \(\mathcal{D}_{2}\) into \(m^{\prime}\) batches \(\mathcal{D}_{2,1},...,\mathcal{D}_{2,m^{\prime}}\).

**1.** For each \(j\in[m^{\prime}]\), input \(\mathcal{D}_{1,j}\cup\mathcal{D}_{2,j}\) in Algorithm 4 to output \(\hat{h}_{n_{1}^{\prime},\theta^{b}}^{j}(x)-\hat{\phi}_{n_{2}^{\prime},\theta^ {b}}^{j}(x)\). Note that \(n_{1}^{\prime}=\frac{n_{1}}{m^{\prime}}\) and \(n_{2}^{\prime}=\frac{n_{2}}{m^{\prime}}\).

**3.** Compute \(\psi_{B}(x)=\frac{1}{m^{\prime}}\sum_{j=1}^{m^{\prime}}\left(\hat{h}_{n_{1}^{ \prime},\theta^{b}}^{j}(x)-\hat{\phi}_{n_{2}^{\prime},\theta^{b}}^{j}(x)\right)\),

and \(S_{B}(x)^{2}=\frac{1}{m^{\prime}-1}\sum_{j=1}^{m^{\prime}}\left(\hat{h}_{n_{1}^ {\prime},\theta^{b}}^{j}(x)-\hat{\phi}_{n_{2}^{\prime},\theta^{b}}^{j}(x)- \psi_{B}(x)\right)^{2}\).

**Output:** At point \(x\), output \(\psi_{B}(x)\) and \(S_{B}(x)^{2}\).

**Algorithm 5** IPNC-Enhanced Batching

**Theorem F.2** (Exact coverage of IPNC-enhanced batching confidence interval).: _Suppose that Assumption C.2 holds. Suppose that Assumptions B.3 and B.4 hold when \(Y\) is replaced by \(Y-s_{\theta^{b}}(X)\) or \(\bar{s}(X)-s_{\theta^{b}}(X)\). Let \(n_{2}=\zeta n_{1}\) where \(\zeta>0\) is a fixed constant. Then an asymptotically exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x)\) is \([\psi_{B}(x)-\frac{S_{B}(x)}{\sqrt{m^{\prime}}}q_{1-\frac{\alpha}{2}},\psi_{B} (x)+\frac{S_{B}(x)}{\sqrt{m^{\prime}}}q_{1-\frac{\alpha}{2}}]\) where \(q_{\alpha}\) is the \(\alpha\)-quantile of the \(t\) distribution \(t_{m^{\prime}-1}\) with degree of freedom \(m^{\prime}-1\)._

The proof is similar to Theorem 3.6.

**IPNC-enhanced cheap bootstrap.** See Algorithm 6.

**Theorem F.3** (Exact coverage of IPNC-enhanced cheap bootstrap confidence interval).: _Suppose that Assumption C.2 holds. Suppose that Assumptions B.3 and B.4 hold when \(Y\) is replaced by \(Y-s_{\theta^{b}}(X)\) or \(\bar{s}(X)-s_{\theta^{b}}(X)\). Let \(n_{2}=\zeta n_{1}\) where \(\zeta>0\) is a fixed constant. Then an asymptotically exact \((1-\alpha)\)-level confidence interval of \(h^{*}(x)\) is \([\psi_{C}(x)-S_{C}(x)q_{1-\frac{\alpha}{2}},\psi_{C}(x)+S_{C}(x)q_{1-\frac{ \alpha}{2}}]\) where \(q_{\alpha}\) is the \(\alpha\)-quantile of the \(t\) distribution \(t_{R}\) with degree of freedom \(R\)._

The proof is similar to Theorem 3.7.
**Input:** Training dataset \(\mathcal{D}_{1,2}=\mathcal{D}_{1}\cup\mathcal{D}_{2}\) with size \(n_{1}+n_{2}\), which is divided into two disjoint (independent) datasets \(\mathcal{D}_{1}\) and \(\mathcal{D}_{2}\) where \(|\mathcal{D}_{1}|=n_{1}\) and \(|\mathcal{D}_{2}|=n_{2}\). \(\mathcal{D}_{1}\) is labeled while \(\mathcal{D}_{1}\) may not be labeled. The number of replications \(R\geq 1\).

**Procedure: 1.** Input \(\mathcal{D}_{1}\cup\mathcal{D}_{2}\) in Algorithm 4 to output \(\hat{h}_{n_{1},\theta^{b}}(x)-\hat{\phi}_{n_{2},\theta^{b}}(x)\).

**2.** For each replication \(b\in[R]\), resample \(\mathcal{D}_{1}\), i.e., independently and uniformly sample with replacement from \(\{(x_{1},y_{1}),...,(x_{n_{1}},y_{n_{1}})\}\)\(n_{1}\) times to obtain \(\mathcal{D}_{1}^{*b}=\{(x_{1}^{*b},y_{1}^{*b}),...,(x_{n_{1}}^{*b},y_{n_{1}}^{* b})\}\). Similarly, resample \(\mathcal{D}_{2}\) without labels to obtain \(\mathcal{D}_{2}^{*b}=\{x_{n_{1}+1}^{*b},...,x_{n_{1}+n_{2}}^{*b}\}\). Input \(\mathcal{D}_{1}^{*b}\cup\mathcal{D}_{2}^{*b}\) in Algorithm 4 to output \(\hat{h}_{n_{1},\theta^{b}}^{*b}(x)-\hat{\phi}_{n_{2},\theta^{b}}^{*b}(x)\).

**3.** Compute \(\psi_{C}(x)=\hat{h}_{n_{1},\theta^{b}}(x)-\hat{\phi}_{n_{2},\theta^{b}}(x)\),

and \(S_{C}(x)^{2}=\frac{1}{R}\sum_{b=1}^{R}\left(\hat{h}_{n_{1},\theta^{b}}^{*b}(x )-\hat{\phi}_{n_{2},\theta^{b}}^{*b}(x)-\psi_{C}(x)\right)^{2}\).

**Output:** At point \(x\), output \(\psi_{C}(x)\) and \(S_{C}(x)^{2}\).

**Algorithm 6** IPNC-enhanced Cheap Bootstrap
