Bootstrapping $\ell_p$-Statistics in High Dimensions

Bootstrapping $_p$-Statistics in High Dimensions

\quad { Keywords: Bootstrap; high-dimensional inference; Berry-Esseen bound; anti-concen-

tration; Gaussian approximation; Gaussian comparison inequality.}

Introduction

Let $X=\{X_i\}_{i=1}^n$ be a random sample of independent and centered random vectors in $\mathbb{R}^d$, where dimension $d = d_n$ may grow with sample size $n$. Consider the re-scaled sum

and define the $\ell_p$-statistic $T_{n,p}$ by

This paper is concerned with developing a bootstrap procedure to estimate the distribution of $\ell_p$-statistics when the dimension $d$ exceeds the sample size $n$. This distribution is of interest in many statistical applications. In particular, $\ell_2$-statistic $T_{n,2}$ and maximum statistic $T_{n,\infty}$ are frequently applied to a broad spectrum of statistical problems such as testing of multiple means, construction of simultaneous confidence regions, and model selection bai1996effect,chen2010two,fan2015power.

In low dimensions, when the dimension $d$ is fixed, the asymptotic properties of $\ell_p$-statistics are well-understood: If the data are i.i.d. with finite second moments, the central limit theorem (CLT) applied to the re-scaled sum $S_n^X$ and the continuous mapping theorem guarantee that $T_{n,p}\overset{d}{\rightarrow} \|Z\|_p$, where $Z \sim N\left(0, \mathrm{E}\left[X_1X_1'\right]\right)$. Thus, the limiting distribution of $T_{n,p}$ depends on the data only through the first two moments. Closed-form expressions of the limiting distribution of $\ell_p$-statistics remain somewhat elusive, but for $T_{n,2}$ and $T_{n,\infty}$ tractable characterizations exist under additional assumptions on the covariance structure.

The situation is very different in high dimensions. If the dimension $d$ grows faster than $\sqrt{n}$, the classical CLT does no longer apply to the re-scaled sum $S_n^X$. So, to approximate the distribution of $T_{n,p}$ one has to target directly the scalar random variable $\|S_n^X\|_p$. Since $\|S_n^X\|_p$ is a highly non-linear function of the random sample $X$, this calls for a non-parametric approach. In this direction, chernozhukov2013GaussianApproxVec, chernozhukov2015ComparisonAnti, chernozhukov2017CLTHighDim have made important progress by developing a non-parametric multiplier bootstrap procedure to approximate the distribution of the maximum statistic $T_{n,\infty}$. In this paper, we further develop this line of research. While there exist specialized results for high-dimensional sum-of-squares type $T_{n,2}^2$-statistics bai1996effect,bentkus2003DependenceBerryEsseen, chen2010two, fan2015power, pouzo2015bootstrap, xu2019L2Asymptotics, a unified investigation on the weak convergence of general $\ell_p$-statistics remains highly challenging due to the lack of smoothness of the $\ell_p$-norm. In this sense, our work solves a long-standing open problem initiated by the aforementioned pioneering work.

The primary methodological contribution of this paper is a bootstrap procedure for $\ell_p$-statistics with $p \in [1, \infty]$. Our bootstrap procedure draws inspiration from above observation that in low dimensions the limiting distribution of $\ell_p$-statistics depends only on the first two moments of the data. Specifically, the proposed algorithm involves sampling bootstrap data from a Gaussian distribution that is parameterized by an estimate of the covariance matrix. The algorithm works with any estimate of the covariance matrix; it is easy to implement and very versatile. In particular, it can be combined with estimates of the covariance matrix that leverage special structures such as low rank, (approximate) sparsity or bandedness. The algorithm is best understood as a hybridization of non-parametric and parametric bootstrap, and we call it the Gaussian parametric bootstrap.

A secondary methodological contribution is a bootstrap hypothesis test for testing many linear restrictions on high-dimensional mean vectors. This hypothesis test is based on the Gaussian parametric bootstrap for $\ell_p$-statistics and it is asymptotically correct and consistent under certain high-dimensional alternatives. We give precise recommendations on how to choose the exponent $p \in [1, \infty]$ based on characteristics of the random sample (tails and covariance structure) and to maximize the power for given alternative hypotheses. For small exponents $p$, the test is useful when the goal is to identify significant subsets from a large collection of means, e.g. sets of genes in micro-array and genetic sequence studies. Whereas for large exponents $p$, the test is powerful when the purpose is to detect significant singletons, e.g. anomaly detection in materials science and medical imaging.

The two main theoretical contributions of this paper are a non-asymptotic characterization of the sampling distribution of $\ell_p$-statistics $T_{n,p}$ in high dimensions and the consistency of the Gaussian parametric bootstrap. The non-asymptotic characterization is based on a Gaussian approximation, i.e. a proxy statistic constructed from Gaussian random vectors. The quality of the Gaussian and bootstrap approximation improves as the sample size $n$ increases and shows a subtle interplay between dimension $d$, exponent $p$, and the tail distribution of the data. Among other things, we demonstrate that if the data has light tails the approximation errors vanish for $ \log d = o(n)$ and all $p \in [1, \infty]$; whereas if the data is heavy-tailed with at most $s \geq 4$ finite moments the approximation errors are negligible for $d \log d= o(n^{s/4})$ and all $p \in [1, s]$. These theoretical results provide a comprehensive view on the asymptotic distribution theory of $\ell_p$-statistics and are relevant in guiding practitioners in choosing between different $\ell_p$-statistics given the properties of the random sample at hand. Qualitatively, our numerical experiments lend further support to these theoretical findings.

Establishing the Gaussian approximation and consistency of the bootstrap is non-trivial and we develop a significant amount of new technical tools. The following three technical results are of interest beyond the scope of this paper: First, we derive an abstract Berry-Esseen-type CLT for $\ell_p$-statistics in high dimensions, which extends and improves the known Berry-Esseen-type CLTs for $p=2$ bentkus2003DependenceBerryEsseen and $p = \infty$ chernozhukov2017CLTHighDim. Second, we establish an anti-concentration inequality for $\ell_p$-norms of random vectors with log-concave probability measure. For $p \in \{2, \infty\}$ this inequality is sharper than related inequalities by goetze2019LargeBall and chernozhukov2017Nazarov. Third, we develop a Gaussian comparison inequality to compare the distributions of $\ell_p$-norms of different Gaussian random vectors in Kolmogorov-Smirnov distance. For $p = \infty$ this inequality improves the corresponding result in chernozhukov2015ComparisonAnti.

Organization. The paper is organized as follows. We introduce the Gaussian parametric bootstrap in Section (ref) and present our main theoretical results on the Gaussian approximation of $\ell_p$-statistics and the consistency of the Gaussian parametric bootstrap in Section (ref). We develop applications to testing high-dimensional mean vectors in Section (ref) and report results from several numerical experiments in Section (ref). In Appendix (ref) we discuss technical results, including the abstract Berry-Esseen-type CLT, the anti-concentration inequalities for $\ell_p$-statistics, and the new Gaussian comparison theorems. Appendix (ref) contains proofs to all our results.

Notation. For non-negative real-valued sequences $\{a_n\}_{n\geq 1}$ and $\{b_n\}_{n \geq 1}$, the relation $a_n \lesssim b_n$ means that there exists an absolute constant $c > 0$ independent of $n, d, p$ and an integer $n_0 \in \mathbb{N}$ such that $a_n \leq c b_n$ for all $n \geq n_0$. We write $a_n \asymp b_n$ if $a_n \lesssim b_n$ and $b_n \lesssim a_n$. We define $a_n \vee b_n = \max\{a_n, b_n\}$ and $a_n \wedge b_n = \min\{a_n, b_n\}$. For a vector $a \in \mathbb{R}^d$ and $p \in [1, \infty)$ we write $\|a\|_p = (\sum_{k=1}^d |a_k|^p)^{1/p}$. Also, we write $\|a\|_\infty = \max_{1 \leq k \leq d} |a_k|$. For a scalar random variable $\xi$ and $\alpha \in (0, 2]$ we define the $\psi_\alpha$-Orlicz norm by $\|\xi \|_{\psi_\alpha} = \inf\{t > 0 : \mathrm{E}[\exp(|\xi|^\alpha/t^\alpha)] \leq 2\}$. For a sequence of scalar random variables $\{\xi_n\}_{n \geq 1}$ we write $\xi_n = O_p(a_n)$ if $\xi_n/a_n$ is stochastically bounded. For any symmetric real-valued matrix $M \in \mathbb{R}^{d \times d}$ we denote its largest and smallest eigenvalue by $\lambda_{\max}(M)$ and $\lambda_{\min}(M)$, respectively. We denote its operator norm by $\|M\|_{op}$ (its largest singular value) and $\|M\|_{2 \rightarrow p} = \sup_{\|u\|_2 \leq 1} \|Mu\|_p$. We write $M \succeq 0$ to indicate that $M$ is positive semi-definite. For any convex body $K \subset \mathbb{R}^d$ we write $\mathrm{Vol}(K) = \int_K d \lambda^d$, where $\lambda^d$ is the Lebesgue measure in $d$ dimensions.

Methodology

We introduce the new Gaussian parametric bootstrap for $\ell_p$-statistics and discuss its relation to the non-parametric Gaussian multiplier bootstrap.

Gaussian parametric bootstrap

Let $X=\{X_i\}_{i=1}^n$ be a random sample of independent and centered random vectors. The Gaussian parametric bootstrap algorithm requires as input a consistent and positive semi-definite estimate $\widehat{\Sigma}_n$ of the (averaged) population covariance matrix,

We will discuss candidates for $\widehat{\Sigma}_n$ in subsequent sections. Let $V^X \mid X \sim N(0, \widehat{\Sigma}_n)$ and define the Gaussian parametric bootstrap estimate of the $\ell_p$-statistic $T_{n,p}$ by

The rationale for this bootstrap statistic is easiest to understand in low dimensions: If dimension $d$ is fixed and the data $X= \{X_i\}_{i=1}^n$ is i.i.d. with finite second moments, the CLT and the continuous mapping theorem imply that $T_{n,p}\overset{d}{\rightarrow} \|Z\|_p$, where $Z \sim N\left(0, \mathrm{E}\left[X_1X_1'\right]\right)$. Hence, in this scenario, the bootstrap statistic $T_{n,p}^*$ is just the parametric bootstrap estimate of the limiting random variable $\|Z\|_p$. Of course, if $d \geq \sqrt{n}$ and the data is non-identically distributed, the CLT does not apply and the limiting random variable $Z$ needs not to exist. The gist of the theoretical results in Sections (ref) and (ref) is that we do not need the CLT to hold for the distributions of $T_{n,p}^*$ and $T_{n,p}$ to be close. For $\ell_p$-statistics this result is new, but it is in line with similar results on linear regression functions, empirical processes in infinite-dimensional Banach spaces, as well as maximum and spectral statistics in high dimensions bickel1983bootstrapping, radulovic1998can, chernozhukov2013GaussianApproxVec, roellin2013stein, lopes2019BootstrappingSpectral.

Relation to the Gaussian multiplier bootstrap

The Gaussian multiplier bootstrap was first proposed by chernozhukov2013GaussianApproxVec in the context of the maximum statistic $T_{n,\infty}$. It is a special case of the wild bootstrap method wu1986jacknife, liu1988bootstrap, mammen1993bootstrap and its adaptation to general $\ell_p$-statistics $T_{n,p}$ is straightforward:

Let $g= \{g_i\}_{i=1}^n$ be a sequence of i.i.d. standard normal random variables independent of the random sample $X = \{X_i\}_{i=1}^n$. The Gaussian multiplier bootstrap algorithm builds on the centered random sample $X_1 - \bar{X}_n, \ldots, X_n - \bar{X}_n$, where $\bar{X}_n := n^{-1}\sum_{i=1}^n X_i$. We set

and define the Gaussian multiplier bootstrap estimate of the $\ell_p$-statistic $T_{n,p}$ by

Since $S_n^{gX} \mid X \sim N(0, \widehat{\Sigma}_{\mathrm{naive}})$ with $\widehat{\Sigma}_{\mathrm{naive}} = n^{-1}\sum_{i=1}^n(X_i - \bar{X}_n)(X_i - \bar{X}_n)'$, the Gaussian multiplier bootstrap statistic is in fact equivalent to a Gaussian parametric bootstrap statistic based on the sample covariance matrix $\widehat{\Sigma}_{\mathrm{naive}}$. The key advantage of the Gaussian parametric over the Gaussian multiplier bootstrap is that it allows for more refined estimates of the population covariance matrix $\Sigma_n$ that leverage additional structure such as low-rank, (approximate) sparsity, and bandedness. This is particularly important in high dimensions where the sample covariance matrix $\widehat{\Sigma}_{\mathrm{naive}}$ is a poor estimate of the population covariance matrix.

Theoretical analysis

We present a non-asymptotic characterization of $\ell_p$-statistics via Gaussian approximation and establish the consistency of the Gaussian parametric bootstrap procedure.

Assumptions

Unless otherwise stated, $X=\{X_i\}_{i=1}^n$ denotes a random sample of independent and centered random vectors in dimension $d$, where $d = d_n$ grows with the sample size $n$. We analyze the theoretical properties of $\ell_p$-statistics and the Gaussian parametric bootstrap under the following three different assumptions on the tails of random vectors.

Assumption (ref) is one of the many equivalent definitions of sub-Gaussian random vectors antonini1997subgaussian, vershynin2018HighDimProb. This specific formulation is useful for applications in high-dimensional statistics because $\mathrm{E}\left[(u'X_i)^2\right]^{1/2} \leq \|\mathrm{E}[X_iX_i']\|_{op}\|u\|_2$. Hence, we can easily incorporate characteristics of the covariance matrix such as sparsity, bandedness, low-rank, etc. Assumptions (ref) and (ref) relax and generalize Assumption (ref) in an obvious way. Most importantly, if $X$ satisfy Assumption (ref) for all $s \geq 1$ and with $K_s = \sqrt{s}$ ($K_s = s$) then $X$ is sub-Gaussian (sub-Exponential) and also satisfy Assumption (ref) (Assumption (ref)).

Gaussian approximation

In this section we show that the distribution of the $\ell_p$-statistic $T_{n,p}$ can be approximated by the distribution of a proxy statistic based on Gaussian random vectors. This result rationalizes the Gaussian parametric bootstrap procedure in high dimensions. It is also relevant for establishing bootstrap consistency in the next section.

Let $Z =\{Z_i\}_{i=1}^n$ be a sequence of independent multivariate Gaussian random vectors $Z_i \sim N(0,\mathrm{E}[X_iX_i'])$ which are independent of $X = \{X_i\}_{i=1}^n$. We define the Gaussian proxy statistic of the $\ell_p$-statistic $T_{n,p}$ as

To state the Gaussian approximation result we need to define the following additional quantities: the rank of the (averaged) covariance matrices of the $X_i$'s,

the smallest and largest variances of the $X_i$'s,

and the largest ratio of the variances of the $X_i$'s,

Our first theorem shows that the distribution of $\widetilde{T}_{n,p}$ can approximate the distribution of $T_{n,p}$ in Kolmogorov-Smirnov distance uniformly over all $p \in [1, \infty]$.

The theorem reveals that even in high dimensions the distribution of $T_{n,p}$ depends on the data mostly through the first and second moments, i.e. mean zero and covariance matrix $\Sigma_n$. This insight significantly simplifies the task of estimating the distribution of $T_{n,p}$ and is the rationale for the Gaussian parametric bootstrap procedure.

Another striking aspect of this result is the dependence on exponent $p\in [1, \infty]$. Namely, as the exponent $p$ crosses the threshold $\log d$, the upper bounds in $(i)-(iv)$ undergo a phase transition from polynomial in $r_n$ to logarithmic in $d$. This phase transition is directly related to similar behavior of the variance of $\ell_p$-norms of Gaussian random vectors paouris2018Dvoretzky. We discuss this technical aspect in greater detail in Appendix (ref).

Since this Gaussian approximation result is non-asymptotic we can take limits (with respect to $n,d, p$) in any order. Given the scope of the paper, we are most interested in the high-dimensional setting with $n,d \rightarrow \infty$ and $p \in [1, \infty)$ fixed. For this asymptotic regime we note the following: The bounds in cases $(i)$, $(ii)$, and $(iv)$ imply that the larger the exponent $p$ and the stronger the moment conditions on the $X_i$'s, the faster $d$ can grow (relative to $n$) while still guaranteeing that the distributions of $T_{n,p}$ and $\widetilde{T}_{n,p}$ are close. Case $(iii)$ (with $p = \infty$) covers the case of the max-statistic considered in chernozhukov2013GaussianApproxVec, chernozhukov2015ComparisonAnti, chernozhukov2017CLTHighDim and improves their bound by removing the dependence on the inverse of $\sigma_{n,\min}^2$. If the $X_i$'s are identically distributed then $\kappa_n^2 = 1$ and the bound is independent of any characteristic of the covariance matrix of the data (rank, eigenvalues, or diagonal values).

Since the $T_{n,2}$ statistic is of particular interest in many statistical applications, we provide the following easy corollary with a short discussion.

The result that is most related to Corollary (ref) is the dimension-dependent Berry-Essen bound by bentkus2003DependenceBerryEsseen. bentkus2003DependenceBerryEsseen addresses a slightly more general problem than we do: He derives a Berry-Esseen-type CLT for $S_n^X$ that holds uniformly over the class of Euclidean balls with arbitrary radii and arbitrary centers. In contrast, our Corollary (ref) corresponds to a Berry-Esseen-type CLT for $S_n^X$ that holds uniformly over the class of localized Euclidean balls with arbitrary radii but center fixed to the origin. The upper bound in Theorem 1.1 in bentkus2003DependenceBerryEsseen is at least of order $d^{3/2} n^{-1/2}$. It appears that part of the reason why we obtain a better dependence on dimension $d$ (relative to $n$) is that we consider only localized Euclidean balls.

There is a rich literature on the closely related problem of Gaussian approximations of quadratic forms bentkus1997uniform, goetze2014ExplicitRates, pouzo2015bootstrap, spokoiny2015bootstrap, goetze2019LargeBall, xu2019L2Asymptotics. The Berry-Esseen-type bounds in this literature often feature a better dependence on the sample size $n$, but either have a worse dependence on dimension $d$ relative to $n$, leave the dependence on $d$ wholly unaddressed, or do not apply to degenerate distributions (i.e. low-rank covariance matrix). In general, the existing bounds appear to be less useful for applications to high-dimensional statistics than our results in this section.

Bootstrap consistency

In this section we provide non-asymptotic bounds on the Kolmogorov-Smirnov distance between the distributions of the $\ell_p$-statistic $T_{n,p}$ and the Gaussian parametric bootstrap statistic $T_{n,p}^*$. As corollary we also show the consistency of the Gaussian parametric bootstrap.

Recall from Section (ref) that the Gaussian parametric bootstrap requires a positive semi-definite estimate $\widehat{\Sigma}_n$ of the (averaged) population covariance matrix $\Sigma_n$. The non-asymptotic bounds in this section depend on the following quantities

Note that $\widehat{\Delta}_p$ corresponds to the entry-wise $\ell_p$-norm of $\widehat{\Sigma}_n - \Sigma_n$ with $\widehat{\Delta}_2$ being the Frobenius norm. To establish the bootstrap consistency, we use

The first term on the right hand side in above display is deterministic and can be bounded by using Theorem (ref). The second term is stochastic and can be handled by the Gaussian comparison inequality in Appendix (ref).

The following theorem shows that the distributions of $T_{n,p}$ and $T_{n,p}^*$ are close in Kolmogorov-Smirnov distance uniformly over all $p \in [1, \infty]$ and for generic estimates $\widehat{\Sigma}_n$.

Theorem (ref) is only practically relevant in combination with estimates $\widehat{\Sigma}_n$ for which the stochastic estimation errors $\widehat{\Delta}_p$ and $\widehat{\Delta}_{op} \wedge \widehat{\Delta}_\infty$ are small. In Appendix (ref) we provide bounds on these quantities for several different estimates $\widehat{\Sigma}_n$. For the remainder of this section we consider the special case $\widehat{\Sigma}_n = \widehat{\Sigma}_{\mathrm{naive}} := n^{-1}\sum_{i=1}^n(X_i - \bar{X}_n)(X_i - \bar{X}_n)'$. We define the naive Gaussian parametric bootstrap estimate based on the sample covariance matrix $\widehat{\Sigma}_{\mathrm{naive}}$ by

Since $T_{n,p, \mathrm{naive}}^*$ is equivalent to the Gaussian multiplier statistic $T_{n,p}^g$, the following result is also a statement about the Gaussian multiplier bootstrap.

The main message of this corollary is that in high dimensions the naive Gaussian parametric and the Gaussian multiplier bootstrap can be consistent for large exponents $p \geq \log d$ but may fail to be consistent for small exponents $p \in [1, \log d)$. More precisely, cases $(i)$ and $(ii)$ imply that the naive Gaussian parametric and the Gaussian multiplier bootstrap are consistent in probability for small $p \in [1, \log d)$ if $d^{2/p} \log d = o(n)$ and for large $p \in [\log d, \infty]$ if $\log^7 d = o(n)$. Using the Borel-Cantelli lemma, one can easily turn this into sufficient conditions for “almost sure” bootstrap consistency.

Bootstrap consistency under structured covariance matrices

We establish two refined consistency results for the Gaussian parametric bootstrap in high dimensions. In particular, we significantly improve the rates of bootstrap consistency for small exponents $p \in [1, \log d)$ (cf. Corollary (ref) $(i)$) by exploiting certain sparsity and bandedness properties of the covariance matrix. We do not present results for large exponents $p \in [\log d, \infty]$ because in this regime sparsity and bandedness properties cannot be leveraged (and are also not needed) to further improve the rates given in Corollary (ref) $(ii)$.

To keep the discussion simple, we now assume that $X = \{X_i\}_{i=1}^n$ is a random sample of i.i.d. random vectors in $\mathbb{R}^d$ with mean zero and covariance matrix $\Sigma = (\sigma_{jk})_{j,k=1}^d$. We will drop the subscript $n$ in $r$, $\sigma_{\min}^2$, and $\sigma_{\max}^2$.

For $\gamma=0$ this assumption is most restrictive and implies that the covariance matrix is sparse with at most $R_{0,\theta}^\theta$ non-zero entries in each row. The covariance matrix of an $\mathrm{AR}$-process is a prominent example satisfying this assumption for some positive $\gamma > 0$.

The larger $\alpha > 0$, the more the covariance matrix $\Sigma$ resembles a diagonal matrix. Covariance matrices of $\mathrm{MA}$-processes satisfies this assumption for some finite $\alpha > 0$.

For $\theta=1$ Assumptions (ref) and (ref) reduce to two frequently adopted assumptions in the literature on high-dimensional covariance estimation bickel2008covariance, bickel2008regularization, cai2010optimal, cai2011adaptive, acella-medina2018robust. The larger $\theta$, the milder are the restrictions imposed on the covariance matrix.

Under Assumption (ref) it is natural to estimate the covariance matrix via thresholding of the naive sample covariance bickel2008covariance, Lam2009sparsistency. For simplicity, here we only consider the hard-thresholding operator; Appendix (ref) contains results for more general thresholding operators. For a matrix $M = (m_{jk})_{j,k=1}^d$ and $\lambda > 0$, we define the hard-thresholding operator by

Under Assumption (ref) it is common to estimate the covariance matrix via banding of the naive sample covariance bickel2008regularization: For a given $\ell > 0$, define

Recall that the Gaussian parametric bootstrap procedure requires a positive semi-definite estimate of the covariance matrix. If $\lambda_{\min}(\Sigma)$ and sample size $n$ are sufficiently large, bickel2008covariance and bickel2008regularization show that $T_{\lambda}(\widehat{\Sigma}_{\mathrm{naive}})$ and $B_{\ell}(\widehat{\Sigma}_{\mathrm{naive}})$ are positive definite with probability one. If the sample size is small we suggest projecting these estimates onto the cone of positive semi-definite matrices. Since the resulting positive semi-definite projections $T_{\lambda}^+(\widehat{\Sigma}_{\mathrm{naive}})$ and $B_{\ell}^+(\widehat{\Sigma}_{\mathrm{naive}})$ maintain the same order of $\ell_p$-error as the original estimates, this projection step does not add any additional theoretical challenges. Indeed, define

and observe that by triangular inequality and contraction property of projections,

The same reasoning applies to $B_{\ell}^+(\widehat{\Sigma}_{\mathrm{naive}})$. In the following, we therefore tacitly assume that this projection step has been applied and drop the superscript “+”.

We define the Gaussian parametric bootstrap statistics based on $ T_{\lambda}(\widehat{\Sigma}_{\mathrm{naive}})$ and $B_{\ell}(\widehat{\Sigma}_{\mathrm{naive}})$, respectively, by

and

where thresholding level $\lambda > 0$ and banding parameter $\ell > 0$ will be specified below. The next two corollaries refine Corollary (ref).

The main takeaway from these two corollaries is that under reasonable assumptions on the covariance structure and the tails of the data there exist Gaussian parametric bootstrap statistics $T_{n,p}^*$ that are consistent in high dimensions for any fixed $p \in [1, \infty)$.

In particular, inequality (ref) (inequality (ref)) implies that if the data is sub-Gaussian and the population covariance matrix is approximately sparse (approximately bandable) the Gaussian parametric bootstrap based on the thresholded covariance matrix (the banded covariance matrix) is consistent in probability for all $p \in [1, \infty)$ provided that $\log d = o(n^3)$.

Application: Testing high-dimensional mean vectors

As an application of the Gaussian parametric bootstrap we now present a bootstrap hypothesis test based on $\ell_p$-statistics for testing linear restrictions on high-dimensional mean vectors. We show that this test is asymptotic correct and consistent. Moreover, we discuss the effect of the exponent $p$ on the size of simultaneous confidence sets and the power of the test. Lastly, we discuss an extension of the generic testing framework to simultaneous inference on high-dimensional linear models.

Setup and test statistic

Given a random sample $X = \{X_i\}_{i=1}^n$ of i.i.d. random vectors in $\mathbb{R}^d$ with unknown mean $\mu$ and unknown covariance matrix $\Sigma$ we are interested in testing the high-dimensional linear restrictions

for some $M \in \mathbb{R}^{d' \times d}$ and $m_0 \in \mathbb{R}^{d'}$ when dimension $d$ and number of restrictions $d'$ may exceed the sample size $n$.

We propose to test hypothesis (ref) on the basis of the $\ell_p$-statistic

and, given a nominal level $\alpha \in (0,1)$, reject the null hypothesis if and only if

where $c^*_{n,p}(\alpha)$ is the $\alpha$-quantile of the Gaussian parametric bootstrap estimate

and $\widehat{\Omega}_n$ is a positive semi-definite estimate of $\Omega = M \Sigma M'$.

A distinguishing feature of this bootstrap hypothesis test is the exponent $p \in [1, \infty]$ and we show that the exponent $p$ has significant impact on the asymptotic correctness and the power of the test. In practice, tests based on $\ell_p$-statistics $S_{n,p}$ with exponents $1$, $2$, and $\infty$ are of particular interest. For one, the $\ell_1$-statistic $S_{n,1}$ and the maximum statistic $S_{n,\infty}$ lie at opposite ends of the spectrum of possible exponents $p$ and therefore have power functions that are complementary in a sense to be made precise below. For another, the maximum statistic $S_{n,\infty}$ can also be applied to the problem of multiple hypothesis testing. Since the bootstrap test based on $S_{n,\infty}$ accounts for the dependence between the multiple tests, it is (asymptotically) less conservative than the Bonferroni adjustment. Lastly, the sum-of-squares type statistic $S_{n,2}$ is essentially a feasible version of Hotelling's $T^2$-statistic in high dimensions and as such interesting in its own right fan2015power.

Let $\mathcal{H}_0 = \{\mu \in \mathbb{R}^d : M\mu = m_0\}$ and $\mathcal{H}_1 = \mathcal{H}_0^c$. Write $\Omega =(\omega_{jk})_{j,k=1}^{d'}$, and

Let $\widehat{\Omega}_n$ be a positive semi-definite estimate of $\Omega$ and define

We also introduce the following high-level assumption.

We emphasize that under rather mild conditions there exist estimates $\widehat{\Omega}_n$ such that $\widehat{\Gamma}_n$ and $\widehat{\Gamma}_{op} \wedge\widehat{\Gamma}_\infty$ satisfy the conditions in Assumption (ref); see Appendix (ref) for details.

Asymptotic correctness

In this section we show that the bootstrap hypothesis test has asymptotic correct size. We state the theorem in a non-asymptotic fashion to match the results from previous sections.

A special feature of this result is the real-valued random variable $\xi$. For now, assume that $\xi \equiv 0$ and let $\eta\downarrow 0 $ arbitrarily fast. In this case, Theorem (ref) provides non-asymptotic error bounds on the type I error of the bootstrap hypothesis test based on $\ell_p$-statistic $S_{n,p}$.

Next, consider the case in which $\xi$ is not identical to zero. Then, Theorem (ref) is a statement about the test statistic $R_{n,p} : = S_{n,p} + \xi$, where $\xi$ may be interpreted as approximation error. This is particularly useful if we want to test hypotheses about a parameter $\beta_0 \in \mathbb{R}^d$ for which there exists an estimator $\hat{\beta}$ that admits the expansion

In this case, the triangle inequality yields $|\xi| \leq \|r_n\|_p$. The primary example that we have in mind is the de-biased lasso estimator for linear models vandegeer2014on, zhang2014confidence. We elaborate on this idea in detail in Section (ref).

Confidence sets for high-dimensional parameters

We can use Theorem (ref) to construct consistent confidence sets $\mathcal{C}_{n,p} \subset \mathbb{R}^d$ for a high-dimensional parameter $\mu_0 \in \mathbb{R}^d$. To this end, set $M = I_d$, $m_0 = \mu_0$, and define

for a given nominal level $\alpha \in (0,1)$. Then, under Assumption (ref), Theorem (ref) guarantees that

Given the collection $\{\mathcal{C}_{n,p}\}_{p \geq 1}$ a practitioner will be most interested in knowing which of these confidence sets is “smallest”. To answer this question, we study how $p \in [1, \infty]$ affects the volume of $\mathcal{C}_{n,p}$ as $d, n \rightarrow \infty$. To simplify matters, we only consider $\Sigma = \sigma^2 I_d$.

Obviously, the confidence sets $\mathcal{C}_{n,p}$ are just $\ell_p$-norm balls with center $\bar{X}_n$ and radii $c^*_{n,p}(1- \alpha)/\sqrt{n}$. Recall that the volume of centered $d$-dimensional $\ell_p$-balls with radius $r > 0$, say $\mathcal{B}_p^d(r)$, is given by

Also, by Lemma (ref), Remark (ref), and Lemma 2 in schechtman1990volume, with probability approaching one, for all $\alpha \in (0, 1/2)$,

Whence, by (ref), (ref), and Sterling's formula we have

where $c_{p}^{1/p} \in (0.8856, 1]$. It is now easy to check that (asymptotically) the volume of $\mathcal{C}_{n,p}$ is a monotonically increasing function of the exponent $p$. In other words, confidence sets based on $\ell_p$-statistics $S_{n,p}$ with small exponents are less conservative than confidence sets based on, say, the maximum statistic $S_{n,\infty}$. Asymptotically, $\mathcal{C}_{n,1}$ is the smallest confidence set.

Consistency under high-dimensional alternatives

We now analyze the consistency of the bootstrap hypothesis test under high-dimensional alternatives. Let $Z \sim N(0, I_{d'})$ and define

and its “complement”

In words, $\mathcal{A}_p$ contains alternatives $(\mu_n)_{n \in \mathbb{N}}$ whose signals $\sqrt{n}\|M\mu_n - m_0\|_p$ asymptotically dominate the mean and standard deviation of the Gaussian proxy statistic $\|\Omega^{1/2}Z\|_p$; whereas $\mathcal{Z}_p$ consists of alternatives whose signals are asymptotically negligible compared to mean and standard deviation of $\|\Omega^{1/2}Z\|_p$.

The following result shows that the bootstrap hypothesis test is consistent for all $(\mu_n)_{n \in \mathbb{N}} \in \mathcal{A}_p$ and inconsistent for all $(\mu_n)_{n \in \mathbb{N}} \in \mathcal{Z}_p$.

Power and the role of the exponent $p$

It is part of statistical folklore that sum-of-squares type statistics have good power against “dense” alternatives, i.e alternatives whose signals in $M\mu$ are spread out over a large number of coordinates, whereas maximum type statistics are more powerful against “sparse” alternatives, i.e. alternatives with only a few strong signals in $M\mu$ fan2015power. Theorem (ref) allows us to verify this statement more formally. Let $M = I_d$, $m_0 = 0$, $\Sigma = \sigma^2I_d$, and define the set of alternatives

where $c \geq 1 $ is an absolute constant, $\delta > 0$ regulates the signal strength, and $s \in \{1, \ldots, d\}$ controls the sparsity.

Given this setup, we ask the following question: What is the minimum signal strength $\delta \equiv \delta(n, d, s, p)$ needed for the bootstrap test based on $S_{n,p}$ to reject the null hypothesis $H_0: \mu = 0$ at significance level $\alpha \in (0,1/2)$ when $\mu \in \mathcal{D}_{\delta, s}$?

By Remark (ref) $\sqrt{\mathrm{Var}\|\Omega^{1/2}Z\|_p} \leq \sigma$ and by Lemma 2 in schechtman1990volume $\mathrm{E}\|\Omega^{1/2}Z\|_p \asymp \sqrt{p} d^{1/p}$ for $p < \log d$ and $\mathrm{E}\|\Omega^{1/2}Z\|_p \asymp \sigma \sqrt{\log d}$ for $p \geq \log d$. Thus, by Theorem (ref) $(ii)$, a necessary condition for correctly rejecting the null hypothesis (with probability approaching one) is

Now, suppose that $s \asymp d$, i.e. $\mathcal{D}_{\delta, s}$ contains only dense alternatives. Then, for $p \in [1, \log d)$ (ref) holds if $\delta \gtrsim \sqrt{p/n}$, whereas for $p \in [\log d, \infty]$ (ref) holds only if $\delta \gtrsim \sqrt{(\log d)/n}$. Thus, bootstrap tests based on $\ell_p$-statistics with small exponents are more powerful in detecting dense alternatives than those based on $\ell_p$-statistics with large exponents.

Next, assume that $s \ll d$, i.e. $\mathcal{D}_{\delta, s}$ contains only sparse alternatives. Then, for $p \in [1, \log d)$ (ref) holds if $\delta \gtrsim \sqrt{p(d/s)^{2/p}/n}$, whereas for $p \in [\log d, \infty]$ (ref) holds already if $\delta \gtrsim \sqrt{(\log d)/n}$. Therefore, tests based on $\ell_p$-statistics with large exponents are more responsive to sparse alternatives than those based on $\ell_p$-statistics with small exponents.

Simultaneous inference on high-dimensional linear models

The bootstrap hypothesis test based on the $\ell_p$-statistic $S_{n,p}$ can be combined with the de-biased Lasso estimator vandegeer2014on, zhang2014confidence to conduct simultaneous inference on high-dimensional linear models. This approach extends the one by zhang2017SimultaneousInf, who propose a bootstrap test for the de-biased lasso estimator based on the Gaussian multiplier bootstrap for the maximum statistic $S_{n, \infty}$.

Consider the high-dimensional sparse model

with response $Y_i \in \mathbb{R}$, i.i.d. predictors $X_i \in \mathbb{R}^d$ with mean $\mu$ and covariance matrix $\Sigma$, i.i.d. errors $\varepsilon_i$ (independent of $X_i$) with mean 0 and variance $\sigma_\varepsilon^2$, and sparse regression vector $\beta_0$. We are interested in testing the linear hypothesis

Write $Y = (Y_1, \ldots,Y_n) \in \mathbb{R}^n$, $\varepsilon = (\varepsilon_1, \ldots, \varepsilon_n) \in \mathbb{R}^n$, and $\mathbf{X} = [X_1, \ldots, X_n]' \in \mathbb{R}^{n \times d}$. For $\lambda > 0$ define the ordinary lasso estimate by

and the de-biased lasso estimate by

where $\widehat{\Theta}$ is a suitable approximation of the inverse of the Gram matrix $\widehat{\Sigma} = \mathbf{X}'\mathbf{X}/n$. Define the $\ell_p$-statistic

and observe that

Further, note that the first term on the right hand side in above display is the re-scaled sum of $n$ i.i.d. random vectors with mean zero and covariance matrix $\sigma_\varepsilon^2 M\Sigma^{-1}M'$. Hence, $R_{n,p} = S_{n,p} + \xi$, where $S_{n,p} = \|n^{-1/2}\sum_{i=}^n M\Sigma^{-1} X_i \varepsilon_i\|_p$ and $|\xi| \leq \|r_1\|_p + \|r_2\|_p$. Under mild assumptions, $\|r_1\|_p \leq \|M(\widehat{\Theta} - \Sigma^{-1})\|_{q \rightarrow p} \|\mathbf{X}'\varepsilon\|_q/\sqrt{n}$ and $\|r_2\|_p \leq \|M(\widehat{\Theta} \widehat{\Sigma} - I_d)\|_{q \rightarrow p}\|\hat{\beta}_\lambda - \beta_0\|_q$, $q \geq 1$, are negligible vandegeer2014on. Thus, based on the expansion (ref) and the discussion in Section (ref) we can approximate the distribution of $R_{n,p}$ under the null hypothesis by the distribution of the Gaussian parametric bootstrap estimate

and $\hat{\sigma}_\varepsilon^2$ is a consistent estimate of the error variance $\sigma_{\varepsilon}^2$ fan2012variance.

We can now use the quantiles of $S^*_{n,p}$ to compute (bootstrap) critical values for the $\ell_p$-statistic $S_{n,p}$ and to construct confidence sets for $M\beta_0$.

Numerical experiments

The purpose of the numerical experiments is in this section is threefold. First, they show that for small exponents $p \in [1, \log d)$ the Gaussian parametric bootstrap outperforms the Gaussian multiplier bootstrap, while for large exponents $p \in [\log d, \infty)$ both bootstrap procedures perform similarly. Second, they confirm the theoretical claims from Section (ref) that for heavy-tailed data the accuracy of the Gaussian parametric and multiplier bootstrap suffers as the exponent $p$ increases. Third, they show that the exponent $p$ affects the power of the bootstrap hypothesis test as described in Section (ref).

Data generation

We generate vectors $X_1, \ldots, X_n \in \mathbb{R}^d$ via a Gaussian copula model

where the random vectors $Y_1, \ldots, Y_n \in \mathbb{R}^d$ are sampled independently and identically from a centered Gaussian distribution with sparse covariance matrix $\Sigma$, $\Phi$ is the cdf of the $N(0,1)$ distribution, and $F$ is the distribution function of either the uniform distribution on $[-1,1]$ (“light-tailed”) or Student's $t$-distribution with $4$ degrees of freedom (“heavy-tailed”). We create the sparse and low-rank covariance matrix $\Sigma$ in two steps: First, define the block diagonal matrix $\widetilde{\Sigma} = \mathrm{diag}(\Lambda, \ldots, \Lambda) \in \mathbb{R}^{d \times d}$, where $\Lambda = (\Lambda_{jk})_{j,k=1}^{d/100}$ with $\Lambda_{jk} = 0.8^{j+k -2}$ for all $1 \leq j,k \leq d/100$ is a rank-one matrix. Then, (randomly) generate a permutation matrix $P$ and set $\Sigma = P\widetilde{\Sigma} P'$. The matrix $\Sigma$ is positive semi-definite, sparse with $d/100$ non-zero elements in each row, and has rank $100$. The permutation matrix $P$ is generated only once and is the same throughout all Monte Carlo simulations.

Specific implementation of a hard-thresholding estimator

The Gaussian parametric bootstrap procedure requires as input a positive semi-definite estimate of the population covariance matrix $\Sigma$. To exploit the sparsity of $\Sigma$ while also ensuring positive semi-definiteness of the estimate, we propose the following two-step procedure:

First, compute a pilot estimate via correlation thresholding fan2011high of the sample covariance matrix $\widehat{\Sigma}_{\mathrm{naive}} = \left(\hat{\sigma}_{jk}\right)_{j,k =1}^d$, i.e. for $\lambda > 0$ compute

Then, project the pilot estimate $\widehat{\Sigma}_n(\lambda)$ onto the cone of positive semi-definite matrices by setting all negative eigenvalues equal to 0. Denote the resulting estimate by $\widehat{\Sigma}_n^+(\lambda)$.

It remains to choose the thresholding level $\lambda > 0$. We proceed as in bickel2008covariance, bickel2008regularization and select $\lambda$ by cross-validation: At each fold $\nu \in \{1, \ldots, N\}$, randomly split the sample $X= \{X_i\}_{i=1}^n$ into two sub-samples $X^1$ and $X^2$ of sizes $n_1 = \lceil n/3 \rceil$ and $n_2 = n - n_1$, respectively. Denote by $\widehat{\Sigma}_{1, \nu}$ and $\widehat{\Sigma}_{2,\nu}$ the sample covariance matrices of the $\nu$th split based on $X^1$ and $X^2$. Let $\widehat{\Sigma}_{1,\nu}^{+}(\lambda)$ be the correlation-thresholded and projected estimate based on $\widehat{\Sigma}_{1, \nu}$. Define the cross-validated risk at level $\lambda > 0$ by

and select the “optimal” thresholding level as

In practice, we set $N = 10$ and minimize the risk $\widehat{R}(\lambda)$ over a grid $\mathcal{G} \subset [0, 1]$ with $|\mathcal{G}| =40$ equally spaced points. The algorithm is sensitive to the number of folds and grid points; increasing $N$ and $|\mathcal{G}|$ beyond 10 and 40, respectively, can improve the accuracy of the bootstrap approximation (at the cost of additional computational complexity).

Performance of Gaussian parametric and multiplier bootstrap

To assess the performance of the Gaussian parametric and the multiplier bootstrap in finite samples, we provide two types of plots:

In the following discussion the Gaussian proxy statistic $\widetilde{T}_{n,p}$ serves as an oracle estimator. It tells us how good the bootstrap procedures could be if we knew the true covariance matrix. Any difference between $\widetilde{T}_{n,p}$ and the other statistics solely arises from the different estimates of the covariance matrix.

Figure (ref) shows that if the data has light tails, the distribution of the Gaussian proxy statistic $\widetilde{T}_{n,p}$ provides an excellent approximation of the distribution of the $\ell_p$-statistic $T_{n,p}$ for all $p \in \{1, 2, \log d , \infty\}$. Moreover, the distribution of the Gaussian parametric bootstrap statistic $T^*_{n,p,\lambda}$ based on $\widehat{\Sigma}_n^+(\hat{\lambda})$ yields a comparably good approximation to the truth. In contrast, the distributions of the Gaussian multiplier statistic and the naive Gaussian parametric bootstrap are significantly poorer approximations to the truth. For large exponents $p \in \{\log d, \infty\}$ all four bootstrap approximations perform similarly. Thus, this plot fully supports every aspect of the theoretical results derived in Sections (ref) and (ref).

Figure (ref) shows that if the data has heavy tails, Gaussian proxy statistic $\widetilde{T}_{n,p}$ and the bootstrap procedures yield poorer approximations to the truth. In particular, we see that the quality of the approximation worsens substantially as the exponent $p$ increases. This further corroborates the theoretical results derived in Sections (ref) and (ref).

Figures (ref) and (ref) tell a similar, but more nuanced, story. From Figure (ref) we infer that if the data has light tails, the 95% quantiles of the distributions of the Gaussian proxy statistic $\widetilde{T}_{n,p}$ and the Gaussian parametric bootstrap statistic $T^*_{n,p,\lambda}$ based on $\widehat{\Sigma}_n^+(\hat{\lambda})$ yield good approximations to the 95% quantile of the true distribution for all exponents $p \in \{1, 2, \log d, \infty\}$. From Figure (ref) we learn that if the data has heavy tails, the approximations are fairly good for small exponents $p \in \{1, 2\}$, but fail spectacularly for large exponents $p \in \{\log d , \infty\}$. Moreover, Gaussian multiplier and naive Gaussian parametric bootstrap yield accurate estimates of the 95% quantiles of the target distribution only for large exponents $p \in \{\log d, \infty\}$ and only when the data has light tails. This again supports the theoretical results from Sections (ref) and (ref).

Power of the bootstrap hypothesis test

To illustrate the effect of the exponent $p$ on size and power of the bootstrap hypothesis test we consider its power function in the following two high-dimensional testing scenarios:

In Figures (ref) and (ref) we plot Monte Carlo estimates of the power function $\beta(\delta) = \mathrm{P}_{\mu(\delta)}\big(T_{n,p} > c^*_{n,p}(0.95)\big)$, where $c^*_{n,p}(0.95) = \inf\left\{t \in \mathbb{R} : \mathrm{P}(S_{n,p}^* \leq t \mid X) \geq 0.05 \right\}$ and $S_{n,p}^*$ is the Gaussian parametric bootstrap test statistic based on $\widehat{\Sigma}_n^+(\hat{\lambda})$. The Monte Carlo estimate of $\beta(\delta)$ is based on 1000 Monte Carlo samples of $X = \{X_i \in \mathbb{R}^d, 1\leq i \leq n\}$ and 1000 bootstrap samples for each observed $X = \{X_i \in \mathbb{R}^d, 1\leq i \leq n\}$. The specific estimation procedure is identical to the one used to compute the lower tail probabilities in Section (ref). We report results for sample size $n = 200$, dimension $d= 400$, exponents $p \in\{1, 2, \log d, \infty\}$, and light-tailed data.

Figure (ref) shows the power function for dense alternatives. We observe that tests based on $S_{n,1}$ and $S_{n,2}$ (they are nearly indistinguishable in the figure) are more powerful than those based on $S_{n,\log d}$ and $S_{n,\infty}$. This fully matches the theoretical predictions from Section (ref). Figure (ref) displays the power function for sparse alternatives. In this case the bootstrap tests based on $S_{n,\log d}$ and $S_{n,\infty}$ are more powerful than those based on $S_{n,1}$ and $S_{n,2}$. The power functions associated with $S_{n,\log d}$ and $S_{n,\infty}$ are essentially the same with $S_{n,\log d}$ being slightly more powerful because of a larger constant (note $\mu$ has $2 \lceil \sqrt{\log d}/2\rceil = 4$ non-zero entries). Again, these findings fit well into the discussion in Section (ref).

Lastly, at $\delta = 0$ the power functions of all four tests are about 0.05 in both Figures (ref) and (ref). Thus, all four tests successfully control the type I error a 5% significance level. For large values of $\delta$ all four tests unanimously reject the null hypothesis with probability (close to) one. This confirms the results from Sections (ref) and (ref).

Conclusion

In this paper we have introduced the Gaussian parametric bootstrap to estimate the distribution of $\ell_p$-statistics of high-dimensional random vectors. The procedure is versatile and user-friendly, since its implementation requires only a positive semi-definite estimate of the population covariance matrix. The main theoretical contributions state the consistency of the Gaussian parametric bootstrap under various conditions on the covariance structure of the data. To showcase the applicability of the Gaussian parametric bootstrap we propose a bootstrap hypothesis test for simultaneous inference on high-dimensional mean vectors. We discuss in detail asymptotic correctness, confidence sets, consistency under high-dimensional alternatives, and power of the test.

One of the current challenges in theoretical statistics is to understand when bootstrap procedures work in high-dimensional problems. At least for bootstrapping $\ell_p$-statistics of high-dimensional random vectors we can now give a definitive answer. The technical results in the appendix to this paper clarify that the success of bootstrapping $\ell_p$-statistics hinges on three factors: (a) $\ell_p$-norms of high-dimensional random vectors satisfy a Berry-Esseen-type central limit theorem under relatively mild moment conditions; (b) $\ell_p$-norms of Gaussian random vectors satisfy powerful anti-concentration inequalities; (c) the distributions of $\ell_p$-statistics under centered Gaussian distributions vary smoothly over their covariance matrices.

\expandafter\def\expandafter\appendixpagename \expandafter{\expandafter\appendixpagename}