Fast, Distribution-free Predictive Inference for Neural Networks with
  Coverage Guarantees
This paper introduces a novel, computationally-efficient algorithm for
predictive inference (PI) that requires no distributional assumptions on the
data and can be computed faster than existing bootstrap-type methods for neural
networks. Specifically, if there are $n$ training samples, bootstrap methods
require training a model on each of the $n$ subsamples of size $n-1$; for large
models like neural networks, this process can be computationally prohibitive.
In contrast, our proposed method trains one neural network on the full dataset
with $(\epsilon, \delta)$-differential privacy (DP) and then approximates each
leave-one-out model efficiently using a linear approximation around the
differentially-private neural network estimate. With exchangeable data, we
prove that our approach has a rigorous coverage guarantee that depends on the
preset privacy parameters and the stability of the neural network, regardless
of the data distribution. Simulations and experiments on real data demonstrate
that our method satisfies the coverage guarantees with substantially reduced
computation compared to bootstrap methods.
Yue Gao, Garvesh Raskutti, Rebecca Willet
# Fast, Distribution-free Predictive Inference for Neural Networks with Coverage Guarantees

###### Abstract

This paper introduces a novel, computationally-efficient algorithm for predictive inference (PI) that requires no distributional assumptions on the data and can be computed faster than existing bootstrap-type methods for neural networks. Specifically, if there are \(n\) training samples, bootstrap methods require training a model on each of the \(n\) subsamples of size \(n-1\); for large models like neural networks, this process can be computationally prohibitive. In contrast, our proposed method trains one neural network on the full dataset with \((\epsilon,\delta)\)-differential privacy (DP) and then approximates each leave-one-out model efficiently using a linear approximation around the differentially-private neural network estimate. With exchangeable data, we prove that our approach has a rigorous coverage guarantee that depends on the preset privacy parameters and the stability of the neural network, regardless of the data distribution. Simulations and experiments on real data demonstrate that our method satisfies the coverage guarantees with substantially reduced computation compared to bootstrap methods.

## 1 Introduction

To assess the accuracy of parameter estimates or predictions without specific distributional knowledge of the data, the idea of re-sampling or sub-sampling on the available data has been long-established to construct prediction intervals, and there is a rich history in the statistics literature on the jackknife and bootstrap methods, see Stine (1985), Efron (1979), Quenouille (1949), Efron and Gong (1983). Among these re-sampling methods, leave-one-out methods (generally referred to as "cross-validation" or "jackknife") are widely used to assess or calibrate predictive accuracy, and can be found in a large line of literature (Stone, 1974, Geisser, 1975).

While it has been demonstrated in a large body of past work with extensive evidence that jackknife-type methods have reliable empirical performance, the theoretical properties of these types of methods are studied relatively little until recently, see Steinberger and Leeb (2018), Bousquet and Elisseeff (2002). One of the most important results among these theoretically guaranteed works is Foygel Barber et al. (2019), which introduces a crucial modification compared to the traditional jackknife method that permits rigorous coverage guarantees of at least \(1-2\alpha\) regardless of the distribution of the data points, for any algorithm that treats the training points symmetrically. We will revisit this work and give more relative details in Section 2.1.

Although theoretically jackknife+ has been proven to have coverage guarantees without distributional assumptions, in practice, this method is computationally costly, since we need to train \(n\) (which is the training sample size) leave-one-out models from scratch to find the predictive interval. Especially for large and complicated models like neural networks, this computational cost is prohibitive. The goal of this paper is to provide a _fast_ algorithm that provides similar theoretical coverage guarantees to those in jackknife+.

To achieve this goal, we develop a new procedure, called Differentially Private Lazy Predictive Inference (DP-Lazy PI), which combines two ideas: lazy training of neural networks and differentially private stochastic gradient descent (DP-SGD). To accelerate the procedure, we introduce a lazy
training scheme inspired by Chizat et al. (2020), Gao et al. (2022) to train the leave-one-out models. The intuition is that with data exchangeability, the leave-one-out models should be quite close or similar to each other and there is no need to train each one from scratch ( _i.e._, random initialization). Instead, we first train a model on the full data and use this model as a good initialization. By using DP-SGD as the initializer for our full model, we are able to provide coverage guarantees since the privacy mechanism prevents information leakage across leave-one-out estimators. In particular, we prove that our DP-Lazy PI procedure has a coverage of at least \(1-2\alpha-3\sqrt{2\eta+2\epsilon+\delta}\) where \(\eta\) represents an out-of-sample stability parameter and \((\epsilon,\delta)\) are the differential privacy parameters. Empirically, we show through simulations and real-data experiments that our method has significant advantage over the existing jackknife+ method in run-time while still maintaining good coverage.

## 2 Preliminaries

We first define key notation. For any values \(v_{i}\) indexed by \(i\in[n]\), define

\[Q_{\alpha}^{+}\left(\{v_{i}\}_{i\in[n]}\right)=\lceil(1-\alpha)(n+1)\rceil^{ \text{th}}\text{ smallest value of }v_{1},\ldots,v_{n}, \tag{1}\]

_i.e._, the \(1-\alpha^{\text{th}}\) quantile of the empirical distribution of these values; Similarly,

\[Q_{\alpha}^{-}\left(\{v_{i}\}_{i\in[n]}\right)=\lfloor(\alpha)(n+1)\rfloor^{ \text{th}}\text{ smallest value of }v_{1},\ldots,v_{n}=-Q_{\alpha}^{+}\left(\{-v_{i}\}_{i\in[n]}\right) \tag{2}\]

is the \(\alpha\) quantile of the empirical distribution.

### Distribution-free Prediction Intervals

Suppose we have training data \((X_{i},Y_{i})\in\mathcal{X}\times\mathcal{Y}\) for \(i=1,\ldots,n\), and a new test point \((X_{n+1},Y_{n+1})\), \(X_{i}\in\mathbb{R}^{p}\), \(Y_{i}\in\mathbb{R}\) drawn independently from the same distribution. We fit a regression model to the training data, _i.e._, a function \(\hat{f}:\mathcal{X}\mapsto\mathcal{Y}\) where \(\hat{f}(x)\) predicts \(Y_{n+1}\) given a new feature vector \(X_{n+1}=x\), and then provide a prediction interval centered around \(\hat{f}(X_{n+1})\) for the test point. Specifically, given some target coverage level \(1-\alpha\), we aim to construct a prediction interval \(\hat{C}_{n,\alpha}\), such that \(\mathbb{P}\left[Y_{n+1}\in\widehat{C}_{n,\alpha}(X_{n+1})\right]\geq 1-\alpha.\) We call the probability \(\mathbb{P}\left[Y_{n+1}\in\widehat{C}_{n,\alpha}(X_{n+1})\right]\) as the coverage, and \(\text{len}[\widehat{C}_{n,\alpha}(X_{n+1})]\) as the interval width of \(\widehat{C}_{n,\alpha}(X_{n+1})\), which is the distance between the left and right endpoints.

Naive Prediction IntervalA naive way to construct a predictive interval at the new test point \((X_{n+1},Y_{n+1})\) is to center the interval at \(\hat{f}(X_{n+1})\) and estimate the margin ( _i.e._, half of the interval width) from the training residuals \(|Y_{i}-\hat{f}(X_{i})|,\ i=1,\ldots,n\). Therefore, a naive prediction interval can be constructed as:

\[\widehat{C}_{n,\alpha}^{\text{naive}}(x)\coloneqq\left[\hat{f}(x)\pm Q_{ \alpha}^{+}\left(\{|Y_{i}-\hat{f}(X_{i})|\}_{i\in[n]}\right)\right].\]

where \(Q_{\alpha}^{+}\left(\{|Y_{i}-\hat{f}(X_{i})|\}_{i\in[n]}\right)\) is defined as (1). Due to the problem of overfitting when the training errors are typically smaller than the test errors, this naive interval \(\widehat{C}_{n,\alpha}(x)\) may undercover-- _i.e._, the probability that \(Y_{n+1}\) falls outside the interval can be larger than \(\alpha\), as discussed in Foygel Barber et al. (2019).

Jackknife Prediction IntervalTo avoid undercoverage due to model overfitting, the _jackknife_ method estimates the margin of errors is estimated by leave-one-out residuals instead of training residuals (Steinberger and Leeb, 2016, 2018). The idea is straightforward: for each \(j=1,\ldots,n\), we fit a regression function \(\hat{f}_{-j}\) using all the training data except the \(j\)-th training sample. Based on these leave-one-out models \(\hat{f}_{-1},\ldots,\hat{f}_{-n}\), we can therefore compute the leave-one-out residuals: \(|Y_{1}-\hat{f}_{-1}(X_{1})|,\ldots,|Y_{n}-\hat{f}_{-n}(X_{n})|\).

With the leave-one-out residuals as well as the regression function \(\hat{f}\) fitted on the full training data, the jackknife prediction interval is constructed as:

\[\widehat{C}_{n,\alpha}^{\text{jackknife}}(x)\coloneqq\left[\hat{f}(x)\pm Q_{ \alpha}^{+}\left(\{|Y_{i}-\hat{f}_{-i}(X_{i})|\}_{i\in[n]}\right)\right].\]
Jackknife+ Prediction IntervalThe jackknife+ is a modification of jackknife, and both of these methods use the leave-one-out residuals when constructing the prediction intervals. The difference is that the jackknife interval is centered at the predicted value \(\hat{f}(X_{n+1})\) (where \(\hat{f}\) is fitted on the full training data), whereas jackknife+ uses the leave-one-out predictions \(\hat{f}_{-i}(X_{n+1}),i=1,\ldots,n\) instead to the build the prediction interval:

\[\widehat{C}_{n,\alpha}^{\text{jackknife+}}(x)\] \[\quad\coloneqq\left[Q_{\alpha}^{-}\left(\left\{\hat{f}_{-i}(x)-|Y _{i}-\hat{f}_{-i}(X_{i})|\right\}_{i\in[n]}\right),\ \ Q_{\alpha}^{+}\left(\left\{\hat{f}_{-i}(x)+|Y_{i}-\hat{f}_{-i}(X_{i})| \right\}_{i\in[n]}\right)\right].\]

Foygel Barber et al. (2019) provides the non-asymptotic coverage guarantee for jackknife+ without assumptions beyond the training and test data being exchangeable. Theoretically, a distribution-free coverage guarantee of at least \(1-2\alpha\) is provided, and it's observed that the method can achieve \(1-\alpha\) coverage in empirical studies.

Summary of coverage guarantees and computational costsThe coverage guarantees and the computational cost of the above methods are summarized in Table 1. Concretely computational cost refers to the computation at two stages, model training and evaluation and we measure cost in terms of number of models trained during the training phase and number of calls to the trained function in the evaluation stage.

Clearly the cost of model training is the largest cost, especially when training a complicated model such as a neural network. Our method attempts to only train a single model (like the naive model) while still providing a coverage guarantee (like jackknife+).

## 3 Method and Algorithm

Our method involves using the lazy estimation framework (Chizat et al., 2020) with an initialization using differential privacy stochastic gradient descent (DP-SGD, Abadi et al. (2016)).

### Lazy Estimation

To calculate the predictive interval fast and reduce the computational cost, we consider using a lazy estimation scheme when estimating the leave-one-out (LOO) models for calculating the predictive interval. The key idea of our approach is given an initialization parameterized by \(\theta_{0}\), for each \(i\), we approximate \(\hat{f}_{-i}\) by linearizing around initial base model and training the linearized model using all but the \(i^{\text{th}}\) training sample. Thus instead of training \(n\) neural networks explicitly, we train a single neural network followed by \(n\) linear models, leading to a substantial computational saving.

For any dataset \(\mathcal{D}=\{(X_{i},Y_{i})\}\) and a large model \(f(x;\ \theta)\), where \(\theta\) denotes the model parameters, we define a lazy estimation operator \(\text{LAZY}_{\lambda}(\theta;\mathcal{D}):\Theta\times\mathcal{D}\mapsto\Theta \subset\mathbb{R}^{M}\) with the ridge regression parameter \(\lambda>0\) and initialization \(\theta_{0}\), we estimate the LOO models by taking a linearization around \(\theta_{0}\) and minimize the penalized training loss on data \(\mathcal{D}_{\setminus j}=\{(X_{i},Y_{i})\}_{i\in[n]\setminus[j]}\),

\[\text{LAZY}_{\lambda}(\theta_{0};\mathcal{D}_{\setminus j})=\theta_{0}+ \operatorname*{argmin}_{\Delta\in\mathbb{R}^{M}}\left\{\sum_{i\in[n]\setminus \{j\}}\mathcal{L}\left(Y_{i},f(X_{i},\theta_{0})+\Delta^{\top}\nabla_{\theta} f(X_{i};\theta)|_{\theta=\theta_{0}}\right)+\lambda\|\Delta\|^{2}\right\}. \tag{3}\]

\begin{table}
\begin{tabular}{l|c c|c c} \hline  & \multicolumn{2}{c|}{Coverage Guarantee} & \multicolumn{2}{c}{Computation} \\ \cline{2-5}  & Distribution-free Theory & Empirical & \# models trained & \# calls to trained function \\ \hline
**Naïve** & No guarantee & \(<1-\alpha\) & 1 & \(n+n_{\text{test}}\) \\ \hline
**jackknife** & No guarantee & \(\approx 1-\alpha\) & \(n\) & \(n+n_{\text{test}}\) \\ \hline
**jackknife+** & \(\geq 1-2\alpha\) & \(\approx 1-\alpha\) & \(n\) & \(n\cdot n_{\text{test}}\) \\ \hline \end{tabular}
\end{table}
Table 1: Coverage guarantees and computation costs of existing methods for estimating prediction interval. Here \(n\) is the training sample size; \(n_{\text{test}}\) is the number of test points for which we seek prediction intervals.
### Differential Privacy (DP) initialization

One of the important questions with lazy estimation is how to choose the initialization \(\theta_{0}\). For this we use the concept of _differential privacy_, which allows us to achieve coverage guarantees with a large reduction in computational costs.

**Definition 1** (Differential Privacy (Dwork and Lei, 2009)).: _A randomized mechanism \(\mathcal{M}:\mathcal{D}\mapsto\mathcal{R}\) with domain \(\mathcal{D}\) and range \(\mathcal{R}\) satisfies \((\epsilon,\delta)\)-differential privacy if for any two adjacent datasets \(d,d^{\prime}\subseteq\mathcal{D}\) and for any subset of outputs \(S\subseteq\mathcal{R}\) it holds that_

\[P\left(\mathcal{M}(d)\in S\right)\leq e^{\epsilon}P\left(\mathcal{M}(d^{ \prime})\in S\right)+\delta. \tag{4}\]

In essence, this differential privacy (DP) condition ensures that the output distribution does not change much when the input data has some small change, which makes it hard to distinguish between the input databases on the basis of the output. Abadi et al. (2016) proposes Algorithm 1, which performs stochastic gradient descent with noisy gradients and showed that it achieves \((\epsilon,\delta)\)-differential privacy, where noise is added to the gradients in the scale of \(O(\frac{\sqrt{\log\delta}}{n\epsilon})\). Algorithms to ensure DP are always designed based on the sensitivity of the original algorithm. More precisely, given an algorithm \(\mathcal{A}\) and a norm function \(\|\cdot\|\) over the range of \(\mathcal{A}\), the sensitivity of \(\mathcal{A}\) is defined as

\[s(\mathcal{A})=\max_{d(D,D^{\prime})=1,\text{card}(D)=n}\|\mathcal{A}(D)- \mathcal{A}(D^{\prime})\| \tag{5}\]

The Laplacian mechanism to construct a differentially private algorithm is as follows (see (Koufogiannis et al., 2015, Holohan et al., 2018)): for an algorithm (or function) \(\mathcal{A}:\mathcal{D}\mapsto\Theta\subset\mathbb{R}^{M}\), the random function \(\mathcal{A}^{\text{DP}}(d)=\mathcal{A}(d)+\xi,\ \forall d\in\mathcal{D}\) satisfies \((\epsilon,0)\)-differential privacy, where the elements \(\xi_{i}\) follows a Laplacian distribution \(\text{Lap}(0,s(\mathcal{A})/\epsilon)\): \(p(\xi)\propto e^{-\epsilon|\xi|/s(\mathcal{A})}\). We can set \(\epsilon\) and \(\delta\) as small constants ( _e.g._, \(\epsilon=0.01\) and \(\delta=10^{-3}\)), and adjust the noise levels in different DP mechanisms.

```
Parameters: Noise scale \(\sigma\); learning rate \(\eta_{t}\); group size \(l\), gradient norm bound \(C\); number of interactions \(T\). Data:\(\mathcal{D}=\{Z_{1},\ldots,Z_{n}\}=\{(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})\}\); Initialize:\(\theta_{0}\) randomly; for\(t\in[T]\)do  Take a random sample \(L_{t}\) with sampling probability \(l/n\);  For each \(i\in L_{t}\), compute the gradient \(g_{t}(Z_{i})\leftarrow\nabla_{\theta_{t}}\mathcal{L}(\theta_{t};Z_{i})\);  Clip gradient: \(\tilde{g}_{t}(Z_{i})=g_{t}(Z_{i})/\max(1,\frac{|g_{t}(Z_{i})|_{2}}{C})\);  Add Noise: \(\tilde{g}_{t}(Z_{i})\leftarrow\frac{1}{|L_{t}|}(\sum_{i}\tilde{g}_{t}(Z_{i})+ \mathcal{N}(0,\sigma^{2}C^{2}I))\);  Descent: \(\theta_{t+1}\leftarrow\theta_{t}-\eta_{t}\tilde{g}_{t}\).  Output \(\theta_{T}\) with the overall privacy cost \((\epsilon,\delta)\) using a privacy accounting method.
```

**Algorithm 1**DP-SGD (Abadi et al., 2016)

### Our method: Lazy Estimation with Differential Privacy

First, we use a noisy SGD method (Abadi et al., 2016) with \((\epsilon,\delta)\) differential privacy to get a full model parameter estimation \(\widehat{\theta}_{n}^{\text{DP}}(\epsilon,\delta)\), which for convenience will be denoted as \(\widehat{\theta}_{n}^{\text{DP}}\). Fig. 3.3 provides an overview of our procedure.

Based on this \(\widehat{\theta}_{n}^{\text{DP}}\), we estimate the LOO model parameters by using lazy training to take a linearization around \(\widehat{\theta}_{n}^{\text{DP}}\) and minimize the penalized training loss on data \(\mathcal{D}_{\setminus j}=\{(X_{i},Y_{i})\}_{i\in[n]\setminus[j]}\),

\[\widetilde{\theta}_{n,\setminus j}=\text{LAZY}_{\lambda}(\widehat{\theta}_{n}^ {\text{DP}};\ \mathcal{D}_{\setminus j}). \tag{6}\]

Here \(\text{LAZY}_{\lambda}(\cdot;\cdot)\) is defined in (3). Based on the LOO model parameters, we can calculate the LOO residuals by

\[R_{j}:=|Y_{j}-f(X_{j};\ \widetilde{\theta}_{n,\setminus j})|.\]

For a test data \(X_{n+1}\), its predictive interval using our method is:

\[\widehat{C}_{\alpha,n}^{\nu}(X_{n+1})\coloneqq\left[Q_{\alpha}^{-}\left\{f(X_ {n+1};\ \widetilde{\theta}_{n,\setminus j})-R_{j}\right\}-\nu,\ Q_{\alpha}^{+}\left\{f(X_ {n+1};\ \widetilde{\theta}_{n,\setminus j})+R_{j}\right\}+\nu\right] \tag{7}\]

where \(\nu\) is a relaxation term that is close to \(0\). Our full method is presented in Algorithm 2.
## 4 Coverage Guarantee

Recall that \(\widehat{\theta}_{n}^{\text{DP}}\) is the NN parameter estimate from \((\epsilon,\delta)\)-DP algorithm using the full training data \(\mathcal{D}=\{Z_{1},\ldots,Z_{n}\}=\{(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})\}\). For any \(j\in[n]\), define \(\widehat{\theta}_{\setminus j}^{\text{DP}}\) as NN parameter estimate from an \((\epsilon,\delta)\)-DP algorithm \(\mathcal{A}^{\text{DP}}\) trained with the \(j\)-th sample deleted, _i.e._, using data \(\mathcal{D}_{\setminus j}=\mathcal{D}\backslash\{Z_{j}\}\). Based on \(\widehat{\theta}_{\setminus j}^{\text{DP}}\), we define another lazy estimate:

\[\widetilde{\theta}_{\setminus j,\setminus j}=\widehat{\theta}_{\setminus j }^{\text{DP}}+\arg\min_{\Delta\theta}\left\{\sum_{i\in[n]\backslash\{j\}} \mathcal{L}\left(Y_{i},f(X_{i},\widehat{\theta}_{\setminus j}^{\text{DP}})+ \Delta\theta^{\top}\nabla_{\theta}f(X_{i};\theta)|_{\theta=\widehat{\theta}_ {\setminus j}^{\text{DP}}}\right)+\lambda\|\Delta\theta\|^{2}\right\}. \tag{8}\]

**Theorem 1**.: _For any \(\nu>0\) and \(\eta>0\), such that_

\[\mathbb{P}\left(|f(X_{n+1};\ \widetilde{\theta}_{n,\setminus j})-f(X_{n+1};\ \widetilde{\theta}_{\setminus j,\setminus j})|>\nu/2\right)\leq\eta, \tag{9}\]

_the coverage of the DP-Lazy prediction interval \(\widehat{C}_{\alpha,n}^{\nu}(x)\) defined in Equation (7) is larger than \(1-2\alpha-3\sqrt{2\eta+2\epsilon+\delta}\), where \((\epsilon,\delta)\) are the DP parameters, i.e.,_

\[\mathbb{P}\left[Y_{n+1}\in\widehat{C}_{\alpha,n}^{\nu}(X_{n+1})\right]\geq 1 -2\alpha-3\sqrt{2\eta+2\epsilon+\delta}. \tag{10}\]

**Remark 1**.: _In Theorem 1, essentially (9) requires an out-of-sample stability in the DP algorithm, which appears elsewhere in the literature Foygel Barber et al. (2019) and is referred to as "hypothesis stability" in Bousquet and Elisseeff (2002). We want to emphasize that this out-of-sample stability only requires that, for a test data point that is independent of the training data, the predicted value does not change much if we remove one point in the training data set, which can hold even for algorithms that suffer from strong overfitting in contrast to the in-sample stability. As an illustrative example, the K-nearest-neighbor algorithm is shown to satisfy the out-of-sample stability with \(\nu=0\) and \(\eta=K/n\) (see Example 5.5 in Foygel Barber et al. (2019))._

_For neural networks, notions of stability are also assumed in the literature (Verma and Zhang, 2019; Forti et al., 1994). In the Supplementary Material, we also give an example showing that if \(x\) is multivariate Gaussian, for a two-layer neural network (one hidden layer) with activation functions like ReLU, sigmoid or tanh and \(s(\mathcal{A})=O(1/n)\) (defined in Equation (5)), (9) holds true with \(\nu=O(1/n)\) up to some log terms and \(\eta=O(e^{-\epsilon n})\) for some constant \(c\)._

Figure 1: The outline of DP-Lazy PI method to calculate prediction intervals
### Proof Overview of Theorem 1

In this section, we'll provide the key ideas and the skeleton of the proof for the coverage guarantee, while the completed proof is postponed in the Appendix. Recall \(R_{j}:=|Y_{j}-f(X_{j};\ \widetilde{\theta}_{n,\setminus j})|\) and define \(R_{j,\setminus j}:=|Y_{j}-f(X_{j};\ \widetilde{\theta}_{\langle j,\setminus j \rangle})|\).

First, for any \(\alpha<\alpha^{\prime}<1\) we show that a jackknife+ type interval \(\widetilde{C}_{\alpha^{\prime},n}(X_{n+1})\) defined as

\[\widetilde{C}_{\alpha^{\prime},n}(x)\coloneqq\left[Q_{\alpha^{\prime}}^{-} \left\{f(x;\ \widetilde{\theta}_{\langle j,\setminus j\rangle}-R_{j,\setminus j}\right\},\ Q_{\alpha^{\prime}}^{+}\left\{f(x;\ \widetilde{\theta}_{\langle j,\setminus j \rangle}+R_{j,\setminus j}\right\}\right] \tag{11}\]

is contained within our rapidly calculated prediction interval \(\widehat{C}_{\alpha,n}^{\nu}(X_{n+1})\) with high probability that depends on the relaxation term \(\eta\) and the DP parameters \((\epsilon,\delta)\):

\[\mathbb{P}\left[\widetilde{C}_{\alpha^{\prime},n}(X_{n+1}) \nsubseteq\widehat{C}_{\alpha,n}^{\nu}(X_{n+1})\right] \tag{12}\] \[\quad=\mathbb{P}\Bigg{[}\left\{Q_{\alpha}^{+}\left\{f(X_{n+1};\ \widetilde{\theta}_{n,\setminus j})+R_{j}+\nu\right\}<Q_{\alpha^{\prime}}^{+} \left\{f(X_{n+1};\ \widetilde{\theta}_{\langle j,\setminus j\rangle})+R_{j,\setminus j}\right\}\right\}\] (13) \[\qquad\qquad\cup\left\{Q_{\alpha}^{-}\left\{f(X_{n+1};\ \widetilde{\theta}_{n,\setminus j})-R_{j}-\nu\right\}>Q_{\alpha^{\prime}}^{-} \left\{f(X_{n+1};\ \widetilde{\theta}_{\langle j,\setminus j\rangle})-R_{j,\setminus j}\right\} \right\}\Bigg{]}\] \[\leq\mathbb{P}\Bigg{[}\Bigg{\{}\sum_{j=1}^{n}\mathbb{1}\Big{(} \Big{|}f(X_{n+1};\ \widetilde{\theta}_{\langle j,\setminus j \rangle})-f(X_{n+1};\ \widetilde{\theta}_{n,\setminus j})\Big{|}\] (14) \[\qquad\qquad+\Big{|}f(X_{j};\ \widetilde{\theta}_{\langle j, \setminus j\rangle})-f(X_{j};\ \widetilde{\theta}_{n,\setminus j})\Big{|}>\nu\Big{)}\geq( \alpha^{\prime}-\alpha)(n+1)\Bigg{\}}\Bigg{]}\] \[\leq\frac{1}{(\alpha^{\prime}-\alpha)}\Bigg{\{}\mathbb{P}\left[ \Big{|}f(X_{n+1};\ \widetilde{\theta}_{\langle j,\setminus j\rangle})-f(X_{n+1};\ \widetilde{\theta}_{n,\setminus j})\Big{|}>\nu/2\right]\] \[\qquad\qquad+\mathbb{P}\left[\Big{|}f(X_{j};\ \widetilde{\theta}_{\langle j, \setminus j\rangle})-f(X_{j};\ \widetilde{\theta}_{n,\setminus j})\Big{|}>\nu/2\right] \Bigg{\}}\] (15) \[\leq\frac{1}{(\alpha^{\prime}-\alpha)}\left\{2\mathbb{P}\left[ \Big{|}f(X_{n+1};\ \widetilde{\theta}_{\langle j,\setminus j \rangle})-f(X_{n+1};\ \widetilde{\theta}_{n,\setminus j})\Big{|}>\nu/2\right]+2 \epsilon+\delta\right\}\] (16) \[\leq\frac{2\eta+2\epsilon+\delta}{\alpha^{\prime}-\alpha}. \tag{17}\]
The probabilities are taken with respect to all the training data \(\mathcal{D}\) as well as the test data \((X_{n+1},Y_{n+1})\). (13) holds by the definitions of the prediction intervals of interest, which is equivalent to (14); by data exchangeability and Markov inequality, (15) holds true; by the property of differential privacy in the Appendix, the in-sample stability term in (15) is relaxed to an out-of-sample stability condition in (16), which is bounded by the condition in Equation (9).

By the jackknife+ coverage guarantee in Foygel Barber et al. (2019) that

\[\mathbb{P}(Y_{n+1}\notin\widetilde{C}_{\alpha^{\prime},n}(X_{n+1}))\leq 2 \alpha^{\prime}, \tag{18}\]

we can bound the miscoverage rate for the prediction interval \(\widehat{C}_{\alpha,n}^{\nu}(X_{n+1})\):

\[\mathbb{P}(Y_{n+1}\notin \widehat{C}_{\alpha,n}^{\nu}(X_{n+1})) \tag{19}\] \[\leq \mathbb{P}(Y_{n+1}\notin\widetilde{C}_{\alpha^{\prime},n}(X_{n+1} ))+\mathbb{P}\left[\widetilde{C}_{\alpha^{\prime},n}(X_{n+1})\nsubseteq\widehat {C}_{\alpha,n}^{\nu}(X_{n+1})\right]\] (20) \[\leq 2\alpha^{\prime}+\frac{2\eta+2\epsilon+\delta}{\alpha^{\prime}-\alpha} \tag{21}\]

for all \(\alpha^{\prime}>\alpha\). Take \(\alpha^{\prime}=\alpha+\sqrt{2\eta+2\epsilon+\delta}\), we therefore have the coverage of \(\widehat{C}_{\alpha,n}^{\nu}(X_{n+1})\):

\[\mathbb{P}(Y_{n+1}\in \widehat{C}_{\alpha,n}^{\nu}(X_{n+1}))=1-\mathbb{P}(Y_{n+1}\notin \widehat{C}_{\alpha,n}^{\nu}(X_{n+1}))\] \[\geq 1-2\alpha-3\sqrt{2\eta+2\epsilon+\delta}.\]

## 5 Experiments

We compare the following methods for estimating prediction intervals: (1) jackknife+: defined in Section 2.1, with the base algorithm being neural networks with random initialization; (2) DP-Lazy PI (labeled _lazy_dp_ in plots): our proposed method with the same NN architecture used in jackknife+, with the relaxation term \(\nu\) set to zero, and differential parameters set as \(\epsilon=0.01,\ \delta=10^{-3}\); (3) Lazy PI without DP (labeled _lazy_finetune_ in plots): removes the privacy mechanism compared to DP-Lazy PI, _i.e._, use \(\widehat{\theta}_{n}\) instead of \(\widehat{\theta}_{n}^{\mathrm{DP}}\) as \(\theta_{0}\) in the lazy estimation step (3).

For evaluation, we consider the following three performance aspects on the test data set \(\mathcal{D}_{\mathrm{test}}\) with size \(n_{\mathrm{test}}\) for a prediction interval \(\widehat{C}_{n}(\cdot)\): (1) Coverage: \(\frac{1}{n_{\mathrm{test}}}\sum_{(x,y)\in\mathcal{D}_{\mathrm{test}}}\mathbbm{1 }\left[y\in\widehat{C}_{n}(x)\right](2)\) Compute time; (3) Average interval width: \(\frac{1}{n_{\mathrm{test}}}\sum_{(x,y)\in\mathcal{D}_{\mathrm{test}}}\mathrm{ len}[\widehat{C}_{n}(x)]\). As we set the level as \(\alpha=0.1\), our target of the coverage is \(1-2\alpha=80\%\), we want to emphasize that a higher coverage of the prediction intervals is not always desirable, since it may suggest that the prediction intervals are overly wide and conservative.

### Simulation

To generate the data, \(\{X_{i},\ i\in[N]\}\) are randomly selected from a Gaussian distribution \(X_{i}\in\mathbb{R}^{p}\sim\mathcal{N}(0,5\cdot\mathbb{I}_{p}),\ i\in N\), where \(N=5000\). The responses \(Y_{i},\ i\in[N]\) are generated by:

\[Y_{i}=\sqrt{\mathrm{ReLU}(X_{i}^{\top}\beta)}+\epsilon_{i},\ \text{where}\ \epsilon_{i}\sim\mathcal{N}(0,0.5),\ \beta\sim\mathrm{Beta}(1.0,\,2.5). \tag{22}\]

We consider a neural network with two hidden layers, each contains 64 nodes. The data are randomly split into the training data with the training size \(n=100\), and evaluation set with the size \(n_{\mathrm{test}}=4900\).

We construct prediction intervals using jackknife+, lazy finetune and DP-Lazy PI respectively on the training data with \(\alpha=0.1\). When training the neural networks, the training batch size is 10 and the max number of epochs is 10. The penalty parameter is taken as \(\lambda=10\) in the lazy type of methods.

As shown in Figure 2, the coverage of our method is closer to the target 80% coverage and in terms of the compute time, lazy methods reduce the average compute time significantly compared to jackknife+. As the feature dimension increases, the coverage is still guaranteed; we can also see in Figure 2 (b) that compared to lazy finetune without noise added in the full model training, lazy DP becomes better in the interval width, which might come from the interpolation (near zero training error) in the high-dimensional case that makes the LOO estimates in lazy finetune stays at the initial base model in the lazy procedure.
### Real Data

Data sets(1) _The BlogFeedback1_ data set (Spiliopoulou et al., 2014, Foygel Barber et al., 2019), contains information on 52397 blog posts with \(p=280\) covariates. The response is the number of comments left on the blog post in the following 24 hours, which we transformed as \(Y=\log(1+\#\text{comments})\). (2) _The Medical Expenditure Panel Survey 2016 data set2_ contains 33005 records on individuals' utilization of medical services such as visits to the doctor, hospital stays etc with feature dimension \(p=107\) with relevant features such as age, race/ethnicity, family income, occupation type, etc. Our goal is the predict the health care system utilization of each individual, which is a composite score reflecting the number of visits to a doctor's office, hospital visits, days in nursing home care, etc.

Footnote 1: [https://archive.ics.uci.edu/ml/datasets/BlogFeedback](https://archive.ics.uci.edu/ml/datasets/BlogFeedback)

Footnote 2: [https://meps.ahrq.gov/mepsweb/data_stats/download_data_files_detail.jsp?cboPufNumber=HC-192](https://meps.ahrq.gov/mepsweb/data_stats/download_data_files_detail.jsp?cboPufNumber=HC-192)

Experiment SettingsThe training size to construct the prediction interval is set \(n=100\). We consider the 3-layer neural network with hidden layers \((64,64)\) as the base NN architecture. To train the model parameters in jackknife+, we initialize the NN randomly for each LOO models. The penalty level is set as \(\lambda=10\), and the DP parameters are \(\epsilon=0.01,\ \delta=10^{-3}\). The batch size is 10 and the maximum number of epochs is 10. We repeat the trials 15 times with different random seeds for the train-test split and the random initialization of NN parameters in jackknife+. The figures show the average coverage, computing time and the interval width across these random trails.

ResultsAs shown in Figure 3 (a) and (b), when we look into the real data sets, DP-Lazy PI decreases the computing time while enjoys a narrower interval width, while achieving the coverage guarantee at 80%. Compared to the data generated in the simulation, for real data with much complicated data structure, the model mis-specification makes the width of the jackknife+ estimator greater and the implicit regularization of the lazy approach reduces this width. The implementation codes can be found in [https://github.com/Vioyue6/DP_LAZY_PI](https://github.com/Vioyue6/DP_LAZY_PI).

Figure 2: The average coverage, average computing time and width of the prediction intervals with \(\alpha=0.1\) in 15 trials, with the error bars to show the \(\pm\) standard errors. The dotted line in the coverage panel indicates the target coverage level 80%.
## 6 Discussion and Limitations

In this paper, we describe a new method, Differential Privacy Lazy Predictive Inference or DP-Lazy PI, which provides a fast method for distribution-free inference for neural networks. Our method involves using differentially private stochastic gradient descent (DP-SGD) as an initial estimate and then computing leave-one-out approximations using lazy training and then re-combining the leave-one-out estimators to form the predictive interval. Importantly, we are able to provide coverage guarantees that closely match the coverage guarantees provided for the jackknife+ procedure with a fraction of the computational cost since only a single neural network model needs to be trained.

The two real data examples also suggest that when we have significant model misspecification, the implicit regularization that lazy training provides means that our DP Lazy-PI method has narrower width prediction intervals than jackknife+.

An alternative method we evaluated in our experiments was lazy finetune, which removes the privacy mechanism prior to lazy training. The simulation results show very similar and at times slightly improved performance compared to our DP Lazy-PI approach; however, a limitation is that we lack theoretical guarantees for this approach. It remains an open question whether we can provide theoretical guarantees for lazy finetune, noting that in the zero training error regime where we can perfectly interpolate, lazy training does not apply since we are at a zero gradient initialization, making linearized models useless. However, in higher dimensions, the implicit regularization through early stopping of neural network training may be the reason lazy training works well. A further limitation of this work is that our method requires the base algorithm ( _e.g._, the neural network estimation) to be stable and satisfy the condition (9).

Figure 3: The coverage, average computing time and width of the prediction intervals on real data sets with \(\alpha=0.1\). The dotted line in the left panels in (a) and (b) is the target 80% coverage.
## Acknowledgements

R. Willett gratefully acknowledges the support of AFOSR grant FA9550-18-1-0166 and NSF grants DMS-2023109 and DMS-1925101. G. Raskutti acknowledges the support of NIH grant R01 GM131381-03.
