Title: THE DISCRIMINATIVE JACKKNIFE: QUANTIFYING PREDICTIVE UNCERTAINTY VIA HIGHER-ORDER INFLUENCE FUNCTIONS Anonymous

Abstract: Deep learning models achieve high predictive accuracy in a broad spectrum of tasks, but rigorously quantifying their predictive uncertainty remains challenging. Usable estimates of predictive uncertainty should (1) cover the true prediction target with a high probability, and (2) discriminate between high-and low-confidence prediction instances. State-of-the-art methods for uncertainty quantification are based predominantly on Bayesian neural networks. However, Bayesian methods may fall short of ( 1) and (2) -i.e., Bayesian credible intervals do not guarantee frequentist coverage, and approximate posterior inference may undermine discriminative accuracy. To this end, this paper tackles the following question: can we devise an alternative frequentist approach for uncertainty quantification that satisfies (1) and ( 2)? To address this question, we develop the discriminative jackknife (DJ), a formal inference procedure that constructs predictive confidence intervals for a wide range of regression models, is easy to implement, and provides rigorous theoretical guarantees on (1) and ( 2). The DJ procedure uses higher-order influence functions (HOIFs) of the trained model parameters to construct a jackknife (leave-one-out) estimator of predictive confidence intervals. DJ computes HOIFs using a recursive formula that requires only oracle access to loss gradients and Hessian-vector products, hence it can be applied in a post-hoc fashion without compromising model accuracy or interfering with model training. Experiments demonstrate that DJ performs competitively compared to existing Bayesian and non-Bayesian baselines. Model Input Model Input Model Input Model Output Model Output Model Output Accurate coverage

1. INTRODUCTION:

Deep learning models have achieved state-of-the-art performance on a variety of supervised learning tasks, and are becoming increasingly popular in various application domains (LeCun et al. (2015) ).

A key question often asked of such models is "Can we trust this particular model prediction?" This question is highly relevant in high-stakes applications where the model is used to guide and inform critical decision-making -examples of such applications include: medical decision support, autonomous vehicle control, and financial forecasting (Amodei et al. (2016) ). Despite their impressive predictive accuracy, rigorously quantifying the predictive uncertainty of deep learning models is a challenging and yet an unresolved problem (Gal (2016) ; Ovadia et al. (2019) ).

Actionable estimates of predictive uncertainty are ones that (1) cover the true prediction target with high probability, and (2) discriminate between high-and low-confidence predictions. (Figure 1 depicts a pictorial visualization for the coverage and discrimination requirements.) The first requirement, frequentist coverage, is especially relevant in applications where predictive uncertainty is incorporated in a decision-theoretic framework (e.g., administering medical treatments (Dusenberry et al. (2019) ), or estimating value functions in model-free reinforcement learning (White & White (2010) )). The second requirement, discrimination, is crucial for auditing model reliability (Schulam & Saria (2019) ), detecting dataset shifts and out-of-distribution samples (Barber et al. (2019a) ), and actively collecting new training examples for which the model is not confident (Cohn et al. (1996) ).

Existing methods for uncertainty quantification are based predominantly on Bayesian neural networks (BNNs), whereby predictive uncertainty is evaluated via posterior credible intervals (Welling & Teh (2011) ; Hernández-Lobato & Adams (2015) ; Ritter et al. (2018) ; Maddox et al. (2019) ). However, Figure 1 : Pictorial depiction for coverage and discrimination in uncertainty estimates. Red dots correspond to training data and dotted black line corresponds to the target function. Confidence intervals are visualized as shaded blue regions, and dotted blue lines are the upper and lower confidence bounds. The left panel shows a confidence interval that perfectly covers the data points, but does not discriminate high-confidence predictions (regions with dense training data) and low-confidence ones (regions with scarce training data). The middle panel shows a confidence interval with a width proportional to the density of training data, but does not cover any data point. The right panel shows a confidence interval that satisfies both coverage and discrimination requirements.

BNNs require significant modifications to the training procedure, and exact Bayesian inference is computationally prohibitive in practice. Approximate dropout-based inference schemes (e.g., Monte Carlo dropout (Gal & Ghahramani (2016) ) and variational dropout (Kingma et al. (2015) )) have been recently proposed as computationally efficient alternatives. However, Bayesian inference in dropout-based models has been shown to be ill-posed, since the induced posterior distributions in such models do not concentrate asymptotically (Osband (2016) ; Hron et al. (2017) ), which jeopardizes both the coverage and discrimination performance of the resulting credible intervals. Moreover, even with exact inference, Bayesian credible intervals generally do not guarantee frequentist coverage (Bayarri & Berger (2004) ). Alternative non-Bayesian methods have been recently developed -mostly based on ad-hoc ensemble approaches (Lakshminarayanan et al. (2017) ) -but formal frequentist methods with rigorous theoretical guarantees are still lacking.

Summary of Contributions.:

In this paper, we develop a formal procedure for constructing frequentist (pointwise) confidence intervals on the predictions of a broad class of deep learning models. Our method -which we call the discriminative jackknife (DJ) -is inspired by the classic jackknife leave-one-out (LOO) re-sampling procedure for estimating variability in statistical models (Miller (1974) ; Efron (1992) ). In order to ensure both frequentist coverage and discrimination, DJ constructs feature-dependent confidence intervals using the LOO local prediction variance at the input feature, and adjusts the interval width (for a given coverage probability) using the model's average LOO error residuals. Whereas the classic jackknife satisfies neither the coverage nor the discrimination requirements (Barber et al. (2019b) ), DJ provides strong theoretical guarantees on both (i.e., DJ generates confidence intervals resembling those in the rightmost panel of Figure 1 with high probability).

Central to the DJ procedure is the usage of influence functions -a key concept in robust statistics and variational calculus (Efron (1992) ; Cook & Weisberg (1982) ; Hampel et al. (2011) ) -in order to estimate the parameters of models trained on LOO versions of the training data, without the need to exhaustively re-train the model for each held-out data point. That is, using the von Mises expansion (Fernholz ( 2012)) -a variant of Taylor series expansion for statistical functionals -we represent the (counter-factual) model parameters that would have been learned on LOO versions of the training dataset in terms of an infinite series of higher-order influence functions (HOIFs) for the model parameters trained on the complete dataset. We show that a truncated von Mises expansion (with a finite number of HOIFs) is sufficient to preserve the guarantees on coverage and discrimination. To enable efficient computation of the von Mises expansion, we derive a recursive formula that links HOIFs with lower-order influence functions and the partial derivatives of the model loss function.

Comprehensive experimental evaluation demonstrates that the DJ procedure performs competitively compared to both Bayesian and non-Bayesian methods. In particular, we found that the DJ procedure outperforms Bayesian neural networks with respect to discriminative performance, and provides better frequentist coverage compared to existing non-Bayesian methods, while being significantly less computationally expensive than both classes of methods.

Learning Setup. We consider a standard supervised learning setup with (x, y) being a feature-label pair, where the feature x belongs to a d-dimensional feature space X ⊆ R d , and y ∈ R. A model is trained to predict y using a dataset D {(x i , y i )} n i=1 of n training examples which are drawn i.i.d from an unknown distribution P. Let the learning model be denoted as f (x; θ) : X → Y, where θ ∈ Θ are the model parameters, and Θ is the parameter space. The learned model parameters θ ∈ Θ are obtained by solving the optimization problem θ = arg min θ∈Θ L(D, θ), with a loss function

L(D, θ) 1 n n i=1 (y i , f (x i ; θ)),

(1) where we fold in any regularization terms into (.). We do not pose any assumptions on how the loss function in (1) is optimized. The model f (x; θ) can be any neural network variant (e.g., feed-forward, convolutional, etc). Throughout the paper, we limit our discussion to the regression task; extending our framework to classification tasks is an interesting direction for future work.

Uncertainty Quantification. The predictions issued by the trained model (with learned parameters θ) are given by f (x; θ). Our main goal is to obtain estimates for the uncertainty in the model predictions, expressed through a pointwise confidence interval C(x; θ) defined as follows:

C(x; θ) [ f L (x; θ), f U (x; θ) ], f L (x; θ) < f (x; θ) < f U (x; θ), ∀x ∈ X .

(2) The extent of uncertainty in the model's prediction (for a data point with feature x) is quantified by the interval width W (.) of C(x; θ), defined as W (C(x; θ)) f U (x; θ) -f L (x; θ). Wider intervals imply less confidence, and vice versa. For C(x; θ) to be usable, it has to satisfy the following requirements:

(i) Frequentist coverage: The confidence interval C(x; θ) covers the true prediction target with an arbitrarily prespecified coverage probability of (1 -α), α ∈ (0, 1), i.e.,

P y ∈ C(x; θ) ≥ 1 -α,

where the probability is taken with respect to a (new) test point (x, y) as well as with respect to the training data D (Lawless & Fredette (2005) ; Barber et al. (2019b) ).

(ii) Discrimination: The confidence interval C(x; θ) is wider (on average) for test points with less accurate predictions (Leonard et al. (1992) ), i.e., for the test points x, x ∈ X , we have

E[W (C(x; θ))] ≥ E[W (C(x ; θ))] ⇔ E[ (y, f (x; θ))] ≥ E[ (y , f (x ; θ))],

where the expectation E[.] is taken with respect to the training data D.

In the next Section, we develop a formal post-hoc procedure that takes the dataset D and the trained model f (x; θ) as inputs, and evaluates an uncertainty estimate C(x; θ) that satisfies both requirements.

3. THE DISCRIMINATIVE JACKKNIFE:

Before presenting our proposed discriminative jackknife (DJ) procedure, we start with a brief recap of the classical jackknife. The jackknife quantifies predictive uncertainty in terms of the (average) prediction error, which is estimated with a leave-one-out (LOO) construction found by systematically leaving out each data point from the training dataset D, and evaluating the error of the re-trained model on the held-out sample. That is, for a target coverage of (1 -α), a naïve jackknife estimates the uncertainty in the prediction for a test point x as (Miller (1974) ; Efron (1992) ):

EQUATION

where

r i = | y i -f (x i ; θ-i ) |

is the error residual on the i th data point, θ-i are the parameters of the model re-trained on the dataset D \ {(x i , y i )} with the i th point removed, and Q α,n is the (1 -α) quantile for the empirical distribution of elements of the set R = {r 1 , . . . , r n } defined as1 

Q α,n (R)

the (1 -α) (n + 1) -th smallest value of the elements of R.

Albeit intuitive, the naïve jackknife estimate in (3) is not guaranteed to achieve the target coverage of (1 -α) (Barber et al. (2019b) ). More crucially, the interval width W ( C α,n (x; θ)) is a constant that does not depend on the feature x, which renders discrimination impossible, i.e., when applied, the naïve jackknife would result in confidence intervals resembling the leftmost panel in Figure 1 .

3.1. EXACT CONSTRUCTION OF THE DJ CONFIDENCE INTERVALS:

We propose an ameliorated jackknife procedure, the DJ, which addresses the shortcomings of the classical jackknife. We start by defining some notation. Let the sets V + n (x) and V - n (x) be defined as:

EQUATION

where

v i (x) = f (x; θ) -f (x; θ-i ), n + = |V + n |, and n -= |V - n |.

Our DJ procedure estimates the predictive confidence interval for test point x as follows:

C α,n (x; θ) = [ f L (x; θ), f U (x; θ) ],

where

f L (x; θ) = f (x; θ) -Q α,n (R n ) -Q 1-α,n -(V - n (x)), f U (x; θ) = f (x; θ) + Q α,n (R n ) Marginal Error + Q α,n + (V + n (x))

Local Variability

.

(5)

Unlike the traditional jackknife, the construction of the DJ confidence intervals comprises two terms: marginal prediction error and local prediction variability. The prediction error term, Q α,n (R n ), uses the LOO residuals to estimate the model's (marginal) generalization error. The prediction variability terms,

Q α,n + (V + n (x)) and Q 1-α,n -(V - n (x)

), quantify the extent to which each training data point impacts the value of the prediction at test point x. The prediction error is constant, i.e., does not depend on x, hence it only contributes to coverage but does not contribute to discrimination. On the contrary, the local variability terms depend on the test point feature x, hence they fully determine the discrimination performance. As we show later in Section 3.4, constructing confidence intervals this way provides strong theoretical guarantees on the achieved frequentist coverage and discrimination.

Figure 2 illustrates the construction of the DJ confidence intervals in (5). The confidence intervals are chosen so that the boundaries of the average error and local variability are exceeded by (n + 1)(1 -α) out of the n LOO samples -these are marked with a star. For the average prediction error term, the width of the resulting boundary is the same for any test data point x ∈ X . For the local prediction variability term, the width of the boundary depends on the specific value of the feature x, and should get wider for less confident prediction instances. Based on the confidence intervals in (5), it follows that the DJ interval width is given by

EQUATION

The marginal error and local variability terms in (6) jointly capture two types of uncertainty: epistemic and aleatoric uncertainties (Gal (2016) ). Epistemic uncertainty measures how well the model is matched to the data, and is reducible as the size of training data increases. On the contrary, aleatoric uncertainty is the irreducible uncertainty arising from the inherent sources of complexity in the data, such as label noise or hidden features (Malinin & Gales (2018) ). Consistency of the DJ confidence estimates requires that W ( C α,n (x; θ)) → 0, i.e., the interval width vanishes, as the size of the training data increases (n → ∞). It follows from (6) that if there is no aleatoric uncertainty, and if the underlying model is stable (i.e., lim n→∞ v i = 0) and consistent (i.e., lim n→∞ r i = 0), then the interval width W ( C α,n (x; θ)) vanishes as n grows asymptotically (more training data is collected).

3.2. EFFICIENT IMPLEMENTATION VIA HIGHER-ORDER INFLUENCE FUNCTIONS:

Exact computation of the DJ confidence intervals in (5) requires re-training the model n times in order to collect the "perturbed" LOO parameters { θ-i } n i=1 . To scale up the DJ procedure, we propose to use influence functions -a classic tool from robust statistics (Huber & Ronchetti (1981); Hampel et al. (2011) ) -in order to recover the parameters { θ-i } n i=1 on the basis of the trained model f (x; θ), without the need for repeated re-training. Through our influence function-based implementation, the DJ can be applied in a post-hoc fashion, requiring only knowledge of the model loss gradients.

Influence functions enable efficient computation of the effect of a training data point (x i , y i ) on θ. This is achieved by evaluating the change in θ if (x i , y i ) was up-weighted2 by some small , resulting in a new parameter θi, arg min θ∈Θ L(D, θ) + • (y i , f (x i ; θ)). The rate of change in θ due to an infinitesimal perturbation in data point i is give by the (first-order) influence function, defined as

EQUATION

)

Note that the model parameter θ is a statistical functional of the data distribution P. Perturbing the i th training point is equivalent to perturbing P to create a new distribution

P i, = (1 -) P + ∆(x i , y i ),

where ∆(x i , y i ) denotes the Dirac distribution in the point (x i , y i ). In this sense, the influence function in ( 7) operationalizes the concept of derivatives to statistical functionals, i.e., it quantifies the derivative of the learned model parameters θ with respect to the data distribution P.

By recognizing that influence functions are the "derivatives" of the model parameters θ with respect to the data distribution P, we can use a Taylor-type expansion to represent the counter-factual model parameter θi, (that would have been learned from a dataset with the i th data point up-weighted) in terms of the true model parameter θ (learned from the true training dataset D) as follows:

EQUATION

where

I (k)

θ (x i , y i ) is the k th order influence function, defined as

I (k) θ (x i , y i ) = ∂ k θi, /∂ k | =0 .

The expansion in (8), known as the von Mises expansion (Fernholz (2012)), is a distributional analog of the Taylor expansion for a statistical functional. If all of the higher-order influence functions (HOIFs) in ( 8) exist, then we can recover θi, without re-training the model on the perturbed training dataset. In Section 3.3, we will derive a recursive formula for efficiently computing the HOIFs in (8). Since exact reconstruction of θi, requires an infinite number of HOIFs, we can only approximate θi, by including a sufficient number of HOIF terms from the von Mises expansion. We discuss the impact of this approximation on the coverage and discrimination performance of the DJ in Section 3.4.

Note that the LOO model parameters { θ-i } n i=1 (required for the construction of the DJ confidence intervals) can be obtained by setting = -1/n, i.e., θ-i = θi, -1 n , since removing a training point is equivalent to up-weighting it by -1/n in the loss function L(D; θ). Thus, by setting = -1/n and selecting a prespecified number of HOIF terms m for obtaining the approximate LOO parameters θ(m) -i , the DJ confidence intervals can be computed using the steps in Algorithm 1.

Algorithm 1: The Discriminative Jackknife

Input :Learned parameter θ, influence order m, coverage α, training data D, test point x Output :Approximate DJ confidence interval C (m) α,n (x; θ) 1: for i = 1 to n do 2: θ(m) -i ← θ -m k=1 (n -k /k!) • I (k) θ (xi, yi) 3: ri ← Yi -f (xi; θ(m) -i ) 4: vi(x) ← f (x; θ) -f (x; θ(m) -i ) 5: end 6: fL(x; θ) ← f (x; θ) -Qα,n(Rn) -Q 1-α,n -(V - n (x)) 7: fU (x; θ) ← f (x; θ) + Qα,n(Rn) + Q α,n + (V + n (x)) 8: return: C (m) α,n (x; θ) ← [ fL(x; θ), fU (x; θ) ]

3.3. COMPUTING HIGHER-ORDER INFLUENCE FUNCTIONS:

How can we compute the influence functions of the von Mises expansion in (8)? The recent work on model interpretability in (Koh & Liang (2017) ) has studied the use of influence functions to quantify the impact of each training data point on a model's predictions. In this work, first-order influence functions were computed -using the classical result in (Cook & Weisberg (1982) ) -as follows:

EQUATION

where

H θ ∇ 2 θ (y, f (x, θ))

is the Hessian of the loss function, which is assumed to be positive definite. However, the expression in ( 9) is only limited to first-order influence functions. In the following Theorem3 , we derive a recursive formula for computing HOIFs of the model parameters θ.

Theorem 1. Let (.) be a differentiable function with a positive definite Hessian matrix H θ . If (.) is locally convex in the neighborhood of θ, then the (k + 1)-th order influence function of θ exists for all (x, y) ∈ X × Y, and for all k ≥ 1 we have

I (k+1) θ (x, y) = -H -1 θ k m=1 a m • g m I (j) θ (x, y), ∇ j θ (y, f (x, θ), ∇ j θ L(D, θ)) m j=1

, where a m is a constant, and g m is a polynomial function, ∀m ∈ {1, . . . , m}.

Similar expressions have been arrived at in (Giordano et al. (2019) ; Debruyne et al. (2008) ). We give the precise definition of the constant term a m (.) and the function g m (.) -in terms of Bell polynomials -in Appendix B. A key consequence of Theorem 1 is that to compute all HOIFs, we need to evaluate the inverse of the Hessian matrix H -1 θ only once, hence evaluating HOIFs scales linearly with the number of training examples. With automatic differentiation tools at our disposal (Paszke et al. (2017) ), computing the loss gradients (arguments of g m (.)) can be done efficiently.

As it is the case in first-order influence functions, inverting H θ is the bottleneck operation in evaluating HOIFs, hence much of the machinery involved in (Koh & Liang (2017)) can be imported into our framework. With p parameters in θ, inverting H θ requires O(np 2 + p 3 ) operations, which can be cumbersome for deep learning models with large number of parameters. Following (Koh & Liang (2017)), we avoid explicitly computing H θ , and instead use implicit Hessian-vector products to efficiently approximate the Hessian products in (9). For models with a very large p, we can use the stochastic estimation approach in (Agarwal et al. (2016) ) to estimate H -1 θ with O(p) complexity.

3.4. THEORETICAL GUARANTEES:

Now that we have developed a practical algorithm for implementing the DJ, we conclude this Section by revisiting the performance requirements in Section 2. In the following Theorem, we show that the approximate DJ confidence intervals satisfy the desired coverage and discrimination requirements.

Theorem 2. Let C (m) α,n (x; θ) be the approximate DJ confidence interval evaluated with m HOIFs using Algorithm 1. Then, for any model f (x; θ), the coverage rate achieved by the DJ is

lim m→∞ P y ∈ C (m) α,n (x; θ) ≥ 1 -α. If the model is stable, i.e., P{ |f (x; θ) -f (x; θ-i )| ≤ ε } ≥ 1 -δ, ∀i,

for some ε and δ ∈ (0, 1) (Bousquet & Elisseeff (2002) ), then ∀x, x ∈ X , y, y ∈ Y, and for a sufficiently large n, we have

E[ W ( C (∞) α,n (x; θ)) ] ≥ E[ W ( C (∞) α,n (x ; θ)) ] ⇔ E[ (y, f (x; θ)) ] ≥ E[ (y , f (x ; θ)) ].

Theorem 2 provides a very strong guarantee on the frequentist coverage of the DJ confidence intervals. That is, without any assumptions, the exact DJ interval C (∞) α,n (x; θ) has a worst case coverage that fulfills the target coverage rate of (1 -α). (The approximate DJ intervals, C (m) α,n (x; θ), achieve the target coverage in the limit of m → ∞.) In Section 5, we show through empirical evaluation that in practice -even with a finite number of HOIFs -the DJ intervals will achieve the (1 -α) coverage. With further assumption on algorithmic stability (Bousquet & Elisseeff (2002) ), Theorem 2 states that if a large number of HOIFs are involved in constructing the approximate DJ intervals, then the discrimination condition will be satisfied if we have a sufficiently large training sample.

4. RELATED WORK:

The DJ is a post-hoc method that operates on deep learning models after being trained. Post-hoc methods are advantageous in that they neither compromise the model accuracy, nor require any modifications in model structure or training. (To implement Algorithm 1 (Section 3.3), we only need access to the model's loss gradients and Hessian-vector products.) It is worth noting that, while our primary motivation is to estimate uncertainty in deep learning models, our method is applicable to any machine learning model for which HOIFs are accessible. Additionally, the exact DJ procedure can be applied to any model via exhaustive leave-one-out re-training, even when gradients are inaccessible.

Post-hoc methods have been traditionally underexplored since existing approaches, such as calibration via temperature scaling (Platt et al. (1999) ), are known to under-perform compared to built-in methods (Ovadia et al. (2019) ). However, recent works have revived post-hoc approaches using ideas based on bootstrapping (Schulam & Saria (2019) ), jackknife resampling (Barber et al. (2019b) ; Giordano et al. (2018) ) and cross-validation (Vovk et al. (2018) ; Barber et al. (2019a) ). An overview of the different classes of post-hoc and built-in methods proposed in recent literature is provided in Table 1 . 2019)). By their very nature, Bayesian models cannot be applied in a post-hoc fashion since they require specifying priors over model parameters, which leads to major modifications in the inference algorithms. While Bayesian models provide a formal framework for uncertainty estimation, posterior credible intervals do not guarantee frequentist coverage, and more crucially, the achieved coverage can be very sensitive to hyper-parameter tuning (Bayarri & Berger (2004) ). Moreover, exact Bayesian inference is often computationally prohibitive, and alternative approximations -e.g., (Gal & Ghahramani (2016) ) and (Kingma et al. ( 2015))) -may induced posterior distributions that do not concentrate asymptotically (Osband (2016) ; Hron et al. (2017) ).

Deep ensembles (Lakshminarayanan et al. ( 2017)) are currently considered the most competitive (non-Bayesian) benchmark for uncertainty quantification (Ovadia et al. (2019) ). This method repeatedly retrains the model on sub-samples of the data (using adversarial training), and then estimates uncertainty through the variance of the aggregate predictions. A similar bootstrapping approach, developed in (Schulam & Saria (2019) ), uses the model's Hessian and loss gradients to create an ensemble without re-training. While these methods may perform favorably in terms of discrimination, they are likely to undercover, since they only consider local variability terms akin to those in (5). Additionally, ensemble methods do not provide theoretical guarantees on coverage and discrimination.

The (infinitesimal) jackknife method was previously applied for quantifying the predictive uncertainty in random forest models (Wager et al. (2014) ; Mentch & Hooker (2016) ; Wager & Athey (2018) ). In these works, however, the developed jackknife estimators are bespokely tailored to bagging predictors, and cannot be straightforwardly extended to deep neural networks. More recently, general-purpose jackknife estimators were developed in (Barber et al. (2019b)) , where two exhaustive leave-oneout procedures: the jackknife+ and the jackknife-minmax where shown to have assumption-free worst-case coverage guarantees of (1 -2α) and (1 -α), respectively. Our work improves on these results in two ways. First, the DJ achieves the target coverage rate of (1 -α) with a strictly narrower interval width compared to the jackknife-minmax (proof is provided in Appendix D). Second, our proposed HOIF-based implementation of the jackknife alleviates the need for exhaustive leave-one-out re-training, which enables scaling up any jackknife-based method to large datasets. 

5. EXPERIMENTS:

Our DJ procedure (Algorithm 1) was implemented on top of PyTorch models, where we computed the loss gradients and Hessian required for evaluating the HOIFs of a model's parameters (Theorem 1) using autograd (Maclaurin et al. (2015) ). The only hyper-parameter involved in our method is the number of HOIFs m -this was tuned by optimizing the evaluation metrics (listed below) on the training samples. Code snippets for the DJ procedure are provided in the supplementary appendix.

Baselines. We compared our DJ procedure with 6 state-of-the-art baseline methods. These included 3 built-in Bayesian methods: (1) Monte Carlo Dropout (MCDP) (Gal & Ghahramani (2016) ), (2) Probabilistic backpropagation (PBP) (Hernández-Lobato & Adams (2015)), and (3) Bayesian neural networks with inference via stochastic gradient Langevin dynamics (BNN-SGLD) (Welling & Teh (2011) ). In addition, we considered deep ensembles (DE) (Lakshminarayanan et al. ( 2017)), which is deemed the most competitive built-in frequentist method (Ovadia et al. (2019) ). We also included 2 post-hoc methods: resampling uncertainty estimation (RUE) (Schulam & Saria (2019) ) and crossconformal predictive distributions (CCP) (Vovk et al. (2018) ). For a target coverage of (1 -α), uncertainty estimates were obtained by setting the posterior quantile functions to (1-α) (for Bayesian methods), or obtaining the (1 -α) percentile point of a normal distribution (for frequentist methods). We implemented all baselines in PyTorch (details are provided in the supplementary appendix).

Evaluation metrics. In all experiments, we used the mean squared error (MSE) as the loss L(D, θ) for the trained model f (x; θ). To ensure a fair comparison, the hyper-parameters of the model f (x; θ) were the same for all baselines. In each experiment, an uncertainty estimate C α,n (x; θ) is obtained from a training sample, and then coverage and discrimination are evaluated on a test sample. To evaluate empirical coverage probability, we compute the fraction of test samples for which y resides in C α,n (x; θ). Discrimination was evaluated as follows. For each baseline method, we evaluate the interval width W ( C α,n (x; θ)) for all test points. For a given error threshold E, we use the interval width to detect whether the test prediction f (x; θ) is a high-confidence, i.e., (y, f (x; θ)) ≤ E, or a low-confidence prediction, i.e., (y, f (x; θ)) > E. We use the area under the ROC curve (AUC-ROC) to evaluate the accuracy of classifying high-and low-confidence predictions using W ( C α,n (x; θ)).

5.1. SYNTHETIC DATA:

Qualitative assessment. First, we qualitatively assess the DJ confidence intervals using the synthetic data model in (Hernández-Lobato & Adams (2015) ). This dataset is based on the one-dimensional model: y = x 3 + , where ∼ N (0, 3 2 ). We fit a feed-forward neural network model (with the same architecture as (Hernández-Lobato & Adams (2015))) on 50 training examples drawn from a uniform distribution on X = [-4, 4], and evaluate confidence intervals on its predictions with a target coverage of 90%. The confidence intervals constructed by all baselines are visualized in Figure 3 .

As we can see in Figure 3 , Bayesian methods tend to undercover the target function -i.e., PBP overestimated model confidence, whereas BNN-SGLD failed to recognize that the region |x| > 4 should be assigned higher uncertainty since it contains no training data. For MCDP, we can see that higher uncertainty is assigned to the region |x| > 4, but the posterior does not concentrate around the training data, and hence credible intervals are poorly calibrated. Post-hoc methods (RUE and CCP) tend to overestimate confidence in regions with dense training data (|x| ≤ 4) since they rely mainly on the (in-sample) error variance to estimate uncertainty. DE produces a favorable confidence interval, but with little variation in interval width and with few training points being uncovered.

The DJ confidence intervals hit the right balance between coverage and discrimination. The addition of the marginal error Q0.1,50(R50) in ( 5) ensures that the training data points and target function are covered everywhere in X . The variability terms, Q0.1,23(V + 50 ) and Q0.9,27(V - 50 ), add the discrimination components, which enable the interval width to clearly distinguish high-confidence regions (|x| ≤ 4) from low-confidence ones (|x| > 4). By contrasting the DJ estimates with m = 1 and m = 2, we see that adding more HOIF terms gives smoother and tighter confidence intervals. Quantitative assessment. Next, we examine coverage and discrimination performance of all baselines on the following synthetic model. For n = 100 and 1 -α = 0.9, we conduct 100 experiments, each for different feature dimensions d = 1, 11, 21, . . . , 101, with i.i.d. samples (x i , y i ) generated as β ∼ U nif orm[ 0, 1 ] ⊗d , x i ∼ N (0, I d ), and y i | x i ∼ N ((x T i β) 3 , 1). In each experiment, we fit a 2-layer neural network with 100 hidden units to the training data, and then generate 100 test data points from the same distribution to evaluate the coverage and discrimination performance for different values of d. (Discrimination is evaluated by setting the error threshold E to be the 95% percentile of the trained model error on the test samples.) As we can see in Figure 4 , DJ is the only method that always manages to achieve the 90% coverage rate for any number of features d (as predicted by Theorem 2). It does so while also achieving the most competitive discrimination (AUC-ROC) performance, together with the RUE and PBP baselines. We also found that, for all values of d, the DJ interval width was strictly tighter than that of the jackknife-minmax (JM) procedure (Barber et al. (2019b) ), which is the only other procedure that guarantees the 1 -α target coverage.

In this experiment, we found that m = 2 HOIFs are sufficient for achieving the best possible performance by DJ. More HOIFs would generally be needed for more complex networks, with more layers or hidden units (see further results in the supplementary appendix). Because of its post-hoc nature, the computational complexity of DJ was significantly less than built-in method (average runtime of DJ was 50 times faster than DE, 30 times faster than BNN-SGLD and 10 times faster than PBP).

5.2. REAL DATA: AUDITING MODEL RELIABILITY:

We conducted a series of experiments on real-world datasets in order to evaluate the accuracy of uncertainty estimates issued by the DJ procedure in various contexts. In this Section, we show how uncertainty estimates can be used to audit the reliability of a predictive model. Further experiments on active learning and out-of-distribution detection are deferred to the supplementary appendix.

In this experiment, we compare how well DJ and the competing baselines are able to detect when a model's prediction will be far from the true target value. We conduct our experiments on 8 common UCI benchmark datasets for regression -those included: Boston housing, concrete compression strength, energy efficiency, robot kinematics (Kin8nm), naval propulsion, combined cycle power plant, red wine quality, and yacht hydrodynamics (Dua & Graff (2017) ). Table 2 : AUC-ROC performance (± 95% confidence intervals) of all baselines on the real-world UCI regression datasets. Numbers between brackets are the achieved (empirical) coverage rates. Bold numbers correspond to the best performing method. Coverage rates marked with an asterisk achieve the 90% target rate.

For each dataset, we sample 20 train-test splits, and use 90% of the data points for training. To control for model complexity, we fix a single feed-forward neural network architecture with two hidden layers and 100 hidden units, logistic activation functions, MSE loss, and a single set of learning hyper-parameters (500 epochs with 100 samples per minibatch, and an Adam optimizer with default settings). We fit a single model for each split of each dataset, and compute the confidence intervals (with 1 -α = 0.9) for all baselines on the test data. We implemented two versions of our method: a DJ with a first-order von Mises approximation (m = 1), and a DJ with m = 10 HOIFs. Similar to the previous experiment, we evaluate the empirical coverage rate and discrimination (AUC-ROC) performance assuming the threshold E to be the 95% percentile of the empirical distribution over test errors. Results are provided in Table 2 .

We observe that, in terms of discriminative performance, DJ outperforms (or is competitive with) the baselines -the AUC-ROC of detecting low-quality predictions by DJ was better than other methods in 6 out of the 8 datasets. Including HOIFs was important for achieving good performance in most datasets. A first-order von Mises approximation was sufficient only in the "Energy" dataset, but under-performed compared to the DJ with m = 10 in all other datasets. Another key observation is that -as predicted by Theorem 2 -the DJ procedure always achieves the 90% target coverage rate. All other methods tend to undercover the test points, and some methods (such as MCDP and PBP) can even achieve close to zero coverage. It is worth mentioning that on the dataset where DJ provided the biggest improvement in discriminative performance (Boston housing), it was also the only method that achieved the target coverage.

A INFLUENCE FUNCTIONS: BACKGROUND & KEY CONCEPTS:

A.1 FORMAL DEFINITION Robust statistics is the branch of statistics concerned with the detection of outlying observations. An estimator is deemed robust if it produces similar results as the majority of observations indicates, regardless of how a minority of other observations is perturbed (Huber & Ronchetti (1981) ). The influence function measures these effects in statistical functionals by analyzing the behavior of a functional not only at the distribution of interest, but also in an entire neighborhood of distributions around it. Lack of model robustness is a clear indicator of model uncertainty, and hence influence functions arise naturally in our method as a (pointwise) surrogate measure of model uncertainty. In this section we formally define influence functions and discuss its properties.

The pioneering works in (Hampel et al. (2011) ) and (Huber & Ronchetti (1981) ) coined the notion of influence functions to assess the robustness of statistical functionals to perturbations in the underlying distributions. Consider a statistical functional T : P → R, defined on a probability space P, and a probability distribution P ∈ P. Consider distributions of the form P ε,z = (1 -ε)P + ε∆z where ∆z denotes the Dirac distribution in the point z = (x, y), representing the contaminated part of the data.

For the functional T to be considered robust, T (P ε,z ) should not be too far away from T (P) for any possible z and any small ε. The limiting case of ε → 0 defines the influence function. That is, Then the influence function of T at P in the point z is defined as

I(z; P) = lim ε→0 T (P ε,z ) -T (P) ε ∂ ∂ε T (P ε,z ) ε=0 , (A.1)

The influence function measures the robustness of T by quantifying the effect on the estimator T when adding an infinitesimally small amount of contamination at the point z. If the supremum of I(.) over z is bounded, then an infinitesimally small amount of perturbation cannot cause arbitrary large changes in the estimate. Then small amounts of perturbation cannot completely change the estimate which ensures the robustness of the estimator.

A.2 THE VON MISES EXPANSION:

The von Mises expansion is a distributional analog of the Taylor expansion applied for a functional instead of a function. For two distributions P and Q, the Von Mises expansion is (Fernholz (2012)):

EQUATION

where I (k) (z; P) is the k th order influence function. By setting Q to be a perturbed version of P, i.e., Q = P ε , the von Mises expansion at point z reduces to: A.3) and so the k th order influence function is operationalized through the derivative

EQUATION

I (k) (z; P) ∂ ∂ k ε T (P ε,z ) ε=0 . (A.4) A.3 INFLUENCE FUNCTION OF MODEL LOSS

Now we apply the mathematical definitions in Sections A.1 and A.2 to our learning setup. In our setting, the functional T (.) corresponds to the (trained) model parameters θ and the distribution P. In this case, influence functions of θ computes how much the model parameters would change if the underlying data distribution was perturbed infinitesimally.

I:

(1)

θ (z) = ∂ θε,z ∂ ε ε=0 , θε,z arg min θ∈Θ 1 n n i=1 (z i ; θ) + ε (z; θ). (A.5)

Recall that in the definition of the influence function P ε,z = (1 -ε)P + ε∆z where ∆z denotes the Dirac distribution in the point z = (x, y). Thus, the (first-order) influence function in (A.5) corresponds to perturbing a training data point z by an infinitesimally small change ε and evaluating the corresponding change in the learned model parameters θ. More generally, the k th order influence function of θ is defined as follows:

I (k) θ (z) = ∂ k θε,z ∂ ε k ε=0 . (A.6)

By applying the von Mises expansion, we can approximate the parameter of a model trained on the training dataset with perturbed data point z as follows:

θε,z ≈ θ + ε I (1) θ (z) + ε 2 2 I (2) θ (z) + . . . + ε m m! I (m) θ (z), (A.7)

where m is the number of terms included in the truncated expansion. When m = ∞, the exact parameter θε,z without the need to re-train the model.

A.4 CONNECTION TO LEAVE-ONE-OUT ESTIMATORS:

Our uncertainty estimator depends on perturbing the model parameters by removing a single training point at a time. Note that removing a point z is the same as perturbing z by ε = -1 n , hence we obtain an (m th order) approximation of the parameter change due to removing the point z as follows:

θ-z -θ ≈ -1 n I (1) θ (z) + 1 2n 2 I (2) θ (z) + . . . + (-1) m n m • m! I (m) θ (z), (A.8)

where θ-z is the model parameter learned by removing the data point z from the training data.

B PROOF OF THEOREM 1:

Recall that the LOO parameter θi, is obtained by solving the optimization problem:

θi, = arg min θ∈Θ L(D, θ) + • (y i , f (x i ; θ)). (B.1)

Let us first derive the first order influence function I (1) (x i , y i ). Let us first define ∆ i, θi, -θ. The first order influence function is given by:

I (1) (x i , y i ) = ∂ θi, ∂ = ∂∆ i, ∂ . (B.2)

Note that, since θi, is the minimizer of (B.1), then the perturbed loss has to satisfy the following (first order) optimality condition:

∇ θ {L(D, θ) + • (y i , f (x i ; θ))} θ= θi, = 0. (B.3)

Since lim →0 θi, = θ, then we can write the following Taylor expansion:

∇ θ ∞ k=0 ∆ k i, k! • ∇ k θ L(D, θ) + • (y i , f (x i ; θ)) = 0. (B.4)

Now by dropping the o( ∆ i, ) terms, we have:

∇ θ L(D, θ) + • (y i , f (x i ; θ)) + ∆ i, • ∇ θ L(D, θ) + • (y i , f (x i ; θ)) = 0. (B.5)

Since θ is a indeed a minimizer of the loss function (.), then we have ∇ θ (.) = 0. Thus, (B.5) reduces to the following condition:

• ∇ θ (y i , f (x i ; θ)) + ∆ i, • ∇ 2 θ L(D, θ) + • ∇ 2 θ (y i , f (x i ; θ)) = 0. (B.6)

By solving for ∇ θ , we have

∆ i, = -∇ 2 θ L(D, θ) + • ∇ 2 θ (y i , f (x i ; θ)) -1 • • ∇ θ (y i , f (x i ; θ)) , (B.7)

which can be approximated by keeping only the O( ) terms as follows:

∆ i, = -∇ 2 θ L(D, θ) -1 • • ∇ θ (y i , f (x i ; θ)) . (B.8)

Noting that ∇ 2 θ L(D, θ) is the Hessian matrix H θ , we have:

∆ i, = -H -1 θ • • ∇ θ (y i , f (x i ; θ)) . (B.9)

By taking the derivative with respect to , we arrive at the expression for first order influence functions:

I (1) (x i , y i ) = ∆ i, =0 = -H -1 θ • ∇ θ (y i , f (x i ; θ)). (B.10)

Now let us examine the second order influence functions. In order to obtain I (2) (x i , y i ), we need to differentiate (B.6) after omitting the O( ) once again as follows:

2∆ i, • • ∇ 2 θ (y i , f (x i ; θ)) + ∆ 2 i, • ∇ 2 θ L(D, θ) + ∆ i, • ∆ i, • ∇ 3 θ L(D, θ) = 0. (B.11)

Where we have applied the chain rule to obtain the above. By substituting ∇ 2 θ L(D, θ) = H θ and dividing both sides of (B.11) by 2 , we have

2 ∆ i, • ∇ 2 θ (y i , f (x i ; θ)) + ∆ 2 i, 2 • H θ + ∆ i, 2 • ∇ 3 θ L(D, θ) = 0. (B.12)

Thus, by re-arranging (B.11), we can obtain I (2) (x i , y i ) in terms of I (1) (x i , y i ) as follows:

I (2) (x i , y i ) = -H -1 θ I (1) (x i , y i ) 2 • ∇ 3 θ L(D, θ) + 2 I (1) (x i , y i ) • ∇ 2 θ (y i , f (x i ; θ)) .

Similarly, we can obtain the k th order influence function, for any k > 1, by repeatedly differentiating equation (B.6) k times, i.e.,

∂ ∂ k • ∇ θ (y i , f (x i ; θ)) + ∆ i, • ∇ 2 θ L(D, θ) = 0. (B.13)

and solving for ∂∆ k i, /∂ k . By applying the higher-order chain rule to the two added terms in (B.13) (or equivalently, take the derivative of I (2) (x i , y i ) for k -2 times), we can represent the formula for I (k) (x i , y i ) using Bell polynomials (Ma (2013) ) as follows:

I (k+1) θ (x, y) = -H -1 θ k m=1 B m,k+1 I (j) θ (x, y), ∇ j θ (y, f (x, θ)) m j=1 + k m=1 B m,k I (j) θ (x, y), ∇ j θ L(D, θ) m j=1 , (B.14)

where B m,k is the Bell polynomial, which is defined as follows:

B m,k (x 1 , x 2 , . . . , x m-k+1 ) = m! j 1 !j 2 ! • • • j m-k+1 ! x 1 1! j1 x 2 2! j2 • • • x m-k+1 (m -k + 1)! j m-k+1 . and Q + α,n ({r i } n i=1 ) + v j = K α n+1 , we have (1 -α)(n + 1) ≤ n j=1 1 K n+1,j > K α n+1 ,

which implies that n + 1 ∈ Q(K), and hence

P y n+1 / ∈ C n,α (x n+1 ) ≤ P {n + 1 ∈ Q(K)} .

which concludes the proof of the Lemma. 

P {n + 1 ∈ Q(K)} = P {j ∈ Q(K)} ,

for all j ∈ {1, . . . , n}. Based on the above, we can then evaluate the coverage probability using the equivalent condition in Lemma 1 as follows:

P {n + 1 ∈ Q(K)} = 1 n + 1 n+1 j=1 P {j ∈ Q(K)} = 1 n + 1 E   n+1 j=1 1 {j ∈ Q(K)}   = 1 n + 1 E [ |Q(K)| ] . (C.3)

Finally, we compute the expected size of the set Q(K) as follows. By combining the original coverage condition in (C.2) together with the equivalence result in Lemma 1, we have

(1 -α)(n + 1) ≤ n+1 j=1 1 {K ji < K α i } ≤ n+1 j=1 1 {j / ∈ Q(K)} = n + 1 -|Q(K)|,

which implies that |Q(K)| ≤ α(n + 1). By combining this with (C.3), we arrive at the following inequality

P y n+1 / ∈ C n,α (x n+1 ) ≤ α.

The result of the Theorem follows from the fact that the approximate DJ confidence interval C (m) n,α (x n+1 ) converges in probability to the exact interval as m → ∞.

C.2 DISCRIMINATION:

Recall that the exact DJ interval width is given by:

W ( C (∞) α,n (x; θ)) = 2 Q α,n (R n ) + 2 Q α,n (V n (x)). (C.4)

Since the term Q α,n (R n ) is constant for any x, discrimination boils down to the following condition:

E[ Q α,n (V n (x)) ] ≥ E[ Q α,n (V n (x )) ] ⇔ E[ (y, f (x; θ)) ] ≥ E[ (y , f (x ; θ)) ].

(C.5) Note that to prove the above, it suffices to prove the following:

E[ v i (x) ] ≥ E[ v i (x ) ] ⇔ E[ (y, f (x; θ)) ] ≥ E[ (y , f (x ; θ)) ]. (C.6)

If the model is stable (based on the definition in (Bousquet & Elisseeff (2002) )), then a classical result by (Devroye & Wagner (1979) ) states that:

E[ | (y, f (x; θ)) -n (y, f (x; θ))| 2 ] ≈ E[ | (y, f (x; θ)) -(y, f (x; θ-i ))

| 2 ] + Const., (C.7) as n → ∞, where n (.) is the empirical risk on the training sample, and the expectation above is taken over y | x. From (C.7), we can see that an increase in the LOO risk (y, f (x; θ-i )) implies an increase in the empirical risk n (y, f (x; θ)), and vice versa. Thus, for any two feature points x and x , if v(x) is greater than v(x ), then on average, the empirical risk at x is greater than that at x .

D EXPERIMENTAL DETAILS:

In this Section, we provide more details for the experimental setups in Section 5, in addition to further experimental results in the context of active learning and out-of-distribution detection.

D.1 CODE SNIPPETS FOR THE DJ PROCEDURE:

The DJ procedure is very simple to implement. The key component of its implementation is the computation of HOIFs. Let model be an instance of the nn module of Pytorch, with a set of parameters model.parameters() that correspond to θ in our formulation. Below is a Python code snippet for computing the first-order influence function of model.

. Computing the HOIFs can be easily implemented by recursively evaluating the gradients and lower order influence function as stated in Theorem 1. Once all the ingredients of the HOIFs are evaluated, then the perturbed parameters can be evaluated as explained in the paper, and a modified LOO model is created for every sample in the training data.

D.2 IMPLEMENTATION OF BASELINES:

In what follows, we provide details for the implementation and hyper-parameter settings for all baseline methods involved in Section 5.

Probabilistic backpropagation (PBP). We implemented the PBP method proposed in (Hernández-Lobato & Adams (2015)) with inference via expectation propagation using the theano code provided by the authors in (github.com/HIPS/Probabilistic-Backpropagation). Training was conducted via 40 epochs.

Monte Carlo Dropout (MCDP). We implemented a Pytorch version of the MCDP method proposed in (Gal & Ghahramani (2016) ). In all experiments, we tuned the dropout probability using Bayesian optimization to optimize the AUC-ROC performance on the training sample. We used 1000 samples at inference time to compute the mean and variance of the predictions. The credible intervals were constructed as the (1 -α) quantile function of a posterior Gaussian distribution defined by the predicted mean and variance estimated through the Monte Carlo outputs. Similar to the other baselines, we conducted training via 40 epochs for the SGD algorithm.

Bayesian neural networks (BNN). We implemented BNNs with inference via stochastic gradient Langevin dynamics (SGLD) (Welling & Teh (2011) ). We initialized the prior weights and biases through a uniform distribution over [-0.01, 0.01]. We run 40 epochs of the SGLD inference procedure and collect the posterior distributions to construct the credible intervals.

Deep ensembles (DE).:

We implemented a Pytorch version of the DE metho (without adversarial training)d proposed in (Lakshminarayanan et al. ( 2017)). We used the number of ensemble members M = 5 as recommended in the recent study in (Ovadia et al. (2019) ). Predictions of the different ensembles were averaged and the confidence interval was estimated as 1.645 multiplied by the empirical variance for a target coverage of 90%. We trained the model through 40 epochs.

Resampling uncertainty estimation (RUE). We implemented RUE (Schulam & Saria (2019) ) in Pytorch using the same Hessian computation routines involved in implementing our method. We tuned the dampening factor λ in the stabilized Hessian H + λ I used for matrix inversion. We used 20 bootstrapped sub-samples of the training data to compute the prediction variance.

D.3 FURTHER RESULTS ON SYNTHETIC DATA:

The minimum number of HOIFs m * required to achieve the best possible AUC-ROC for the quantitative experiment (d = 20) in Section 5.1 is shown in Figure 5 . As expected, more complex models will generally require more high order terms of the von Mises expansion in order to recover the parameters of the LOO re-trained model. 

D.4 FURTHER EXPERIMENTS WITH REAL DATA:

We performed another series of experiments to evaluate the accuracy of the estimates of the uncertainty estimates produced by different methods. We study two particular applications of uncertainty quantification: active learning and out-of-distribution detection (OOD).

Active Learning. For this setup, we replicate the active learning experiment in (Hernández-Lobato & Adams (2015) ). In active learning, it is necessary to produce accurate estimates of uncertainty in order to decide which data samples we should actively collect. In this experiment, we used a neural network with a single hidden layer and 50 hidden units. We split each data set into training and test sets with 20 and 100 data instances, respectively, and pool sets with all the remaining data. After this, one data point is collected from the pool set and then moved into the training set. The process repeats until 9 of these active additions to the training set have been completed, that is, until we have performed 10 evaluations on the test set. The entire process, including the random data set splitting, is repeated 40 times. The pool data is initially lacking the target variables and these become available only once the data is moved to the training set.

We compare a baseline neural network model trained on random selections from the pool data with version of the model trained by data collected actively (PBP-A, MCDP-A, DE-A, BNN-A and DJ-A). Results are provided in Table 3 . As we can see, DJ-based active collection of data is significantly better than the random baseline model on all datasets. Moreover, data actively collected using DJ-based uncertainty estimates were superior to all other methods on 6 out of the 8 datasets. Out-of-distribution detection. Next, we evaluate the accuracy of the DJ uncertainty estimates in detecting whether a new test example came from the same distribution as the training data, or from another unknown distribution. This task is known as Out-of-distribution detection (OOD), and is commonly used to benchmark uncertainty estimation methods (Ovadia et al. (2019) ).

In this experiment, we use data from the Infant Health and Development Program (IHDP) (Alaa & van der Schaar (2017)). The IHDP is an interventional program intended to enhance the health of premature infants. This dataset comprises observational data with treatment assignments from a real-world clinical trial, and introduced selection bias to the data artificially by removing a subset of the patients. Every sample in the dataset is of the form (x, y, w), where w ∈ {0, 1} is a binary treatment assignment variable which is equal to 1 if the patient was treated, and is 0 otherwise. Since patients are not treated at random, i.e., the likelihood of a patient getting treated depends on their feature x, the treated and untreated patient populations have different distributions, i.e., we have P(x | w = 0) = P(x | w = 1). The entire dataset comprised 747 subjects, with 25 features for each subject, of which 165 were treated and the rest were untreated. We conducted our experiment as follows. We train the model on 100 samples for treated patients, and leave out 65 samples for testing. In testing, we use the interval width to predict whether the test samples comes from the 65 treated patient (whose features resemble the training data), or from the untreated patients who were not involved in training. We repeated the experiment 20 times with random train/test splitting. The model we used is a feed-forward neural network with 2 hidden layers and 100 hidden units. AUC-ROC performance is reported in Table 4 .

We found that DJ outperforms other methods in detecting the OOD samples, with DE and RUE being the most competitive baselines. This result is of practical significance since it implies that we can accurately determine which patients should be treated based on the recommendations of a machine learning-based decision-support system, and which patients should be examined in a case by case basis through experts. Moreover, by determining which samples in an observational dataset have the most unreliable predictions, and consequently determine what features are required in a prospective design of a clinical trial.

:

Note that, from the definition of Qα,n, it follows that Qα,n(R) = Q1-α,n(-R).

Estimators that perturb parameters instead of a LOO construct are known as the infinitesimal jackknife (IJ)(Giordano et al. (2018)). A detailed technical background on influence functions is provided in Appendix A.

All proofs are provided in the supplementary appendix.

CONCLUSIONUncertainty quantification is a crucial requirement in various high-stakes applications where deep learning can inform critical decision-making. In this paper, we developed a rigorous frequentist procedure for quantifying the uncertainty in predictions issued by deep learning models in a post-hoc fashion. Our procedure, which is inspired by classical jackknife estimators, does not require any modifications in the underlying deep learning model, and provides theoretical guarantees on its achieved performance. Implementation of our procedure using higher-order influence functions entails minimal computational costs, and allows our method to scale for big datasets, without the need for repeated re-training as in the classical jackknife procedure. While this paper focuses on uncertainty estimation in deep learning, the methods therein can be applied to any machine learning model as long as we are granted access to its loss gradients and Hessian.

Figures and Tables:

Figure 2: Illustration of the discriminative jackknife.

Figure 3: Uncertainty estimates by all baselines on the synthetic dataset. Training data is shown as red dots and the true function is displayed as a black dotted line. Confidence intervals are shown as a shaded blue region.

Figure 4: Coverage and discrimination performance of all baselines on high-dimensional synthetic data.

Given the result of Lemma 1, we use the exchangeability of the training and test points to further bound the probability P y n+1 / ∈ C n,α (x n+1 ) as follows. Since the training and testing points are exchangeable, we have K d = Π K Π T for any (n + 1) × (n + 1) permutation matrix Π, where d = denotes equality in distribution. Thus, we have

Figure 5: Impact of model complexity on number of HOIFs required for accurate recovery of LOO parameters.

Overview of existing methods for uncertainty quantification in deep learning.

Table content: <table><tr><td>Method</td><td>Bayesian vs. Frequentist</td><td>Post-hoc vs. Built-in</td><td>Coverage</td></tr><tr><td>Bayesian neural nets (Ritter et al. (2018))</td><td>Bayesian</td><td>Built-in</td><td>No guarantees</td></tr><tr><td>Probabilistic backprop. (Blundell et al. (2015))</td><td>Bayesian</td><td>Built-in</td><td>No guarantees</td></tr><tr><td>Monte Carlo dropout (Gal &amp; Ghahramani (2016))</td><td>Bayesian</td><td>Built-in</td><td>No guarantees</td></tr><tr><td>Deep ensembles (Lakshminarayanan et al. (2017))</td><td>Frequentist</td><td>Built-in</td><td>No guarantees</td></tr><tr><td>Resampling uncertainty (Schulam &amp; Saria (2019))</td><td>Frequentist</td><td>Post-hoc</td><td>No guarantees</td></tr><tr><td>Cross-conformal (Vovk et al. (2018))</td><td>Frequentist</td><td>Post-hoc</td><td>1 -2α</td></tr><tr><td>Jackknife+ (Barber et al. (2019b))</td><td>Frequentist</td><td>Post-hoc</td><td>1 -2α</td></tr><tr><td>Jackknife-minmax (Barber et al. (2019b))</td><td>Frequentist</td><td>Post-hoc</td><td>1 -α</td></tr></table>

......

Table content: <table><tr><td>def influence_function(model, training_data_indexes):</td></tr><tr><td>....... IF.... = []</td></tr><tr><td>....... H...... = hessian(model.loss, model.parameters(), create_graph=True)</td></tr><tr><td>....... Lamb = 1e-3</td></tr><tr><td>....... Hinv = torch.inverse(H + Lamb * torch.eye(H.shape[0]))</td></tr><tr><td>.......</td></tr></table>

Average test RMSE and standard errors in active learning.

Table content: <table/>

± 0.01 0.69 ± 0.03 0.72 ± 0.01 0.75 ± 0.02 0.65 ± 0.01 0.75 ± 0.01 0.72 ± 0.02 AUC-ROC accuracy of all baselines on the OOD task.

Table content: <table><tr><td>DJ</td><td>PBP</td><td>MCDP</td><td>DE</td><td>BNN</td><td>RUE</td><td>CCP</td></tr><tr><td>0.78</td><td/><td/><td/><td/><td/><td/></tr></table>

Peer review (iclr2020):
This paper studies how to construct confidence intervals for deep neural networks with guaranteed coverage. The authors propose an algorithm, “discriminative jackknife”, based on the standard jackknife confidence interval estimate which they augment by a “local uncertainty estimate” based on the variability of the n leave-one-out fitted versions of the underlying algorithm (n = # data points). The whole study is concluded with toy and real-world examples showing the proposed algorithm is competitive with existing methods while also achieving the desired coverage. I am currently leaning towards recommending rejection of the paper. The main reasons for this are: (i) A potential failure to cite and acknowledge the prior contribution of [1] which seems to have non-trivial overlap with this paper (on arxiv since end of July which is more than 30 days before the ICLR submission deadline and thus should be treated as prior work); (ii) The claims of guaranteed frequentist coverage are not backed up as, according to thm.2, they only hold when n >> 0 and the number of influence functions used goes to infinity (ideally, the authors would provide non-asymptotic bounds as in [1], but at the very least, these limitations should have been clearly pointed out and their practical implications discussed). Finally, I would like to say that I am holding the paper to a higher standard due to its 10 page length as instructed by the guidelines. Major comments: - Can you please explain the relation of this work to [1]? It seems that [1] already proposes use of the higher order expansions, provides an efficient implementation based on forward mode autodiff (do you plan to release code?), and moreover provides non-asymptotic bounds which are not present in your work?! - On a related note, Giordano et al. provide a careful analysis and discussion of the assumptions in their sect.4. Can you please clarify which of the assumptions you also make, and why you don’t need the others (if any)? - Throughout the paper (e.g., in and around eq.1), you seem to assume that there exists a unique minimiser of the objective which generally won’t be true (especially in your application to deep neural networks). In [1], the technical assumptions ensure this is true but I don’t see how this is handled in your case?! Can you please help me understand how to interpret your work in case there are multiple (possibly local) minima, and whether this has any effect on the results in thm.2? - I think it would be beneficial to the reader if you could please provide a discussion and/or formula for how large n has to be for Theorem 2 to apply. - On p.8, you say “To ensure a fair comparison, the hyper-parameters of the model f(x; ?) were the same for all baselines.” I am not sure I agree this is a fair comparison. Commonly, one has the opportunity to select hyperparameters for their algorithm (e.g., using a validation set, cross-validation, etc.) so it seems it would have been fair to run each algorithm with its best hyperparameters. Can you please provide a discussion of how this would affect the reported results? Minor comments: - On p.3, you say “We do not pose any assumptions on how the loss function in (1) is optimized.” Do you mean to say that you do not assume anything **but** that the chosen optimiser reaches a (global?!) minimum of the loss function? It seems like you would want to exclude pathological optimisers (e.g., one that always outputs zero) but also more realistically think about the known pathologies of commonly used optimisers (see, e.g., [2]). - Can you please clarify if and how the use of the algorithm from (Agarwal et al., 2016) for approximation of the Hessian products affects accuracy of your confidence intervals? - In fig.3, it seems like some of the methods are not properly tuned. For example, MC-dropout should not have zero uncertainty around zero (did you by any chance set bias variance to zero?!), and BNN-SGLD does not seem to have converged (can you please provide plots providing some evidence that the MCMC sampler has mixed + information about how the hyperparameters were selected?). - I am somewhat confused by the statement “The only hyper-parameter involved in our method is the number of HOIFs m — this was tuned by optimizing the evaluation metrics ...” on p.8. Wouldn’t thm.2 suggest that you should use as high m as possible? - Also on p.8, can you please clarify how you selected the threshold for the evaluation of “discriminative power”? In particular, thm.2 seems to suggest that what one would desire is that the ranking based on width of predictive intervals is equivalent to the ranking based on the actual prediction error. This would suggest, for example, comparing these two rankings using Kendall’s tau coefficient (perhaps a binned modification where one would count only the number of times true error puts a point into a different bin than the width of its confidence interval). Note that I am not suggesting the above method is perfect either (it completely ignores the actual sizes of the intervals), but I’m currently having trouble interpreting the results you report, so it would be very helpful to understand how you selected this particular measure of “discriminative power” and why alternatives like the example above were discarded please. References: [1] Ryan Giordano, Michael I. Jordan, Tamara Broderick. A Higher-Order Swiss Army Infinitesimal Jackknife. https://arxiv.org/abs/1907.12116 [2] Ashia C. Wilson, Rebecca Roelofs, Mitchell Stern, Nathan Srebro, Benjamin Recht. The Marginal Value of Adaptive Gradient Methods in Machine Learning. https://arxiv.org/abs/1705.08292