A Scalable Bootstrap for Massive Data

Introduction

The development of the bootstrap and related resampling-based methods in the 1960s and 1970s heralded an era in statistics in which inference and computation became increasingly intertwined Efron (1979); Diaconis and Efron (1983). By exploiting the basic capabilities of the classical von Neumann computer to simulate and iterate, the bootstrap made it possible to use computers not only to compute estimates but also to assess the quality of estimators, yielding results that are quite generally consistent Bickel and Freedman (1981); Giné and Zinn (1990); van der Vaart and Wellner (1996) and often more accurate than those based upon asymptotic approximation Hall (1992). Moreover, the bootstrap aligned statistics to computing technology, such that advances in speed and storage capacity of computers could immediately allow statistical methods to scale to larger datasets.

Two recent trends are worthy of attention in this regard. First, the growth in size of datasets is accelerating, with “massive” datasets becoming increasingly prevalent. Second, computational resources are shifting toward parallel and distributed architectures, with multicore and cloud computing platforms providing access to hundreds or thousands of processors. The second trend is seen as a mitigating factor with respect to the first, in that parallel and distributed architectures present new capabilities for storage and manipulation of data. However, from an inferential point of view, it is not yet clear how statistical methodology will transport to a world involving massive data on parallel and distributed computing platforms.

While massive data bring many statistical issues to the fore, including issues in exploratory data analysis and data visualization, there remains the core inferential need to assess the quality of estimators. Indeed, the uncertainty and biases in estimates based on large data can remain quite significant, as large datasets are often high dimensional, are frequently used to fit complex models with large numbers of parameters, and can have many potential sources of bias. Furthermore, even if sufficient data are available to allow highly accurate estimation, the ability to efficiently assess estimator quality remains essential to allow efficient use of available resources by processing only as much data as is necessary to achieve a desired accuracy or confidence.

The bootstrap brings to bear various desirable features in the massive data setting, notably its relatively automatic nature and its applicability to a wide variety of inferential problems. It can be used to assess bias, to quantify the uncertainty in an estimate (e.g., via a standard error or a confidence interval), or to assess risk. However, these virtues are realized at the expense of a substantial computational burden. Bootstrap-based quantities typically must be computed via a form of Monte Carlo approximation in which the estimator in question is repeatedly applied to resamples of the entire original observed dataset.

Because these resamples have size on the order of that of the original data, with approximately 63% of data points appearing at least once in each resample, the usefulness of the bootstrap is severely blunted by the large datasets increasingly encountered in practice. In the massive data setting, computation of even a single point estimate on the full dataset can be quite computationally demanding, and so repeated computation of an estimator on comparably sized resamples can be prohibitively costly. To mitigate this problem, one might naturally attempt to exploit the modern trend toward parallel and distributed computing. Indeed, at first glance, the bootstrap would seem ideally suited to straightforwardly leveraging parallel and distributed computing architectures: one might imagine using different processors or compute nodes to process different bootstrap resamples independently in parallel. However, the large size of bootstrap resamples in the massive data setting renders this approach problematic, as the cost of transferring data to independent processors or compute nodes can be overly high, as is the cost of operating on even a single resample using an independent set of computing resources.

While the literature does contain some discussion of techniques for improving the computational efficiency of the bootstrap, that work is largely devoted to reducing the number of resamples required Efron (1988); Efron and Tibshirani (1993). These techniques in general introduce significant additional complexity of implementation and do not eliminate the crippling need for repeated computation of the estimator on resamples having size comparable to that of the original dataset.

Another landmark in the development of simulation-based inference is subsampling Politis et al. (1999) and the closely related $m$ out of $n$ bootstrap Bickel et al. (1997). These methods (which were introduced to achieve statistical consistency in edge cases in which the bootstrap fails) initially appear to remedy the bootstrap’s key computational shortcoming, as they only require repeated computation of the estimator under consideration on resamples (or subsamples) that can be significantly smaller than the original dataset. However, these procedures also have drawbacks. As we show in our simulation study, their success is sensitive to the choice of resample (or subsample) size (i.e., $m$ in the $m$ out of $n$ bootstrap). Additionally, because the variability of an estimator on a subsample differs from its variability on the full dataset, these procedures must perform a rescaling of their output, and this rescaling requires knowledge and explicit use of the convergence rate of the estimator in question; these methods are thus less automatic and easily deployable than the bootstrap. While schemes have been proposed for data-driven selection of an optimal resample size Bickel and Sakov (2008), they require significantly greater computation which would eliminate any computational gains. Also, there has been work on the $m$ out of $n$ bootstrap that has sought to reduce computational costs using two different values of $m$ in conjunction with extrapolation Bickel and Yahav (1988); Bickel and Sakov (2002). However, these approaches explicitly utilize series expansions of the estimator’s sampling distribution and hence are less automatically usable; they also require execution of the $m$ out of $n$ bootstrap for multiple values of $m$.

Motivated by the need for an automatic, accurate means of assessing estimator quality that is scalable to large datasets, we introduce a new procedure, the Bag of Little Bootstraps (BLB), which functions by combining the results of bootstrapping multiple small subsets of a larger original dataset. Instead of applying an estimator directly to each small subset, as in the $m$ out of $n$ bootstrap and subsampling, the Bag of Little Bootstraps (BLB) applies the bootstrap to each small subset, where in the resampling process of each individual bootstrap run, weighted samples are formed such that the effect is that of sampling the small subset $n$ times with replacement, but the computational cost is that associated with the size of the small subset. This has the effect that, despite operating only on subsets of the original dataset, BLB does not require analytical rescaling of its output. Overall, BLB has a significantly more favorable computational profile than the bootstrap, as it only requires repeated computation of the estimator under consideration on quantities of data that can be much smaller than the original dataset. As a result, BLB is well suited to implementation on modern distributed and parallel computing architectures which are often used to process large datasets. Also, our procedure maintains the bootstrap’s generic applicability, favorable statistical properties (i.e., consistency and higher-order correctness), and simplicity of implementation. Finally, as we show in experiments, BLB is consistently more robust than alternatives such as the $m$ out of $n$ bootstrap and subsampling.

The remainder of our presentation is organized as follows. In Section 2, we formalize our statistical setting and notation, present BLB in detail, and discuss the procedure’s computational characteristics. Subsequently, in Section 3, we elucidate BLB’s statistical properties via a theoretical analysis (Section 3.1) showing that BLB shares the bootstrap’s consistency and higher-order correctness, as well as a simulation study (Section 3.2) which compares BLB to the bootstrap, the $m$ out of $n$ bootstrap, and subsampling. Section 4 discusses a large-scale implementation of BLB on a distributed computing system and presents results illustrating the procedure’s superior computational performance in the massive data setting. We present a method for adaptively selecting BLB’s hyperparameters in Section 5. Finally, we apply BLB (as well as the bootstrap and the $m$ out of $n$ bootstrap, for comparison) to several real datasets in Section 6, we present an extension of BLB to time series data in Section 7, and we conclude in Section 8.

Bag of Little Bootstraps (BLB)

Note that we allow $\xi$ to depend directly on $P$ in addition to $Q_{n}(P)$ because $\xi$ might operate on the distribution of a centered and normalized version of $\hat{\theta}_{n}$. For example, if $\xi$ computes a confidence region, it might manipulate the distribution of the statistic $\sqrt{n}(\hat{\theta}_{n}-\theta)$, which is determined by both $Q_{n}(P)$ and $\theta$; because $\theta$ cannot in general be obtained directly from $Q_{n}(P)$, a direct dependence on $P$ is required in this case. Nonetheless, given knowledge of $Q_{n}(P)$, any direct dependence of $\xi$ on $P$ generally has a simple form, often only involving the parameter $\theta$. Additionally, rather than restricting $Q_{n}(P)$ to be the distribution of $\hat{\theta}_{n}$, we could instead allow it to be the distribution of a more general statistic, such as $(\hat{\theta}_{n},\hat{\sigma}_{n})$, where $\hat{\sigma}_{n}$ is an estimate of the standard deviation of $\hat{\theta}_{n}$ (e.g., this would apply when constructing confidence intervals based on the distribution of the studentized statistic $(\hat{\theta}_{n}-\theta)/\hat{\sigma}_{n}$). Our subsequent development generalizes straightforwardly to this setting, but to simplify the exposition, we will largely assume that $Q_{n}(P)$ is the distribution of $\hat{\theta}_{n}$.

We use $\mathbf{1}_{d}$ to denote the $d$-dimensional vector of ones, and we let $I_{d}$ denote the $d\times d$ identity matrix.

2 Bag of Little Bootstraps

Now, to realize the substantial computational benefits afforded by BLB, we utilize the following crucial fact: each BLB resample, despite having nominal size $n$, contains at most $b$ distinct data points. In particular, to generate each resample, it suffices to draw a vector of counts from an $n$-trial uniform multinomial distribution over $b$ objects. We can then represent each resample by simply maintaining the at most $b$ distinct points present within it, accompanied by corresponding sampled counts (i.e., each resample requires only storage space in $O(b)$). In turn, if the estimator can work directly with this weighted data representation, then the computational requirements of the estimator—with respect to both time and storage space—scale only in $b$, rather than $n$. Fortunately, this property does indeed hold for many if not most commonly used estimators, such as general M-estimators. The resulting BLB algorithm, including Monte Carlo resampling, is shown in Algorithm 1.

Thus, BLB only requires repeated computation on small subsets of the original dataset and avoids the bootstrap’s problematic need for repeated computation of the estimate on resamples having size comparable to that of the original dataset. A simple and standard calculation Efron and Tibshirani (1993) shows that each bootstrap resample contains approximately $0.632n$ distinct points, which is large if $n$ is large. In contrast, as discussed above, each BLB resample contains at most $b$ distinct points, and $b$ can be chosen to be much smaller than $n$ or $0.632n$. For example, we might take $b=n^{\gamma}$ where $\gamma\in[0.5,1]$. More concretely, if $n=1,000,000$, then each bootstrap resample would contain approximately $632,000$ distinct points, whereas with $b=n^{0.6}$ each BLB subsample and resample would contain at most $3,981$ distinct points. If each data point occupies 1 MB of storage space, then the original dataset would occupy 1 TB, a bootstrap resample would occupy approximately 632 GB, and each BLB subsample or resample would occupy at most 4 GB. As a result, the cost of computing the estimate on each BLB resample is generally substantially lower than the cost of computing the estimate on each bootstrap resample, or on the full dataset. Furthermore, as we show in our simulation study and scalability experiments below, BLB typically requires less total computation (across multiple data subsets and resamples) than the bootstrap to reach comparably high accuracy; fairly modest values of $s$ and $r$ suffice.

Due to its much smaller subsample and resample sizes, BLB is also significantly more amenable than the bootstrap to distribution of different subsamples and resamples and their associated computations to independent compute nodes; therefore, BLB allows for simple distributed and parallel implementations, enabling additional large computational gains. In the large data setting, computing a single full-data point estimate often requires simultaneous distributed computation across multiple compute nodes, among which the observed dataset is partitioned. Given the large size of each bootstrap resample, computing the estimate on even a single such resample in turn also requires the use of a comparably large cluster of compute nodes; the bootstrap requires repetition of this computation for multiple resamples. Each computation of the estimate is thus quite costly, and the aggregate computational costs of this repeated distributed computation are quite high (indeed, the computation for each bootstrap resample requires use of an entire cluster of compute nodes and incurs the associated overhead).

In contrast, BLB straightforwardly permits computation on multiple (or even all) subsamples and resamples simultaneously in parallel: because BLB subsamples and resamples can be significantly smaller than the original dataset, they can be transferred to, stored by, and processed on individual (or very small sets of) compute nodes. For example, we could naturally leverage modern hierarchical distributed architectures by distributing subsamples to different compute nodes and subsequently using intra-node parallelism to compute across different resamples generated from the same subsample. Thus, relative to the bootstrap, BLB both decreases the total computational cost of assessing estimator quality and allows more natural use of parallel and distributed computational resources. Moreover, even if only a single compute node is available, BLB allows the following somewhat counterintuitive possibility: even if it is prohibitive to actually compute a point estimate for the full observed data using a single compute node (because the full dataset is large), it may still be possible to efficiently assess such a point estimate’s quality using only a single compute node by processing one subsample (and the associated resamples) at a time.

Statistical Performance

as $n\rightarrow\infty$, for any sequence $b\rightarrow\infty$ and for any fixed $s$.

Moving beyond analysis of the asymptotic consistency of BLB, we now characterize its higher-order correctness (i.e., the rate of convergence of its output to $\xi(Q_{n}(P))$). A great deal of prior work has been devoted to showing that the bootstrap is higher-order correct in many cases (e.g., see the seminal book by Hall (1992)), meaning that it converges to the true value $\xi(Q_{n}(P))$ at a rate of $O_{P}(1/n)$ or faster. In contrast, methods based on analytical asymptotic approximation are generally correct only at order $O_{P}(1/\sqrt{n})$. The bootstrap converges more quickly due to its more data-driven nature, which allows it to better capture finite-sample deviations of the distribution of $\hat{\theta}_{n}$ from its asymptotic limiting distribution.

As shown by the following theorem, BLB shares the same degree of higher-order correctness as the bootstrap, assuming that $s$ and $b$ are chosen to be sufficiently large. Importantly, sufficiently large values of $b$ here can still be significantly smaller than $n$, with $b/n\rightarrow 0$ as $n\rightarrow\infty$. Following prior analyses of the bootstrap, we now make the standard assumption that $\xi$ can be represented via an asymptotic series expansion in powers of $1/\sqrt{n}$. In fact, prior work provides such expansions in a variety of settings. When $\xi$ computes a cdf value, these expansions are termed Edgeworth expansions; if $\xi$ computes a quantile, then the relevant expansions are Cornish-Fisher expansions. See Hall (1992) for a full development of such expansions both in generality as well as for specific forms of the estimator, including smooth functions of mean-like statistics and curve estimators.

Suppose that $\xi(Q_{n}(P))$ admits an expansion as an asymptotic series

where $z$ is a constant independent of $P$ and the $p_{k}$ are polynomials in the moments of $P$. Additionally, assume that the empirical version of $\xi(Q_{n}(P))$ for any $j$ admits a similar expansion

in which case BLB enjoys the same level of higher-order correctness as the bootstrap.

The following result, which applies to the alternative variant of BLB that constrains the $s$ randomly sampled subsets to be disjoint, also highlights the fact that $s$ can grow substantially more slowly than $n$:

Under the assumptions of Theorem 2, and assuming that BLB uses disjoint random subsets of the observed data (rather than simple random subsamples), we have

Therefore, if $s\sim(n/b)$ and $b=\Omega(\sqrt{n})$, then

in which case BLB enjoys the same level of higher-order correctness as the bootstrap.

Finally, while the assumptions of the two preceding theorems generally require that $\xi$ studentizes the estimator under consideration (which involves dividing by an estimate of standard error), similar results hold even if the estimator is not studentized. In particular, not studentizing slows the convergence rate of both the bootstrap and BLB by the same factor, generally causing the loss of a factor of $O_{P}(1/\sqrt{n})$ van der Vaart (1998).

2 Simulation Study

We investigate empirically the statistical performance characteristics of BLB and compare to the statistical performance of existing methods via experiments on simulated data. Use of simulated data is necessary here because it allows knowledge of $P$, $Q_{n}(P)$, and hence $\xi(Q_{n}(P))$; this ground truth is required for evaluation of statistical correctness. For different datasets and estimation tasks, we study the convergence properties of BLB as well as the bootstrap, the $m$ out of $n$ bootstrap, and subsampling.

To evaluate the various quality assessment procedures on a given estimation task and true underlying data distribution $P$, we first compute the ground truth $\xi(Q_{n}(P))$ by generating $2,000$ realizations of datasets of size $n$ from $P$, computing $\hat{\theta}_{n}$ on each, and using this collection of $\hat{\theta}_{n}$’s to form a high-fidelity approximation to $Q_{n}(P)$. Then, for an independent dataset realization of size $n$ from the true underlying distribution, we run each quality assessment procedure (without parallelization) until it converges and record the estimate of $\xi(Q_{n}(P))$ produced after each iteration (e.g., after each bootstrap resample or BLB subsample is processed), as well as the cumulative processing time required to produce that estimate. Every such estimate is evaluated based on the average (across dimensions) relative deviation of its component-wise confidence intervals’ widths from the corresponding true widths; given an estimated confidence interval width $c$ and a true width $c_{o}$, the relative deviation of $c$ from $c_{o}$ is defined as $|c-c_{o}|/c_{o}$. We repeat this process on five independent dataset realizations of size $n$ and average the resulting relative errors and corresponding processing times across these five datasets to obtain a trajectory of relative error versus time for each quality assessment procedure. The relative errors’ variances are small relative to the relevant differences between their means, and so these variances are not shown in our plots. Note that we evaluate based on confidence interval widths, rather than coverage probabilities, to control the running times of our experiments: in our experimental setting, even a single run of a quality assessment procedure requires non-trivial time, and computing coverage probabilities would require a large number of such runs. All experiments in this section were implemented and executed using MATLAB on a single processor. To maintain consistency of notation, we refer to the $m$ out of $n$ bootstrap as the $b$ out of $n$ bootstrap throughout the remainder of this section. For BLB, the $b$ out of $n$ bootstrap, and subsampling, we consider $b=n^{\gamma}$ with $\gamma\in\{0.5,0.6,0.7,0.8,0.9\}$; we use $r=100$ in all runs of BLB.

Computational Scalability

The experiments of the preceding section, though primarily intended to investigate statistical performance, also provide some insight into computational performance: as seen in Figures 1 and 2, when computing on a single processor, BLB generally requires less time, and hence less total computation, than the bootstrap to attain comparably high accuracy. Those results only hint at BLB’s superior ability to scale computationally to large datasets, which we now demonstrate in full in the following discussion and via large-scale experiments on a distributed computing platform.

As discussed in Section 2, modern massive datasets often exceed both the processing and storage capabilities of individual processors or compute nodes, thus necessitating the use of parallel and distributed computing architectures. As a result, the scalability of a quality assessment method is closely tied to its ability to effectively utilize such computing resources.

Recall from our exposition in preceding sections that, due to the large size of bootstrap resamples, the following is the most natural avenue for applying the bootstrap to large-scale data using distributed computing: given data partitioned across a cluster of compute nodes, parallelize the estimate computation on each resample across the cluster, and compute on one resample at a time. This approach, while at least potentially feasible, remains quite problematic. Each computation of the estimate will require the use of an entire cluster of compute nodes, and the bootstrap repeatedly incurs the associated overhead, such as the cost of repeatedly communicating intermediate data among nodes. Additionally, many cluster computing systems currently in widespread use (e.g., Hadoop MapReduce Hadoop (2012)) store data only on disk, rather than in memory, due to physical size constraints (if the dataset size exceeds the amount of available memory) or architectural constraints (e.g., the need for fault tolerance). In that case, the bootstrap incurs the extreme costs associated with repeatedly reading a very large dataset from disk—reads from disk are orders of magnitude slower than reads from memory. Though disk read costs may be acceptable when (slowly) computing only a single full-data point estimate, they easily become prohibitive when computing many estimates on one hundred or more resamples. Furthermore, as we have seen, executing the bootstrap at scale requires implementing the estimator such that it can be run on data distributed over a cluster of compute nodes.

In contrast, BLB permits computation on multiple (or even all) subsamples and resamples simultaneously in parallel, allowing for straightforward distributed and parallel implementations which enable effective scalability and large computational gains. Because BLB subsamples and resamples can be significantly smaller than the original dataset, they can be transferred to, stored by, and processed independently on individual (or very small sets of) compute nodes. For instance, we can distribute subsamples to different compute nodes and subsequently use intra-node parallelism to compute across different resamples generated from the same subsample. Note that generation and distribution of the subsamples requires only a single pass over the full dataset (i.e., only a single read of the full dataset from disk, if it is stored only on disk), after which all required data (i.e., the subsamples) can potentially be stored in memory. Beyond this significant architectural benefit, we also achieve implementation and algorithmic benefits: we do not need to parallelize the estimator internally to take advantage of the available parallelism, as BLB uses this available parallelism to compute on multiple resamples simultaneously, and exposing the estimator to only $b$ rather than $n$ distinct points significantly reduces the computational cost of estimation, particularly if the estimator computation scales super-linearly.

Given the shortcomings of the $m$ out of $n$ bootstrap and subsampling illustrated in the preceding section, we do not include these methods in the scalability experiments of this section. However, it is worth noting that these procedures have a significant computational shortcoming in the setting of large-scale data: the $m$ out of $n$ bootstrap and subsampling require repeated access to many different random subsets of the original dataset (in contrast to the relatively few, potentially disjoint, subsamples required by BLB), and this access can be quite costly when the data is distributed across a cluster of compute nodes.

We compare the performance of BLB and the bootstrap, both implemented as described above. That is, our implementation of BLB processes all subsamples simultaneously in parallel on independent compute nodes; we use $r=50$, $s=5$, and $b=n^{0.7}$. Our implementation of the bootstrap uses all available processors to compute on one resample at a time, with computation of the logistic regression parameter estimates parallelized across the available compute nodes by simply distributing the relevant gradient computations among the different nodes upon which the data is partitioned. We utilize Poisson resampling van der Vaart and Wellner (1996) to generate bootstrap resamples, thereby avoiding the complexity of generating a random multinomial vector of length $n$ in a distributed fashion. Due to high running times, we show results for a single trial of each method, though we have observed little variability in qualitative outcomes during development of these experiments. All experiments in this section are run on Amazon EC2 and implemented in the Scala programming language using the Spark cluster computing framework Zaharia et al. (2012), which provides the ability to either read data from disk (in which case performance is similar to that of Hadoop MapReduce) or cache it in memory across a cluster of compute nodes (provided that sufficient memory is available) for faster repeated access.

In the left plot of Figure 4, we show results obtained using a cluster of 10 worker nodes, each having 6 GB of memory and 8 compute cores; thus, the total memory of the cluster is 60 GB, and the full dataset (150 GB) can only be stored on disk (the available disk space is ample and far exceeds the dataset size). As expected, the time required by the bootstrap to produce even a low-accuracy output is prohibitively high, while BLB provides a high-accuracy output quite quickly, in less than the time required to process even a single bootstrap resample. In the right plot of Figure 4, we show results obtained using a cluster of 20 worker nodes, each having 12 GB of memory and 4 compute cores; thus, the total memory of the cluster is 240 GB, and we cache the full dataset in memory for faster repeated access. Unsurprisingly, the bootstrap’s performance improves significantly with respect to the previous disk-bound experiment. However, the performance of BLB (which also improves), remains substantially better than that of the bootstrap.

Hyperparameter Selection

Like existing resampling-based procedures such as the bootstrap, BLB requires the specification of hyperparameters controlling the number of subsamples and resamples processed. Setting such hyperparameters to be sufficiently large is necessary to ensure good statistical performance; however, setting them to be unnecessarily large results in wasted computation. Prior work on the bootstrap and related procedures—which largely does not address computational issues—generally assumes that a procedure’s user will simply select a priori a large, constant number of resamples to be processed (with the exception of Tibshirani (1985), who does not provide a general solution for this issue). However, this approach reduces the level of automation of these methods and can be quite inefficient in the large data setting, in which each subsample or resample can require a substantial amount of computation.

Thus, we now examine the dependence of BLB’s performance on the choice of $r$ and $s$, with the goal of better understanding their influence and providing guidance toward achieving adaptive methods for their selection. For any particular application of BLB, we seek to select the minimal values of $r$ and $s$ which are sufficiently large to yield good statistical performance.

While these results are useful and provide some guidance for hyperparameter selection, we expect the sufficient values of $r$ and $s$ to change based on the identity of $\xi$ (e.g., we expect a confidence interval to be harder to compute and hence to require larger $r$ than a standard error) and the properties of the underlying data. Thus, to help avoid the need to set $r$ and $s$ to be conservatively and inefficiently large, we now provide a means for adaptive hyperparameter selection, which we validate empirically.

Concretely, to select $r$ adaptively in the inner loop of Algorithm 1, we propose an iterative scheme whereby, for any given subsample $j$, we continue to process resamples and update ${\xi}^{*}_{n,j}$ until it has ceased to change significantly. Noting that the values ${\hat{\theta}}^{*}_{n,k}$ used to compute ${\xi}^{*}_{n,j}$ are conditionally i.i.d. given a subsample, for most forms of $\xi$ the series of computed ${\xi}^{*}_{n,j}$ values will be well behaved and will converge (in many cases at rate $O(1/\sqrt{r})$, though with unknown constant) to a constant target value as more resamples are processed. Therefore, it suffices to process resamples (i.e., to increase $r$) until we are satisfied that ${\xi}^{*}_{n,j}$ has ceased to fluctuate significantly; we propose using Algorithm 2 to assess this convergence. The same scheme can be used to select $s$ adaptively by processing more subsamples (i.e., increasing $s$) until BLB’s output value $s^{-1}\sum_{j=1}^{s}{\xi}^{*}_{n,j}$ has stabilized; in this case, one can simultaneously also choose $r$ adaptively and independently for each subsample. When parallelizing across subsamples and resamples, one can simply process batches of subsamples and resamples (with batch size determined by the available parallelism) until the output stabilizes.

The right plot of Figure 5 shows the results of applying such adaptive hyperparameter selection in a representative empirical setting from our earlier simulation study (without parallelization). For selection of $r$ we use $\epsilon=0.05$ and $w=20$, and for selection of $s$ we use $\epsilon=0.05$ and $w=3$. As illustrated in the plot, the adaptive hyperparameter selection allows BLB to cease computing shortly after it has converged (to low relative error), limiting the amount of unnecessary computation that is performed without degradation of statistical performance. Though selected a priori, $\epsilon$ and $w$ are more intuitively interpretable and less dependent on the details of $\xi$ and the underlying data generating distribution than $r$ and $s$. Indeed, the aforementioned specific values of $\epsilon$ and $w$ yield results of comparably good quality when also used for the other data generation settings considered in Section 3.2, when applied to a variety of real datasets in Section 6 below, and when used in conjunction with different forms of $\xi$ (see the table in Figure 5, which shows that smaller values of $r$ are selected when $\xi$ is easier to compute). Thus, our scheme significantly helps to alleviate the burden of a priori hyperparameter selection.

Automatic selection of a value of $b$ in a computationally efficient manner would also be desirable but is more difficult due to the inability to easily reuse computations performed for different values of $b$. One could consider similarly increasing $b$ from some small value until the output of BLB stabilizes (an approach reminiscent of the method proposed in Bickel and Sakov (2008) for the $m$ out of $n$ bootstrap); devising a means of doing so efficiently is the subject of future work. Nonetheless, based on our fairly extensive empirical investigation, it seems that $b=n^{0.7}$ is a reasonable and effective choice in many situations.

Real Data

In this section, we present the results of applying BLB to several different real datasets. In this case, given the absence of ground truth, it is not possible to objectively evaluate the statistical correctness of any particular estimator quality assessment method; rather, we are reduced to comparing the outputs of various methods (in this case, BLB, the bootstrap, and the $b$ out of $n$ bootstrap) to each other. Because we cannot determine the relative error of each procedure’s output without knowledge of ground truth, we now instead report the average (across dimensions) absolute confidence interval width yielded by each procedure.

Figure 6 shows results for BLB, the bootstrap, and the $b$ out of $n$ bootstrap on the UCI connect4 dataset Frank and Asuncion (2010), where the model is logistic regression (as in the classification setting of our simulation study above), $d=42$, and $n=67,557$. We select the BLB hyperparameters $r$ and $s$ using the adaptive method described in the preceding section. Notably, the outputs of BLB for all values of $b$ considered, and the output of the bootstrap, are tightly clustered around the same value; additionally, as expected, BLB converges more quickly than the bootstrap. However, the values produced by the $b$ out of $n$ bootstrap vary significantly as $b$ changes, thus further highlighting this procedure’s lack of robustness. We have obtained qualitatively similar results on six additional datasets from the UCI dataset repository (ct-slice, magic, millionsong, parkinsons, poker, shuttle) Frank and Asuncion (2010) with different estimators (linear regression and logistic regression) and a range of different values of $n$ and $d$ (see the appendix for plots of these results).

Time Series

To extend BLB in this manner, we must simply alter both the subsample selection mechanism and the resample generation mechanism such that both of these processes respect the underlying data generating process. In particular, for stationary time series data it suffices to select each subsample as a (uniformly) randomly positioned block of length $b$ within the observed time series of length $n$. Given a subsample of size $b$, we generate each resample by applying the stationary bootstrap to the subsample to obtain a series of length $n$. That is, given $p\in$ (a hyperparameter of the stationary bootstrap), we first select uniformly at random a data point in the subsample series and then repeat the following process until we have amassed a new series of length $n$: with probability $1-p$ we append to our resample the next point in the subsample series (wrapping around to the beginning if we reach the end of the subsample series), and with probability $p$ we (uniformly at random) select and append a new point in the subsample series. Given subsamples and resamples generated in this manner, we execute the remainder of the BLB procedure as described in Algorithm 1.

Conclusion

We have presented a new procedure, BLB, which provides a powerful new alternative for automatic, accurate assessment of estimator quality that is well suited to large-scale data and modern parallel and distributed computing architectures. BLB shares the favorable statistical properties (i.e., consistency and higher-order correctness) and generic applicability of the bootstrap, while typically having a markedly better computational profile, as we have demonstrated via large-scale experiments on a distributed computing platform. Additionally, BLB is consistently more robust than the $m$ out of $n$ bootstrap and subsampling to the choice of subset size and does not require the use of analytical corrections. To enhance our procedure’s computational efficiency and render it more automatically usable, we have introduced a means of adaptively selecting its hyperparameters. We have also applied BLB to several real datasets and presented an extension to non-i.i.d. time series data.

A number of open questions and possible extensions remain. Though we have constructed an adaptive hyperparameter selection method based on the properties of the subsampling and resampling processes used in BLB, as well as empirically validated the method, it would be useful to develop a more precise theoretical characterization of its behavior. Additionally, as discussed in Section 5, it would be beneficial to develop a computationally efficient means of adaptively selecting $b$. It may also be possible to further reduce $r$ by using methods that have been proposed for reducing the number of resamples required by the bootstrap Efron (1988); Efron and Tibshirani (1993).

References

Appendix A Appendix: Proofs

We provide here full proofs of the theoretical results included in Section 3.1 above.

Under the assumptions of Theorem 1, for any $j$,

as $n\rightarrow\infty$, for any sequence $b\rightarrow\infty$.

Lemma 2 in conjunction with the continuous mapping theorem van der Vaart (1998) implies the desired result. ∎

A.2 Higher-Order Correctness

Assume that $X_{1},\ldots,X_{b}\sim P$ are i.i.d., and let $\hat{p}_{k}(X_{1},\ldots,X_{b})$ be the sample version of $p_{k}$ based on $X_{1},\ldots,X_{b}$, as defined in Theorem 2. Then, assuming that $E[\hat{p}_{k}(X_{1},\ldots,X_{b})^{2}]<\infty$,

By definition, the $\hat{p}_{k}$ are simply polynomials in sample moments. Thus, we can write

where each $g^{(\beta)}_{\alpha}$ raises its argument to some power. Now, observe that for any $\beta$,

is a V-statistic of order $A_{\beta}$ applied to the $b$ observations $X_{1},\ldots,X_{b}$. Let $h_{\beta}(x_{1},\ldots,x_{A_{\beta}})$ denote the kernel of this V-statistic, which is a symmetrized version of $\prod_{\alpha=1}^{A_{\beta}}g^{(\beta)}_{\alpha}(x_{\alpha})$. It follows that $\hat{p}_{k}=\sum_{\beta=1}^{B}c_{\beta}V_{\beta}$ is itself a V-statistic of order $A=\max_{\beta}A_{\beta}$ with kernel $h(x_{1},\ldots,x_{A})=\sum_{\beta=1}^{B}c_{\beta}h_{\beta}(x_{1},\ldots,x_{A_{\beta}})$, applied to the $b$ observations $X_{1},\ldots,X_{b}$. Let $U$ denote the corresponding U-statistic having kernel $h$. Then, using Proposition 3.5(ii) and Corollary 3.2(i) of Shao (2003), we have

Assume that $X_{1},\ldots,X_{b}\sim P$ are i.i.d., and let $\hat{p}_{k}(X_{1},\ldots,X_{b})$ be the sample version of $p_{k}$ based on $X_{1},\ldots,X_{b}$, as defined in Theorem 2. Then, assuming that $E|\hat{p}_{k}(X_{1},\ldots,X_{b})|<\infty$,

As noted in the proof of Lemma 3, we can write

where each $g^{(\beta)}_{\alpha}$ raises its argument to some power. Similarly,

Given that the number of terms in the outer sum on the right-hand side is constant with respect to $b$, to prove the desired result it is sufficient to show that, for any $\beta$,

If $i_{1},\ldots,i_{A_{\beta}}$ are all distinct, then $E\prod_{\alpha=1}^{A_{\beta}}g^{(\beta)}_{\alpha}(X_{i_{\alpha}})=\prod_{\alpha=1}^{A_{\beta}}Eg^{(\beta)}_{\alpha}(X_{1})$ because $X_{1},\ldots,X_{b}$ are i.i.d.. Additionally, the right-hand summation in (7) has $b!/(b-A_{\beta})!$ terms in which $i_{1},\ldots,i_{A_{\beta}}$ are all distinct; correspondingly, there are $b^{A_{\beta}}-b!/(b-A_{\beta})!$ terms in which $\exists\alpha,\alpha^{\prime}\text{ s.t.\ }i_{\alpha}=i_{\alpha^{\prime}}$. Therefore, it follows that

Note that $C$ is a constant with respect to $b$. Also, simple algebraic manipulation shows that

for some $\kappa>0$. Thus, plugging into equation (8), we obtain the desired result:

We now provide full proofs of Theorem 2, Remark 1, and Theorem 3.

Summing the expansion (3) over $j$, we have

Subtracting the corresponding expansion (2) for $\xi(Q_{n}(P))$, we then obtain

We now further analyze the first two terms on the right-hand side of the above expression; for the remainder of this proof, we assume that $k\in\{1,2\}$. Observe that, for fixed $k$, the $\hat{p}^{(j)}_{k}$ are conditionally i.i.d. given $X_{1},\ldots,X_{n}$ for all $j$, and so

Combining the expressions in the previous three panels, we find that

Finally, plugging into equation (9) with $k=1$ and $k=2$, we obtain the desired result. ∎

By Lemmas 3 and 4, $\operatorname{Var}(\hat{p}^{(1)}_{k}-p_{k})=O(1/b)$ and $E[\hat{p}^{(1)}_{k}-p_{k}]^{2}=O(1/b^{2})$. Combining with the expressions in the previous three panels, we obtain the desired result:

Additionally, from the result of Lemma 4, we have

Combining the expressions in the previous two panels, we find that

Finally, plugging into equation (10) with $k=1$ and $k=2$, we obtain the desired result. ∎

Appendix B Appendix: Additional Real Data Results