High-dimensional Data Bootstrap

High-dimensional Data Bootstrap

\address[V. Chernozhukov]{ Department of Economics and Center for Statistics and Data Science, MIT. } \email{[email removed]}

\address[D. Chetverikov]{ Department of Economics, UCLA.} \email{[email removed]}

\address[K. Kato]{ Department of Statistics and Data Science, Cornell University. } \email{[email removed]}

\address[Y. Koike]{ Mathematics and Informatics Center and Graduate School of Mathematical Sciences, The University of Tokyo.} \email{[email removed]}

Introduction

Overview

The bootstrap is a generic method to estimate the sampling distribution of a statistic, typically by resampling one's own data. Since the seminal work of Efron1979, there has been a substantial amount of research that explores the theoretical properties of the bootstrap. In classical settings where the data dimension is fixed, the bootstrap often yields more accurate confidence intervals or tests than those based on first-order asymptotic approximations. Also, the bootstrap provides a practical method of inference for models where the analytical estimation of the asymptotic distribution of an estimator or a test statistic is difficult, such as in quantile regression or semiparametric models.

Recently, there has been growing interest in extending the scope of the bootstrap to high-dimensional (a.k.a. “$p \gg n$") settings. Due to the advancement of science and technology, data sets in which the number of features (e.g., genes) well exceeds the sample size (e.g., the number of patients) have become common in many application domains. Analysis of such high-dimensional data has been a major focus in statistics in the last two or three decades; see, e.g., BuhlmannvandeGeer2011, Giraud2014, HastieTibshiraniWainwright2015, Wainwright2018 as textbook references on high-dimensional statistics. Yet, developing high-dimensional inference methods with provable accuracy under weak conditions is challenging, since classical statistical theory (which presumes fixed data dimensions) does not carry over to high dimensions, at least directly.

CCK2013AoS contributed to the above goal and established the consistency of the Gaussian multiplier (or wild) and empirical bootstraps for the maximum of the sum of independent high-dimensional random vectors. Notably, CCK2013AoS allow the data dimension $p$ to be much larger than the sample size $n$. CCK2017AoP extended the results of CCK2013AoS and proved that, for $S_{n} = \sum_{i=1}^{n}X_{i}/\sqrt{n}$ the scaled sum of independent centered random vectors in $\mathbb{R}^{p}$ with $p = p_{n} \to \infty$, it holds that

under moment conditions, where $\mathcal{R}$ is the class of rectangles in $\mathbb{R}^p$ and $S_n^B$ is either the Gaussian multiplier or empirical bootstrap statistic. Importantly, the error bound ((ref)) only requires $p$ to be $\log p = o(n^{1/7})$, which allows $p \gg n$, for the bootstrap to be consistent. Subsequent works have explored refinements and extensions of these results, including deng2020beyond,lopes2020bootstrapping,kuchibhotla2021high,koike2021notes,chernozhukov2019improved,fang2021high,lopes2020central,kuchibhotla2020high,chernozhukov2020nearly.

These extensive theoretical developments stimulated many new applications of bootstrap to high-dimensional or nonparametric statistics, including simultaneous inference for high-dimensional (or infinite-dimensional) parameters CCK2013AoS,CCK2014AoS,CCK2014AoSdensity,wasserman2014,Belloni2015,chen2015,chen2016,ZhangCheng2017JASA,Dezeure2017,chang2017simulation, ning2017general,belloni2018uniformly,rinaldo2019, kuchibhotla2020valid, testing for shape restrictions chetverikov2019testing, detection of spurious correlations fan2018, comparison of large covariance matrices chang2017, inference for partially identified models chetverikov2018adaptive,CCK2018RES, goodness-of-fit testing jankova2020goodness, error estimation of matrix randomized algorithms Lopes2019JMLR, testing the mean function for functional data lopes2020bootstrapping, and many more.

This article aims to provide a brief overview of the current literature on high-dimensional bootstrap. In Section (ref), we review classical asymptotics for the empirical bootstrap and discuss challenges in inference for high-dimensional data. In Section (ref), we review basic high-dimensional central limit theorems (CLTs) and bootstrap consistency results in high dimensions. In Section (ref), we discuss selected applications of high-dimensional bootstrap. Section (ref) leaves some concluding remarks.

Finally, while this review aims to convey the main ideas and techniques in high-dimensional bootstrap, the literature is now quite broad and has seen rapid expansion, and thus this review is by no means exhaustive. For instance, due to space limitation, we omit several topics such as the empirical process extension CCK2014AoS,chernozhukov2016empirical, high-dimensional CLTs and bootstrap for dependent data ZhangWu2017, ZhangCheng2017, CCK2018RES, chang2021central,chiang2021inference,kurisu2021gaussian and $U$-statistics Chen2018,chen2019randomized,chen2020jackknife,song2019approximating,song2020stratified.

Notation

We use $N(\mu,\Sigma)$ to denote the normal distribution with mean vector $\mu$ and covariance matrix $\Sigma$. Abusing the notation, we also use $N(\mu,\Sigma)$ to denote a random vector following $N(\mu,\Sigma)$. Let $\stackrel{d}{\to}$ and $\stackrel{P}{\to}$ denote convergence in distribution and convergence in probability, respectively. Let $\| \cdot \|$ and $\| \cdot \|_{\infty}$ denote the Euclidean and $\ell^\infty$-norms for vectors, respectively, i.e., for $x = (x_1,\dots,x_p)^T$, $\| x \| = \sqrt{\sum_{j=1}^p x_j^2}$ and $\| x \|_{\infty} = \max_{1 \leqslant j \leqslant p} |x_j|$. The notation $\delta_x$ denotes the Dirac delta at point $x$. For two numbers $a,b \in \mathbb{R}$, we use the notation $a \vee b = \max \{ a,b \}$. We use $\lesssim$ to denote inequalities up to numerical constants.

Challenges in Inference for High-dimensional Data

Classical Asymptotics

We first review how the bootstrap works in the classical setting where the data dimension is fixed while the sample size $n$ tends to infinity. Consider first inference on the mean parameter $\mu$ from independent and identically distributed (i.i.d.) univariate ($p=1$) random variables $X_1,\dots,X_n$. As long as the common distribution has finite second moment, the CLT yields that, for the sample mean $\overline{X}_n = n^{-1}\sum_{i=1}^n X_i$, \[ \sqrt{n}(\overline{X}_n - \mu) \stackrel{d}{\to} N(0,\sigma^2), \] where $\sigma^2 > 0$ is the population variance. We shall apply the empirical bootstrap to estimate the limit distribution $N(0,\sigma^2)$ and hence the sampling distribution of $\sqrt{n}(\overline{X}_n-\mu)$. To this end, conditionally on the data $X_1,\dots,X_n$, generate an independent sample $X_1^B,\dots,X_n^B$ from the empirical distribution $P_n = n^{-1}\sum_{i=1}^n \delta_{X_i}$ of $X_1,\dots,X_n$. Since, conditionally on the data, $P_n$ has mean $\overline{X}_n$ and variance $\widehat{\sigma}^2 = n^{-1}\sum_{i=1}^n (X_i-\overline{X}_n)^2$, one can expect that, conditionally on the data, \[ \sqrt{n}(\overline{X}_n^B - \overline{X}_n) \stackrel{d}{\approx} N(0,\widehat{\sigma}^2) \stackrel{d}{\approx} N(0,\sigma^2), \] where $\overline{X}_n^B = n^{-1}\sum_{i=1}^n X_i^B$. This suggests that the conditional distribution of $\sqrt{n}(\overline{X}_n^B - \overline{X}_n)$ given the data approximates the $N(0,\sigma^2)$ distribution. Indeed, since $N(0,\sigma^2)$ has a continuous distribution function, it is not difficult to show that \[ \sup_{t \in \mathbb{R}} \Big | \mathbb{P}^B \big(\sqrt{n}(\overline{X}_n^B - \overline{X}_n) \leqslant t\big) - \mathbb{P} \big(N(0,\sigma^2) \leqslant t\big) \Big | \stackrel{P}{\to} 0. \] Here and in what follows, $\mathbb{P}^B$ denotes the conditional probability given the data. This result implies that the conditional distribution of $\sqrt{n}(\overline{X}_n^B - \overline{X}_n)$ given the data consistently estimates the sampling distribution of $\sqrt{n}(\overline{X}_n-\mu)$.

Suppose, more generally, that we want to make inference on a scalar parameter $\theta^{\star}$, for which a suitable estimator $\widehat{\theta}_n$ (based on i.i.d. data $X_1,\dots,X_n$) is available. Suppose also that $\widehat{\theta}_n$ admits an asymptotically linear expansion of the form

where $\psi$ is an influence function such that $\psi (X_i)$ has mean zero. This expansion implies that $ \sqrt{n}(\widehat{\theta}_n - \theta^{\star}) \stackrel{d}{\to} N(0,\sigma^2_\psi), $ where $\sigma^2_{\psi}$ is the variance of $\psi (X_1)$ (assume $\sigma^2_{\psi} > 0$). To estimate the $N(0,\sigma^2_\psi)$ distribution or the sampling distribution of $\sqrt{n}(\widehat{\theta}_n - \theta^{\star})$, we may apply bootstrap to $\widehat{\theta}_n$ by replacing the data $X_1,\dots,X_n$ with the bootstrap sample $X_1^B,\dots,X_n^B$, which gives rise to the bootstrap version of the estimator $\widehat{\theta}_n^B$. Then, under regularity conditions, we can expect that the bootstrap estimator $\widehat{\theta}_n^B$ admits a similar asymptotic linear expansion as

Subtracting ((ref)) from ((ref)), we have \[ \sqrt{n}(\widehat{\theta}_n^B - \widehat{\theta}_n) = \frac{1}{\sqrt{n}} \sum_{i=1}^n (\psi (X_i^B) - \overline{\psi}_n) + o_{P}(1), \] where $\overline{\psi}_n = n^{-1}\sum_{i=1}^n \psi (X_i)$. Thus, as in the sample mean case, we obtain the bootstrap consistency \[ \sup_{t \in \mathbb{R}} \Big | \mathbb{P}^B \big(\sqrt{n}(\widehat{\theta}_n^B - \widehat{\theta}_n) \leqslant t\big) - \mathbb{P} \big(N(0,\sigma_{\psi}^2) \leqslant t\big) \Big | \stackrel{P}{\to} 0. \] See, e.g., wellner1996bootstrapping.

The bootstrap consistency result yields methods to construct (asymptotically valid) confidence intervals for $\theta^{\star}$. For example, let $\widehat{q}_\alpha$ denote the conditional $\alpha$-quantile of $\widehat{\theta}_n^B$ given the data, i.e., $\widehat{q}_\alpha = \inf\{ t : \mathbb{P}^B\big( \widehat{\theta}_n^B \leqslant t\big) \geqslant \alpha \}$. Then, from the bootstrap consistency result, it is not difficult to verify that the following percentile confidence interval, $[ 2\widehat{\theta}_n - \widehat{q}_{1-\alpha/2}, 2\widehat{\theta}_n - \widehat{q}_{\alpha/2} ]$, contains $\theta^{\star}$ with probability $1-\alpha+o(1)$ (cf. Lemma 23.3 in van2000asymptotic). We refer the reader to Hall1992,VW1996,horowitz2001bootstrap for several other classical asymptotic results on the bootstrap, including higher-order refinements and empirical process extension.

Challenges in High Dimensions

The discussion in the preceding section immediately extends to the case where $\overline{X}_n$ or $\widehat{\theta}_n$ are multivariate, as long as the dimension is fixed. However, in modern statistical applications, we are often interested in situations where the data dimension $p$ is comparable to or even much larger than the sample size $n$. Mathematically, such scenarios can be captured by allowing $p$ to depend on $n$ and considering the case where $p_n \to \infty$ as $n \to \infty$.

The main issue in high-dimensional inference problems is the lack of explicit limit distributions. To fix ideas, let $X_1,\dots,X_n$ be i.i.d. $p$-dimensional random vector with mean $\mu$ and covariance matrix $\Sigma$. Even when $p = p_n \to \infty$, one can still expect that $S_n = \sqrt{n}(\overline{X}_n - \mu)$ be approximated by $N(0,\Sigma)$, but since the dimension of $N(0,\Sigma)$ is $p = p_n \to \infty$, there will be no obvious “final” limit distribution for $S_n$. Of course, we can think of $S_n$ as random elements in $\mathbb{R}^\mathbb{N}$ (the countable product of $\mathbb{R}$), but weak convergence in $\mathbb{R}^{\mathbb{N}}$ is equivalent to finite dimensional convergence, which is too weak a result to develop genuinely high-dimensional inference methods.

Also, often, we are interested in approximating the sampling distribution of a certain functional of $S_n$, such as the $\ell^\infty$-norm, i.e., $\| S_n \|_{\infty} = \max_{1 \leqslant j \leqslant p} | S_{n,j}|$. In some situations, it could be possible to derive (after a proper normalization) limit distributions to such one-dimensional functionals; however, required regularity conditions tend to be restrictive. For instance, the $\ell^\infty$-functional $\| S_n\|_{\infty}$ may converge to the type I extreme value distribution after a proper normalization, but the derivation of such results typically requires the coordinates of the random vectors to be weakly dependent; see, e.g., Section 2.2 in chang2017 for relevant discussion. Further, in practical situations, we often deal with high-dimensional statistics that are approximately linear (like ((ref))), but how fast the linearization error should vanish is a nontrivial question because of the lack of explicit limit distributions.

Another challenge is that a large number of features induce a potentially complex dependence structure among the features. Failing to take into account the dependence structure among the features would result in (sometimes severely) conservative inference methods, especially in high dimensions. For instance, consider testing the global hypothesis $H_0: \mu_1=\cdots=\mu_p$ using $\| S_n \|_{\infty}$. Assume each coordinate of $S_n$ has unit variance (for simplicity) and apply a \v{S}idak-type correction, which yields a critical value $t_\alpha$ given by $\mathbb{P} (|Z| > t_\alpha) = 1-(1-\alpha)^{1/p}$ for $Z \sim N(0,1)$. This is the critical value computed as if the coordinates were independent and each coordinate of $S_n$ were distributed as $N(0,1)$ fan2007many. However, in the extreme situation where the coordinates are perfectly correlated, then the actual rejection probability under the null is, using the classical Berry-Esseen theorem, $\mathbb{P}(\| S \|_{\infty} > t_\alpha) = \mathbb{P}(|Z| > t_{\alpha}) + O(n^{-1/2}) = 1-(1-\alpha)^{1/p} + O(n^{-1/2}) = -p^{-1}\log(1-\alpha) + O(p^{-2} + n^{-1/2})$, which is very close to $0$ when $p$ is large, thus yielding conservative tests. Of course, the perfect correlation is an extreme case, but the above discussion indicates that if $S_n$ contains a large number of strongly correlated components, then Bonferoni or \v{S}idak-type corrections would yield overly conservative tests.

We will review below several techniques (high-dimensional CLTs, Gaussian anticoncentration inequality, Gaussian comparison) that help overcome the lack of explicit limit distributions and justify bootstrap methods for high-dimensional data. Also, we will demonstrate that the bootstrap is able to automatically take into account the dependence among the coordinates and yields asymptotically exact inference methods even in high dimensions. It should be noted that the bootstrap allows for arbitrary correlations between the coordinates, thereby accommodating a broader class of high-dimensional inference tasks.

Bootstrap Methods in High Dimensions

A key step toward establishing the validity of bootstrap methods in high dimensions is a high-dimensional CLT. The high-dimensional CLT amounts to bounding a certain discrepancy measure between the sampling distribution of the scaled sample mean and the corresponding Gaussian distribution in such a way that the bound permits the dimension to increase with $n$. One natural approach in this direction is to compare the distributions over a (sufficiently large) collection of subsets $\mathcal{A}$ in $\mathbb{R}^p$ and try to derive as explicit bounds as possible on $\sup_{A \in \mathcal{A}} | \mathbb{P} (S_n \in A) - \mathbb{P} (N(0,\Sigma) \in A)|$, where $S_n$ is the scaled sample mean of mean-zero independent random vectors in $\mathbb{R}^p$ and $N(0,\Sigma)$ is the corresponding Gaussian distribution. As CCK2013AoS,CCK2017AoP observed, choosing $\mathcal{A}$ to be the collection of rectangles allows us to establish error bounds that depend on the dimension $p$ only through $\log p$, thereby permitting $p \gg n$. Such high-dimensional CLTs and accompanying bootstrap consistency results, together with techniques developed therein, paved the way for constructing inference methods with theoretical guarantees for many high-dimensional inference problems under relatively mild regularity conditions. This section reviews high-dimensional CLTs over the rectangles and bootstrap consistency results in high dimensions, following the recent work of chernozhukov2019improved that improves earlier results of CCK2013AoS,CCK2017AoP. We also briefly review relevant developments in the literature. Throughout this and the following sections, we always assume $n\geq2$ and $p\geq2$.

High-Dimensional CLTs

Let $X_1,\dots,X_n$ be independent (but not necessarily identically distributed) random vectors with dimension $p$. Assume that $X_i$'s have mean zero (otherwise, work with $X_i - \mathbb{E}[X_i]$ instead of $X_i$). Consider the scaled sample mean \[ S_n = \frac{1}{\sqrt{n}} \sum_{i=1}^n X_i. \] Let $\overline{\sigma}, \underline{\sigma}$ be given positive constants such that $\underline{\sigma} \leqslant \overline{\sigma}$, and let $B_n \geqslant 1$ be a sequence of constants that may diverge as $n \to \infty$. Let $\Sigma = \mathbb{E}[S_nS_n^T] = n^{-1}\sum_{i=1}^n \mathbb{E}[X_iX_i^T]$. Also, let $\mathcal{R}$ denote the collection of closed rectangles in $\mathbb{R}^p$, \[ \mathcal{R} = \Big \{ \prod_{j=1}^p [a_j,b_j] : -\infty \leqslant a_j \leqslant b_j \leqslant \infty, \ j=1,\dots,p \Big \}. \] Finally, for notational convenience, let \[ \delta_{1,n} = \left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4} \quad \text{and} \quad \delta_{2,n}^{[q]} =\sqrt{ \frac{B_n^2(\log (pn))^{3-2/q}}{n^{1-2/q}}} \quad \text{for $q>2$}. \] We first present a high-dimensional CLT over the rectangles under a sub-exponential condition on the coordinates.

The theorem follows from Lemma 4.3 in chernozhukov2019improved. The assumption of the theorem is satisfied if, for example, $|X_{ij}| \leqslant B_n$ almost surely for all $(i,j)$ and $\underline{\sigma}^2 \leqslant n^{-1}\sum_{i=1}^n \mathbb{E}[X_{ij}^2] \leqslant \overline{\sigma}^2$ for all $j$. Notably, the above theorem does not impose any restrictions on the correlation structure between the coordinates of the random vectors, so $\Sigma$ is permitted to be singular. This aspect is relevant in applications to high-dimensional data since, in the presence of many features, some of them are likely to have strong correlations, causing $\Sigma$ to be (or close to being) singular. For example, data from microarray or transcriptome experiments contain genes that are divided into groups with varying sizes according to functionalities, and genes from the same group tend to have relatively strong (sometimes very strong) within-correlations (cf. chang2017).

Theorem (ref) shows that, under the assumption of the theorem,

provided that $B_n^2 \log^5 (pn) = o(n)$, which allows $p$ to be much larger than $n$.

Theorem (ref) requires that the coordinates of $X_i$ be sub-exponential. The following theorem, adapted from Theorem 2.5 in chernozhukov2019improved, presents a version of the high-dimensional CLT under polynomial moment conditions.

\if0

\fi

Theorem (ref) covers the following scenario relevant to regression applications: $X_i = \epsilon_i z_i$ where $\epsilon_i$ is a univariate “error” term while $z_i \in \mathbb{R}^p$ is a vector of fixed “covariates”. In this case, $\mathbb{E}[\| X_i \|_{\infty}^q] \leqslant \| z_i \|_{\infty}^q \mathbb{E}[|\epsilon_i|^q]$, so if the covariates are uniformly bounded and the $q$-th moments of the error terms are bounded, then $B_n = O(1)$. Notably this only requires $\epsilon_i$ to have $q=2+ \delta$ bounded moments.

In practice, statistics of interest may not be exactly sample means but be only approximately linear. Since, in high dimensions, the approximating Gaussian distribution changes with $n$, the linearization error has to vanish at a sufficiently fast rate for such approximate sample means to satisfy the high-dimensional CLT. Fortunately, for the rectangle case, the required condition on the linearization error is rather mild, as shown next.

The lemma follows from the following Gaussian anticoncentration inequality for the rectangles, due (essentially) to nazarov2003maximal; see also klivans2008learning. A self-contained proof can be found in chernozhukov2017detailed.

Nazarov's inequality states that the Gaussian maximum $W^{\vee}$ possesses a mass at most $O(\delta \sqrt{\log p})$ around the $\delta$-neighborhood of any point $t$. This is an opposite statement to the Gaussian concentration inequality, which asserts that the Gaussian maximum possesses a large mass around its expectation or median (cf. boucheron2013concentration). As such, Nazarov's inequality is an instance of an anticoncentration inequality for the Gaussian distribution. The $\sqrt{\log p}$ dependence in Proposition (ref) is sharp in the sense that the opposite inequality (up to multiplicative constants) holds when $W_1,\dots,W_p$ are i.i.d. as $N(0,1)$; see CCK2015PTRF. Nazarov's inequality also plays a crucial role in establishing high-dimensional CLTs over the rectangles.

To see how Lemma (ref) follows from Nazarov's inequality, let us pick a particular rectangle $R = \{ x \in \mathbb{R}^p: \max_{1 \leqslant j \leqslant p} x_j \leqslant t \}$ for some $t \in \mathbb{R}$ (the general case follows by translation and scaling; cf. the proof of Corollary 5.1 in CCK2017AoP). Then, for $\widehat{S}_n^{\vee} = \max_{1 \leqslant j \leqslant p}\widehat{S}_{n,j}$ and $W \sim N(0,\Sigma)$, \[

\] where $o(1)$ is uniform in $t$. Likewise, we have \[ \mathbb{P} \big(\widehat{S}_n \in R\big) = \mathbb{P} \big(\widehat{S}_n^{\vee} \leqslant t\big) \geqslant \mathbb{P} \big(W^{\vee} \leqslant t\big) - O\big(\epsilon \sqrt{\log p}\big) - \mathbb{P} (\| \mathsf{R}_n \|_{\infty} > \epsilon) - o(1). \] Since $\| \mathsf{R}_n \|_{\infty} = o_P(1/\sqrt{\log p})$, we can choose $\epsilon = \epsilon_n = o(1/\sqrt{\log p})$ in such a way that $\mathbb{P} (\| \mathsf{R}_n \|_{\infty} > \epsilon) = o(1)$, which yields that $\sup_{t \in \mathbb{R}} |\mathbb{P} (S_n^{\vee} \leqslant t) - \mathbb{P} (W^{\vee} \leqslant t)| = o(1)$.

Multiplier and Empirical Bootstraps

The high-dimensional CLTs discussed in the preceding section show that the sampling distribution of $S_n$ can be approximated by the Gaussian distribution $N(0,\Sigma)$ uniformly over the rectangles, even when $p \gg n$. However, in practice, the approximating Gaussian distribution $N(0,\Sigma)$ is infeasible since the covariance matrix $\Sigma$ is unknown. We shall consider further estimating the $N(0,\Sigma)$ distribution by using bootstrap methods. In this section, for simplicity, we shall focus our attention on the Gaussian multiplier and empirical bootstraps. We will briefly discuss other bootstrap methods at the end of this section. We first formally define the Gaussian multiplier and empirical bootstraps. Recall $\overline{X}_n = n^{-1}\sum_{i=1}^nX_i$.

In either case, we approximate $N(0,\Sigma)$ or the sampling distribution of $S_n$ by the conditional distribution of $S_{n}^{B}$.

Consider first the Gaussian multiplier bootstrap. Observe that, conditionally on the data, we have $S_{n}^{B} \sim N(0,\widehat{\Sigma})$ with $\widehat{\Sigma} = n^{-1} \sum_{i=1}^n (X_i-\overline{X}_n)(X_i - \overline{X}_n)^T$. Thus, for the Gaussian multiplier bootstrap, the problem reduces to comparing two Gaussian distributions with covariance matrices $\widehat{\Sigma}$ and $\Sigma$. The following Gaussian comparison inequality, taken from Proposition 2.1 in chernozhukov2019improved, implies that the Gaussian multiplier bootstrap is consistent over the rectangles provided that $$\max_{1 \leqslant j ,k \leqslant p}| \widehat{\Sigma}_{jk} - \Sigma_{jk} | = o_{P}(1/\log^2 p),$$ which can hold under mild moment conditions.

As for the empirical bootstrap, in view of the fact that the empirical distribution $P_n$ has mean $\overline{X}_n$ and covariance matrix $\widehat{\Sigma}$, we may first approximate the conditional distribution of $S_n^{B}$ by $N(0,\widehat{\Sigma})$ over the rectangles by applying the high-dimensional CLT, and then use the Gaussian comparison inequality to further approximate $N(0,\widehat{\Sigma})$ by $N(0,\Sigma)$. The following theorem presents finite sample error bounds for the Gaussian multiplier and empirical bootstraps under a sub-exponential condition. Recall that $\mathbb{P}^B$ denotes the conditional probability given $X_1,\dots,X_n$.

The theorem follows from Lemmas 4.5 and 4.6 in chernozhukov2019improved. Precisely speaking, the actual proof of Theorem (ref) for the empirical bootstrap exploits the fact that the empirical bootstrap matches approximately higher-order moments of $S_n$. See chernozhukov2019improved for details.

As in the high-dimensional CLT, the error bound for the bootstrap depends on $p$ only through $\log p$, which allows $p$ to be much larger than $n$ for the bootstrap to be consistent. Also, the theorem imposes no restrictions on the correlation structure between the coordinates of the random vectors, allowing $\Sigma$ to be singular.

Like Theorem (ref), the following theorem, adapted from Theorem 2.6 in chernozhukov2019improved, covers polynomial moment conditions.

Figure (ref) illustrates the finite sample performance of the Gaussian, Gaussian multiplier bootstrap, and empirical bootstrap approximations for $\| S_n \|_{\infty}$ under a regression setup, namely compares $\mathbb{P}(\|S_n\|_\infty \leqslant x)$, $\mathbb{P}(\| N(0,\Sigma) \|_{\infty} \leqslant x)$, and $\mathbb{P}^B(\| S_n^B \|_\infty \leqslant x)$. The figure indicates that both Gaussian and bootstrap approximations are reasonably good, especially in the tails.

In applications, we often normalize the coordinates of the sample mean by estimates of the standard deviations, so that each coordinate is approximately distributed as $N(0,1)$. We may estimate the variance $\Sigma_{jj}$ of the $j$-th coordinate of $S_n$ by the sample variance $\widehat{\Sigma}_{jj} = n^{-1}\sum_{i=1}^n (X_{ij} - \overline{X}_{n,j})^2$. We shall then consider approximating the sampling distribution of the normalized sample mean $\widehat{\Lambda}^{-1/2}S_n$ by the conditional distribution of $\widehat{\Lambda}^{-1/2}S_n^B$, where $\widehat{\Lambda} = \mathrm{diag} \{ \widehat{\Sigma}_{11},\dots,\widehat{\Sigma}_{pp} \}$. Let $\Sigma_0 = \Lambda^{-1/2}\Sigma \Lambda^{-1/2}$ denote the correlation matrix of $S_n$. Consider the assumption of Theorem (ref) with $B_n^2 \log^5 (pn) = o(n)$. Then, we have $\max_{1 \leqslant j \leqslant p} |\Sigma_{jj}/\widehat{\Sigma}_{jj}-1| = o_{P}(1/\log^2p)$, so that combining the high-dimensional CLT and Gaussian concentration, it holds that $\| (\widehat{\Lambda}^{-1/2}-\Lambda^{-1/2})S_n\|_{\infty} = o_{P}(1/\sqrt{\log p})$. Thus, arguing as in the proof of Lemma (ref), we have \[ \sup_{R \in \mathcal{R}} \left|\mathbb{P}\big(\widehat{\Lambda}^{-1/2}S_n \in R \big) - \mathbb{P} \big(N(0,\Sigma_0) \in R\big) \right| \to 0. \] Likewise, for either the Gaussian multiplier or empirical bootstrap, it holds that \[ \sup_{R \in \mathcal{R}} \left|\mathbb{P}^B\big(\widehat{\Lambda}^{-1/2}S_n^B \in R \big) - \mathbb{P} \big(N(0,\Sigma_0) \in R\big) \right| \stackrel{P}{\to} 0. \] Similar results hold under polynomial moment conditions. See Appendix A.2 in chen2019randomized for relevant arguments.

The above theoretical results demonstrate that the bootstrap can adequately capture the dependence between the coordinates, thereby yielding asymptotically exact coverage or size control even in high dimensions. Together with the fact that the bootstrap allows for arbitrary correlations between the coordinates, the bootstrap is a particularly powerful inferential tool in high dimensions.

Near $n^{-1/2}$ Rates under Structural Assumptions on Covariance Matrix

Recall that, in $p=1$, the classical Berry-Esseen theorem shows that, for $X,X_1,\dots,X_n$ i.i.d. univariate $p=1$ random variables with mean zero and unit variance (for simplicity), $\sup_{t \in \mathbb{R}} |\mathbb{P} (S_n \leqslant t) - \mathbb{P}(N(0,1) \leqslant t)| \lesssim \mathbb{E}[|X|^3]/\sqrt{n}$, and this bound is known to be sharp in terms of dependence on $n$. Given this, there has been great interest in deriving high-dimensional CLTs and bootstraps over the rectangles that achieve (near) $n^{-1/2}$ rates while allowing for $p \gg n$. Recent progress shows that such (near) $n^{-1/2}$ rates are possible under structural assumptions on $\Sigma$. In the following discussion, we assume that $X,X_1,\dots,X_n$ are i.i.d. with mean zero and covariance matrix $\Sigma$. Recall $S_n = n^{-1/2}\sum_{i=1}^n X_i$.

lopes2020bootstrapping derived a first such result under the variance decay condition. Namely, assuming $\max_{1 \leqslant j \leqslant p}\operatorname{Var}(X_{j}) = O(j^{-a})$ for arbitrarily small $a >0$ and other technical conditions, they derive error bounds of order $n^{-1/2+\delta}$ with arbitrarily small $\delta > 0$ for $\sup_{A \in \mathcal{A}}|\mathbb{P} (S_n \in A) - \mathbb{P} (N(0,\Sigma) \in A)|$ for a certain subclass $\mathcal{A}$ of the rectangles. lopes2020bootstrapping also derive similar error bounds for the Gaussian multiplier bootstrap.

In a different direction, if we assume that the smallest eigenvalue of $\Sigma$ is bounded away from zero, as in fang2021high, then it was shown by chernozhukov2020nearly that under regularity conditions: $$\sup_{R \in \mathcal{R}} |\mathbb{P} (S_n \in R) - \mathbb{P}(N(0,\Sigma) \in R)| \leq C \left ( \frac{B_n^2 (\log p)^3}{n}\right )^{1/2} \log n, $$ and a similar result was obtained for bootstrap approximation. This result builds on and refines a sequence of other important results in this direction obtained by fang2021high, lopes2020central, kuchibhotla2020high. Under the nondegeneracy of the covariance matrix $\Sigma$, the above bound implies that the high-dimensional CLT holds if $\log p = o(n^{1/3})$ up to logarithmic factors, which weakens the previous requirement that $\log p = o(n^{1/5})$. See also das2021central who investigate necessary conditions for high-dimensional CLTs.

Selcted Applications

This section reviews applications of high-dimensional bootstrap to several inference tasks. Specifically, we discuss penalty choice for the Lasso, simultaneous confidence intervals for high-dimensional parameters, estimation and inference for maximum/minumum effects, comparing large covariance matrices, and large-scale multiple testing. There are many other applications of high-dimensional bootstrap, some of which are delineated in the introduction, and the reader is also encouraged to look at the references therein.

Penalty Choice for High-Dimensional Regression

Consider a high-dimensional regression with non-Gaussian errors \[ y_i =x_i^T\beta^{\star} + \epsilon_{i}, \quad \mathbb{E}[\epsilon_i]=0, \quad i=1,\dots,n, \] where $y_i$ is a scalar response variable, $x_i$ is a $p$-dimensional vector of (fixed) covariates, and $\epsilon_i$ is an error term. Here the dimension $p$ of the covariate vector $x_i$ can be much larger than the sample size $n$, $p \gg n$, but we assume that the model is sparse in the sense that the number of nonzero components of $\beta^{\star}$ is $s = \# \{ j : \beta_j^\star \ne 0 \} \ll n$.

Arguably, one of the most popular estimates for such high-dimensional linear regression is the Lasso tibshirani1996regression, which is defined by \[ \widehat{\beta} \in \operatorname*{arg\,min}_{\beta \in \mathbb{R}^p} \left [ \frac{1}{n} \sum_{i=1}^n (y_i - x_i^{T}\beta )^2 + \lambda \sum_{j=1}^p | \beta_{j} | \right], \] where $\lambda \geqslant 0$ is a tuning parameter. It is well known that the statistical performance of the Lasso crucially relies on the choice of the tuning parameter $\lambda$. From a seminal work of bickel2009simultaneous, for a given confidence level $\alpha \in (0,1)$, if we choose $\lambda$ in such a way that \[ \lambda = (1-\alpha)\text{-quantile of } S(\beta^\star), \quad S(\beta^\star) = 2\cdot \max_{1\leq j \leq p} \left | n^{-1} {\textstyle \sum}_{i=1}^n x_{ij} \epsilon_i \right |, \] then it holds that $\| \widehat{\beta} - \beta^{\star} \|_{2,n} \lesssim \lambda \sqrt{s}$ with probability at least $1-\alpha$, provided that the restricted eigenvalue condition is satisfied; see bickel2009simultaneous and also belloni2013least. Here $\| \delta \|_{2,n} =\sqrt{ n^{-1}\sum_{i=1}^n (x_i^T \delta)^2 }$.

If $\epsilon_1,\dots,\epsilon_n$ are i.i.d. sub-Gaussian, then the above choice of $\lambda$ is of order $O(\sqrt{\log p/n})$ (which follows by using a standard concentration inequality) and thus $\| \widehat{\beta} - \beta^{\star} \|_{2,n} \lesssim \sqrt{(s\log p)/n}$. However, if the error distribution has heavier tails than sub-Gaussian, then bounding $\lambda$ by concentration inequalities would lead to sub-optimal rates, and the bound itself contains distribution-dependent parameters that are unknown in practice. Instead, we may approximate or estimate $\lambda$ by applying the high-dimensional CLT (homoscedastic case) or using the multiplier bootstrap (heteroscedastic case).

Consider first the case where the error terms are homoscedastic, $ \sigma^2 = \mathbb{E}[\epsilon_1^2]=\dots=\mathbb{E}[\epsilon_n^2]. $ Here we assume that $\epsilon_1,\dots,\epsilon_n$ are independent. In this case, we can directly apply the high-dimensional CLT to approximate the quantile of $S(\beta^\star)$ by that of $2\sigma \cdot \max_{1\leq j \leq p} | n^{-1}\sum_{i=1}^n x_{ij} \xi_i |$, where $\xi_1,\dots,\xi_n$ are i.i.d. $N(0,1)$. Thus, under homoscedasticity, we can approximate $\lambda$ by

In practice, $\sigma$ is unknown but can be consistently estimated by pre-estimating $\beta^{\star}$ by the Lasso with a crude-choice of the $\lambda$-parameter.

In general, if the error terms are heteroscedastic, then we first pre-estimate $\beta^{\star}$ (again by using the Lasso with a crude-choice of the $\lambda$-parameter) to construct estimates $\widehat{\epsilon}_i$ of $\epsilon_i$, and apply the multiplier bootstrap to estimate the quantile of $\max_{1\leq j \leq p} | n^{-1} \sum_{i=1}^n x_{ij} \epsilon_i |$, namely,

Under regularity conditions, it is shown that the Lasso with this data-driven choice of $\lambda$ satisfies $\| \widehat{\beta} - \beta^{\star} \|_{2,n} \lesssim \sqrt{(s \log p)/n}$ with high probability. See Section 4 of CCK2013AoS for related results. Note that for the justification of both methods ((ref)) and ((ref)), we may apply Theorems (ref) and (ref), which only require that the error terms have finite polynomial moments.

The function rlasso in the R package hdm (chernozhukov2016high) implements the above methods of choosing the penalty parameter.\footnote{ The options homoscedastic=TRUE and X.dependent.lambda=TRUE in the argument penalty implement ((ref)), while homoscedastic=FALSE (default) and \texttt{X.dependent.lambda=TRUE} implement ((ref)), where the default value of $\alpha$ is set to $\alpha=0.1$} The function \texttt{rlasso} also offers a joint test of significance of variables -- the test of the null hypothesis that $\beta^\star = 0$ -- based on the sup-score statistic $ S(\beta^\star)$ constrained under the null hypothesis $\beta^\star = 0$; the critical value for such test is $\widehat \lambda$.

Simultaneous Confidence Intervals for High-Dimensional Parameters

Modern statistical and machine learning problems often entail estimation and inference for a large number of parameters, the number of which may exceed the sample size. In such cases, researchers are interested in conducting inference for not only individual parameters but also groups of parameters simultaneously. Methods of uniform inference for high-dimensional parameters are also a basis of large-scale multiple testing; cf. Section (ref).

In this section, we consider constructing simultaneous confidence intervals (rectangles) for a high-dimensional parameter $\theta^{\star} \in \mathbb{R}^p$. In many settings, e.g., those in examples discussed in below, we have an estimator $\widehat{\theta}_n$ for $\theta^{\star}$ such that it admits an asymptotically linear expansion of the form

where $\psi_1,\dots,\psi_n$ are independent random vectors (influence functions) with mean zero and $\mathsf{R}_n$ is a remainder term such that $\| \sqrt{n} \mathsf{R}_n \|_{\infty} = o_{P}(1/\sqrt{\log p})$. Applying the high-dimensional CLT to the leading term $n^{-1}\sum_{i=1}^n \psi_i$ yields that (cf. Lemma (ref)), for $\Sigma = n^{-1} \sum_{i=1}^n \mathbb{E}[\psi_i \psi_i^T]$, \[ \sup_{R \in \mathcal{R}} \Big |\mathbb{P} \big(\sqrt{n}(\widehat{\theta}_n - \theta^{\star}) \in R\big) - \mathbb{P} \big(N(0,\Sigma) \in R\big) \Big | \to 0, \] even when $p \gg n$, provided that certain moment conditions on $\psi_1,\dots,\psi_n$ are satisfied.

We shall estimate the $N(0,\Sigma)$ distribution by using the multiplier bootstrap. In practice, influence functions $\psi_1,\dots,\psi_n$ may be unknown, so we assume that there are suitable estimates $\widehat{\psi}_1,\dots,\widehat{\psi}_n$ for the influence functions. We then apply the multiplier bootstrap to the estimated influence functions, i.e., \[ S_n^B = \frac{1}{\sqrt{n}} \sum_{i=1}^n \xi_i (\widehat{\psi}_i - \overline{\widehat{\psi}}), \] where $ \overline{\widehat{\psi}} = n^{-1}\sum_{i=1}^n \widehat{\psi}_i$. Let $\widehat{\Sigma} = n^{-1}\sum_{i=1}^n (\widehat{\psi}_{i} - \overline{\widehat{\psi}})(\widehat{\psi}_{i} - \overline{\widehat{\psi}})^T$. In view of the Gaussian comparison inequality (Proposition (ref)), the Gaussian multiplier bootstrap is consistent over the rectangles, $ \sup_{R \in \mathcal{R}} | \mathbb{P}^B (S_n^B \in R) - \mathbb{P} (N(0,\Sigma) \in R) | \stackrel{P}{\to} 0, $ provided that $\max_{1 \leqslant j,k \leqslant p} | \widehat{\Sigma}_{jk} - \Sigma_{jk}| = o_{P}(1/\log^2 p)$, which can hold even when $p \gg n$. Alternatively, we may apply the empirical bootstrap by generating an independent sample $\widehat{\psi}_1^B,\dots,\widehat{\psi}_n^B$ from the empirical distribution $n^{-1}\sum_{i=1}^n \delta_{\widehat{\psi}_i}$ and construct $S_n^B = n^{-1/2}\sum_{i=1}^n (\widehat{\psi}_i^B - \overline{\widehat{\psi}})$.

Also, the high-dimensional CLT and bootstrap consistency hold for the normalized statistics $\widehat{\Lambda}^{-1/2}\sqrt{n}(\widehat{\theta}_n - \theta^{\star})$ and $\widehat{\Lambda}^{-1/2}S_n^B$, under regularity conditions, where $\widehat{\Lambda} = \mathrm{diag} \{ \widehat{\Sigma}_{11},\dots,\widehat{\Sigma}_{pp} \}$ with $\widehat{\Sigma}_{jj} = n^{-1}\sum_{i=1}^n (\widehat{\psi}_{ij} - \overline{\widehat{\psi}}_{j})^2$ (see the discussion at the end of Section (ref)), namely,

Here $\Sigma_0 = \Lambda^{-1/2} \Sigma \Lambda^{-1/2}$ with $\Lambda = \mathrm{diag} \{ \Sigma_{11},\dots,\Sigma_{pp} \}$, and $p$ is allowed to increase with $n$, $p=p_n \to \infty$ (we assume that the diagonal elements $\Sigma_{jj}$ are bounded away from zero). For a given $\alpha \in (0,1)$, let \[ \widehat{q}_{1-\alpha} = \text{conditional $(1-\alpha)$-quantile of $\| \widehat{\Lambda}^{-1/2}S_n^{B} \|_{\infty}$}. \] We claim that the rectangle of the form \[ \widehat{R}_{1-\alpha}= \prod_{j=1}^p \left [\widehat{\theta}_{n,j} \pm \widehat{\Sigma}_{jj}^{1/2}\widehat{q}_{1-\alpha}/\sqrt{n} \right ] =: \prod_{j=1}^p \widehat{\mathrm{CI}}_{j,1-\alpha} \] contains $\theta^{\star}$ with probability $1-\alpha + o(1)$. Indeed, $\theta^{\star} \in \widehat{R}_{1-\alpha}$ if and only if $\| \widehat{\Lambda}^{-1/2}\sqrt{n}(\widehat{\theta}_n - \theta^{\star}) \|_{\infty} \leqslant \widehat{q}_{1-\alpha}$. Let $q_{1-\alpha}$ denote the $(1-\alpha)$-quantile of $\| N(0,\Sigma_0) \|_{\infty}$. Then, from the bootstrap consistency result ((ref)), there exists a sequence $\epsilon_n \to 0$ such that $q_{1-\alpha-\epsilon_n} \leqslant \widehat{q}_{1-\alpha} \leqslant q_{1-\alpha+\epsilon_n}$ with probability approaching one (cf. Theorem 2.5 in belloni2018high). Combining the high-dimensional CLT from ((ref)) and the fact that $\| N(0,\Sigma_0) \|_{\infty}$ has a continuous distribution function (cf. Proposition (ref)), we see that \[ \mathbb{P} (\| \widehat{\Lambda}^{-1/2}\sqrt{n}(\widehat{\theta}_n - \theta^{\star}) \|_{\infty} \leqslant \widehat{q}_{1-\alpha}) \leqslant 1-\alpha+ \epsilon_n + o(1) = 1-\alpha + o(1). \] The reverse inequality follows similarly. Conclude that $\widehat{R}_{1-\alpha}$ is a valid simultaneous confidence interval for $\theta^{\star}$ with level $1-\alpha+o(1)$, i.e., $\mathbb{P} (\theta^{\star} \in \widehat{R}_{1-\alpha}) = 1-\alpha+o(1)$.

Estimation and Inference for Maximum/Minimum Effects

In certain applications, researchers are interested in the maximum or minimum of high-dimensional parameters. We start with discussing two such examples.

Consider the setting of Section (ref) and estimation of $\vartheta^{\star} = \max_{1 \leqslant j \leqslant p}\theta_j^\star$ (the minimum effect can be dealt with analogously, as $\min_{1 \leqslant j \leqslant p}\theta_j^\star = - \max_{1 \leqslant j \leqslant p} (-\theta_j^\star)$). As observed in chernozhukov2013intersection and guo2021inference, the plug-in estimator $\max_{1 \leqslant j \leqslant p}\widehat{\theta}_{n,j}$ tends to be upward biased in the finite sample. To tackle this issue, chernozhukov2013intersection proposed the following precision corrected estimator \[ \widehat{\vartheta}_n (\alpha) = \max_{1 \leqslant j \leqslant p} \big [ \widehat{\theta}_{n,j} - \widehat k_{1-\alpha} \widehat{\Sigma}^{1/2}_{jj}/\sqrt{n} \big ], \] which can make the estimator upward $(1-\alpha)$-quantile unbiased, for example median upward unbiased for $\alpha=1/2$, and $\widehat k_{1-\alpha}$ is the bootstrap estimate of $(1-\alpha)$-quantile of the studentized maximum estimation error $ T= \sqrt{n} \max_{1 \leqslant j \leqslant p} (\widehat \theta_{n,j} - \theta_j^{\star})/\widehat{\Sigma}^{1/2}_{jj}.$ Indeed, \[

\] under regularity conditions like those in Section (ref). Thus, the precision correction guarantees that the estimator is upward biased with probability at most $\alpha+o(1)$. Also, $[\widehat{\vartheta}_n(\alpha),\infty)$ is an asymptotically valid one-sided confidence interval for $\vartheta^*$ with level $1-\alpha+o(1)$. The approach can be refined by a conservative pre-estimation of the argmax set $J_0 = \arg \max_{1 \leq j \leq p}\theta^*_j$, and then working with such set in place of $\{1,\dots, p\}$.

Comparing Large Covariance Matrices

In modern genomics, understanding how the dependencies among many genes vary between different biological states (e.g., healthy or with disease) has received significant attention. Formally, the problem amounts to comparing large covariance matrices across different populations. Let $X$ and $Y$ be two random vectors in $\mathbb{R}^p$ with covariance matrices $$ \Sigma_1 = \Big (\sigma_{1,jk} \Big )_{1 \leqslant j,k \leqslant p} \ \ \mathrm{ and } \ \ \Sigma_2 = \Big (\sigma_{2,jk} \Big )_{1 \leqslant j,k \leqslant p}, $$ respectively, and consider testing the hypothesis $H_0: \Sigma_1 = \Sigma_2$ against the alternative $H_1: \Sigma_1 \ne \Sigma_2$ in high dimensions.

Let $X_1,\dots,X_n$ and $Y_1,\dots,Y_m$ be independent observations from $X$ and $Y$, respectively, and let $\widehat{\Sigma}_1 = (\widehat{\sigma}_{1,jk})_{1 \leqslant j,k \leqslant p}$ and $\widehat{\Sigma}_2 = (\widehat{\sigma}_{2,jk})_{1 \leqslant j,k \leqslant p}$ be the corresponding empirical covariance matrices, respectively (e.g., $\widehat{\Sigma}_1 = n^{-1}\sum_{i=1}^n (X_i - \overline{X}_n)(X_i - \overline{X}_n)^T$ with $\overline{X}_n = n^{-1}\sum_{i=1}^n X_i$). chang2017 consider the max-type test statistic $\widehat{T}_{\max} = \max_{1 \leqslant j \leqslant k \leqslant p} |\widehat{t}_{jk}|$, where \[ \widehat{t}_{jk} = \frac{\widehat{\sigma}_{1,jk}-\widehat{\sigma}_{2,jk}}{\sqrt{n^{-1}\widehat{s}_{1,jk}+m^{-1}\widehat{s}_{2,jk}}}. \] Here $\widehat{s}_{1,jk} = n^{-1}\sum_{i=1}^n \{ (X_{ij} -\overline{X}_{n,j})(X_{ik} -\overline{X}_{n,k}) - \widehat{\sigma}_{1,jk} \}^2$ and $\widehat{s}_{2,jk}$ is defined analogously.

To calibrate critical values for the test, chang2017 propose the following Gaussian multiplier bootstrap procedure:

Assuming that $m$ is of comparative size as $n$, chang2017 show that the test has an asymptotic size $\alpha$ even if $p \gg n$, namely, they show that $\mathbb{P} (\widehat{T}_{\max}^\dagger > \widehat{q}_{1-\alpha}) = \alpha+o(1)$ if $H_0$ holds true. (From the results reviewed in this article, the improved condition on $p$ is that $\log p = o(n^{1/5})$.) Further, chang2017 combine this test with the Benjamini–Hochberg procedure to develop a gene clustering algorithm.

Large-Scale Multiple Testing

Suppose that we have a large collection of null $H_j$ and alternative $H_j'$ hypotheses for $j=1,\dots,p$. Such large-scale multiple testing problems commonly appear in biological applications. Here, we are interested in testing these hypotheses simultaneously for all $j=1,\dots,p$, so we aim to construct a testing procedure that would reject at least one true null hypothesis with probability at most $\alpha+o(1)$, uniformly over the set of true null hypotheses. Procedures with this property are said to control the Family-Wise Error Rate (FWER).

To this end, one can adopt the step-down procedure of romano2005. Specifically, consider the setting of Section (ref) and the testing problems $H_j: \theta_j^\star \leqslant 0$ against $H_j': \theta_j^\star > 0$ for $j=1,\dots,p$ (the equality case $\theta_j^\star = 0$ can be dealt with by splitting $\theta_j^\star = 0$ into two hypotheses $\theta_j^\star \leqslant 0$ and $-\theta_j^\star \leqslant 0$). For each $j=1,\dots,p$, let $\widehat{t}_j = \sqrt{n}\widehat{\theta}_{n,j}/\widehat{\Sigma}_{jj}^{1/2}$ be the studentized test statistic for $H_j$ against $H_j'$.

For a subset $w\subset \{1,\dots,p\}$, let $c_{1-\alpha,w}$ be the bootstrap estimate of the $(1-\alpha)$-quantile of $\max_{j\in w}\sqrt{n}(\widehat{\theta}_{n,j}-\theta_j^\star)/\widehat{\Sigma}_{jj}^{1/2}$ by either applying the Gaussian multiplier or empirical bootstrap (cf. Section (ref)). On the first step, let $w(1)=\{1,\dots,p\}$. Reject all hypotheses $H_{j}$ satisfying $\widehat{t}_{j}>c_{1-\alpha,w(1)}$. If no null hypothesis is rejected, then stop. If some $H_{j}$ are rejected, let $w(2)$ be the set of all null hypotheses that were not rejected in the first step. In step $l\geq2$, let $w(l)\subset \{1,\dots,p \}$ be the subset of null hypotheses that were not rejected up to step $l$. Reject all hypotheses $H_{j}$, $j\in w(l)$, satisfying $\widehat{t}_{j}>c_{1-\alpha,w(l)}$. If no null hypothesis is rejected, then stop. If some $H_{j}$ are rejected, let $w(l+1)$ be the subset of all null hypotheses among $j\in w(l)$ that were not rejected. Proceed in this way until the algorithm stops.

CCK2013AoS show that this step-down procedure can achieve the FWER control under regularity conditions, even when $p \gg n$, extending the analysis of romano2005 to high dimensions. See Section 5 in CCK2013AoS and also Section 2.4 in belloni2018high for more details.

The following R code illustrates implementation of Romano--Wolf's step-down procedure for multiple one-sample $t$-tests. We compute the adjusted $p$-values following RoWo16. The dataset Fund is taken from the ISRL2 package associated with the textbook JWHT21. The dataset contains $n=50$ rows and $p=2000$ columns, and each column corresponds to the returns of a hedge fund manager. Writing $\mu_j$ for the $j$-th fund manager's mean return, we test for $H_j:\mu_j\leq0$ against $H_j':\mu_j>0$ for all $j=1,\dots,2000$. See Section 13.3 of JWHT21 for more illustration of this dataset.

Finally, the R package hdm offers functionality on multiple hypothesis testing in high-dimensional approximately sparse linear regression models based upon Romano--Wolf step down procedures (see bach2018valid for documentation).

Concluding Remarks

The field of high-dimensional bootstrap has seen a rapid development, and this article reviewed the main ideas and key techniques used there. Notably, the bootstrap offers the following advantages to inference for high-dimensional data:

The theoretical development stimulated many new applications of bootstrap to high-dimensional inference tasks, some of which we reviewed in this article.

We end this article with commenting on a couple of future research topics. First, Section (ref) presents error bounds for the Gaussian multiplier and empirical bootstraps, and the same error bounds hold for other wild bootstraps such as Mammen's bootstrap. As observed numerically in deng2020beyond and chernozhukov2019improved, however, the empirical and Mammen's bootstraps perform (slightly) better than the Gaussian multiplier bootstrap in practice. Deriving error bounds that certify a specific bootstrap to be preferred over others is an interesting direction in future research. Also, while there is considerable progress on extending high-dimensional bootstrap for (temporally or graph-) dependent data, the derived error bounds in the dependent case are substantially slower than the independent case. More research is needed to explore sharp error bounds for high-dimensional bootstrap for dependent data. Finally, the literature on high-dimensional bootstrap has been focused on the case where the approximating distribution is Gaussian. However, several important statistics such as degenerate $U$-statistics have non-Gaussian approximating distributions. With an exception of koike2019mixed, extending the scope of high-dimensional bootstrap to non-Gaussian approximating distributions is open.