EconBase
← Back to paper

Orthogonal Bootstrap: Efficient Simulation of Input Uncertainty

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

50,286 characters · 14 sections · 65 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Orthogonal Bootstrap: Efficient Simulation of Input Uncertainty

\twocolumn[ \icmltitle{Orthogonal Bootstrap: Efficient Simulation of Input Uncertainty}

icmlauthorlist\icmlauthor{Kaizhao Liu}{yyy} \icmlauthor{Jose Blanchet}{comp1} \icmlauthor{Lexing Ying}{comp} \icmlauthor{Yiping Lu}{sch1,sch2}

\icmlaffiliation{yyy}{Department of Mathematics, Peking University, Beijing, China} \icmlaffiliation{comp1}{Department of Management Science and Engineering, Stanford University} \icmlaffiliation{comp}{Department of Mathematics, Stanford University} \icmlaffiliation{sch1}{Courant Institute of Mathematical Sciences, New York University} \icmlaffiliation{sch2}{Department of Industrial Engineering and Management Sciences, Northwestern University}

\icmlcorrespondingauthor{Kaizhao Liu}{[email removed]} \icmlcorrespondingauthor{Yiping Lu}{[email removed]}

\icmlkeywords{Machine Learning, ICML}

\vskip 0.3in ]

\printAffiliationsAndNotice

abstractBootstrap is a popular methodology for simulating input uncertainty. However, it can be computationally expensive when the number of samples is large. We propose a new approach called Orthogonal Bootstrap that reduces the number of required Monte Carlo replications. We decomposes the target being simulated into two parts: the non-orthogonal part which has a closed-form result known as Infinitesimal Jackknife and the orthogonal part which is easier to be simulated. We theoretically and numerically show that Orthogonal Bootstrap significantly reduces the computational cost of Bootstrap while improving empirical accuracy and maintaining the same width of the constructed interval.

Introduction

The input uncertainty problem lam2018subsampling manifests as the propagation of statistical noise from input models, typically derived and calibrated from data, into the subsequent output analysis. This noise can significantly impact the accuracy and reliability of the resulting conclusions, necessitating robust strategies for quantification and mitigation. Bootstrap stine1985bootstrap,efron1992bootstrap,efron1994introduction is a non-parametric method that uses random resampling with replacement to estimate this uncertainty. In this paper, we consider using Bootstrap to estimate the mean and variance of a function with input uncertainty, which has wide application in debiasing the functional estimation quenouille1949approximate,efron1982jackknife,adams1971asymptotic,cordeiro2014introduction,etter2020operator,jiao2020bias,koltchinskii2021estimation,koltchinskii2022estimation,zhou2021high,etter2021operator,ma2022correcting, improving worst group generalization sagawa2020investigation,nguyen2022improved and, most influential, constructing the confidence interval tukey1958bias,stine1985bootstrap,efron1992bootstrap,efron1994introduction.

Despite its benefits, Bootstrap is computationally demanding. Ideally, to simulate the Bootstrap estimation perfectly, an infinite number of Monte Carlo replications of resamples is required. In practice, only a finite number of Monte Carlo replications can be actually performed, and this introduce an error in the Bootstrap estimation (that can be characterized by its variance), which we call the simulation error. Another source of error in the Bootstrap estimation arises from the randomness in the data. The former can be controlled by increasing the number of Monte Carlo replications, while the latter can not be controlled if the samples are given. To obtain a reasonable Bootstrap estimation, the number of Monte Carlo replications required should at least ensure that the simulation error is smaller than the error arising from the data. In most applications, the number of Monte Carlo replications required scales up with the number of data points lam2018subsampling,lam2022subsampling, where each replication can involve expensive optimization procedures to refit models, making it extremely challenging to scale Bootstrap resampling to large datasets. This drawback also occurs in bagging wager2014confidence.

To address these issues, we propose an alternative method called Orthogonal Bootstrap. We use semiparametric techniques to reduce the number of required Monte Carlo replications for stochastic simulation under input uncertainty. Our method is inspired by recently proposed double/orthogonal machine learning foster2019orthogonal,chernozhukov2018double,chernozhukov2022automatic,chernozhukov2022locally which utilize the influence function cook1982residuals,efron1992jackknife to debias parameter estimates in the presence of nuisance parameters. Suppose we want to estimate the uncertainty of the output using Bootstrap. We can regard the Bootstrap method as a two-stage estimator. In the first stage, we generate resample distributions by resampling the original data with replacement and calculate the output for each resampled distribution. In the second stage, we simply calculate the variance of the outputs from the first stage. We consider the first stage to be “nuisance" because our primary interest is in the final uncertainty, and we do not need to know the simulation output for each resampled distribution. Motivated by this, we show in this paper that we can reduce the simulation error of Bootstrap resampling by dealing with the non-orthogonal and orthogonal parts separately. The non-orthogonal part enjoys a closed form result using the influence function rousseeuw2011robust,cook1980characterizations,koh2017understanding, also known as Infinitesimal Jackknife (IJ) jaeckel1972infinitesimal,giordano2019swiss,giordano2019higher,lu2020uncertainty,alaa2020discriminative,abad2022approximating. For modern machine learning, the influence function can be calculated efficiently using implicit Hessian-vector products and is much faster than retraining the model cook1980characterizations,koh2017understanding,giordano2019swiss. Note that when only calculating the non-orthogonal part, i.e. using the IJ method, the final variance estimate tends to be conservative efron1981jackknife,efron1992jackknife. In Orthogonal Bootstrap, we further simulate the orthogonal part to correct the IJ estimation. Thus, our method provides the same result as the Bootstrap method in contrast to the biased Jackknife estimator. This enables our method to enjoy the accuracy and effectiveness of the Bootstrap method such as higher order coverage for confidence interval construction hall1986bootstrap,efron1992jackknife,hall2013bootstrap while enjoying similar computational cost as the IJ method. We also remark here that, interestingly, our method can also be understood as using IJ as a control variate in the original Bootstrap method.

Related Work

Our paper uses a similar but not identical setting as in lam2022subsampling for the input uncertainty problem. The authors also investigated the simulation effort required by performing Bootstrap. They showed how the total required simulation effort can be reduced from an order bigger than the data size in the conventional approach to an order independent of the data size via subsampling. However, their setting involves two layer of simulation in Bootstrap, while we consider simple plug-in estimator so we consider only one layer of simulation. Interestingly, although their subsampling techniques do not involve influence functions, they leveraged influence functions when proving theoretical results.

A relevant baseline for our paper is the recently proposed Cheap Bootstrap lam2022cheap. The method also aims to maintain desirable statistical guarantees of confidence interval coverage with minimal resampling effort as low as one Monte Carlo replication. However, Cheap Bootstrap enlarges the confidence interval $t_{B,1-\alpha/2}/z_{1-\alpha/2}>1$ times, where $t_{B,1-\alpha/2}$ is the $(1-\alpha/2)-$quantile of $t_B$, the student $t-$distribution with degree of freedom $B$ where $B$ is the number of Monte Carlo replications, and $z_{1-\alpha/2}$ is the $(1-\alpha/2)-$quantile of standard normal. When the number of Monte Carlo replications is constrained, the Cheap Bootstrap will provide a confidence interval with a huge width mean. By contrast, our proposed Orthogonal Bootstrap provides a confidence interval with the same expected width mean as the Standard Bootstrap, even when the number of Monte Carlo replications is constrained. At the same time, the Orthogonal Bootstrap can be used beyond confidence interval construction. It can be used for all settings when the Bootstrap technique is applied, for example, bias reduction for functional estimation quenouille1949approximate,efron1982jackknife,efron1992bootstrap,jiao2020bias,koltchinskii2021estimation,ma2022correcting.

Another relevant paper to ours is kline2012score, which proposed a higher order approximation (up to order $O(n^{-1})$) to the Wild Bootstrap wu1986jackknife,liu1988bootstrap . However, kline2012score's methodology can only be used to construct confidence intervals for linear estimators while orthogonal bootstrap works without assuming any structure of the output functional. Recently, zhou2021high proposed Higher-Order Statistical Expansion (HODSE) based on the closed-form representation of the Jackknife estimator of ideal degenerate expansion of the target statistical functional. However, these types of constructions can only be used for certain statistical models, such as smooth function-of-mean model or linear regression. In contrast, our method functions as a general procedure for Bootstrap methods that doesn't depend on either the linear structure of the estimator or the smooth function-of-mean model. We achieve this by considering the Taylor expansion over the input distribution as a control variate for the Bootstrap estimator.

Contribution

Our contributions are summarized as follows:

itemize{0pt} {0pt} • We propose Orthogonal Bootstrap, a brand new Bootstrap method that provides the same expected simulation result as the original Standard Bootstrap but with reduced simulation effort via separately treating the non-orthogonal part and orthogonal part. • Theoretically, we show that the Orthogonal Bootstrap can provably reduce the simulation error so that the required number of Monte Carlo replications required decreased from $\Omega(n)$ to $O(1)$, under the mild assumption that the performance measure has a continuous Fr\'echet derivative under the Kernel Maximum Mean Discrepancy (MMD) distance. As far as the authors know, we are the first paper to link the Bootstrap simulation error with differentiability in Kernel MMD. • Empirically, we show that Orthogonal Bootstrap can significantly improve the result on both simulated and real datasets when the number of Monte Carlo replications is limited.

Organization of the Paper

We organize our paper as follows. In Section (ref), we demonstrate how Orthogonal Bootstrap can be used to reduce the bias of a statistical estimator. In Section (ref), we demonstrate how Orthogonal Bootstrap can be used to quantify the uncertainty of a statistical estimator. In Section (ref), we present some numerical examples on both simulated and real datasets.

Notation

Throughout the paper, we adopt the following notation conventions lam2022subsampling. The symbol $n$ represents the sample size. The notation $F$ refers to a distribution, $\hat{F}$ stands for the empirical distribution created through sampling from $F$, and $\hat{F}^b$ signifies the empirical distribution created through bootstrap resampling from $\hat{F}$. In particular, superscripts are used to distinguish bootstrapped distributions and subscripts to denote distinct input distributions; for instance, $\hat{F}^b_i$ denotes the empirical distribution constructed from the $b$-th bootstrap resample of the $i$-th empirical distribution $\hat{F}_i$. For a random variable $X_{i,j}$ with double subscripts, the first subscript denotes distinct input distributions, and the second subscript denotes diverse samples drawn from the $i$-th input distribution. We use $\mathbb{E}_*$ and $\text{Var}_*$ to denote the expectation and variance over the bootstrap resamples from the data, conditional on the original data $\hat{F}$. That is, $\mathbb{E}_*$ and $\text{Var}_*$ only accounts for the randomness in simulation, i.e. the simulation error. We use $\mathbb{E}$ and $\text{Var}$ to denote the expectation and variance over the data distribution $F$. We also use the standard $O_p$ notations: $\Theta_p(\cdot)$, $O_p(\cdot)$, $\Omega_p(\cdot)$ in probability statement about the data distribution $F$, to only hide constants that do not depend on $n$.

Debiasing via Orthogonal Bootstrap

In this section, we study the problem of estimating function/functional values when uncertainty exists on the input value quenouille1949approximate,efron1982jackknife,efron1992bootstrap,jiao2020bias,koltchinskii2021estimation,koltchinskii2022estimation,zhou2021high,etter2020operator,etter2021operator,ma2022correcting. We first formulate the problem of simulating with input uncertainty song2014advanced,lam2018subsampling, then describe the motivation of our Orthogonal Bootstrap method.

Suppose one aim to estimate a generic performance measure $\phi(F_1,\cdots,F_m)$ depend on $m$ independent input distributions $F_1,\cdots,F_m$. We consider the setting when the input distribution is unknown and only $n_i$ i.i.d. generated data $\{X_{i,1},\cdots, X_{i,n_i}\}$ from distribution $F_i$ is available, forming the empirical distributions $\hat F_i:=\sum_{j=1}^{n_i}\delta_{X_{i,j}}(x)/n_i$. The Bootstrap efron1994introduction methods simulates $B\in\mathbb{Z}^+$ Monte Carlo replications. In each Monte Carlo replication, resampling with replacement is performed independently and uniformly from $\{X_{i,1},\cdots,X_{i,n_i}\}$, repeated $n_i$ times. This process yields sets $\{X_{i,1}^b,\cdots,X_{i,n_i}^b\} $ for each replications $b=1,\cdots,B$. Then according to the bootstrap principle, we can estimate $\mathbb{E}\phi(\hat{F}_1,\cdots,\hat{F}_m)$ by the bootstrap mean $$ \mathbb{E}_* \phi(\hat F_{1}^{b},\cdots,\hat F_{m}^{b})=\lim_{B\rightarrow \infty}\frac{1}{B}\sum_{i=1}^B \phi(\hat F_{1}^{b},\cdots,\hat F_{m}^{b}), $$where we recall from the notation section that $\mathbb E_*$ is the expectation with respect to resampling conditional on the original data and $\hat F_{i}^{b}=\sum_{j=1}^{n_i}\delta_{X_{i,j}^b}(x)/n_i$ is the resampled empirical distribution.

The simulated bootstrap mean can be used for debiasing the naive estimation of function/functional values when uncertainty exists on the input value quenouille1949approximate,efron1982jackknife,efron1992bootstrap,jiao2020bias,koltchinskii2021estimation,koltchinskii2022estimation,zhou2021high,etter2020operator,etter2021operator,ma2022correcting. The debiasing procedure appears in many applications, including inverse a noisy elliptic system etter2020operator,etter2021operator, optimal stopping zhou2021unbiased, online learning chen2022debiasing and stochastic optimization ma2022correcting,li2020debiasing. ma2022correcting showed that the estimation can be improved from the naive estimator with $\Omega(n)$ times of bootstrap simulation. Bootstrap function/functional estimates in real applications is computationally expensive, as it requires retraining a machine learning model. Therefore, bootstrapping a large number of times (proportional to the number of data points) is impractical, and Bootstrap a limited number of times can lead to high variance and bad performance lam2018subsampling.

In this paper, we propose methods to reduce resampling effort while maintaining the expected accuracy. To achieve this, suppose $\nabla \phi$ is the von Mises derivative (defined in Section (ref) in the appendix) of the statistical functional $\phi$. Then we decompose the bootstrap simulation target into the non-orthogonal part and the orthogonal part

equation*[equation* omitted — 656 chars of source]

\newline To compute $\nabla \phi(\hat F_1,\cdots,\hat F_m)[\hat F_{1}^{b}-\hat F_1,\cdots,\hat F_{m}^{b}-\hat F_m]$, we use the influence function technique cook1980characterizations,koh2017understanding,giordano2019higher. Noting that typically $\nabla \phi(\hat F_1,\cdots,\hat F_m)[\tilde{F}_{1}-\hat F_{1},\cdots,\tilde{F}_{m}-\hat F_m]$ is a linear functional in $(\tilde{F}_{1},\cdots,\tilde{F}_{m})$, and the Riesz representation theorem implies the existence of a mean zero, finite variance functions $(\mathcal{I}_1(X_1),\cdots, \mathcal{I}_m(X_m))$ such that $\nabla \phi(\hat F_1,\cdots,\hat F_m)[\tilde{F}_{1}-\hat F_{1},\cdots,\tilde{F}_{m}-\hat F_m]=\sum_{i=1}^m \mathbb{E}_{X_i\sim \tilde{F}_i}\mathcal{I}_i(X_i)$ cook1980characterizations,serfling2009approximation,van2000asymptotic. As the influence function (the non-orthogonal part) is mean zero, i.e. $\mathbb{E}_* \nabla \phi(\hat F_1,\cdots,\hat F_m)[\hat F_{1}^{b}-\hat F_1,\cdots,\hat F_{m}^{b}-\hat F_m]=0$, it doesn't need to be simulated. We can also understand the non-orthogonal part as a control variate asmussen2007stochastic to reduce the simulation error of bootstrap simulation. The whole idea is summarized in Figure (ref).

Applying the above idea to bias correction, we can simulate the bias $\phi(F_1,\cdots,F_m)-\mathbb{E}\phi(\hat{F}_1,\cdots,\hat{F}_m)$ by \[\phi(\hat{F}_1,\cdots,\hat{F}_m)-\frac{1}{B}\sum_{b=1}^B(\hat{\phi}^b-\sum_{i=1}^m\frac{1}{n_i}\sum_{j=1}^{n_i}\mathcal{I}_i^\phi(x_{i,j}^b))\]

The whole procedure is summarized in Algorithm (ref).

algorithm[algorithm omitted — 1,015 chars of source]

Provable Improvement of Orthogonal Bootstrap

In this section, we show our Orthogonal Bootstrap method can provably reduce the required Monte Carlo replications under mild assumption. For simplicity, we assume $n_1=n_2=\cdots=n_m=n$ in our theoretical investigation. We assume our simulation functional has a continuous Fr\'echet gradient in tht under the Kernel Maximum Mean Discrepancy (MMD) distance muandet2017kernel. For background in kernel mean embeddings, we refer to Section (ref) in the appendix. We first rewrite the simulation functional in terms of the kernel mean embeddings, i.e.

equation[equation omitted — 101 chars of source]

where $\mu_i$ are kernel mean embeddings using kernel $k_i$ muandet2017kernel, $\mathcal{H}:=\mathcal{H}_1\times\cdots\times\mathcal{H}_m$ where $\mathcal{H}_i$ is the reproducing kernel Hilbert space respect to the kernel $k_i$ and $h:\mathcal{H}\to\mathbb{R}$ is a functional on $\mathcal{H}$. Denote $\mu:=\mu_1\times\cdots\mu_m$ and $F:=F_1\times\cdots F_m$ for simplicity.

restatable{assumption}{mainassumptionone} There exists kernel mean embeddings $\mu_i:\mathcal{F}_i\to\mathcal{H}_i$ which maps $F_i$ into $\mathcal{H}_i$ for $i=1,\cdots,m$. Moreover, for all $i=1,\cdots,m$, $\mathbb{E}k_i(X_i,X_i)^4<\infty$ and $\mathbb{E}k_i(X_i,Y_i)^4<\infty$ where $X_i,Y_i$ are independent samples from $F_i$.

For a functional on $\mathcal{H}$, we say that it is of class $C^1$ if its Fr\'echet derivative exists and is continuous.

restatable{assumption}{mainassumptiontwo} The non-constant functional $h:\mathcal{H}_1\times\cdots\times\mathcal{H}_m\to\mathbb{R}$ is of class $C^1$ and its derivative is Lipschitz in the sense that $|Dh(x_1)(v)-Dh(x_2)(v)|\le L\|x_1-x_2\|_{\mathcal{H}}\|v\|_{\mathcal{H}}$. Moreover, $\|\partial_ih(\mu(F))^4\|_{\mathcal{H}_i}<\infty$.

Under the above two assumption on the performance measure (ref), we have

theoremLet $X_{ob}$ be the Orthogonal Bootstrap estimator defined in Equation (ref) and $X_{sb}$ be the Standard Bootstrap estimator defined by $X_{sb}:=2\hat \phi-\frac{1}{B}\sum_{b=1}^B \hat\phi^b$. Under Assumption (ref) and Assumption (ref), if the number of Monte Carlo replications $B\geq Cn^\alpha$ for some absolute constant $C>0$ and $\alpha\geq 0$, then the simulation error for the Orthogonal Bootstrap estimator satisfies $\text{Var}_*(X_{ob})=O_p(\frac{1}{n^{2+\alpha}})$ and the simulation error for the Standard Bootstrap estimator satisfies $\text{Var}_*(X_{sb})=\Theta_p(\frac{1}{n^{1+\alpha}})$.

The theorem is proved by combining Theorem (ref) and Theorem (ref) with Theorem (ref) in the appendix. Here we provide an informal proof to illustrate our idea.

iproofThe complexity of bootstrap simulation is related to the variance of the following random variable {\tiny\begin{equation*} \xi=\phi(\hat F_{1}^{b},\cdots,\hat F_{m}^{b})-{\nabla \phi(\hat F_1,\cdots,\hat F_m)\left[\hat F_{1}^{b}-\hat F_1,\cdots,\hat F_{m}^{b}-\hat F_m\right]}. \end{equation*}} The random variable is at the scale of $O_p(\hat F_{1}^{b}-\hat F_1,\cdots,\hat F_{m}^{b}-\hat F_m)^2$, i.e. $O_p(n^{-1})$. Thus its variance is at scale $O_p(n^{-2})$. Thus to achieve simulation error of order $O_p(n^{-\alpha})$, one needs $O(n^{\alpha-2})$ Monte Carlo replications.
remarkTo debias a functional estimation, one typically needs the simulation error to be $O_p(n^{-2})$ (see Theorem 2 and its proof in ma2022correcting). According to our theorem, Standard Bootstrap needs $\Omega(n)$ Monte Carlo replications but Orthogonal Bootstrap only needs $O(1)$ Monte Carlo replications.

Variance Estimation via Orthogonal Bootstrap

In this section, we discuss how the idea of Orthogonal Bootstrap can be used to accelerate bootstrap simulation for estimating the variance, which provides a versatile nonparametric method for constructing confidence intervals and prediction intervals without detailed model knowledge. Following the same setting of the previous section, suppose one aims to estimate a generic performance measure $\phi(F_1,\cdots,F_m)$. To quantify the uncertainty of the plug-in estimator $\phi(\hat F_1,\cdots,\hat F_m)$, one needs to know its variance, which can be simulated by $\text{Var}_*(\phi(\hat F_{1}^{b},\cdots,\hat F_{m}^{b}))$ according to the bootstrap principle (recall that $\text{Var}$ without subscript denotes the simulation variance conditioned on the original data). Similar to the previous section, we decompose the variance into the variance of the non-orthogonal part, the variance of the orthogonal part, and their cross-covariance as in Equation (ref) as follows.

equation[equation omitted — 789 chars of source]

The variance of the non-orthogonal part enjoys a closed-form representation using the influence function as follows.

lemma{\scriptsize $\text{Var}_* \left(\nabla \phi(\hat F_1,\cdots,\hat F_m)\left[\hat F_{1}^{b}-\hat F_1,\cdots,\hat F_{m}^{b}-\hat F_m\right]\right)\\= \sum_{i=1}^m\frac{1}{n_i^2}\sum_{j=1}^{n_i}(\mathcal{I}_i^\phi (X_{i,j}))^2$}

For the proof, see Lemma (ref) in the appendix. Based on this observation and the orthogonal and non-orthogonal decomposition, we propose Algorithm (ref). We calculate the non-orthogonal part (influence function) using the closed-form representation and simulate the remainder part using the Monte-Carlo method.

Under the same assumption as debiasing, Orthogonal Bootstrap provably improves the required Monte Carlo replications when simulating the variance.

theoremFor simplicity, we assume $n_1=n_2=\cdots=n_m=n$. Let $X_{ob}$ be the Orthogonal Bootstrap estimator defined in Equation (ref) and $X_{sb}$ be the Standard Bootstrap estimator defined by $X_{sb}:=\frac{1}{B}\sum_{b=1}^{B}(\hat{\phi}^b-\overline{\phi})^2$. Under Assumption (ref) and Assumption (ref), if the number of Monte Carlo replications $B\geq Cn^\alpha$ for some absolute constant $C>0$ and $\alpha\geq 0$, then the simulation error for the Orthogonal Bootstrap estimator satisfies $\text{Var}_*(X_{ob})=O_p(\frac{1}{n^{3+\alpha}})$ and the simulation error for the Standard Bootstrap estimator satisfies $\text{Var}_*(X_{sb})=\Theta_p(\frac{1}{n^{2+\alpha}})$.
remarkTo achieve relative consistency in variance estimation, one typically needs the simulation error to be $O_p(n^{-3})$ (see Theorem 1 in lam2022subsampling). According to Theorem (ref), the Standard Bootstrap needs $\Omega(n)$ Monte Carlo replications but the Orthogonal Bootstrap only needs $O(1)$ Monte Carlo replications.

The theorem is proved by combining Theorem (ref) and Theorem (ref) with Theorem (ref) in the appendix.

When the sample size $n$ is small, our Orthogonal Bootstrap estimator for variance can obtain negative values sometimes. An improved Orthogonal Bootstrap estimator (Algorithm (ref) in the appendix) should be used in this case. We note here that the improved Orthogonal Bootstrap estimator enjoys the same theoretical properties (Theorem (ref) in the appendix) as the original Orthogonal Bootstrap estimator. See Section (ref) in the appendix for details.

algorithm*\caption{Variance Estimation via Orthogonal Bootstrap} Input: A generic performance measure $\phi(F_1,\cdots,F_m)$, i.i.d samples $\{X_{i,1},\cdots,X_{i,n_i}\}\in\mathbb{R}^{d_1}$ of $F_i$, and influence function $\mathcal{I}_i^\phi$ of $\phi$ respect to $\hat{F}_i$. \\ Output: Estimation of {$\text{Var}\phi(\hat F_1,\cdots,\hat F_m)$}. \begin{algorithmic} \STATE $\hat \phi\leftarrow \phi(\hat F_1,\cdots,\hat F_m)$, where $\hat F_i=\frac{1}{n_i}\sum_{j=1}^{n_i}\delta_{X_{i,j}}$ \FOR{b=1:B} \FOR{i=1:m} \STATE Sample $\{x_{i,1}^b,\cdots,x_{i,n_i}^b\}$ i.i.d from $\hat{F_i}$ \ENDFOR \STATE $\hat\phi^b\leftarrow\phi(\hat F_1^b,\cdots,\hat F_m^b)$, where $\hat{F}_i^b=\frac{1}{n_i}\sum_{j=1}^{n_i} \delta_{x_{i,j}^b}$ \STATE Calculate $\mathcal{\hat I}^b=\sum_{i=1}^m\frac{1}{n_i}\sum_{j=1}^{n_i}\mathcal{I}_i^\phi(\tilde{x}_{i,j})$ \ENDFOR \STATE Calculate $\overline{\phi-\mathcal{I}}\leftarrow\frac{1}{B}\sum_{b=1}^{B}(\hat\phi^b-\mathcal{\hat I}^b)$, $\overline{\mathcal{I}}\leftarrow\frac{1}{B}\sum_{b=1}^{B}\mathcal{\hat I}^b$. \STATE Estimate $\text{Var}_{F_1,\cdots,F_m}\phi(\hat F_1,\cdots,\hat F_m)$ by $S^2$, where \begin{equation} \begin{aligned} S^2&=\sum_{i=1}^m\frac{1}{n_i^2}\sum_{j=1}^{n_i}(\mathcal{I}_i^\phi (X_{i,j}))^2+\frac{1}{B}\sum_{b=1}^{B}\left(\hat\phi^b-\mathcal{\hat I}^b-\overline{\phi-\mathcal{I}}\right)^2+\frac{2}{B}\sum_{b=1}^{B}\left(\hat\phi^b-\mathcal{\hat I}^b-\overline{\phi-\mathcal{I}}\right)\left(\mathcal{\hat I}^b-\overline{\mathcal{I}}\right). \end{aligned} \end{equation} \end{algorithmic}
figure*[figure* omitted — 501 chars of source]
table*[table* omitted — 5,551 chars of source]

Numerical Examples

In this section, we test the numerical performances of our Orthogonal Bootstrap and compare it with the Standard Bootstrap. For confidence and prediction interval examples, we also compared our method with the recently proposed Cheap Bootstrap method lam2022cheap. We demonstrated that our Orthogonal Bootstrap achieved significantly smaller bias and higher empirical coverage probability over the original Bootstrap method when the number of Monte Carlo replications is small. Although Cheap Bootstrap lam2022cheap can also provide comparable empirical coverage probability to our method, the mean width of the interval constructed by Cheap Bootstrap is much longer. Our Orthogonal Bootstrap method provides interval with higher empirical coverage probability but {achieves the same expected width as the original Bootstrap.} The experiment details are left in Section (ref) in the appendix.

Debiasing

We consider four numerical examples following ma2022correcting. For all of the examples, we compute the the average bias (BIAS) of the estimates across 1000 experiments. We run the naive estimator without debiasing, Standard Bootstrap, and our Orthogonal Bootstrap for a small number of Monte Carlo replications $B=2,3,4,5,6,7,8,9,10$. To ensure robustness and reliability, we repeat the bootstrap procedure ten times and provide insights through the reporting of quantiles at the 5th, 50th, and 95th percentiles. The results is shown in Figure (ref).

\paragraph{Function of Mean} We simulate the function-of-mean model, namely estimating $\phi=g(\boldsymbol{\mu})$ where $\boldsymbol{\mu}=\mathbb{E}\mathbf{X}$ for a $d$-dimensional random vector $\mathbf{X}$ and $g:\mathbb{R}^d\rightarrow\mathbb{R}$ is a function. We have i.i.d random sample $\{\mathbf{X}_1,\cdots,\mathbf{X}_n\}$ of random vector $\mathbf{X}$. We consider an ellipsoidal estimation problem $g(\boldsymbol{\mu})=\|\boldsymbol{\mu}\|_2^2$ and a fourth-order polynomial estimation problem $g(\boldsymbol{\mu})=\|\boldsymbol{\mu}\|_2^4$. The underlying distribution is set to be $\mathcal{N}(0.2\mathbf{1_d},I_{d})$. We use $d=25$ and a sample size $n=100$. The influence function of $g(\boldsymbol{\mu})=\|\boldsymbol{\mu}\|^2$ is $\mathcal{I}^g(\mathbf{x})=2(\mathbf{x}^T\hat{\boldsymbol{\mu}}-\|\hat{\boldsymbol{\mu}}\|^2)$, and the influence of $g(\boldsymbol{\mu})=\|\boldsymbol{\mu}\|^4$ follows by the chain rule.

\paragraph{Entropy} We consider estimating entropy for discrete probability distributions, namely estimating $\phi=-\sum_{i=1}^d p_i\ln p_i$ where $p=(p_1,\cdots,p_d)$ satisfies $p_i>0$ and $\sum_{i=1}^d p_i =1$. We generate the groundtruth distribution $p$ from symmetric Dirichlet distribution with parameter $\alpha_i=1$ for all $i=1,\cdots,d$. Noisy observations of $p$ are the empirical distributions of single samples from $p$. We use a sample size $n=1000$ and dimension $d=100$. The influence function is $\mathcal{I}^\phi(x)= \sum_{i=1}^d (x_i-\hat{p}_i)(\log \hat{p}_i +1)$.

\paragraph{Constrained Optimization Problem} The debiasing method can also be applied to optimization problems with randomness. Here we consider estimating $\phi = \arg\min_{\mathbf{x}} \mathbf{x}^T B \mathbf{x}$ under the constrain $A\mathbf{x}=\mathbf{b}$, where $\mathbf{x}\in\mathbb{R}^p$, $B\in\mathbb{R}^{p\times p}$ is a positive definite matrix, $A \in \mathbb{R}^{d\times p}$, and $\mathbf{b} \in \mathbb{R}^d$ while we only have access to its noisy observations, i.e. $\mathbf{b}\sim \mathcal{N}(\mathbf{b}^*,I_{d\times d})$. We sample $\mathbf{b}^*$ uniformly from the unit sphere $\mathbb{S}^{d-1}$. We use a sample size $n=10$, dimension $d=100$ and $p=200$.

To compute the influence function, we do not solve the quadratic programming problem directly and differentiate. Instead, we utilize a general procedure that can be applied to optimization problems lacking an explicit expression for their solutions by differentiating the KKT conditions of the constrained optimization problem. The details can be found Section (ref) in the appendix.

figure[figure omitted — 752 chars of source]
table*[table* omitted — 1,365 chars of source]

Confidence Interval Construction

In this section, we aim to construct the confidence interval hall1986bootstrap,hall1988bootstrap,efron1992bootstrap with minimal resampling effort using our Orthogonal Bootstrap technique. The confidence interval constructed by the Standard Bootstrap method is

equation*[equation* omitted — 268 chars of source]

where $z_{1-\alpha/2}$ being the $(1-\alpha/2)-$quantile of the standard normal, $\hat\phi=\phi(\hat F_1,\cdots,\hat F_m)$ is the plug-in estimator of $\phi(F_1,\cdots,F_m)$ and $\hat F_i=\frac{1}{n}\sum_{j=1}^n \delta_{X_{i,j}}$. We use Orthogonal Bootstrap to estimate $\text{Var}_*(\phi(\hat F_{1}^{b},\cdots,\hat F_{m}^{b}))$. The details of the algorithm for confidence interval construction is provided in Algorithm (ref) in the appendix.

\paragraph{Elementary Examples} Following lam2022cheap, we consider a folded standard normal (i.e.,$|N(0,1)|$) and double exponential with rate 1 (i.e. $\text{Sgn}\times \text{Exp}(1)$, where $\text{Sgn}=+1$ or $-1$ with equal probability and is independent with $\text{Exp}(1)$). We aim to provide confidence intervals for their variance. The two setups have explicit ground truth, namely $1-2/\pi$ and $2$, respectively. The influence function of variance is $\mathcal{I}(x)=(x-\hat \mu)^2-\hat \sigma^2$, where $\hat \mu$ is the empirical mean and $\hat\sigma^2$ is the empirical variance. The third example, also follows lam2022cheap, is estimating the correlation of bivariate lognormal (i.e.$(e^{Z_1},e^{Z_2})$, where $(Z_1,Z_2)$ is the bivariate normal with mean zero, unit variance and correlation 0.5). The ground truth of this example is also known as $(e^{3/2}-e)/(e^2-e)$. The influence function of correlation can be determined by the influence of covariance and the chain rules. We use a sample size $n=1000$ for all three examples. We run 1000 independent simulations and report the empirical coverage and the mean width of the confidence interval. We run the Standard Bootstrap, Cheap Bootstrap lam2022cheap and Orthogonal Bootstrap using a small number of Monte Carlo replications $B=2,5,10$. We set the confidence level $(1-\alpha)$ as $95\%$. For each setting, we repeat our experiments 1000 times and report the empirical coverage, interval width mean and standard deviation. The result is shown in Table (ref).

\paragraph{Regression Problem} We apply our method to a linear regression problem in sengupta2016subsampled,lam2022cheap. We aim to fit a model $Y=\beta_1X_1+\cdots+\beta_dX_d+\epsilon$ where we set dimension $d=50$ and use data $\{(x_{1,i},\cdots,x_{d,i},Y_i)\}_{i=1}^n$ of size $n=5000$ to fit the model. We set $\log(X_i)\sim N(0,1)$ and $\epsilon\sim 10*N(0,1)$ as the data generating process. The influence function is easy to calculate in this case as the estimator is a $M$-estimator. For example, cook1980characterizations,koh2017understanding calculate the influence function for M-estimator. Specifically, consider the influence function of the M-estimator $\hat \theta = \arg\min_\theta \frac{1}{n}\sum_{i=1}^n L(z_i,\theta)$ , where $z_i$ represents training data and labels, $L(z,\theta)$ is the loss function for data $z$ and parameter $\theta$. The influence function for the parameter $\theta$ at point $z$ is

equation[equation omitted — 263 chars of source]

where $H_{\hat \theta}=\frac{1}{n}\sum_{i=1}^n \nabla_\theta^2 L(z_i,\hat\theta)$ is the empirical hessian of the loss function and is positive definite by assumption. As shown in koh2017understanding,alaa2020discriminative, the influence function can be solved efficiently via inverse hessian vector product combined with the modern auto-grad systems like TensorFlow abadi2016tensorflow and Pytorch paszke2019pytorch. The calculation of the influence function is much faster than retraining the model. In this example, we run the Standard Bootstrap, Cheap Bootstrap lam2022cheap, and our Orthogonal Bootstrap using a small number of Monte Carlo replications $B=2,5,10$ and report the empirical coverage, interval width mean, and standard deviation for the first coefficient in Table (ref).

Prediction Interval Construction for Real Data

In this section, we use Orthogonal Bootstrap to accelerate the pivot Bootstrap method contarino2022constructing to construct prediction intervals for neural networks. Following alaa2020discriminative, we conduct our experiments on 3 UCI benchmark datasets for regression: yacht hydrodynamics, energy efficiency Dua:2019 and kin8nm. We partition these datasets into distinct training and testing subsets. For the training datasets, we employ two-layer neural networks with a hidden dimension of 100 and a hyperbolic tangent (tanh) activation function. Subsequently, we apply Standard Bootstrap, Cheap Bootstrap, and our proposed Orthogonal Bootstrap to simulate the variance of the resampling step in the pivot Bootstrap method contarino2022constructing. The construction details for prediction intervals can be found in Algorithm (ref) in the appendix. and the training procedures specific to each dataset We set the target coverage to be $(1-\alpha)=0.95$. We report the empirical coverage and interval width mean in Table (ref), where the empirical coverage is the percentage of test data points which fall into the corresponding prediction interval. Our reported results are averaged over three random seeds to ensure robustness and reliability.

Conclusion and Discussion

In summary, our paper introduces the concept of Orthogonal Bootstrap, a novel technique that streamlines the process of bootstrap resampling. By separately treating the non-orthogonal (influence function) and orthogonal parts, akin to utilizing Infinitesimal Jackknife estimator as a control variate, we effectively reduce the number of required Monte Carlo replications. This innovation allows our method to maintain the higher-order coverage properties associated with traditional Bootstrap methods while simultaneously decreasing the computational burden. It's important to note that the control variate significantly reduces the computational costs involved in the simulation process while it doesn't alter the expected length of the constructed confidence interval. In essence, Orthogonal Bootstrap presents a practical and efficient solution for improving the accuracy and speed of Bootstrap resampling, making it a valuable tool for statistical analysis and inference.

Impact Session

Our paper contributes to fast and accurate quantification of the uncertainty for larger scale machine learning problems, which giving policy maker an idea of reliability of the AI prediction and is vital for risk management and help to be responsible in contexts where AI decisions have significant ethical implications. The theoretical results presented in the paper will have no ethical impact.