Year: 2021

Abstract:
Bootstrap is a principled and powerful frequentist statistical tool for uncertainty quantification. Unfortunately, standard bootstrap methods are computationally intensive due to the need of drawing a large i.i.d. bootstrap sample to approximate the ideal bootstrap distribution; this largely hinders their application in large-scale machine learning, especially deep learning problems. In this work, we propose an efficient method to explicitly optimize a small set of high quality “centroid” points to better approximate the ideal bootstrap distribution. We achieve this by minimizing a simple objective function that is asymptotically equivalent to the Wasserstein distance to the ideal bootstrap distribution. This allows us to provide an accurate estimation of uncertainty with a small number of bootstrap centroids, outperforming the naive i.i.d. sampling approach. Empirically, we show that our method can boost the performance of bootstrap in a variety of applications.

text: Bootstrap is a principled and powerful frequentist statistical tool for uncertainty quantification. Unfortunately, standard bootstrap methods are computationally intensive due to the need of drawing a large i.i.d. bootstrap sample to approximate the ideal bootstrap distribution; this largely hinders their application in large-scale machine learning, especially deep learning problems. In this work, we propose an efficient method to explicitly optimize a small set of high quality “centroid” points to better approximate the ideal bootstrap distribution. We achieve this by minimizing a simple objective function that is asymptotically equivalent to the Wasserstein distance to the ideal bootstrap distribution. This allows us to provide an accurate estimation of uncertainty with a small number of bootstrap centroids, outperforming the naive i.i.d. sampling approach. Empirically, we show that our method can boost the performance of bootstrap in a variety of applications.

heading: 1 INTRODUCTION
text: Bootstrap is a simple and principled frequentist uncertainty quantification tool and can be flexibly applied to obtain data uncertainty estimation with strong theoretical guarantees (Hall et al., 1988; Austern & Syrgkanis, 2020; Chatterjee et al., 2005; Cheng et al., 2010). In particular, when combined with the maximum likelihood estimator or more general M-estimators, bootstrap provides a generalpurpose, plug-and-play non-parametric inference framework for general probabilistic models without case-by-case derivations; this makes it a promising frequentist alternative to Bayesian inference.
However, the standard bootstrap inference is highly expensive in both computation and memory as it typically requires drawing a large number1 of i.i.d. bootstrap particles (samples) to obtain an accurate uncertainty estimation. In the context of this paper, as each bootstrap particle/sample/centroid is a machine learning model, we might directly call a model as particle/sample/centroid. With a small number of particles, bootstrap may perform poorly. As a consequence, when applied to deep learning, we need to store a large number of neural networks and feed the input into a tremendous number of networks every time we make inference, which can be quite expensive and even unaffordable for deep learning problems with huge models2. For example, in autonomous driving applications, our device can only store a limited number of models and we need to make decisions within a short time, which makes the standard bootstrap with a large number of models no more feasible. Typical ensemble methods in deep learning, such as Lakshminarayanan et al. (2017); Huang et al. (2017); Vyas et al. (2018); Maddox et al. (2019); Liu & Wang (2016), can only afford to use a small number (e.g., less than 20) of models.
Therefore, to make bootstrap more accessible in modern machine learning, it is essential to develop new approaches that break the key computation and memory barriers mentioned above. This paper aims to improve the bootstrap when the resource at inference is limited. We are motivated to consider the following problem:
How to improve the accuracy of bootstrap inference when the number of particles is limited?
We attack this challenge by presenting an efficient centroid approximation for bootstrap. Our method replaces the i.i.d. bootstrap particles with a set of carefully optimized centroid particles that are
1For example, thousands of, as suggested by Statistics textbooks such as Wasserman (2013). 2While training cost is an extra burden, it is small compared with the cost of making prediction as we only
need to train the model once but make countless predictions at deployment.
guaranteed to provide an accurate and compact approximation to the ideal bootstrap distribution so that only a smaller number of particles is needed to obtain good performance.
Our method is based on minimizing a specially designed objective function that is asymptotically equivalent to the Wasserstein distance between the ideal bootstrap distribution and the particle distribution formed by the learned centroids. During the training, each centroid adjusts its location being aware of the locations of the others so that centroids are diversified and well distributed on the domain. Our method is similar to doing K-means on the ideal bootstrap distribution, finding K representative centroids that well represent K separate parts of the target distribution’s domain in an optimal way. As centroids are optimized to better approximate the distribution, our approach naturally improves over the vanilla bootstrap with i.i.d. particles. See Figure 1 for illustration.
Empirically, we apply the centroid approximation method to various applications, including confidence interval estimation (DiCiccio et al., 1996), bootstrap method for contextual bandit (Riquelme et al., 2018), bootstrap deep Q-network (Osband et al., 2016) and bagging method (Breiman, 1996) for neural networks. We find that our method consistently improves over the standard bootstrap.

heading: 2 BACKGROUND
text: Suppose we have a model fθ parameterized by θ in a parameter space Θ ⊆ Rd. Let {xi}ni=1 ⊂ X be a training set with n data points on X . Assume `(x, fθ) is the negative log-likelihood of data point x with model fθ. A standard approach to estimate θ is maximum likelihood estimator (MLE), which minimizes the negative log-likelihood function (loss) over the training set
θ̂ = arg min θ∈Θ L(θ), L(θ) =
∑n i=1`(xi, fθ)/n.
Here the MLE θ̂ provides a point estimation without any information on the data uncertainty. Bootstrap is a simple and effective frequentist method to quantify the uncertainty. The bootstrap loss is a randomly perturbed loss defined as
Lw(θ) = ∑n i=1wi`(xi, fθ)/n,
where w = [w1, ..., wn]> is a set of random weights of data points drawn from some distribution π. A typical choice of π is the multinomial distribution with uniform probability, which corresponds to resampling on the training set with replacement. Given w, one can calculate its associated bootstrap particle by minimizing the bootstrap loss:
θ̂w = arg min θ∈Θ Lw(θ). (1)
Let ρπ be the distribution of θ̂w when w ∼ π. Bootstrap theory indicates that we can quantify the data uncertainty of θ or any function g(θ) using ρπ . We call ρπ the ideal bootstrap distribution and it is the main object we want to approximate.
Denote δθ as the delta measure centered at θ. Standard bootstrap method approximates ρπ by the particle distribution ρ̂π(·) = ∑m j=1δθ̂wj (·)/m formed by m i.i.d. particles {θ̂wj}mj=1, which can be obtained by drawing m i.i.d. weights {wj}mj=1 from π and calculating each θ̂wj based on (1). However, for deep learning applications, as discussed in the introduction, storing and making inference using a large number m of bootstrap particles can be quite expensive. On the other hand, if m is small, ρ̂π tends to be a poor approximation of ρπ. In this paper, we aim to improve the approximation of the particle distribution when m is small.

heading: 3 METHOD
text: Our idea is simple. Instead of using i.i.d. particles, in which the location of each particle is independent from that of the others, we try to actively optimize the location of each particle so that particles are diversified, better distributed and eventually providing a particle distribution with improved approximation accuracy. A natural way to achieve this goal is to explicitly optimize a set of points {θj}mj=1 (called centroids) jointly such that the Wasserstein distance between ρπ and the induced particle distribution is minimized:
{θ∗j , v∗j }mj=1 = arg min θ1,...,θm∈Θ, [v1,...,vm]∈Cm
W2 [∑m j=1vjδθj , ρπ ] . (2)
Here we consider a Wasserstein distanceW2 equipped with a special data-dependent distance metric || · ||D that will be introduced later in (5). We note that here we also optimize the probability weights {vj}mj=1 of the centroids. Finding the optimal centroids and probability weights can be decomposed into two steps: the centroid learning phase and the probability weights learning phase, based on the facts in (3,4).
W22 [∑ j∈[m]v ∗ j δθ∗j , ρπ ] = Jπ({θ∗j }mj=1), where Jπ({θj}mj=1):=Ew∼π [ minj∈[m]||θj − θ̂w||2D ] . (3)
Here (3) implies that, to find the optimal particle distribution in (2), we can start with the centroid learning phase where we only need to optimize the centroids. It can be achieved by minimizing Jπ({θj}mj=1), which is the averaged distance of bootstrap particles to their closest centroid. After we obtain the optimal centroids, the optimal probability weights can be learned by (4):
v∗j = ṽ ∗ j / ∑ s∈[m] ṽ ∗ s , where ṽ ∗ j = Pw∼π ( j = arg minj∈[m]||θ∗j − θ̂w||2D ) . (4)
Here v∗j is the proportion of bootstrap particles that are closest to the centroid j. We emphasize that the optimal solution to two-stage learning is guaranteed to be the global minimizer of the loss in (2) (see Lemma 3.1 and 3.2 in Canas & Rosasco (2012)).
However, the key issue is that the losses in both (2, 3) can not be computed in practice, as they require us to access ρπ (i.e., obtain θ̂w first in order to calculate the loss). To handle this issue, we seek an easy-to-compute surrogate loss. Our idea is based on the following observation. Assuming the size of training data is large, which is usually the case in deep learning, we can expect that θw will be centered around a small region3. It implies that we should search the centroid in this small region. Notice that when θ is close to θ̂w, based on Taylor approximation, we have
Lw(θ) ≈ Lw(θ̂w) +∇>θ Lw(θ̂w)(θ − θ̂w) + 1/2(θ − θ̂w)>∇2θLw(θ̂w)(θ − θ̂w) ≈ L∞(θ0) + 1/2||θ − θ̂w||2D, where ‖V ‖ 2 D := V
>∇2θL∞(θ0)V. (5) Here L∞(θ) := Ex`(x, fθ) denotes the population loss; θ0 is the minimizer of L∞(θ). In (5), we use the facts4 that ∇>θ Lw(θ̂w) = 0; and with large training set, the empirical distribution ∑n i=1δxi/n
well approximates the whole data population, and hence the bootstrap resampling distribution, i.e.,∑n i=1wiδxi/n on the empirical distribution also well approximates the whole data population. This implies that Lw(·) ≈ L∞(·) and ∇2θLw(·) ≈ ∇2θL∞(·). As the loss are close to each other, their minimizers are also close θ̂w ≈ θ0. Since L∞(θ0) is some (unknown) constant independent with θ, we can replace the ||θj − θ̂w||2D in (3) by Lw(θj) as it only adds some constant into the loss.
3This can be formally characterized by central limit theorem as discussed in Section 4. 4We defer the detailed analysis to Section 4.
Intuitively, we can expect that the centroid closest to θ̂w is the one that gives the smallest loss on Lw. It motivates us to learn the centroids via the modified centroid learning phase:
{θ∗j }mj=1 = arg min θ1,...,θm∈Θ
Ew∼π [ minj∈[m]Lw(θj) ] . (6)
Similarly, the optimal probability weights can be learned via the modified weight learning phase: v∗j = ṽ ∗ j / ∑ s∈[m] ṽ ∗ s , where ṽ ∗ j = Pw∼π ( j ∈ arg minj∈[m]Lw(θ∗j ) ) . (7)
We note that here we slightly abuse the notation of θ∗j and v ∗ j in (3,4) and (6,7) for simplification. In the later context, θ∗j and v ∗ j are used based on their definitions in (6,7).
Connection to K-means By viewing the target distribution as a set of particles that we want to cluster, in K-means clustering, each centroid (i.e., K-means center) represents one of the K disjoint groups5 of particles, which is formed by assigning each particle in the whole set to the closest centroid among all the K centroids. K-means learns the optimal K centroids in the way that they can best approximate the whole set. The ‘closeness’ for assigning the particles is measured by the distance between the two points. As pointed out by Canas & Rosasco (2012), K-means essentially searches the optimal particle distribution formed by the K centroids that minimizes its Wasserstein distance to the target distribution. Our centroid approximation idea follows the same fashion of clustering but our key innovation is to measure the ‘closeness’ by examining the bootstrap loss of the centroids so that we can still learn the optimal centroids without obtaining the i.i.d. bootstrap particles first. We also point out that, while we share the same objective as K-means, the optimization algorithms differ. The Expectation-Maximization type of algorithm used by K-means is not applicable to our scenario.
Comparing with Other Particle Improving Approach Intuitively, from a high level abstracted perspective, we provide an approach to use K-means type of idea to improve the particle quality without accessing to the true target distribution. This is the key differentiator of this work to other approaches that improve the particle quality, as they all require to access the target distribution. For example Claici et al. (2018) requires that sampling from target distribution is cheap and easy. Chen et al. (2012; 2018a); Campbell & Beronov (2019) need to access the logarithm of the probability density function of the target distribution. In our problem, neither sampling from the target distribution is cheap nor the logarithm of the probability density function is available, making those approaches not more applicable.

heading: 3.1 TRAINING
text: The optimization of (6) can be solved by gradient descent. Suppose θ∗j (t) is the j-th centroid at iteration t. We initialize {θ∗j (0)}mj=1 by sampling from ρπ and at iteration t, we update θ∗t by applying the gradient descent on the loss in (6), which yields
θ∗j (t+ 1)← θ∗j (t)− tg(θ∗j (t)), g(θ∗j (t)) = ∇θEw∼π [ I{j ∈ uw(t)}Lw(θ∗j (t)) ] /v∗j (t),
(8)
where we define the index of the closest centroid to particle θ̂w as uw(t) = arg minj∈[m] Lw(θ∗j (t)) and v∗j (t) = Pw∼π (j ∈ uw(t)) denotes the probability that centroid j is the one that gives the lowest bootstrap loss. The denominator v∗j (t) in g(θ ∗ j (t)) is optional. However, notice that the magnitude of numerator in g(θ∗k(t)) decays with larger m, which might require an adjustment of the learning rate when m changes. This adjustment can be avoided by rescaling with v∗k(t).
We note that {θ∗j (0)}mj=1 is just m i.i.d. bootstrap particles which is not optimal for approximation and our algorithm can be viewed as an approach for refining the m particles by solving (6). In practice, we find that we can simply use random initialization (e.g., draw θ from some Gaussian distribution) instead.
Centroid Degeneration Phenomenon Naively applying the updating rule (8) may cause a degeneration phenomenon: When a centroid happens to give considerably worse performance than others, which can be caused by the stochasticity of gradient or worse initialization, the performance of this centroid will remain considerably worse throughout the optimization. The reason is simple. As this
5i.e. the regions separated by the dashed lines in the right plot of Figure 1.
centroid (e.g. θ∗j (t)) gives a considerably worse performance, the probability that it gives the lowest bootstrap loss, i.e., v∗j (t), is small. As a consequence, the gradient that updates this centroid is only based on aggregating information from a small low-density region of π and hence can be unstable and further degrades this centroid. Note that this mechanism is self-reinforced since when this centroid cannot be effectively improved in the current iteration, it faces the same issue in the next one. As a result, this centroid is always significantly worse than the others.
We call this undesirable phenomenon centroid degeneration and we want to prevent this phenomenon because when it happens, we have a centroid that is not representative and contributes less to approximating ρπ. We solve this issue with a simple solution and here is the intuition. The reason that a centroid degenerates lies in that this centroid is far from the good region where it gives a good performance. And when this happens, we should push the centroid to move towards this good region, which can be achieved by using the common gradient over the whole training data. Specifically, we define a threshold γ, indicating centroid j is degenerated if v∗j (t) ≤ γ. And when it happens, we update using the common gradient over the whole data:
θ∗j (t+ 1)← θ∗j (t)− t∇θL(θ∗j (t)). (9) In section 4, we give a theoretical analysis on why this modification is important and is able to solve the centroid degeneration issue.
Practical Algorithm In practice, we estimate the gradient by replacing the expectation over w ∼ π in (8) with averaging over M i.i.d. samples {wh}Mh=1 drawn from π:
ĝ(θ∗j (t)) =
∑M h=1 [ I{j ∈ uwh(t)}∇θLwh(θ∗j (t)) ]∑M h=1 I{j ∈ uwh(t)} . (10)
We emphasize that here uwh and Lwh for all wh can be computed very cheaply, enabling us to use a very large M to reduce the error of gradient estimation. Specifically, at iteration t, for each j ∈ [m], we first calculate
L(θ∗j (t)) = [`(x1, fθ∗j (t)), ...`(xn, fθ∗j (t))] > ∈ Rn, (11)
which is the vector encodes the loss of centroid j at each data point. This procedure does not introduce any extra computational overhead compared with standard gradient descent. After that, the bootstrap loss of centroid j can be computed cheaply by Lw(θ∗j (t)) = w>L(θ∗j (t)). Similarly, it is cheap to obtain uwh by
uwh(t) = arg min j∈[m]
w>h L(θ ∗ j (t)). (12)
Taking the modified updating rule introduced to prevent the centroid degeneration phenomenon into account, we update θ∗j (t) by θ ∗ j (t+ 1)← θ∗j (t)− tφ(θ∗j (t)), where
φ(θ∗j ) =
{ ĝ(θ∗j (t)) if ∑ h∈[M ] I{uwh(t) = j}/M > γ
∇θL(θ∗j (t)) otherwise. (13)
Notice that as as L(θ∗j (t)) is pre-computed, calculating Lwh(θ∗j (t)) for many (e.g., M ) different wh is almost free, since it only requires a simple matrix multiplication with O(nM) complexity. Similarly, calculating uwh(t) is also very cheap. In practical implementation, as uwh(t) do not change much within a few iterations, we can update uwh(t) every a few iterations (e.g., every epoch). We can also replace the ∇θLwh(θ∗j (t)) or∇θL(θ∗j (t)) in (13) using a mini-batch of data instead of the whole data, which leads to a stochastic gradient version of our algorithm. We refer readers to Algorithm 1 for the ideal updating and Algorithm 2 for the practical implementation in Appendix B for more details.

heading: 4 THEORY
text: Recall that, as discussed in (5), our approach relies on the intuition that bootstrap particles are nested in a small region so that we can approximate the distance between the centroid and a bootstrap particle by the bootstrap loss of that centroid. The main goal of this section is to give a formal theoretical justification of this intuition.
Before we proceed, we clarify several important setups for establishing and interpreting the theoretical result. As discussed in the introduction, we are mainly interested in the scenerio that the number of available particles/centroids m is small while the number of training data n is large, which motivates us to establish theoretical result in the region of small m and large n. This is significantly different from conventional asymptotic analysis in which we aim to show the behavior when m → ∞. We consider the setting that the parameter dimension d is fixed and does not scale with n.
We are mainly interested in characterizing the approximation of the proposed loss in (6) to the ideal loss in (3), given any small and fixed number m of centroids when n → ∞. This justifies why the proposed centroid approximation method can be viewed as minimizing the Wasserstein distance between the particle distribution ρ∗π and the target bootstrap distribution ρπ .
For simplicity, we build our analysis assuming the ideal update rule (8,9) is used. We start with the following main assumptions.
Assumption 1 (Smoothness and boundedness) Assume that the following quantities are upper bounded by some constant c <∞:
1. max i,j,k∈[d] sup θ∈Θ,x∈X
∂3`(x, fθ)
∂iθi∂θj∂θk ; 2. sup θ1,θ2∈Θ sup x∈X ||∇2θ`(x, fθ1)−∇2θ`(x, fθ2)|| ||θ1 − θ2|| ;
3. sup x∈X ,θ∈Θ ∥∥∇2θ`(x, fθ)∥∥ ; 4. sup θ∈Θ ‖θ‖ .
Assumption 1 is a standard regularity condition on the boundness and smoothness of the problem. Assumption 2 (Asymptotic normality) Assume √ n ( θ̂w − θ̂ ) d→ N (0, A) and √ n ( θ̂ − θ0 ) d→
N (0, A) as n→∞, where A is a positive-definite matrix with the largest eigenvalue bounded.
Assumption 2 is a higher level assumption on the asymptotic normality of the estimators. Such result is classic and can be derived with some weak and technical regularity conditions. See examples in Chatterjee et al. (2005); Cheng et al. (2010).
Assumption 3 (On the global minimizer) Suppose that λmin ( ∇2θL∞(θ0) ) > 0.
Assumption 3 is also standard showing the locally strongly convexity of the loss around the truth θ0.
Assumption 4 (On the learning rate) Suppose that maxt t = O(n−1).
Assumption 4 assumes that the learning rate of the algorithm is sufficiently small such that its induced discretization error is not the dominating term.
The key challenge of our analysis is to show that our dynamics is B(θ0, r)-stable (defined below in Definition 1) for some small r, saying that {θ∗j (t)}mj=1 stay in a small region that is close to θ0 for any iteration t. Combined with the property6 that θ̂w are also close to θ0, the centroids and the bootstrap particles are close to each other and thus our approximation in (5) holds for all t ≥ 0. In this way, optimizing the centroids by minimizing our loss is almost equivalent to optimizing the centroids by minimizing the Wasserstein distance.
Definition 1 (B(θ, r)-stable) Given some θ ∈ Θ and r ≥ 0, we say our dynamics is B(θ, r)-stable if ∀t ≥ 0 and ∀j ∈ [m], θ∗j (t) ∈ B(θ, r), where B(θ, r) := {θ′ : ‖θ′ − θ‖ ≤ r, θ′ ∈ Θ} is the ball with radius r centered at θ.
The key intuition to establish such B(θ0, r)-stable result is to characterize that our optimization dynamics is implicitly self-controlled: when some centroid approaches the boundary of B(θ0, r), the updating mechanism automatically start to push the centroid to move towards the center of the region. Thus, if all the centroids are within B(θ0, r) at initialization, they will alway stay in this region. Thanks to assumption 2, 3, when the dataset is large, the landscape of our loss is locally strongly convex around θ0. When a centroid j is at the boundary of B(θ0, r), it has v∗j (t) < γ and thus the
6This is implied by the asymptotic normality in assumption 2.
updating direction is the gradient of loss L. By the convexity, such gradient will push the centroid move towards the center of B(θ0, r) where the empirical minimizer locates at. On the other hand, for centroid j with v∗j (t) ≥ γ, its updating direction aggregates information from sufficient data point and thus behaves similarly to that of the common gradient, pushing centroid to move towards the center with the centroid is not close to the center.
Theorem 1 Under Assumptions 1-4 and suppose that we initialize θ∗j (0), j ∈ [m] by sampling from ρπ , given any m <∞ and γ > 0, when n is sufficiently large, we have
max j∈[m] sup t≥0 ∥∥θ∗j (t)− θ0∥∥ = Op(√(log n)/n). Here the probability is taken w.r.t. training data.
Theorem 1 implies our dynamics is B(θ0, rn)-stable with rn = O( √
log n/n). The condition that θ∗j (0) ∼ ρπ i.i.d. can be replaced by the condition that θ∗j (0) is sufficiently close to θ0. We need such condition as we uniformly bound the distance between θ∗j (t) and θ0 at any iteration including the first one. Theorem 1 implies that the approximation stated in (5) holds with high probability and hence the proposed loss in (6) is ‘almost as good as’ the ideal loss in (3).
Theorem 2 Under the same assumptions as Theorem 1, given any m < ∞ and γ > 0, when n is sufficiently large, we have
sup t≥0 ∣∣∣∣Ew∼π[ minj∈[m]Lw(θ∗j (t))]−B − Ew∼π[ minj∈[m] ||θ∗j (t)− θ̂w||2D]/2 ∣∣∣∣ = Op(√(log n)/n3/2).
Here the probability is taken w.r.t. training data and B is some constant independent from θ∗j (t) for any t ≥ 0 and j ∈ [m].
Asymptotics when m also grows Although our main interest is the asymptotics with a small, fixed m and growing n, we discuss here on asymptotics when m also grows. As shown in Section 3 and introduction, our method can be viewed as an ‘approximated’ K-means on the target distribution. From Theorem 5.2 in Canas & Rosasco (2012), the particle distribution formed by the optimal centroids learned by K-means gives improved O(m−1/d) convergence to any general target distribution in terms of Wasserstein distance, where d is data dimension. In comparison, the particle distribution of i.i.d. sample only gives O(m−1/(2d+4)) from Theorem 5.1 in Canas & Rosasco (2012). This implies that our approach potentially also has such a rate improvement. Note that the results in Canas & Rosasco (2012) are for general target distribution without any n involves. To rigorously establish the large m asymptotic result for our problem, we need to study the joint limit of n and m. This is indeed very non-trivial: as discussed in Weed et al. (2019) (i.e. Proposition 14), when n→∞, the target distribution ρπ becomes a sharp Gaussian and the convergence rate of i.i.d. bootstrap particles will gradually improve to O(m−1/2) (in a way that depends on n). It implies that when n m→∞, our improvement may become only constant level. We find establishing such a theory is out the scope of this conference paper and leave it as future work.

heading: 5 EXPERIMENT
text: As discussed in the introduction, our main goal is to improve the quality of the particle distribution when only a limited number of particles/centroids is allowed, so that we can use less particles at
deployment, which reduces the memory consumption and the computational cost for making prediction. Thus, our experiment design will be focusing on comparing the testing performance of vanilla bootstrap and our centroid based approach when the same and a small number of particles/centroids is used. We apply our method to four applications: confidence interval construction, bootstrap method in contextual bandit, bootstrapped deep Q-network and bagging. Due to space limit, we refer to Appendix C.4 for the bagging experiment, Appendix C.5 for ablation study on the importance of modifying the gradient to overcome the centroid degeneration phenomenon. Although we are less interested in the computational cost of training, as discussed in Section 3, our method actually only introduces a little training computation overhead, which is another advantage of our method. We draw analysis on this aspect in Appendix C.6.
5.1 BOOTSTRAP CONFIDENCE INTERVAL
We start with a classic application of bootstrap: confidence interval estimation for linear model with parameter θ. Fix confidence level α, we consider three ways to construct (two-sided) bootstrap confidence interval of θ: the Normal interval, the percentile interval and the pivotal interval. And we test m = 20, 50, 100, 200. For all experiments, we repeat with 1000 independent random trials. We consider the standard bootstrap as baseline. Detailed experimental setup are included in Appendix C.1.
Figure 2 shows the Wasserstein distance between the true target distribution ρπ and the empirical distributions obtained by (a) i.i.d. sampling ρ̂π, (b) the proposed centroid approximation ρ∗π. The centroid approximation significantly reduces the Wasserstein distance by a large margin. We then compare the quality of obtained confidence intervals, which is measured by the difference between the estimated coverage probability and the true confidence level, i.e., |α̂− α| (the lower the better). Here we only consider confidence intervals of the first coordinate of θ: θ1. Table 1 summarizes the result with α = 0.9. We see that using more particles is generally able to improve the constructed confidence intervals. We also compare with two variants of standard bootstrap: Bayesian bootstrap (Rubin, 1981) and residual bootstrap (Efron, 1992). And we consider varying α = 0.8, 0.95. These results are included in Appendix C.1.

heading: 5.2 CENTROID APPROXIMATION FOR BOOTSTRAP METHOD IN CONTEXTUAL BANDIT
text: Contextual bandit is a classic task in sequential decision making, in which accurately quantifying the model uncertainty is important in order to achieve good exploration-exploitation trade-off. As shown in Riquelme et al. (2018), tracking the model uncertainty using bootstrap is a strong method for contextual bandit. However, it is costly to maintain a large number of bootstrap models and thus the number of models is typically within 10 (Osband et al., 2016). We find that applying the proposed centroid approximation here can significantly improve the performance. Riquelme et al. (2018) uses m = 3 bootstrap models and we give a more comprehensive evaluation with m = 3, 4, 5, 10. We consider three datasets: Mushroom, Statlog and Financial. We set γ = 0.5/m. We randomly generate 20 different context sequences, apply all the methods and report the averaged cumulative reward and its standard deviation. Table 2 summarizes the result and note that a large part of variance
can be explained by different context sequences. All results in Table 2 are statistically significant under significant level 5% using matched pair t-test. Table 2 shows that using more bootstrap models generally improves the accumulated reward. And when using the same number of models, the proposed centroid approximation method consistently improves over standard bootstrap method. We refer readers to appendix C.2 for more information on the background and experiment.

heading: 5.3 CENTROID APPROXIMATION FOR BOOTSTRAP DQN
text: “Efficient exploration is a major challenge for reinforcement learning (RL). Common dithering strategies such as -greedy do not carry out temporally-extended exploration, which leads to exponentially larger data requirements” (Osband et al., 2016). To tackle this issue, Osband et al. (2016) proposes the Bootstrapped Deep Q-Network (DQN). We apply our centroid approximation to improve Bootstrapped DQN. We consider m = 2, 5 and similar to the experimental setting in contextual bandit, we set γ = 0.5/m. We consider two benchmark environments: LunarLander-v2 and Catcher-v0 from GYM (Brockman et al., 2016) and PyGame learning environment (Tasfi, 2016). We conduct the experiment with 5 independent random trails and report the averaged result with its standard deviation. We refer readers to Appendix C.3 for more background and other experiment details. Figure 3 summarizes the result. For LunarLander-v2, Bootstrap DQN with 2 and 5 heads give similar performance but both converge to a less optimal model compared with the centroid approximation method. Centroid approximation method with 2 and 5 heads performs similarly at convergence but centroid approximation method with 5 heads is able to converge faster than 2-head model and thus has lower regret. For Catcher-v0, adding more heads to the model is able to improve the performance for both methods. The proposed centroid approximation consistently improves over baselines.

heading: 6 RELATED WORK
text: Bootstrap is an classical statistical inference method, which was developed by Efron (1992) and generalized by, i.e., Mammen (1993); Shao (2010); Efron (2012) (just to name a few). Bootstrap can be widely applied to various statistical inference problem, such as confidence interval estimation (DiCiccio et al., 1996), model selection (Shao, 1996), high-dimensional inference (Chen et al., 2018b; El Karoui & Purdom, 2018; Nie & Ročková, 2020), off-policy evaluation (Hanna et al., 2017), distributed inference (Yu et al., 2020) and inference for ensemble model (Kim et al., 2020), etc.
Despite its wide applications and nice theoretical properties, there has been very few works on discussing and improving the approximation efficiency in the region of small bootstrap sample size, beyond the i.i.d. sampling paradigm. While the m-out-of-n bootstrap (Bickel et al., 2012) and the bag of little bootstrap (Kleiner et al., 2014) are designed to reduce the computational cost with the subsampling techniques in the big data settings (large n), they still require a large bootstrap sample size and thus are still not scalable to large deep learning applications.
Bayesian Inference is a different approach to quantify the model uncertainty. Different from frequentists’ method, Bayesian assumes a prior over the model and the uncertainty can be captured by the posterior. Bayesian inference have been largely popularized in machine learning, largely thanks to the recent development in scalable sampling method (Welling & Teh, 2011; Chen et al., 2014; Seita et al., 2018; Wu et al., 2020), variational inference (Blei et al., 2017; Liu & Wang, 2016), and other approximation methods such as Gal & Ghahramani (2016); Lee et al. (2018). In comparison,
bootstrap has been much less widely used in modern machine learning and deep learning. We believe this is largely attributed to the lack of similarly efficient computational methods in the small sample size m region, which is the very problem that we aim to address with our new centroid approximation method.
Uncertainty in Deep Learning In additional to the applications considered in this paper, uncertainty in deep learning model can also be applied to problems including calibration (Guo et al., 2017) and out-of-distribution detection (Nguyen et al., 2015). The definition of uncertainty of neural network is quite generalized (e.g., Gal & Ghahramani (2016); Ovadia et al. (2019); Maddox et al. (2019); Van Amersfoort et al. (2020)) and can be quite different from the uncertainty that bootstrap inference want to quantify and can be approached with various methods including drop out (Gal & Ghahramani, 2016; Durasov et al., 2020), label smoothing (Qin et al., 2020), designing new modules in the model (Kivaranovic et al., 2020), adversarial training (Lakshminarayanan et al., 2017) and Bayesian modeling (Blundell et al., 2015), etc. This paper focuses on improving the bootstrap method and thus is orthogonal to those previous works. Pearce et al. (2018); Salem et al. (2020) also try to refine the ensemble models to improve the quality of prediction interval. Compare with our method, their method can only be applied to prediction interval and does not have theoretical guarantee.

heading: 7 CONCLUSION
text: We propose a centroid approximation method to learn an improved particle distribution that better approximates the target bootstrap distribution, especially in the region with small particle size. Theoretically, when the size of training data is large, our objective function is surrogate to the Wasserstein distance between the particle distribution and target distribution. Thus, compared with standard bootstrap, the proposed centroid approximation method actively optimizes the distance between particle distribution and target distribution. The proposed method is simple and can be flexibly used for applications of bootstrap with negligible extra computational cost.

heading: A PROOF
text: We also show Theorem 3, which gives a formal characterization of the Taylor approximation intuition introduced in (5).
Theorem 3 Under Assumption 1 and 2, when n is sufficiently large, we have
Lw(θ) = Lw(θ̂w) + 12 ( θ − θ̂w )> ∇2θL∞(θ0) ( θ − θ̂w ) +Op ( ||θ − θ̂w||2(n−1/2 + ||θ − θ̂w||) ) .
Here the stochastic boundedness is taken w.r.t. the training data and w.
In the proof, we may use c to represent some absolute constant, which may vary in different lines.

heading: A.1 PROOF OF THEOREM 3
text: With the fact that∇wL(θ̂w) = 0 and under assumption 1, using Taylor expansion, we have
Lw(θ) = Lw(θ̂w) + 1
2
( θ − θ̂w )> ∇2θLw(θ̂w) ( θ − θ̂w ) +O (∥∥∥θ − θ̂w∥∥∥3) . Notice that ∣∣∣∣(θ − θ̂w)> (∇2θLw(θ̂w)−∇2θL∞(θ0))(θ − θ̂w)∣∣∣∣
≤ ∥∥∥θ − θ̂w∥∥∥2 ∥∥∥∇2θLw(θ̂w)−∇2θL∞(θ0)∥∥∥ ≤ ∥∥∥θ − θ̂w∥∥∥2 (∥∥∥∇2θLw(θ̂w)−∇2θLw(θ0)∥∥∥+ ∥∥∇2θLw(θ0)−∇2θL∞(θ0)∥∥) ≤ ∥∥∥θ − θ̂w∥∥∥2 (C ∥∥∥θ̂w − θ0∥∥∥+ ∥∥∇2θLw(θ0)−∇2θL∞(θ0)∥∥F) ,
where we denote the Frobenius norm as ‖·‖F . With assumption 2, we have ∥∥∥θ̂w − θ0∥∥∥ = Op(n−1/2).
By applying centroid limit theorem and delta method to ∣∣∣∇2θijLw(θ0)−∇2θijL∞(θ0)∣∣∣ for every pair
i, j ∈ [d], we have ∥∥∇2θLw(θ0)−∇2θL∞(θ0)∥∥F = Op(n−1/2). Thus we obtained the desired result.

heading: A.2 PROOF OF THEOREM 1
text: Given any radius r and > 0, with sufficiently large n, we have P (∥∥∥θ̂w − θ0∥∥∥ ≥ r) ≤ exp(−λnr2)+ /m,
for λ = 14λmax(A) −1. Here the probability is the jointly probability of bootstrap weight and training data. Thus, given any r, under the assumption that θ∗j (0) is initialized via sampling θ̂w, then we have
P ( ∪j∈[m] {∥∥θ∗j (0)− θ0∥∥ ≥ r}) ≤ ∑ j∈[m] P (∥∥θ∗j (0)− θ0∥∥ ≥ r) ≤ m exp(−λnr2) + .
We proof by induction. Given any {θj}mj=1, define Rk,r = I { w ∈ supp(π) : arg min
j∈[m] Lw(θj) = k and ∥∥∥θ̂w − θ0∥∥∥ ≤ r} . Suppose at iteration t, we have ‖θ∗k(t)− θ0‖ ≤ cα √ logn n λ0γ for some constant c and λ0, which we denote as the minimum eigenvalue of∇2θL∞(fθ0). Now at iteration t, we have two cases.
Case 1: EπRk,∞ ≥ γ Suppose that at iteration t, for k such that EπRk,∞ ≥ γ, and ‖θ∗k(t)− θ0‖ = qk, we have the following property:
‖θ∗k(t+ 1)− θ0‖ 2 = ∥∥∥∥θ∗k(t)− tEπRk,∞Eπ [∇θLw(θ∗k(t))Rk,∞]− θ0 ∥∥∥∥2
= ‖θ∗k(t)− θ0‖ 2 − 2 t
EπRk,∞ 〈θ∗k(t)− θ0,Eπ [∇θLw(θ∗k(t))Rk,∞]〉+ 2t ‖Eπ [∇θLw(θ∗k(t))Rk,∞]‖ 2 .
Notice that
Eπ [∇θLw(θ∗k(t))Rk,qk ] (1) = Eπ [ ∇2θLw(θ̂w)(θ∗k(t)− θ̂w)Rk,qk ] + o ( q2k )
(2) = Eπ [ ∇2θLw(θ̂w)(θ∗k(t)− θ0)Rk,qk ] +O ( q2k )
(3) = Eπ [ ∇2θLw(θ0)(θ∗k(t)− θ0)Rk,qk ] +O ( q2k )
Here (1) is obtained via applying Taylor expansion on∇θLw(θ∗k(t)) at θ̂w. (2) is by assumption 1 and 2. (3) is by assumption 1. We thus have
− 〈θ∗k(t)− θ0,Eπ [∇θLw(θ∗k(t))Rk,∞]〉 ≤ − 〈θ∗k(t)− θ0,Eπ [∇θLw(θ∗k(t))Rk,qk ]〉+ ‖θ∗k(t)− θ0‖ ‖Eπ∇θLw(θ∗k(t))(1−Rk,qk)‖ ≤ − 〈θ∗k(t)− θ0,Eπ [∇θLw(θ∗k(t))Rk,qk ]〉+ cqk exp(−λnq2k) ≤− EπRk,qk(θ∗k(t)− θ0)>∇2θLw(θ0)(θ∗k(t)− θ0) + cqk exp(−λnq2k) +O ( q3k ) .
Notice that with sufficiently large n, with central limit theorem, we have
− EπRk,qk(θ∗k(t)− θ0)>∇2θLw(θ0)(θ∗k(t)− θ0) ≤‖θ∗k(t)− θ0‖ 2 Eπ ∥∥∇2θLw(θ0)−∇2θLw(θ0)∥∥− EπRk,qk(θ∗k(t)− θ0)>∇2θL∞(θ0)(θ∗k(t)− θ0) =− EπRk,qk(θ∗k(t)− θ0)>∇2θL∞(θ0)(θ∗k(t)− θ0) +O(n−1/2).
This gives that
− 〈θ∗k(t)− θ0,Eπ [∇θLw(θ∗k(t))Rk,∞]〉 ≤ − EπRk,qk(θ∗k(t)− θ0)>∇2θL∞(θ0)(θ∗k(t)− θ0) + cqk exp(−λnq2k) +O ( q3k + qkn −1/2 )
≤− λ0EπRk,qk ‖θ∗k(t)− θ0‖ 2 + cqk exp(−λnq2k) +O ( q3k + qkn −1/2 ) .
Use the above estimation, we have
‖θ∗k(t+ 1)− θ0‖ 2
≤‖θ∗k(t)− θ0‖ 2 − 2 tλmin EπRk,qk EπRk,∞ ‖θ∗k(t)− θ0‖ 2
+2c tqk exp(−λnq2k)/EπRk,∞ +O ( t(q 3 k + qkn −1/2)/EπRk,∞ + 2t )
≤‖θ∗k(t)− θ0‖ 2 + t
EπRk,∞
( −2λ0EπRk,qk ‖θ∗k(t)− θ0‖ 2 − 2cqk exp(−λnq2k) +O ( q3k + t + qkn −1/2 )) .
Notice that by choosing α > √ 1/(2λ) and t = O(n−1), with sufficiently large n, when
‖θ∗k(t)− θ0‖ ≥ cα √ logn n
λ0EπRk,qk for some constant c, we have
‖θ∗k(t+ 1)− θ0‖ ≤ ‖θ∗k(t)− θ0‖ .
Thus ‖θ∗k(t+ 1)− θ0‖ ≤ cα √ logn n λ0EπRk,∞ ≤ cα √ logn n λ0γ for some constant c.
Case 2: EπRk,∞ ≤ γ In this case, we have
‖θ∗k(t+ 1)− θ0‖ 2 = ∥∥∥θ∗k(t)− t∇θL(fθ∗k(t))− θ0∥∥∥2
= ‖θ∗k(t)− θ0‖ 2 − 2 t 〈θ∗k(t)− θ0,∇θL(θ∗k(t))〉+ 2t ∥∥∥∇θL(fθ∗k(t))∥∥∥2 . Notice that −〈θ∗k(t)− θ0,∇θL(θ∗k(t))〉 ≤ − 〈 θ∗k(t)− θ0,∇2θL(θ0) (θ∗k(t)− θ0) 〉 − 〈θ∗k(t)− θ0,∇θL(fθ0)〉+ o(||θ∗k(t)− θ0||3)
= − (θ∗k(t)− θ0) >∇2θL∞(θ0) (θ∗k(t)− θ0) + o(||θ∗k(t)− θ0||3) +Op(n−1/2)||θ∗k(t)− θ0||.
This gives that
‖θ∗k(t+ 1)− θ0‖ 2 ≤ ‖θ∗k(t)− θ0‖ 2−2 tλ0 ‖θ∗k(t)− θ0‖ 2 +o( t||θ∗k(t)−θ0||3+ 2t )+Op(n−1/2) t||θ∗k(t)−θ0||. With t = O(n−1) and sufficiently large n, when
‖θ∗k(t)− θ0‖ ≥ cα √ logn n
λ0γ ,
we have ‖θ∗k(t+ 1)− θ0‖ 2 ≤ ‖θ∗k(t)− θ0‖ 2. Combine this two cases, we conclude that ‖θ∗k(t+ 1)− θ0‖ ≤ cα √ logn n λ0γ for any t, when
‖θ∗k(0)− θ0‖ ≤ cα √ logn n λ0γ . We thus conclude that, for any α > √ 1/(2λ) and > 0, when n
is sufficiently large, with probability at least 1−m exp ( −λ cα
2 logn λ20γ 2
) − , we have
max j∈[m] sup t
∥∥θ∗j (t)− θ0∥∥ ≤ cα √ logn n
λ0γ .

heading: A.3 PROOF FOR THEOREM 2
text: Notice that
Lw(θ∗j (t))− Lw(θ̂w) = 1
2
( θ∗j (t)− θ̂w )> ∇2θLw(θ̂w) ( θ∗j (t)− θ̂w ) +O( ∥∥∥θ∗j (t)− θ̂w∥∥∥3) = 1
2
( θ∗j (t)− θ̂w )> ∇2θLw(θ0) ( θ∗j (t)− θ̂w ) +O( ∥∥∥θ∗j (t)− θ̂w∥∥∥3) +O(∥∥∥θ∗j (t)− θ̂w∥∥∥2 ∥∥∥θ̂w − θ0∥∥∥) = 1
2 ∥∥∥θ∗j (t)− θ̂w∥∥∥2 D + ∥∥∇2θLw(θ0)−∇2θL∞(θ0)∥∥ ∥∥∥θ∗j (t)− θ̂w∥∥∥2 +O(∥∥∥θ∗j (t)− θ̂w∥∥∥3)
+O( ∥∥∥θ∗j (t)− θ̂w∥∥∥2 ∥∥∥θ̂w − θ0∥∥∥)
Given w, we define uw = arg minj∈[m] ∥∥∥θ∗j (t)− θ̂w∥∥∥2
D . For any α >
√ 1/(2λ) and > 0, when n
is sufficiently large, with probability at least 1−m exp ( −λ cα
2 logn λ20γ 2
) − , we have
1 2 Ew∼π [ min j∈[m] ∥∥∥θ∗j (t)− θ̂w∥∥∥2 D ] = 1
2 Ew∼π [∥∥∥θ∗uw − θ̂w∥∥∥2 D ] ≥Ew∼π [ Lw(θ∗uw)− Lw(θ̂w)− c
∥∥∥θ∗uw − θ̂w∥∥∥2 (∥∥∥θ̂w − θ0∥∥∥+ ∥∥∇2θLw(θ0)−∇2θL∞(θ0)∥∥+ ∥∥∥θ∗uw − θ̂w∥∥∥)] ≥Ew∼π [ Lw(θ∗uw) ] − Ew∼π [ Lw(θ̂w) ] − c ( α √ log n
λ0γn3/2 ) =Ew∼π [ min j∈[m] Lw(θ∗j (t)) ] − Ew∼π [ Lw(θ̂w) ] − c ( α √ log n
λ0γn3/2
) .
Algorithm 1 Ideal algorithm for centroid approximation with full-batch gradient used and wh updated every iteration.
1: Initialize θ∗j (0), j ∈ [m] by i.i.d. sampling from ρπ or other distribution such as Gaussian. 2: for t ∈ iterations do 3: ∀j ∈ [m], calculate L(θ∗j (t)) defined in (11) 4: Sample {wh}Mh=1, i.i.d. from π. 5: ∀h ∈ [M ] and j ∈ [m], calculate Lwh(θ∗j (t)) = wThL(θ∗j (t)). 6: ∀h ∈ [M ], calculate uwh defined in (12) for each h. 7: ∀j ∈ [m], update θ∗j by (13). 8: end for
Similarly, we also have, with probability at least 1−m exp ( −λ cα
2 logn λ20γ 2
) ,
Ew∼π [
min j∈[m]
Lw(θ∗j (t)) ] −Ew∼π [ Lw(θ̂w) ] ≥ 1
2 Ew∼π [ min j∈[m] ∥∥∥θ∗j (t)− θ̂w∥∥∥2 D ] −c ( α √ log n λ0γn3/2 ) .
Notice that the above bound holds uniformly for all j ∈ [m] and any iteration t, which implies that with probability at least 1− 2m exp ( −λ cα
2 logn λ20γ 2
) − 2 , we have
sup t≥0 ∣∣∣∣Ew∼π[ minj∈[m]Lw(θ∗j (t))]−B − Ew∼π[ minj∈[m] ||θ∗j (t)− θ̂w||2D]/2 ∣∣∣∣ ≤ c(α√log nλ0γn3/2 ) .

heading: B ALGORITHM BOX
text: We provide pseudo algorithm for the ideal centroid approximation algorithm in Algorithm 1. In practical implementation, we do not need to update wh every iteration and can also replace the full-batch gradient by stochastic gradient. Specifically, notice that ˆg(θ∗j ) defined in (10) can be alternative represented as
ĝ(θ∗j (t)) =
∑M h=1 ∑n i=1 [I{j ∈ uwh(t)}]wh,i∇θ`(xi, fθ∗j (t))/n∑M
h=1 [I{j ∈ uwh}] =
1
n n∑ i=1 qi,j∇θ`(xi, fθ∗j (t)),
(14)
where qi,j is defined by
qij :=
∑M h=1 ∑n i=1 [I{j ∈ uwh(t)}]wh,i∑M h=1 [I{j ∈ uwh(t)}] . (15)
This allows us to use a stochastic gradient version of gradient
ĝsgd(θ ∗ j ) =
1 |B| ∑ i∈[B] qi,j∇θ`(xi, fθ∗j ), (16)
where B is the set of mini-batch data. The detailed algorithm is summarized in Algorithm

heading: C ADDITIONAL EXPERIMENT DETAILS

heading: C.1 BOOTSTRAP CONFIDENCE INTERVAL
text: Given a model fθ parameterized by θ and a training set with n data points i.i.d. sampled from population, our goal is to construct confidence interval for θ. Let ρ̃π be an empirical distribution approximating ρπ, which could be obtained by i.i.d. sampling, or by our centroid method. Denote by Q[α, ρ̃π] the α-quantile function of ρ̃π with some α ∈ [0, 1]. We consider the following three ways to construct (two-sided) bootstrap confidence interval of θ with confidence level α: the Normal interval, the percentile interval and the pivotal interval which are defined below.
Algorithm 2 Practical implementation of centroid approximation with less frequent updating of wh and stochastic gradient enabled.
1: Initialize θ∗j (0), j ∈ [m] by i.i.d. sampling from ρπ or other distribution such as Gaussian. 2: for t ∈ iterations do 3: // Update wh only every a few iterations. 4: if t mod freq == 0 then 5: ∀j ∈ [m], calculate L(θ∗j (t)) defined in (11) 6: Sample {wh}Mh=1, i.i.d. from π. 7: ∀h ∈ [M ] and j ∈ [m], calculate Lwh(θ∗j (t)) = wThL(θ∗j (t)). 8: ∀h ∈ [M ], calculate uwh(t) defined in (12) for each h. 9: else
10: uwh(t) = uwh(t− 1) 11: end if 12: ∀j ∈ [m], update θ∗j (t) by (13). (May use mini-batch gradient defined in (16)). 13: end for
Methods to construct confidence interval The methods we used to construct confidence interval are – The Normal interval:
[θ̂ − z((1 + α)/2)ŝeboot, θ̂ + z((1 + α)/2)ŝeboot],
where z(·) is the inverse cumulative distribution function of standard Normal distribution. And ŝeboot is the standard deviation estimated from ρ̃π .
– The percentile intervals:
[Q[(1− α)/2, ρ̃π], Q[(1 + α)/2, ρ̃π]].
– The pivotal interval:
[2θ̂ −Q[(1 + α)/2, ρ̃π], 2θ̂ −Q[(1− α)/2, ρ̃π]].
We consider the following simple linear regression: x ∼ N (0, I), y | x ∼ N (θ>x, I), where the features x ∈ R4 and we set the true parameter to be θ0 = [1,−1, 1,−1]. We consider n = 50 and the number of particles m = 20, 50, 100, 200. We compare the coverage probability and the confidence level α to measure the quality:
Measuring the quality of confidence interval With a large number N of independently generated training data (we use N = 1000), we are able to obtain the corresponding confidence intervals {CI(α)s}Ns=1and thus obtain the probability that the true parameter falls into the confidence intervals, which is the estimated coverage probability
α̂ = 1
N N∑ s=1 I{θ0 ∈ CI(α)s}.
A good confidence interval should have α̂ close to α. Thus we measure the performance by calculating the difference |α− α̂|.
As θ̂w is the least square estimator of the bootstrapped dataset, it has analytic solution and thus can be obtained via some matrix multiplications. θ∗w is initialized using θ̂w and then updated for 2000 steps. For this experiment, we find that adding the threshold γ does not gives further improvement for this experiment and thus we simply set γ = 0 and use M = 1. We approximate the true bootstrap distribution by sampling 10000 i.i.d. samples.
More experiment result Table 3 all the result we have varying α = 0.8, 0.9, 0.95, m = 20, 50, 100, 200 and three different approaches for constructing confidence interval. As we can see, centroid approximation gives the best performance in most cases compared with the other three baselines.
Algorithm 3 Algorithm for Centroid Approximation Applied to Contextual Bandit. 1: Obtain a randomly initialized θ∗j (0), j ∈ [m]. 2: Initialize a common replay buffer Rc = ∅ recording all the observed contexts. 3: For each model, initialize its own replay buffer Rj = ∅ that is used for training. 4: for t ∈ number of total steps do 5: Obtain the t-th context xt. 6: Sampling one model based on probability {v∗j (t)}mj=1 to make action at and get reward rt. 7: Update the common replay buffer by Rc ← Rc ∪ {(xt, at, rt)} 8: // Update wh and Rj and model every a few iterations. 9: if t mod freq == 0 then 10: ∀j ∈ [m], calculate L(θ∗j (t)) defined in (11) for all the contexts in Rc. // L(θ∗j (t)) ∈ R|Rc|. 11: Generate M sets of random weights {wh}Mh=1 of contexts in Rc from π. 12: ∀h ∈ [M ] and j ∈ [m], calculate Lwh(θ∗j (t)) = wThL(θ∗j (t)). 13: ∀h ∈ [M ], calculate uwh(t) defined in (12) for each h. 14: ∀i ∈ [|Rc|] and j ∈ [m], calculate qi,j by (15) 15: ∀j ∈ [m], update v∗j (t) based on (7). 16: ∀j ∈ [m], if v∗j (t) ≤ γ, construct Rj = Rc, else, construct Rj by sample |Rc| contexts
in Rc. The probability that context i is being sampled is qi,j/ ∑|Rc| i=1 qi,j .
17: ∀j ∈ [m], train model j using the data in Rj for several iterations. 18: end if 19: end for

heading: C.2 CENTROID APPROXIMATION FOR BOOTSTRAP METHOD IN CONTEXTUAL BANDIT

heading: C.2.1 MORE BACKGROUND
text: Contextual bandit is a classic task in sequential decision making problem in which at time t = 1, ..., n, a new context xt arrives and is observed by an algorithm. Based on its internal model, the algorithm selects an actions at and receives a reward rt(xt, at) related to the context and action. During this process, the algorithm may update its internal model with the newly received data. At the end of this process, the cumulative reward of the algorithm is calculated by r = ∑n t=1 rt and the goal for the algorithm is to improve the cumulative reward r. The exploration-exploitation dilemma is a fundamental aspect in sequential decision making problem such as contextual bandit: the algorithm needs to trade-off between the best expected action returned by the internal model at the moment (i.e., exploitation) with potentially sub-optimal exploratory actions. Thompson sampling (Thompson, 1933; Wyatt, 1998; May et al., 2012) is an elegant and effective approach to tackle the explorationexploitation dilemma using the model uncertainty, which can be approached with various methods including Bayesian posterior (Graves, 2011; Welling & Teh, 2011), dropout uncertainty (Srivastava et al., 2014; Hron et al., 2017) and Bootstrap (Osband et al., 2016; Hao et al., 2019). The ability to accurately assess the uncertainty is a key to improve the cumulative reward. Bootstrap method for contextual bandit maintains m bootstrap models trained with different bootstrapped training set. When conducting an action, the algorithm uniformly samples a model and then selects the best action returned by the sampled model.

heading: C.2.2 MORE EXPERIMENT SETUP DETAILS
text: We set all the experimental setting including data preprocessing, network architecture and training pipeline exactly the same as the one used in Riquelme et al. (2018) and adopt its open source code repository.
Network architecture Following Riquelme et al. (2018), we consider a fully connected feed forward network with two hidden layers with 50 hidden units and ReLU activations. The input and output dimensions depends on the dimension of context and number of possible actions.
Training For each dataset, we randomly generate 2000 contexts, and for each algorithm, we update the replay memory buffer for each model every 50 contexts, and each model is also updated every 50 contexts. For the standard bootstrap, when updating the replay buffer of each model, we sample 50 i.i.d. contexts with uniform probability from the latest 50 contexts (each model have different
realizations of the samples) and add the newly sampled contexts to each model’s replay buffer. For the centroid approximation, we update the replay buffer of each model by applying resampling on all the observed contexts up to the current steps. The resampling probability of each context for each model is different and determined by the algorithm. We refer readers to Algorithm 3 for the detailed procedures. Here we choose freq = 50 andM = 100. When at model updating, each model is trained for 100 iterations with batch size 512 using the data from its replay buffer. Following Riquelme et al. (2018), we use RMSprop optimizer with learning rate 0.1 for optimizing. When making actions, we sample the prediction head according to v∗k(t) obtained using the examples between the last two model updates.
Notice that in the implementation, we only need to maintain one common replay buffer and the replay buffer for each model can be implemented by maintaining the number of each context. Thus when sampling batches of context, we simply need to sample the index of the context and refer to the common replay buffer to get the actual data.

heading: C.3 CENTROID APPROXIMATION FOR BOOTSTRAP DQN

heading: C.3.1 MORE BACKGROUND
text: Similar to the bootstrap method in contextual bandit problem, Bootstrap DQN explores using the model uncertainty, which can be assessed via maintaining several models trained with bootstrapped training set. Maintaining several independent models can be very expensive in RL and to reduce the computational cost, Bootstrap DQN uses a multi-head network with a shared base. Each head in the network corresponds to a bootstrap model and the common shared base is thus trained via the union of the bootstrap training set of each head. We train the Bootstrap DQN with standard updating rule for DQN and use Double-DQN (Van Hasselt et al., 2016) to reduce the overestimate issue. Notice that our centroid approximation method only changes the memory buffer for each head and thus introduces no conflict to other possible techniques that can be applied to Bootstrap DQN.

heading: C.3.2 MORE EXPERIMENT SETUP DETAILS
text: Network Architecture Following Osband et al. (2016), we considered multi-head network structure with a shared base layer to save the memory. Specifically, we use a fully connected layer with 256 hidden neurons as the shared base and stack two fully connected layers each with 256 hidden neurons as head. Each head in the model can be viewed as one bootstrap particles and in computation, all the bootstrap particles use the same base layer.
Training and Evaluation For LunarLander-v2, we train the model for 450 episodes with the first 50 episodes used to initialize the common memory buffer. The maximum number of steps within each episode is set to 1000 and we report the moving average reward with window width 100. For Catcher-v0, we train the model for 100 episodes with the first 10 episodes used to initialize the common memory buffer. We set the maximum number of steps within each episodes 2000 and report the moving average reward with window width 25.
For training the Bootstrap DQN, given the current state xt, we sample one particle based on {v∗j }mj=1 and use its policy network to make an action at and get the reward rt and next state xt+1. The Q-value of the state action pair Q(xt, at) is estimated by rt + λ ∗ Q̂(xt+1), where Q̂(xt+1) is the predicted state value by the target network of the sampled particle and λ is the discount factor set to be 0.99. At each step, the policy network of all particles are updated using one step gradient descent with Adam optimizer (β = (0.9, 0.999) and learning rate 0.001) and mini-batch data (size 64) sampled from its replay buffer. We update target model, each particle’s replay buffer and v∗j s every 1000 steps for LunarLander-v2 and every 200 steps for Catcher-v0. The update scheme for replay buffers of each particles and v∗j s is the same as the one in contextual bandit experiment. As the model see significantly larger number of contexts than that in the contextual bandit experiment, to reduce the memory consumption, we set the max capacity of the common replay buffer to 50000 (the oldest data point will be pop out when the size reaches maximum and new data comes in). For training the shared base, following Osband et al. (2016), we adds up all the gradient comes from each head and normalizes it by the number of heads. Algorithm 4 summarizes the whole training pipeline.
Algorithm 4 Algorithm for Centroid Approximation Applied to DQN. 1: Obtain a randomly initialized θ∗j (0), j ∈ [m]. (For the j-th particle, both of its target and policy
network use the same initialization.) 2: Initialize a common replay buffer Rc = ∅ recording all the observed contexts. 3: For each head, initialize its own replay buffer Rj = ∅ that is used for training. 4: for t ∈ number of total episodes do 5: while not at terminal state and the number of steps does not exceed the threshold do 6: Obtain the t-th context xt. 7: Sample an head based on {v∗j } to make action at and get Q(xt, at) using the reward rt
and the prediction of the corresponding target network. 8: Update the common replay buffer by Rc ← Rc ∪ {(xt, Q(xt, at))} 9: ∀j ∈ [m], update its policy network by one step gradient descent using the data from its
reply buffer. 10: // Update wh and Rj and target network every a few iterations. 11: if t mod freq == 0 then 12: ∀j ∈ [m], calculate L(θ∗j (t)) defined in (11) for all the contexts in Rc. 13: Generate M sets of random weights {wh}Mh=1 of contexts in Rc from π. 14: ∀h ∈ [M ] and j ∈ [m], calculate Lwh(θ∗j (t)) = wThL(θ∗j (t)). 15: ∀h ∈ [M ], calculate uwh(t) defined in (12) for each h. 16: ∀i ∈ [|Rc|] and j ∈ [m], calculate qi,j by (15) 17: ∀j ∈ [m], update v∗j (t) based on (7). 18: ∀j ∈ [m], if v∗j (t) ≤ γ, construct Rj = Rc, else, construct Rj by sample |Rc|
contexts in Rc. The probability that context i is being sampled is qi,j/ ∑|Rc| i=1 qi,j .
19: ∀j ∈ [m], update the j-th target network by loading the weights of the j-th policy network. 20: end if 21: end while 22: end for

heading: C.4 BOOTSTRAP ENSEMBLE MODEL
text: Ensemble of deep neural networks have been successfully used to boost predictive performance (Lakshminarayanan et al., 2017). In this experiment, we consider using an ensemble of deep neural network trained on different bootstrapped training set, which is also known as a popular strategy called bagging.
We consider image classification task on CIFAR-100 and use standard VGG-16 (Simonyan & Zisserman, 2014) with batch normalization. We apply a standard training pipeline. We train the bootstrap model for 160 epochs using SGD optimizer with 0.9 momentum and batchsize 128. The learning rate is initialized to be 0.1 and is decayed by a factor of 10 at epoch 80 and 120. We start to apply the centroid approximation at epoch 120 (thus the centroid is initialized with 120 epochs’ training). We generate the bootstrap training set for each centroid every epoch using the proposed centroid approximation method. We consider m = 3, 4, 5, 10 ensembles and use γ = 0.5/m. We
repeat the experiment for 3 random trials and report the averaged top1 and top5 accuracy with the standard deviation of the mean estimator. Algorithm
Figure 4 summarizes the result. Overall, increasing m is able to improve the predictive performance and with the same number of models, our centroid approximation consistently improves over standard bootstrap ensembles.

heading: C.5 ABLATION STUDY
text: We study the effectiveness of using (8) to modify the gradient of centroid with v∗k(t) ≤ γ. We consider the setting γ = 0 (no modification) and γ = m (always modify, equivalent to no bootstrap uncertainty) and applied the method on the mushroom dataset in the contextual bandit problem. Table 4 shows that (i) modifying the gradient of centroid with small v∗k(t) using do improve the overall result; (2) bootstrap uncertainty is important for exploration.

heading: C.6 COMPUTATION OVERHEAD
text: Our main goal is not to decrease the training cost but improve the quality of bootstrap partical distribution so that we can use less models at deployment and hence reduce the memory cost and the computational cost for inference. Actually, as discussed in Section 3, our method actually only introduces a little computation overhead while much improves the quality of the particles, which is another advantage of our method. For example, in bandit problem on mushroom dataset, when m = 3, 10, vanilla bootstrap takes 33s, 101s while our approach takes 35s, 108s per run. For the bagging, when m = 3, 10, vanilla bootstrap takes 11200s, 34240s while ours take 12000s, 36200s. Results are based on the average of 3 runs. Our method only introduces about 7% computational overhead even with an naive implementation.

Review summaries:
This paper aims to improve the quality of the parametric bootstrap when the number of particles is limited. It does so by constructing a small collection of "centroids" that well-approximate the ideal bootstrap distribution on parameter space in a certain Wasserstein distance. Practically, the centroids are constructed by gradient descent with respect to a loss function involving the negative log-likelihood. Theoretically, this paper shows that the surrogate loss converges to the Wasserstein distance under some assumptions. Empirically, this paper applies the centroid methodology to confidence interval construction, contextual bandits, the bootstrap deep q-network, and ensemble learning .

Bootstrap is a versatile technique for uncertainty quantification and a viable alternative to Bayesian inference, but it can be expensive in terms of memory and computations. With a small number of particles/samples, bootstrap may perform poorly. As a consequence, bootstrap can be quite expensive and even unaffordable for deep learning problems with huge models. The paper proposes to address the question of improving the accuracy of bootstrap inference when the number of particles is limited.   The approach is to minimize the Wasserstein distance between the ideal bootstrap distribution (denoted by $\rho_{\pi}$ in the paper) and a probability mass mass function supported on a small number of particles (called "centroids" in the paper). The exposition given around equations (3) and (4) plausible decomposes this problem into a coupled two-stage optimization problem. Unfortunately, $\rho_{\pi}$ is not observe, which is a major challenge to this conceptual development. The paper then relies on the fact that in standard problems when training data size is large, the bootstrap obtains properties similar (or identical) to the MLE or Bayes estimators, and the algorithmic steps are base don this framework. 

The paper proposes a way to optimise the choice of Bootstrap samples so as to obtain efficient inference with smaller number of such samples. While traditional Bootstrap requires the need to use a large number of samples to approximate the target distribution, the authors propose a combination of k-means type algorithm is association with Bootstrapping to achieve the goal. In addition, they provide a theoretical justification for the method along with providing some experimental evaluation .

This paper proposes an acceleration technique for bootstrap. Here a bootstrap particle corresponds to a machine learning model and for a large number of bootstrap particles, a huge computational cost is incurred. The authors propose using a k-means-like centroid approximation which optimizes a set of centroids (each of which corresponds to a machine learning model) minimizing a data-dependent Wasserstein distance.