#  Gaussian and Bootstrap Approximations for Suprema of Empirical Processes

Alexander Giessing Department of Statistics, University of Washington,
Seattle, WA. E-mail<EMAIL_ADDRESS>
###### Abstract

In this paper we develop non-asymptotic Gaussian approximation results for the
sampling distribution of suprema of empirical processes when the indexing
function class $\mathcal{F}_{n}$ varies with the sample size $n$ and may not
be Donsker. Prior approximations of this type required upper bounds on the
metric entropy of $\mathcal{F}_{n}$ and uniform lower bounds on the variance
of $f\in\mathcal{F}_{n}$ which, both, limited their applicability to high-
dimensional inference problems. In contrast, the results in this paper hold
under simpler conditions on boundedness, continuity, and the strong variance
of the approximating Gaussian process. The results are broadly applicable and
yield a novel procedure for bootstrapping the distribution of empirical
process suprema based on the truncated Karhunen-Loève decomposition of the
approximating Gaussian process. We demonstrate the flexibility of this new
bootstrap procedure by applying it to three fundamental problems in high-
dimensional statistics: simultaneous inference on parameter vectors, inference
on the spectral norm of covariance matrices, and construction of simultaneous
confidence bands for functions in reproducing kernel Hilbert spaces.

Keywords: Gaussian Approximation; Gaussian Comparison; Gaussian Process
Bootstrap; High-Dimensional Inference.

##  1 Introduction

###  1.1 Approximating the sampling distribution of suprema of empirical
processes

This paper is concerned with non-asymptotic bounds on the Kolmogorov distance
between the sampling distribution of suprema of empirical processes and the
distribution of suprema of Gaussian proxy processes. Consider a simple random
sample of independent and identically distributed random variables
$X_{1},\ldots,X_{n}$ with common law $P$ taking values in the measurable space
$(S,\mathcal{S})$. Let $\mathcal{F}_{n}$ be a class of measurable functions
$f:S\rightarrow\mathbb{R}$ and define the empirical process

$\displaystyle\mathbb{G}_{n}(f)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\big{(}f(X_{i})-Pf\big{)},\quad\quad
f\in\mathcal{F}_{n}.$

Also, denote by $\\{G(f):f\in\mathcal{F}_{n}\\}$ a centered Gaussian process
with some positive semi-definite covariance function. Within this minimal
setup (and a few additional technical assumptions), the paper addresses the
problem of deriving non-asymptotic bounds on

$\displaystyle\varrho_{n}:=\sup_{s\geq
0}\left|\mathbb{P}\Big{\\{}\sup_{f\in\mathcal{F}_{n}}|\mathbb{G}_{n}(f)|\leq
s\Big{\\}}-\mathbb{P}\Big{\\{}\sup_{f\in\mathcal{F}_{n}}|G(f)|\leq
s\Big{\\}}\right|,$ (1)

and, consequently, conditions under which an empirical process indexed by
$\mathcal{F}_{n}$ is _Gaussian approximable_ , i.e. $\varrho_{n}\rightarrow 0$
as $n\rightarrow\infty$. Evidently, if $\mathcal{F}_{n}$ is a Donsker class
then the empirical process is trivially Gaussian approximable by the centered
Gaussian process with covariance function $(f,g)\mapsto P(fg)-(Pf)(Pg)$. Thus,
in this paper we are concerned with conditions that are strictly weaker than
those that guarantee a central limit theorem.

Gaussian approximation results of this type have gained significant interest
recently, as they prove to be effective in tackling high- and infinite-
dimensional inference problems. Applications include simultaneous inference on
high-dimensional parameters (Dezeure et al., 2017; Zhang and Cheng, 2017),
inference on nonparametric models (Chernozhukov et al., 2014a; Chen et al.,
2015, 2016), testing for spurious correlations (Fan et al., 2018), testing for
shape restrictions (Chetverikov, 2019), inference on large covariance matrices
(Chen, 2018; Han et al., 2018; Lopes et al., 2019), goodness-of-fit tests for
high-dimensional linear models (Janková et al., 2020), simultaneous confidence
bands in functional data analysis (Lopes et al., 2020; Singh and Vijaykumar,
2023), and error quantification in randomized algorithms (Lopes et al., 2023).

The theoretical groundwork on Gaussian approximation was laid in three seminal
papers by Chernozhukov et al. (2013, 2014b, 2015). Since then, subsequent
theoretical works have developed numerous refinements and extensions (e.g.
Chernozhukov et al., 2016, 2019, 2020; Deng and Zhang, 2020; Fang and Koike,
2021; Kuchibhotla et al., 2021; Cattaneo et al., 2022; Lopes, 2022a, b; Bong
et al., 2023). In this paper, we resolve two limitations of Chernozhukov et
al.’s (2014b) original results that have not been previously addressed:

  * •

_Entropy conditions._ The original bounds on the Kolmogorov distance
$\varrho_{n}$ depend on the metric entropy of the function class
$\mathcal{F}_{n}$. These upper bounds are non-trivial only if the metric
entropy grows at a much slower rate than the sample size $n$. As a result,
most applications only address sparse inference problems involving function
classes with either a finite number of functions, discretized functions, or a
small VC-index. The new upper bounds in this paper no longer depend on the
metric entropy of the function class and therefore open the possibility to
tackle non-sparse, high-dimensional inference problems. We call these bounds
_dimension-/entropy-free_.

  * •

_Lower bounds on weak variances._ The original bounds on $\varrho_{n}$ require
a strictly positive lower bound on the _weak variance_ of the Gaussian proxy
process, i.e. $\inf_{f\in\mathcal{F}_{n}}\mathrm{Var}\big{(}G(f)\big{)}>0$.
This condition limits the scope of the original results to problems with
standardized function classes (studentized statistics) and excludes situations
with variance decay, typically observed in non-sparse, high-dimensional
problems. The new Gaussian approximation results in this paper only depend on
the _strong variance_
$\mathrm{Var}\big{(}\sup_{f\in\mathcal{F}_{n}}|G(f)|\big{)}$ and are therefore
applicable to a broad range of problems including those with degenerate
distributions. We say that these bounds are _weak variance-free_.

Even though these limitations (and our solution) are quite technical, the
resulting new Gaussian and bootstrap approximations have immediate practical
consequences. We present three substantial applications to problems in high-
dimensional inference in Section 4 and two toy examples that cast light on why
and when the original approximation results by Chernozhukov et al. (2014b)
fail and the new results succeed in Appendix A.

Notably, in the special case of inference on spectral norms of covariance
matrices, entropy-free bounds on $\varrho_{n}$ already exist. Namely, for
$\mathcal{F}_{n}=\\{x\mapsto f(x)=(x^{\prime}u)(x^{\prime}v):u,v\in
S^{d-1}\\}$ with $S^{d-1}=\\{u\in\mathbb{R}^{d}:\|u\|_{2}=1\\}$ Lopes et al.
(2019; 2020; 2022a; 2022b) have devised an approach to bounding $\varrho_{n}$
that combines a specific variance decay assumption with a truncation argument.
For a carefully chosen truncation level, the resulting bound on $\varrho_{n}$
depends only on the sample size $n$ and the parameter characterizing the
variance decay. Since their approach is intimately related to the bilinearity
of the functions in $\mathcal{F}_{n}$ and the specific variance decay
assumption, it does not easily extend to arbitrary function classes. We
therefore develop a different strategy in this paper.

Furthermore, for finite function classes $|\mathcal{F}_{n}|<\infty$,
Chernozhukov et al. (2020) and Deng and Zhang (2020) have been able to
slightly relax the requirements on the weak variance of the Gaussian proxy
process. However, their results do not generalize to arbitrary function
classes with $|\mathcal{F}_{n}|=\infty$. Our strategy for replacing the weak
variance with the strong variance in the upper bounds on $\varrho_{n}$ is
therefore conceptually completely different from theirs. For details, we refer
to our companion paper on Gaussian anti-concentration inequalities Giessing
(2023).

###  1.2 Contributions and overview of the results

This paper consists of three parts which contribute to probability theory
(Section 2), mathematical statistics and bootstrap methodology (Section 3),
and high-dimensional inference (Section 4). The Appendices A–F contain
additional supporting results and all proofs.

Section 2 contains the main mathematical innovations of this paper. We
establish dimension- and weak variance-free Gaussian and bootstrap
approximations for maxima of sums of independent and identically distributed
high-dimensional random vectors (Section 2.1). Specifically, we derive a
Gaussian approximation inequality (Proposition 1), a Gaussian comparison
inequality (Proposition 3), and two bootstrap approximation inequalities
(Propositions 4 and 5). At the core of these four theoretical results is a new
proof of a multivariate Berry-Esseen-type bound which leverages two new
technical auxiliary results: an anti-concentration inequality for suprema of
separable Gaussian processes (Lemma 6 in Appendix B.2) and a smoothing
inequality for partial derivatives of the Ornstein-Uhlenbeck semigroup
associated to a multivariate normal measure (Lemma 2 in Appendix B.1). We
conclude this section with a comparison to results from the literature
(Section 2.2).

Section 3 contains our contributions to mathematical statistics and bootstrap
methodology. The results in this section generalize the four basic results of
Section 2 from finite dimensional vectors to empirical processes indexed by
totally bounded (and hence separable) function classes. As one would expect,
dimension- (a.k.a entropy-) and weak variance-freeness of the finite
dimensional results carry over to empirical processes. We establish Gaussian
approximation inequalities (Section 3.2), Gaussian comparison inequalities
(Section 3.3), and an abstract bootstrap approximation inequality (Section
3.4). The latter result motivates a new procedure for bootstrapping the
sampling distribution of suprema of empirical processes. We call this
procedure the _Gaussian process bootstrap_ (Algorithm 1) and discuss practical
aspects of its implementations via the truncated Karhunen-Loève decomposition
of a Gaussian proxy process (Section 3.5). We include a selective comparison
with results from the literature (Section 3.6).

In Section 4 we showcase the flexibility of the Gaussian process bootstrap by
applying it to three fundamental problems in high-dimensional statistics:
simultaneous inference on parameter vectors (Section 4.1), inference on the
spectral norm of covariance matrices (Section 4.2), and construction of
simultaneous confidence bands for functions in reproducing kernel Hilbert
spaces (Section 4.3). For each of these three examples we include a brief
comparison with alternative bootstrap methods from the literature to explain
in what sense our method improves over (or matches with) existing results.

##  2 Dimension- and weak variance-free results for the maximum norm of sums
of i.i.d. random vectors

###  2.1 Four basic results

We begin with four propositions on non-asymptotic Gaussian and bootstrap
approximations for finite-dimensional random vectors. These propositions are
our main mathematical contribution and essential for the statistical results
to come. In Section 3 we lift these propositions to general empirical
processes which are widely applicable to problems in mathematical statistics.
Throughout this section, $\|x\|_{\infty}=\max_{1\leq k\leq d}|x_{k}|$ denotes
the maximum norm of a vector $x\in\mathbb{R}^{d}$.

The first result is a Gaussian approximation inequality. This inequality
provides a non-asymptotic bound on the Kolmogorov distance between the laws of
the maximum norm of sums of independent and identically distributed random
vectors and a Gaussian proxy statistic.

######  Proposition 1 (Gaussian approximation).

Let $X,X_{1},\ldots X_{n}\in\mathbb{R}^{d}$ be i.i.d. random vectors with mean
zero and positive semi-definite covariance matrix
$\Sigma\neq\mathbf{0}\in\mathbb{R}^{d\times d}$. Set
$S_{n}=n^{-1/2}\sum_{i=1}^{n}X_{i}$ and $Z\sim N(0,\Sigma)$. Then, for $M\geq
0$, $n\geq 1$,

$\displaystyle\sup_{s\geq 0}\Big{|}\mathbb{P}\left\\{\|S_{n}\|_{\infty}\leq
s\right\\}-\mathbb{P}\left\\{\|Z\|_{\infty}\leq s\right\\}\Big{|}$
$\displaystyle\quad{}\quad{}\quad{}\lesssim\frac{(\mathrm{E}[\|X\|_{\infty}^{3}])^{1/3}}{\sqrt{n^{1/3}\mathrm{Var}(\|Z\|_{\infty})}}+\frac{\mathrm{E}\left[\|X\|_{\infty}^{3}\mathbf{1}\\{\|X\|_{\infty}>M\\}\right]}{\mathrm{E}\left[\|X\|_{\infty}^{3}\right]}+\frac{\mathrm{E}[\|Z\|_{\infty}]+M}{\sqrt{n\mathrm{Var}(\|Z\|_{\infty})}},$

where $\lesssim$ hides an absolute constant independent of $n,d,M$, and the
distribution of the $X_{i}$’s.

######  Remark 1 (Extension to independent, non-identically distributed random
vectors).

The proof of Proposition 1 involves an inductive argument inspired by Theorem
3.7.1 in Nourdin and Peccati (2012) who attribute it to Bolthausen (1984).
This argument relies crucially on independent and identically distributed
data. Generalizing it to independent but non-identically distributed data
would require additional assumptions on the variances similar to the uniform
asymptotic negligibility condition in the classical Lindeberg-Feller CLT for
triangular arrays. For empirical processes such a generalization is less
relevant. We therefore leave this to future research.

######  Remark 2 (Extension to non-centered random vectors).

If the covariance matrix $\Sigma$ is strictly positive definite, then
Proposition 1 also holds for i.i.d. random vectors $X,X_{1},\ldots
X_{n}\in\mathbb{R}^{d}$ with non-zero mean $\mu\in\mathbb{R}^{d}$ and $Z\sim
N(\mu,\Sigma)$. Nevertheless, in the broader context of empirical processes a
strictly positive definite covariance (function) is a very strong assumption.
Therefore, in this paper, we do not pursue this refinement. Instead, the
interested reader may consult the companion paper Giessing (2023).

######  Remark 3 (Gaussian approximation versus CLT).

Proposition 1 is strictly weaker than a multivariate CLT because the
Kolmogorov distance between the maximum norm of two sequences of random
vectors induces a topology on the set of probability measures on
$\mathbb{R}^{d}$ which is coarser than the topology of convergence in
distribution. Indeed, consider $U_{n}=(X,Y)^{\prime}\in\mathbb{R}^{2}$ and
$V_{n}=(Y,X)^{\prime}\in\mathbb{R}^{2}$ for $n\geq 1$ and $X,Y$ arbitrary
random variables. Then, the Kolmogorov distance between $\|U_{n}\|_{\infty}$
and $\|V_{n}\|_{\infty}$ is zero for all $n\geq 1$, but
$U_{n}\overset{d}{=}V_{n}$ if and only if $X$ and $Y$ are exchangeable. For
another perspective on the same issue, see Section 2.2.

Above Gaussian approximation result differs in several ways from related
results in the literature. The two most striking difference are the following:
First, the non-asymptotic upper bound in Proposition 1 does not explicitly
depend on the dimension $d$, but only on the (truncated) third moments of
$\|X\|_{\infty}$ and $\|Z\|_{\infty}$ and the variance of $\|Z\|_{\infty}$.
Under suitable conditions on the marginal and/ or joint distribution of the
coordinates of $X$ these three quantities grow substantially slower than the
dimension $d$. Among other things, this opens up the possibility of applying
this bound in the context of high-dimensional random vectors when $d\gg n$.
Second, Proposition 1 holds without stringent assumptions on the distribution
of $X$ and applies even to degenerate distributions that do not have a
strictly positive definite covariance matrix $\Sigma$. In particular, the
lemma does not require a lower bound on the variances of the coordinates of
$X$ or the minimum eigenvalue of $\Sigma$. This is fundamental for an
effortless extension to general empirical processes. For lower bounds on the
variance of $\|Z\|_{\infty}$ the reader may consult Giessing (2023). For a
comprehensive comparison of Proposition 1 with previous work, we refer to
Section 2.2.

If we are willing to impose a lower bound on the variances of the coordinates
of $X$, we can deduce the following useful inequality:

######  Corollary 2.

Recall the setup of Proposition 1. In addition, suppose that
$\sigma_{(1)}^{2}=\min_{1\leq j\leq d}\mathrm{Var}(X^{(j)})>0$ and that $X$
has coordinate-wise finite $3+\delta$ moments, $\delta>0$. Define
$\widetilde{X}=(X^{(j)}/\sigma_{(1)})_{j=1}^{d}$,
$\widetilde{Z}=(Z^{(j)}/\sigma_{(1)})_{j=1}^{d}$. Then, for all $n\geq 1$,

$\displaystyle\sup_{s\geq 0}\Big{|}\mathbb{P}\left\\{\|S_{n}\|_{\infty}\leq
s\right\\}-\mathbb{P}\left\\{\|Z\|_{\infty}\leq s\right\\}\Big{|}$
$\displaystyle\quad\quad\quad\lesssim\frac{1}{n^{1/6}}\left(\mathrm{E}[\|\widetilde{X}\|_{\infty}^{3+\delta}]+\mathrm{E}[\|\widetilde{Z}\|_{\infty}^{3+\delta}]\right)^{\frac{1}{3+\delta}}\mathrm{E}[\|\widetilde{Z}\|_{\infty}]+\frac{1}{n^{\delta/3}}\left(\frac{\mathrm{E}[\|\widetilde{X}\|_{\infty}^{3+\delta}]^{\frac{1}{3+\delta}}}{\mathrm{E}[\|\widetilde{X}\|_{\infty}^{3}]^{1/3}}\right)^{3},$

where $\lesssim$ hides an absolute constant independent of
$n,d,\sigma_{(1)}^{2}$, and the distribution of the $X_{i}$’s. If the
coordinates in $X$ are equicorrelated with correlation coefficient
$\rho\in(0,1]$, then above inequality holds with
$\mathrm{E}[\|\widetilde{Z}\|_{\infty}]$ replaced by $1/\rho$.

Since $\mathrm{E}[\|\widetilde{X}\|_{\infty}^{3}]^{1/3}\geq\max_{1\leq j\leq
d}\mathrm{E}[|\widetilde{X}^{(j)}|^{2}]^{1/2}=\frac{\sigma_{(n)}}{\sigma_{(1)}}$
for $\sigma_{(n)}^{2}=\max_{1\leq j\leq d}\mathrm{Var}(X^{(j)})$, the right
hand side of the inequality in Corollary 2 can be easily upper bounded under a
variety of moment assumptions on the marginal distributions of the coordinates
of $X$. For a few concrete examples relevant in high-dimensional statistics we
refer to Lemmas 1–3 in Giessing and Fan (2023).

The second basic result is a Gaussian comparison inequality. The novelty of
this result is (again) that it holds even for degenerate Gaussian laws with
singular covariance matrices and does not explicitly depend on the dimension
of the random vectors.

######  Proposition 3 (Gaussian comparison).

Let $Y,Z\in\mathbb{R}^{d}$ be Gaussian random vectors with mean zero and
positive semi-definite covariance matrices $\Sigma$ and
$\Omega\neq\mathbf{0}\in\mathbb{R}^{d\times d}$, respectively. Then,

$\displaystyle\sup_{s\geq 0}\Big{|}\mathbb{P}\left\\{\|Y\|_{\infty}\leq
s\right\\}-\mathbb{P}\left\\{\|Z\|_{\infty}\leq
s\right\\}\Big{|}\lesssim\left(\frac{\max_{j,k}|\Omega_{jk}-\Sigma_{jk}|}{\mathrm{Var}(\|Y\|_{\infty})\vee\mathrm{Var}(\|Z\|_{\infty})}\right)^{1/3},$

where $\lesssim$ hides an absolute constant independent of $d,\Sigma,\Omega$.

Since the Gaussian distribution is fully characterized by its first two
moments, Proposition 1 instantly suggests that it should be possible to
approximate the sampling distribution of $\|S_{n}\|_{\infty}$ with the
sampling distribution of $\|Z_{n}\|_{\infty}$, where $Z_{n}\mid
X_{1},\ldots,X_{n}\sim N(0,\widehat{\Sigma}_{n})$ and $\widehat{\Sigma}_{n}$
is a positive semi-definite estimate of $\Sigma$. The third basic result
formalizes this idea; it is a simple consequence of the triangle inequality
combined with Propositions 1 and 3:

######  Proposition 4 (Bootstrap approximation of the sampling distribution).

Let $X,X_{1},\ldots X_{n}\in\mathbb{R}^{d}$ be i.i.d. random vectors with mean
zero and positive semi-definite covariance matrix
$\Sigma\neq\mathbf{0}\in\mathbb{R}^{d\times d}$. Let
$\widehat{\Sigma}_{n}\equiv\widehat{\Sigma}_{n}(X_{1},\ldots,X_{n})$ be any
positive semi-definite estimate of $\Sigma$. Set
$S_{n}=n^{-1/2}\sum_{i=1}^{n}X_{i}$, $Z_{n}\mid X_{1},\ldots,X_{n}\sim
N(0,\widehat{\Sigma}_{n})$, and $Z\sim N(0,\Sigma)$. Then, for $M\geq 0$,
$n\geq 1$,

$\displaystyle\sup_{s\geq 0}\Big{|}\mathbb{P}\left\\{\|S_{n}\|_{\infty}\leq
s\right\\}-\mathbb{P}\left\\{\|Z_{n}\|_{\infty}\leq s\mid
X_{1},\ldots,X_{n}\right\\}\Big{|}$
$\displaystyle\quad{}\quad{}\lesssim\frac{(\mathrm{E}[\|X\|_{\infty}^{3}])^{1/3}}{\sqrt{n^{1/3}\mathrm{Var}(\|Z\|_{\infty})}}+\frac{\mathrm{E}\left[\|X\|_{\infty}^{3}\mathbf{1}\\{\|X\|_{\infty}>M\\}\right]}{\mathrm{E}\left[\|X\|_{\infty}^{3}\right]}$
$\displaystyle\quad{}\quad{}\quad{}\quad{}+\frac{\mathrm{E}[\|Z\|_{\infty}]+M}{\sqrt{n\mathrm{Var}(\|Z\|_{\infty})}}+\left(\frac{\max_{j,k}|\widehat{\Sigma}_{n,jk}-\Sigma_{jk}|}{\mathrm{Var}(\|Z\|_{\infty})}\right)^{1/3},$

where $\lesssim$ hides an absolute constant independent of $n,d,M$, and the
distribution of the $X_{i}$’s.

Typically, statistical applications require estimates of the quantiles of the
sampling distribution. Since the covariance matrix $\Sigma$ is unknown, the
quantiles of $\|Z\|_{p}$ with $Z\sim N(0,\Sigma)$ are infeasible. Hence, for
$\alpha\in(0,1)$ arbitrary, we define the feasible Gaussian bootstrap
quantiles as

$\displaystyle c_{n}(\alpha;\widehat{\Sigma}_{n})$
$\displaystyle:=\inf\left\\{s\geq 0:\mathbb{P}\left\\{\|Z_{n}\|_{\infty}\leq
s\mid X_{1},\ldots,X_{n}\right\\}\geq\alpha\right\\},\quad\text{where}\quad
Z_{n}\mid X_{1},\ldots,X_{n}\sim N(0,\widehat{\Sigma}_{n}).$

Since this quantity is random, it is not immediately obvious that it is a
valid approximation of the $\alpha$-quantile of the sampling distribution of
$\|S_{n}\|_{p}$. However, combing Proposition 4 with standard arguments (e.g.
Chernozhukov et al., 2013) we obtain the fourth basic result:

######  Proposition 5 (Bootstrap approximation of quantiles).

Consider the setup of Lemma 4. Let $(\Theta_{n})_{n\geq 1}\in\mathbb{R}$ be a
sequence of arbitrary random variables, not necessarily independent of
$X,X_{1},\ldots,X_{n}$. Then, for $M\geq 0$, $n\geq 1$,

$\displaystyle\sup_{\alpha\in(0,1)}\Big{|}\mathbb{P}\left\\{\|S_{n}\|_{\infty}+\Theta_{n}\leq
c_{n}(\alpha;\widehat{\Sigma}_{n})\right\\}-\alpha\Big{|}$
$\displaystyle\quad{}\quad{}\lesssim\frac{(\mathrm{E}[\|X\|_{\infty}^{3}])^{1/3}}{\sqrt{n^{1/3}\mathrm{Var}(\|Z\|_{\infty})}}+\frac{\mathrm{E}\left[\|X\|_{\infty}^{3}\mathbf{1}\\{\|X\|_{\infty}>M\\}\right]}{\mathrm{E}\left[\|X\|_{\infty}^{3}\right]}+\frac{\mathrm{E}[\|Z\|_{\infty}]+M}{\sqrt{n\mathrm{Var}(\|Z\|_{\infty})}}$
$\displaystyle\quad{}\quad{}\quad{}\quad+\inf_{\delta>0}\left\\{\left(\frac{\delta}{\mathrm{Var}(\|Z\|_{\infty})}\right)^{1/3}+\mathrm{P}\left(\max_{j,k}|\widehat{\Sigma}_{n,jk}-\Sigma_{jk}|>\delta\right)\right\\}$
$\displaystyle\quad{}\quad{}\quad{}\quad{}\quad\quad+\inf_{\eta>0}\left\\{\frac{\eta}{\sqrt{\mathrm{Var}(\|Z\|_{\infty})}}+\mathrm{P}\left(|\Theta_{n}|>\eta\right)\right\\},$

where $\lesssim$ hides an absolute constant independent of $n,d,M$, and the
distribution of the $X_{i}$’s.

######  Remark 4 (On the purpose of the random variables $(\Theta_{n})_{n\geq
1}$).

In applications $(\Theta_{n})_{n\geq 1}$ may be a higher-order approximation
error such as the remainder of a first-order Taylor approximation. For
concrete examples we refer to Corollary 2 and Theorem 10 in Giessing and Fan
(2023). The random variables $(\Theta_{n})_{n\geq 1}$ may also be taken
identical to zero. In this case, the expression in the last line in Lemma 5
vanishes.

Propositions 4 and 5 inherit the dimension- and weak variance-freeness from
Propositions 1 and 3. Since these propositions do not require lower bounds the
minimum eigenvalue of the estimate of $\Sigma$, we can always use the naive
sample covariance matrix
$\widehat{\Sigma}_{n}=n^{-1}\sum_{i=1}^{n}(X_{i}-\bar{X}_{n})(X_{i}-\bar{X}_{n})^{\prime}$
with $\bar{X}_{n}=n^{-1}\sum_{i=1}^{n}X_{i}$ to estimate $\Sigma$ even if
$d\gg n$. If additional information about $\Sigma$ is available (viz. low-
rank, bandedness, or approximate sparsity), we can of course use more
sophisticated estimators $\widehat{\Sigma}_{n}$ to improve the non-asymptotic
bounds. This will prove particularly effective when we lift Propositions 4 and
5 to general empirical processes (see Section 3.4). The reader can find
several more examples in Section 4.2 of the companion paper Giessing and Fan
(2023).

###  2.2 Relation to previous work

Here, we compare Propositions 1 and 3 with the relevant results in the
literature. From a mathematical point of view Propositions 4 and 5 are just an
afterthought. For a comprehensive review of the entire literature on Gaussian
and bootstrap approximations of maxima of sums of random vectors see
Chernozhukov et al. (2023).

  * •

_On the dependence on the dimension._ As emphasized above, the upper bounds in
Propositions 1 and 3 are dimension-free. This is a significant improvement
over Theorems 2.2, 2.1, 3.2, 2.1 in Chernozhukov et al. (2013, 2017, 2019,
2020), Theorem 1.1 and Corollary 1.3 in Fang and Koike (2021), and Theorems
2.1 and 2.2 in Lopes (2022a), which all feature logarithmic factors of the
dimension. Such bounds generalize poorly to empirical processes since the
resulting upper bounds necessarily depend on the $\varepsilon$-entropy of the
function class. This precludes (or, at the very least, substantially
complicates) applications to objects as basic as the operator norm of a high-
dimensional covariance matrix.

  * •

_On the moment assumptions._ The above mentioned results by Chernozhukov et
al. (2013, 2017, 2019, 2020), Fang and Koike (2021), and Lopes (2022a) all
require strictly positive lower bounds on the variances of the components of
the random vector $X$ and/ or on the minimum eigenvalue of the covariance
matrix $\Sigma$. The strictly positive lower bounds are especially awkward if
we try to extend their finite dimensional results to general empirical
processes: While is often sensible to impose an upper bound on the variances
of the increments of an empirical process, lower bounds on the variances are
much harder to justify and in general only achievable via a discretization (or
thinning) of the function class. The weak variance-free upper bounds in
Propositions 1 and 3 give us greater leeway and increases the scope of our
approximation results considerably (see also Appendix A).

  * •

_On the sharpness of the upper bounds._ The results in Chernozhukov et al.
(2020), Fang and Koike (2021), and Lopes (2022a) show that the upper bound in
Proposition 1 is sub-optimal and that its dependence on sample size $n$ can be
improved to $n^{-1/2}$. The proof techniques in Chernozhukov et al. (2020) and
Fang and Koike (2021) are both based on delicate estimates of Hermite
polynomials and thus inherently dimension dependent. Extending their
approaches to the coordinate-free Wiener chaos decomposition (which would
yield dimension-free results) is a formidable research task. The proof
strategy in Lopes (2022a) is very different and requires sub-Gaussian data.
The results in the aforementioned papers also show that under additional
distributional assumptions the exponent on the upper bound in Proposition 3
can be improved to $1/2$.

  * •

_Proper generalizations of classical CLTs._ In a certain regard, Propositions
1 and 3 are strictly weaker than the results in Chernozhukov et al. (2017,
2020), Fang and Koike (2021), and Lopes (2022a). Proposition 1 (and similarly
Proposition 3) provides a bound on

$\displaystyle\sup_{s\geq 0}\Big{|}\mathbb{P}\left\\{\|S_{n}\|_{\infty}\leq
s\right\\}-\mathbb{P}\left\\{\|Z\|_{\infty}\leq
s\right\\}\Big{|}=\sup_{A\in\mathcal{Q}_{d}}\Big{|}\mathbb{P}\left\\{S_{n}\in
A\right\\}-\mathbb{P}\left\\{Z\in A\right\\}\Big{|},$

where $\mathcal{Q}_{d}$ be the collection of all hypercubes in
$\mathbb{R}^{d}$ with center at the origin. Chernozhukov et al. (2017), Fang
and Koike (2021), and Lopes (2022a) provide bounds on above quantity not for
the supremum over $\mathcal{Q}_{d}$ but the supremum over the much larger
collection of all hyper-rectangles in $\mathbb{R}^{d}$. When the dimension $d$
is fixed, their results imply convergence in distribution of $S_{n}$ to $Z$ as
$n\rightarrow\infty$ and are therefore stronger than our results (see also
Remark 3). In particular, their results can be considered proper
generalizations of classical CLTs to high dimensions.

The results in this paper depend on a dimension-free anti-concentration
inequality for $\|Z\|_{\infty}$ (or, in other words, for the supremum over
$\mathcal{Q}_{d}$) and a smoothing inequality specific to the map
$x\mapsto\|x\|_{\infty}$ (i.e. Lemmas 2 and 6). Since these inequalities do
not apply to the class of all hyper-rectangles in $\mathbb{R}^{d}$, the
arguments in this paper cannot be easily modified to yield dimension-free
generalizations of the classical CLTs to high dimensions.

##  3 Approximation results for suprema of empirical processes

###  3.1 Empirical process notation and definitions

In the previous section we got away with intuitive standard notation. We now
introduce precise empirical process notation for the remainder of the paper:
We denote by $X,X_{1},X_{2},\ldots$ a sequence of i.i.d. random variables
taking values in a measurable space $(S,\mathcal{S})$ with common distribution
$P$, i.e. $X_{i}:S^{\infty}\rightarrow S$, $i\geq 1$, are the coordinate
projections of the infinite product probability space
$(\Omega,\mathcal{A},\mathbb{P})=\left(S^{\infty},\mathcal{S}^{\infty},P^{\infty}\right)$
with law $\mathbb{P}_{X_{i}}=P$. If auxiliary variables independent of the
$X_{i}$’s are involved, the underlying probability space is assumed to be of
the form
$(\Omega,\mathcal{A},\mathbb{P})=\left(S^{\infty},\mathcal{S}^{\infty},P^{\infty}\right)\times(Z,\mathcal{Z},Q)$.
We define the empirical measures $P_{n}$ associated with observations
$X_{1},\ldots,X_{n}$ as random measures on
$\left(S^{\infty},\mathcal{S}^{\infty}\right)$ given by
$P_{n}(\omega):=n^{-1}\sum_{i=1}^{n}\delta_{X_{i}(\omega_{i})}$ for all
$\omega\in S^{\infty}$, where $\delta_{x}$ is the Dirac measure at $x$.

For a class $\mathcal{F}$ of measurable functions from $S$ onto the real line
$\mathbb{R}$ we define the _empirical process indexed by $\mathcal{F}$_ as

$\displaystyle\mathbb{G}_{n}(f):=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(f(X_{i})-Pf\big{)},\quad{}f\in\mathcal{F}.$

Further, we denote by $\\{G_{P}(f):f\in\mathcal{F}\\}$ the _Gaussian
$P$-bridge process_ with mean zero and the same covariance function
$\mathcal{C}_{P}:\mathcal{F}\times\mathcal{F}\rightarrow\mathbb{R}$ as the
process $\\{f(X):f\in\mathcal{F}\\}$, i.e.

$\displaystyle(f,g)\mapsto\mathcal{C}_{P}(f,g):=\mathrm{E}[G_{P}(f)G_{P}(g)]=Pfg-(Pf)(Pg).$
(2)

Moreover, we denote by $\\{Z_{Q}(f):f\in\mathcal{F}\\}$ the _Gaussian
$Q$-motion_ with mean zero and covariance function covariance function
$\mathcal{E}_{Q}:\mathcal{F}\times\mathcal{F}\rightarrow\mathbb{R}$ given by

$\displaystyle(f,g)\mapsto\mathcal{E}_{Q}(f,g):=\mathrm{E}[Z_{Q}(f)Z_{Q}(g)]=Qfg.$
(3)

For probability measures $Q$ on $(S,\mathcal{S})$ we define the
$L_{q}(Q)$-norm, $q\geq 1$, for function $f\in\mathcal{F}$ by
$\|f\|_{Q,q}=(Q|f|^{q})^{1/q}$, the $L_{2}(Q)$-semimetric by
$e_{Q}(f,g):=\|f-g\|_{Q,2}$ and the intrinsic standard deviation metric by
$d_{P}(f,g):=e_{P}(f-Pf,g-Pg)=\sqrt{P(f-g)^{2}-(P(f-g))^{2}}$,
$f,g\in\mathcal{F}$. We denote by $L_{q}(S,\mathcal{S},Q)$, $q\geq 1$, the
space of all real-valued measurable functions $f$ on $(S,\mathcal{S})$ with
finite $L_{q}(Q)$-norm. A function $F:S\rightarrow\mathbb{R}$ is an _envelop_
for the class $\mathcal{F}$ satisfies $|f(x)|\leq F(x)$ for all
$f\in\mathcal{F}$ and $x\in S$. For a semimetric space $(T,d)$ and any
$\varepsilon>0$ we write $N(T,d,\varepsilon)$ to denote the
$\varepsilon$_-covering numbers_ of $T$ with respect to $d$. For two
deterministic sequences $(a_{n})_{n\geq 1}$ and $(b_{n})_{n\geq 1}$, we write
$a_{n}\lesssim b_{n}$ if $a_{n}=o(b_{n})$ and $a_{n}\asymp b_{n}$ if there
exist absolute constants $C_{1},C_{2}>0$ such that $C_{1}b_{n}\leq a_{n}\leq
C_{2}b_{n}$ for all $n\geq 1$. For $a,b\in\mathbb{R}$ we write $a\wedge b$ for
$\min\\{a,b\\}$ and $a\vee b$ for $\max\\{a,b\\}$.

###  3.2 Gaussian approximation inequalities

We present the first main theoretical result of this paper: a Gaussian
approximation inequality for empirical processes indexed by function classes
that need not be Donsker. This result generalizes Proposition 1 from finite
dimensional index sets to general function classes.

######  Theorem 1 (Gaussian approximation).

Let $(\mathcal{F}_{n},\rho)$ be a totally bounded pseudo-metric space.
Further, let $\mathcal{F}_{n}\subset L_{2}(S,\mathcal{S},P)$ have envelope
function $F_{n}\in L_{3}(S,\mathcal{S},P)$. Suppose that there exist functions
$\psi_{n},\phi_{n}$ such that for $\delta>0$,

$\displaystyle\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n,\delta}^{\prime}}\lesssim\psi_{n}(\delta)\sqrt{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}\quad{}\quad{}\text{and}\quad{}\quad{}\big{\|}\|\mathbb{G}_{n}\|_{\mathcal{F}_{n,\delta}^{\prime}}\big{\|}_{P,1}\lesssim\phi_{n}(\delta)\sqrt{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})},$
(4)

where
$\mathcal{F}_{n,\delta}^{\prime}=\\{f-g:f,g\in\mathcal{F}_{n},\>\rho(f,g)<\delta\|F_{n}\|_{P,2}\\}$.
Let
$r_{n}=\inf\big{\\{}\psi_{n}(\delta)\vee\phi_{n}(\delta):\delta>0\big{\\}}$.
Then, for each $M\geq 0$,

$\displaystyle\sup_{s\geq
0}\Big{|}\mathbb{P}\left\\{\|\mathbb{G}_{n}\|_{\mathcal{F}_{n}}\leq
s\right\\}-\mathbb{P}\left\\{\|G_{P}\|_{\mathcal{F}_{n}}\leq
s\right\\}\Big{|}$
$\displaystyle\quad{}\quad{}\quad\lesssim\frac{\|F_{n}\|_{P,3}}{\sqrt{n^{1/3}\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}+\frac{\|F_{n}\mathbf{1}\\{F_{n}>M\\}\|_{P,3}^{3}}{\|F_{n}\|_{P,3}^{3}}+\frac{\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}+M}{\sqrt{n\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}+\sqrt{r_{n}},$

where $\lesssim$ hides an absolute constant independent of
$n,r_{n},M,\mathcal{F}_{n},F_{n},\psi_{n},\phi_{n}$, and $P$.

######  Remark 1 (Lower and upper bounds on
$\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}}$).

Since $\mathcal{F}_{n}$ may change with the sample size, so may
$\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}}$. In a few special cases it is
possible to exactly compute the variance (Giessing and Fan, 2023), but in most
cases one can only derive lower and upper bounds. In our companion work
Giessing (2023) we derive such bounds under mild conditions on the Gaussian
$P$-bridge process. For the reader’s convenience we include the relevant
results from this paper in Appendix B.2.

The total boundedness of $\mathcal{F}_{n}$ in above theorem is a standard
assumption which allows us to reduce the proof of Theorem 1 to an application
of Proposition 1 combined with a discretization argument. Since
$\mathcal{F}_{n}$ is totally bounded whenever
$\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}<\infty$ (i.e. Lemma 16 in Appendix
B.4), this is a rather mild technical assumption. For a detailed discussion of
this result and a comparison with the literature we refer to Section 3.6.

Under more specific smoothness assumptions on the Gaussian $P$-bridge process
(beyond the abstract control of the moduli of continuity in eq. (4)) Theorem 1
simplifies as follows:

######  Corollary 2.

Let $\mathcal{F}_{n}\subset L_{2}(S,\mathcal{S},P)$ and have envelope
$F_{n}\in L_{3}(S,\mathcal{S},P)$. If the Gaussian $P$-Bridge process
$\\{G_{P}(f):f\in\mathcal{F}_{n}\\}$ has almost surely uniformly
$d_{P}$-continuous sample paths and
$\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}<\infty$, then for each $M\geq 0$,

$\displaystyle\sup_{s\geq
0}\Big{|}\mathbb{P}\left\\{\|\mathbb{G}_{n}\|_{\mathcal{F}_{n}}\leq
s\right\\}-\mathbb{P}\left\\{\|G_{P}\|_{\mathcal{F}_{n}}\leq
s\right\\}\Big{|}$
$\displaystyle\quad{}\quad{}\quad{}\lesssim\frac{\|F_{n}\|_{P,3}}{\sqrt{n^{1/3}\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}+\frac{\|F_{n}\mathbf{1}\\{F_{n}>M\\}\|_{P,3}^{3}}{\|F_{n}\|_{P,3}^{3}}+\frac{\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}+M}{\sqrt{n\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}},$

where $\lesssim$ hides an absolute constant independent of
$n,M,\mathcal{F}_{n},F_{n}$, and $P$.

######  Remark 2 (Further simplification with lower bound on second moment).

Recall the setup of Corollary 2. If in addition there exists $\kappa_{n}>0$
such that $Pf^{2}\geq\kappa_{n}^{2}>0$ for all $f\in\mathcal{F}_{n}$ and
$F_{n}\in L_{3+\delta}(S,\mathcal{S},P)$, $\delta>0$, then

$\displaystyle\sup_{s\geq
0}\Big{|}\mathbb{P}\left\\{\|\mathbb{G}_{n}\|_{\mathcal{F}_{n}}\leq
s\right\\}-\mathbb{P}\left\\{\|G_{P}\|_{\mathcal{F}_{n}}\leq
s\right\\}\Big{|}$
$\displaystyle\quad\quad\quad\lesssim\frac{1}{n^{1/6}\kappa^{2}_{n}}\left(\|F_{n}\|_{P,3+\delta}^{3+\delta}+\mathrm{E}[\|G_{P}\|_{\mathcal{F}_{n}}^{3+\delta}]\right)^{\frac{1}{3+\delta}}\mathrm{E}[\|G_{P}\|_{\mathcal{F}_{n}}]+\frac{1}{n^{\delta/3}}\left(\frac{\|F_{n}\|_{P,3+\delta}}{\|F_{n}\|_{P,3}}\right)^{3},$

where $\lesssim$ hides an absolute constant independent of
$n,d,\kappa_{n}^{2}$, and the distribution of the $X_{i}$’s. This is the
empirical process analogue to Corollary 2. Since the proof of this inequality
is is identical to the one of Corollary 2 and only differs in notation, we
omit the details.

Since centered Gaussian processes are either almost surely uniformly
continuous w.r.t. their intrinsic standard deviation metric or almost surely
discontinuous, the smoothness condition in Corollary 2 is natural. (For the
proofs to go through uniform $d_{P}$-continuity in probability would be
sufficient.) In Lemmas 17 and 18 in Appendix B.4 we provide simple sufficient
and necessary conditions for almost sure uniform $d_{P}$-continuity of $G_{P}$
on $\mathcal{F}_{n}$. Importantly, $n\geq 1$ is fixed; the uniform
$d_{P}$-continuity is a point-wise requirement and does not need to hold when
taking the limit $n\rightarrow\infty$. This is especially relevant in high-
dimensional statistics where Gaussian processes are typically unbounded and
have diverging sample paths as $n\rightarrow\infty$, and, more generally,
whenever the function classes are non-Donsker. Of course, just as with
Proposition 1 the conclusion of Corollary 2 is weaker than a weaker
convergence (see Section 2.2).

###  3.3 Gaussian comparison inequalities

Our second main result is a comparison inequality which bounds the Kolmogorov
distance between a Gaussian $P$-bridge process and a Gaussian $Q$-motion both
indexed by (possibly) different function classes $\mathcal{F}_{n}$ and
$\mathcal{G}_{n}$. This generalizes Proposition 3 to Gaussian processes.

To state the theorem, we introduce the following notation: Given the pseudo-
metric spaces $(\mathcal{F}_{n},\rho_{i})$ and $(\mathcal{G}_{n},\rho_{i})$,
$i\in\\{1,2\\}$, we write
$\pi^{1}_{\mathcal{G}_{n}}:\mathcal{F}_{n}\cup\mathcal{G}_{n}\rightarrow\mathcal{G}_{n}$
to denote the projection from $\mathcal{F}_{n}\cup\mathcal{G}_{n}$ onto
$\mathcal{G}_{n}$ defined by
$\rho_{1}(h,\pi^{1}_{\mathcal{G}_{n}}h)=\inf_{g\in\mathcal{G}_{n}}\rho_{1}(h,g)$
for all $h\in\mathcal{F}_{n}\cup\mathcal{G}_{n}$. In the case that the image
$\pi^{1}_{\mathcal{G}_{n}}h$ is not a singleton for some
$h\in\mathcal{F}_{n}\cup\mathcal{G}_{n}$, we choose any one of the equivalent
points in $\mathcal{G}_{n}$. Similarly, we write
$\pi^{2}_{\mathcal{F}_{n}}:\mathcal{F}_{n}\cup\mathcal{G}_{n}\rightarrow\mathcal{F}_{n}$
for the projection from $\mathcal{F}_{n}\cup\mathcal{G}_{n}$ onto
$\mathcal{F}_{n}$ defined via $\rho_{2}$.

######  Theorem 3 (Gaussian comparison).

Let $(\mathcal{F}_{n},\rho_{1})$ and $(\mathcal{G}_{n},\rho_{2})$ be totally
bounded pseudo-metric spaces. Further, let $\mathcal{F}_{n}\subset
L_{2}(S,\mathcal{S},P)$ and $\mathcal{G}_{n}\subset L_{2}(S,\mathcal{S},Q)$
have envelope functions $F_{n}\in L_{2}(S,\mathcal{S},P)$ and $G_{n}\in
L_{2}(S,\mathcal{S},Q)$, respectively. Suppose that there exist functions
$\psi_{n},\phi_{n}$ such that

$\displaystyle\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n,\delta}^{\prime}}\lesssim\psi_{n}(\delta)\sqrt{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}\quad{}\quad{}\text{and}\quad{}\quad{}\mathrm{E}\|Z_{Q}\|_{\mathcal{G}_{n,\delta}^{\prime}}\lesssim\phi_{n}(\delta)\sqrt{\mathrm{Var}(\|Z_{Q}\|_{\mathcal{G}_{n}})},$

where
$\mathcal{F}_{n,\delta}^{\prime}=\\{f-g:f,g\in\mathcal{F}_{n},\>\rho_{1}(f,g)<\delta\|F_{n}\|_{P,2}\\}$
and
$\mathcal{G}_{n,\delta}^{\prime}=\\{f-g:f,g\in\mathcal{G}_{n},\>\rho_{2}(f,g)<\delta\|G_{n}\|_{Q,2}\\}$.
Let
$r_{n}=\inf\big{\\{}\psi_{n}(\delta)\vee\phi_{n}(\delta):\delta>0\big{\\}}$.
Then,

$\displaystyle\sup_{s\geq
0}\Big{|}\mathbb{P}\left\\{\|G_{P}\|_{\mathcal{F}_{n}}\leq
s\right\\}-\mathbb{P}\left\\{\|Z_{Q}\|_{\mathcal{G}_{n}}\leq
s\right\\}\Big{|}\lesssim\left(\frac{\Delta_{P,Q}(\mathcal{F}_{n},\mathcal{G}_{n})}{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})\vee\mathrm{Var}(\|Z_{Q}\|_{\mathcal{G}_{n}})}\right)^{1/3}+\sqrt{r_{n}},$

where

$\displaystyle\begin{split}\Delta_{P,Q}(\mathcal{F}_{n},\mathcal{G}_{n})&:=\sup_{f_{1},f_{2}\in\mathcal{F}_{n}}\big{|}\mathcal{C}_{P}(f_{1},f_{2})-\mathcal{E}_{Q}(\pi^{1}_{\mathcal{G}_{n}}f_{1}),\pi^{1}_{\mathcal{G}_{n}}f_{2})\big{|}\\\
&\quad\quad\quad\bigvee\sup_{g_{1},g_{2}\in\mathcal{G}_{n}}\big{|}\mathcal{E}_{Q}(g_{1},g_{2})-\mathcal{C}_{P}(\pi^{2}_{\mathcal{F}_{n}}g_{1},\pi^{2}_{\mathcal{F}_{n}}g_{2})\big{|},\end{split}$
(5)

where $\lesssim$ hides an absolute constant independent of
$n,\mathcal{F}_{n},\mathcal{G}_{n},F_{n},G_{n},\psi_{n},\phi_{n},P$ and $Q$.

######  Remark 3 (Comparison inequality for Gaussian $Q$-bridge processes).

Theorem 3 holds without changes (and identical proof) for Gaussian $Q$-bridge
processes $\\{G_{Q}(f):f\in\mathcal{F}_{n}\\}$. This easily follows from the
almost sure representation $G_{Q}(f)=Z_{Q}(f)+(Pf)Z$ for all
$f\in\mathcal{F}_{n}$, where $Z\sim N(0,1)$ is independent of $Z_{Q}(f)$ for
all $f\in\mathcal{F}_{n}$. We state the comparison inequality for $Q$-motions
because this is the version that we use in the context of the Gaussian process
bootstrap in Section 3.4.

Informally speaking, Theorem 3 states that the distributions of the suprema of
any two Gaussian processes are close whenever their covariance functions are
not too different. Note that we compare two Gaussian processes with different
covariance functions $\mathcal{C}_{P}$ and $\mathcal{E}_{Q}$ that differ not
only in their measures ($P$ and $Q$) but also their functional form (viz. eq.
(2) and (3) in Section 3.1). This turns out to be essential for developing
alternatives to the classical Gaussian multiplier bootstrap procedure for
suprema of empirical processes proposed by Chernozhukov et al. (2014a) and
Chernozhukov et al. (2014b). We discuss this in detail in Section 3.4.

Under additional smoothness conditions on the Gaussian processes and
assumptions on the function classes Theorem 3 simplifies. The special cases
$\mathcal{G}_{n}=\mathcal{F}_{n}$ and
$\mathcal{G}_{n}\subseteq\mathcal{F}_{n}$ are of particular interest in the
context of the bootstrap. They are the content of the next two corollaries.

######  Corollary 4.

Let $\mathcal{F}_{n}\subset L_{2}(S,\mathcal{S},P)\cap
L_{2}(S,\mathcal{S},Q)$. If the Gaussian $P$-Bridge process
$\\{G_{P}(f):f\in\mathcal{F}_{n}\\}$ and the Gaussian $Q$-motion
$\\{Z_{Q}(f):f\in\mathcal{F}_{n}\\}$ have almost surely uniformly continuous
sample paths w.r.t. their respective standard deviation metrics $d_{P}$ and
$d_{Q}$, and
$\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}\vee\mathrm{E}\|Z_{Q}\|_{\mathcal{F}_{n}}<\infty$,
then

$\displaystyle\sup_{s\geq
0}\Big{|}\mathbb{P}\left\\{\|G_{P}\|_{\mathcal{F}_{n}}\leq
s\right\\}-\mathbb{P}\left\\{\|Z_{Q}\|_{\mathcal{F}_{n}}\leq
s\right\\}\Big{|}\lesssim\left(\frac{\sup_{f,g\in\mathcal{F}_{n}}\big{|}\mathcal{C}_{P}(f,g)-\mathcal{E}_{Q}(f,g)\big{|}}{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})\vee\mathrm{Var}(\|Z_{Q}\|_{\mathcal{F}_{n}})}\right)^{1/3},$

where $\lesssim$ hides an absolute constant independent of
$n,\mathcal{F}_{n},P$, and $Q$.

######  Corollary 5.

Let $\mathcal{F}_{n}\subset L_{2}(S,\mathcal{S},P)$ and
$\mathcal{G}_{n}\subseteq\mathcal{F}_{n}$ be a $\delta\|F\|_{P,2}$-net of
$\mathcal{F}_{n}$ with respect to $d_{P}$, $\delta>0$. If the Gaussian
$P$-Bridge process $\\{G_{P}(f):f\in\mathcal{F}_{n}\\}$ has almost surely
uniformly continuous sample paths w.r.t. metric $d_{P}$ and
$\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}<\infty$, then

$\displaystyle\sup_{s\geq
0}\Big{|}\mathbb{P}\left\\{\|G_{P}\|_{\mathcal{F}_{n}}\leq
s\right\\}-\mathbb{P}\left\\{\|G_{P}\|_{\mathcal{G}_{n}}\leq
s\right\\}\Big{|}\lesssim\left(\frac{\delta\|F\|_{P,2}\sup_{f\in\mathcal{F}_{n}}\sqrt{Pf^{2}}}{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})\vee\mathrm{Var}(\|G_{P}\|_{\mathcal{G}_{n}})}\right)^{1/3},$

where $\lesssim$ hides an absolute constant independent of
$\delta,n,\mathcal{F}_{n}$, and $P$.

###  3.4 Gaussian process bootstrap

In this section we develop a general framework for bootstrapping the
distribution of suprema of empirical processes. We recover the classical
Gaussian multiplier bootstrap as a special case of this more general
framework. For concrete applications to statistical problems we refer to
Section 4.

The following abstract approximation result generalizes Proposition 4 to
empirical processes.

######  Theorem 6 (Abstract bootstrap approximation).

Let $(\mathcal{F}_{n},\rho_{1})$ and $(\mathcal{G}_{n},\rho_{2})$ be totally
bounded pseudo-metric spaces. Further, let $\mathcal{F}_{n}\subset
L_{2}(S,\mathcal{S},P)$ and $\mathcal{G}_{n}\subset L_{2}(S,\mathcal{S},Q)$
have envelope functions $F_{n}\in L_{3}(S,\mathcal{S},P)$ and $G_{n}\in
L_{2}(S,\mathcal{S},Q)$, respectively. Suppose that there exist functions
$\psi_{n},\phi_{n}$ such that for $\delta>0$,

$\displaystyle\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n,\delta}^{\prime}}\vee\big{\|}\|\mathbb{G}_{n}\|_{\mathcal{F}_{n,\delta}^{\prime}}\big{\|}_{P,1}\lesssim\psi_{n}(\delta)\sqrt{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})},\quad\mathrm{E}\|Z_{Q}\|_{\mathcal{G}_{n,\delta}^{\prime}}\lesssim\phi_{n}(\delta)\sqrt{\mathrm{Var}(\|Z_{Q}\|_{\mathcal{G}_{n}})},$
(6)

where
$\mathcal{F}_{n,\delta}^{\prime}=\\{f-g:f,g\in\mathcal{F}_{n},\>\rho_{1}(f,g)<\delta\|F_{n}\|_{P,2}\\}$
and
$\mathcal{G}_{n,\delta}^{\prime}=\\{f-g:f,g\in\mathcal{G}_{n},\>\rho_{2}(f,g)<\delta\|G_{n}\|_{Q,2}\\}$.
Let
$r_{n}=\inf\big{\\{}\psi_{n}(\delta)\vee\phi_{n}(\delta):\delta>0\big{\\}}$.
Then, for each $M\geq 0$,

$\displaystyle\sup_{s\geq
0}\Big{|}\mathbb{P}\left\\{\|\mathbb{G}_{n}\|_{\mathcal{F}_{n}}\leq
s\right\\}-\mathbb{P}\left\\{\|Z_{Q}\|_{\mathcal{G}_{n}}\leq
s\right\\}\Big{|}$
$\displaystyle\quad{}\quad{}\quad{}\lesssim\frac{\|F_{n}\|_{P,3}}{\sqrt{n^{1/3}\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}+\frac{\|F_{n}\mathbf{1}\\{F_{n}>M\\}\|_{P,3}^{3}}{\|F_{n}\|_{P,3}^{3}}+\frac{\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}+M}{\sqrt{n\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}$
$\displaystyle\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}+\left(\frac{\Delta_{P,Q}(\mathcal{F}_{n},\mathcal{G}_{n})}{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})\vee\mathrm{Var}(\|Z_{Q}\|_{\mathcal{G}_{n}})}\right)^{1/3}+\sqrt{r_{n}},$

where $\Delta_{P,Q}(\mathcal{F}_{n},\mathcal{G}_{n})$ is given in eq. (5) and
$\lesssim$ hides an absolute constant independent of $n,r_{n},M,$
$\mathcal{F}_{n},$ $\mathcal{G}_{n},F_{n},G_{n},\psi_{n},\phi_{n},P$, and $Q$.

######  Remark 4 (Bootstrap approximation with Gaussian $Q$-bridge process).

Theorem 6 holds without changes (and identical proof) with the Gaussian
$Q$-bridge process $\\{G_{Q}(f):f\in\mathcal{F}_{n}\\}$ substituted for the
Gaussian $Q$-motion $\\{Z_{Q}(f):f\in\mathcal{F}_{n}\\}$ (see also Remark 3).

Theorem 6 is an immediate consequence of Theorems 1 and 3 and the triangle
inequality, and thus a mathematical triviality. What matters is its
statistical interpretation: It implies that we can approximate the sampling
distribution of the supremum of an empirical process
$\\{\mathbb{G}_{n}(f):f\in\mathcal{F}_{n}\\}$ with the distribution of the
supremum of a Gaussian $Q$-motion $\\{Z_{Q}(g):g\in\mathcal{G}_{n}\\}$
provided that their covariance functions $\mathcal{C}_{P}$ and
$\mathcal{E}_{Q}$ do not differ by too much. Importantly, up to the smoothness
condition in (6), we are completely free to choose whatever $Q$-measure and
function class $\mathcal{G}_{n}$ suit us best.

This interpretation of Theorem 6 motivates the following bootstrap procedure
for estimating the distribution of the supremum of an empirical process
$\\{\mathbb{G}_{n}(f):f\in\mathcal{F}_{n}\\}$ when the $P$-measure is unknown.

######  Algorithm 1 (Gaussian process bootstrap).

Let $X_{1},\ldots,X_{n}$ be a simple random sample drawn from distribution
$P$.

  * Step 1:

Construct a positive semi-definite estimate $\widehat{\mathcal{C}}_{n}$ of the
covariance function $\mathcal{C}_{P}$ of the empirical process
$\\{\mathbb{G}_{n}(f):f\in\mathcal{F}_{n}\\}$.

  * Step 2:

Construct a Gaussian process
$\widehat{Z}_{n}=\\{\widehat{Z}_{n}(f):f\in\mathcal{F}_{n}\\}$ such that for
all $f,g\in\mathcal{F}_{n}$,

$\displaystyle\mathrm{E}\big{[}\widehat{Z}_{n}(f)\mid
X_{1},\ldots,X_{n}\big{]}=0\quad\mathrm{and}\quad\mathrm{E}\big{[}\widehat{Z}_{n}(f)\widehat{Z}_{n}(g)\mid
X_{1},\ldots,X_{n}\big{]}=\widehat{\mathcal{C}}_{n}(f,g).$ (7)

  * Step 3:

Approximate the distribution of $\|\mathbb{G}_{n}\|_{\mathcal{F}_{n}}$ by
drawing Monte Carlo samples from $\|\widehat{Z}_{n}\|_{\mathcal{F}_{n}}$.

By Kolmogorov’s consistency theorem the Gaussian process $\widehat{Z}_{n}$
will always exist. If the estimate $\widehat{\mathcal{C}}_{n}$ is uniformly
consistent and if the envelope function $F_{n}$ and the supremum of the
Gaussian $P$-bridge process satisfy certain moment conditions, Theorem 6
readily implies uniform consistency of the Gaussian process bootstrap as
$n\rightarrow\infty$. The challenge is to actually construct a (measurable)
version of $\widehat{Z}_{n}$ from which we can draw Monte Carlo samples. This
is the content the next section.

######  Remark 5 (Relation to the _Gaussian multiplier bootstrap_ by
Chernozhukov et al. (2014b)).

The Gaussian multiplier bootstrap is a special case of the Gaussian process
bootstrap. To see this, note that if we estimate the covariance function
$\mathcal{C}_{P}$ nonparametrically via the sample covariance function
$\mathcal{C}_{P_{n}}(f,g)=P_{n}fg-(P_{n}f)(P_{n}g)$, $f,g\in\mathcal{F}_{n}$,
then the corresponding Gaussian $P_{n}$-Bridge process
$\\{G_{P_{n}}(f):f\in\mathcal{F}_{n}\\}$ can be expanded into a finite sum of
non-orthogonal functions

$\displaystyle
G_{P_{n}}(f)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\xi_{i}\big{(}f(X_{i})-P_{n}f\big{)},\quad{}f\in\mathcal{F}_{n},$

where $\xi_{1},\ldots,\xi_{n}$ are i.i.d. standard normal random variables
independent of $X_{1},\ldots,X_{n}$. While this representation is extremely
simple (in theory and practice!), unlike the Gaussian process bootstrap it
does not allow integrating additional structural information about the
covariance function. Empirically, this often results in less accurate
bootstrap estimates (Giessing and Fan, 2023).

###  3.5 A practical guide to the Gaussian process bootstrap

To implement the Gaussian process bootstrap we need two things: a uniformly
consistent estimate of the covariance function and a systematic way of
constructing Gaussian processes for given covariance functions.

Finding consistent estimates of the covariance function is relatively
straightforward. For example, a natural estimate is the nonparametric sample
covariance function defined as
$\mathcal{C}_{P_{n}}(f,g):=P_{n}(fg)-(P_{n}f)(P_{n}g)$,
$f,g\in\mathcal{F}_{n}$. Under suitable conditions on the probability measure
$P$ and the function class $\mathcal{F}_{n}$ this estimate is uniformly
consistent in high dimensions (e.g. Koltchinskii and Lounici, 2017). Of
course, an important feature of the Gaussian process bootstrap is that it can
be combined with any positive semi-definite estimate of the covariance
function. In particular, we can use (semi-)parametric estimators to exploit
structural constraints induced by $P$ and $\mathcal{F}_{n}$. For concrete
examples we refer to Section 4 and also Giessing and Fan (2023).

Constructing Gaussian processes for given covariance functions and defined on
potentially arbitrary index sets is a more challenging problem. We propose to
base the construction on an almost sure version of the classical Karhunen-
Loève decomposition.

To develop this decomposition in the framework of this paper, we need new
notation and concepts. In the following,
$(\mathcal{F}_{n},\mathcal{B}_{n},\mu)$ denotes a measurable space for some
finite measure $\mu$ with support $\mathcal{F}_{n}$. Typically,
$\mathcal{B}_{n}$ is the Borel $\sigma$-algebra on $\mathcal{F}_{n}$ and the
measure $\mu$ is chosen for convenience. For example, $\mu$ can be set to the
Lebesgue measure when $\mathcal{F}$ an interval, or the counting measure when
$\mathcal{F}_{n}$ is a finite (discrete) set. The space
$L_{2}(\mathcal{F}_{n},\mathcal{B}_{n},\mu)$ equipped with the inner product
$\langle\psi,\phi\rangle:=\int_{\mathcal{F}_{n}}\psi(f)\phi(f)d\mu(f)$ is a
Hilbert space. Given a positive semi-definite and continuous kernel
$K:\mathcal{F}_{n}\times\mathcal{F}_{n}\rightarrow\mathbb{R}$ we define the
linear operator $T_{K}$ on $L_{2}(\mathcal{F}_{n},\mathcal{B}_{n},\mu)$ via

$\displaystyle(T_{K}\psi)(g):=\langle
K(\cdot,g),\psi\rangle=\int_{\mathcal{F}_{n}}K(f,g)\psi(f)d\mu(f),\quad\psi\in
L^{2}(\mathcal{F}_{n},\mathcal{B}_{n},\mu),\quad g\in\mathcal{F}_{n}.$ (8)

If
$\int_{\mathcal{F}_{n}\times\mathcal{F}_{n}}K^{2}(f,g)d\mu(f)d\mu(g)<\infty$,
then $T_{K}$ is a bounded linear operator. If $(\mathcal{F}_{n},d_{Q})$ is a
compact metric space, then $T_{K}$ is a compact operator. In this case,
$T_{K}$ has a spectral decomposition which depends on the kernel $K$ alone;
the measure $\mu$ is exogenous. For further details we refer the reader to
Chapters 4 and 7 in Hsing and Eubank (2015).

The following result is well-known (Jain and Kallianpur, 1970); since it is
slightly adapted to our setting we have included its proof in the appendix.

######  Proposition 7 (Karhunen-Loève decomposition of Gaussian processes).

Let $(\mathcal{F}_{n},e_{Q})$ be a compact pseudo-metric space and
$Z_{Q}=\\{Z_{Q}(f):f\in\mathcal{F}_{n}\\}$ be a Gaussian $Q$-motion. Suppose
that $Z_{Q}$ has a continuous covariance function $\mathcal{E}_{Q}$ and
$\mathrm{E}\|Z_{Q}\|_{\mathcal{F}_{n}}<\infty$. Let
$\\{(\lambda_{k},\varphi_{k})\\}_{k=1}^{\infty}$ be the eigenvalue and
eigenfunction pairs of the linear operator $T_{\mathcal{E}_{Q}}$ corresponding
to $\mathcal{E}_{Q}$, and $\\{\xi_{k}\\}_{k=1}^{\infty}$ be a sequence of
i.i.d. standard normal random variables. Then, with probability one,

$\displaystyle\lim_{m\rightarrow\infty}\|Z_{Q}^{m}-Z_{Q}\|_{\mathcal{F}_{n}}=0,$

where $Z_{Q}^{m}(f):=\sum_{k=1}^{m}\xi_{k}\sqrt{\lambda_{k}}\varphi_{k}(f)$,
$f\in\mathcal{F}_{n}$, is an almost surely bounded and continuous Gaussian
process on $\mathcal{F}_{n}$ for all $m\in\mathbb{N}$.

This proposition justifies the following constructive approximation of the
Gaussian proxy process $\widehat{Z}_{n}$ defined in (7): Let
$\widehat{\mathcal{C}}_{n}$ be a generic estimate of the covariance function
$\mathcal{C}_{P}$ of the empirical process
$\\{\mathbb{G}_{n}(f):f\in\mathcal{F}_{n}\\}$. Suppose that
$\widehat{\mathcal{C}}_{n}$ is such that the associated integral operator
$T_{\widehat{\mathcal{C}}_{n}}$ defined in (8) admits a spectral
decomposition. Denote by
$\big{\\{}(\widehat{\lambda}_{k},\widehat{\varphi}_{k})\big{\\}}_{k=1}^{\infty}$
its eigenvalue and eigenfunction pairs, and, for concreteness, assume that the
$\widehat{\lambda}_{k}$’s are sorted in nonincreasing order. Given a sequence
of i.i.d. standard normal random variables $\\{\xi_{k}\\}_{k=1}^{\infty}$ and
truncation level $m\geq 1$, define the Gaussian process and associated
covariance function

$\displaystyle\widehat{Z}_{n}^{m}(f):=\sum_{k=1}^{m}\xi_{k}\sqrt{\widehat{\lambda}_{k}}\widehat{\varphi}_{k}(f),\quad\quad\widehat{\mathcal{C}}_{n}^{m}(f,g):=\sum_{k=1}^{m}\widehat{\lambda}_{k}\widehat{\varphi}_{k}(f)\widehat{\varphi}_{k}(g),\quad
f,g\in\mathcal{F}_{n}.$ (9)

While there are only a few situations in which the eigenvalues and
eigenfunctions of $T_{\widehat{\mathcal{C}}_{n}}$ can be found analytically,
from a computational perspective this is a standard problem and efficient
numerical solvers exist (Ghanem and Spanos, 2003; Berlinet and Thomas-Agnan,
2004). Thus, constructing the process $\widehat{Z}_{n}^{m}$ in (9) does not
pose any practical challenges. Proposition 7 now guarantees (under appropriate
boundedness and smoothness conditions on $\widehat{Z}_{n}$) that the
approximation error between $\widehat{Z}_{n}$ and $\widehat{Z}_{n}^{m}$ can be
made arbitrarily small by choosing $m\geq 1$ sufficiently large. In fact, by
Mercer’s theorem the error can be quantified in terms of the operator norm of
the difference between the covariance functions $\widehat{\mathcal{C}}_{n}$
and $\widehat{\mathcal{C}}_{n}^{m}$ (e.g. Hsing and Eubank, 2015, Lemma 4.6.6
and Corollary 8).

In above discussion we have implicitly imposed several assumptions on the
estimate $\widehat{\mathcal{C}}_{n}$. For future reference, we summarize these
assumptions in a single definition.

######  Definition 1 (Admissibility of $\widehat{\mathcal{C}}_{n}$).

We say that an estimate
$\widehat{\mathcal{C}}_{n}:\mathcal{F}_{n}\times\mathcal{F}_{n}\rightarrow\mathbb{R}$
of the covariance function $\mathcal{C}_{P}$ is admissible if it is
continuous, symmetric, positive semi-definite, and its associated integral
operator $T_{\widehat{\mathcal{C}}_{n}}$ is a bounded linear operator on
$L_{2}(\mathcal{F}_{n},\mathcal{B}_{n},\mu)$ for some finite measure $\mu$
with support $\mathcal{F}_{n}$.

######  Remark 6 (Admissible estimates exist).

Under the assumption that $\mathcal{F}_{n}$ has an envelope function, the
nonparametric sample covariance
$\mathcal{C}_{P_{n}}(f,g)=P_{n}(fg)-(P_{n}f)(P_{n}g)$ is admissible. Indeed,
by definition, $\mathcal{C}_{P_{n}}$ is continuous (w.r.t. $e_{P_{n}}$),
symmetric, positive semi-definite, and, by existence of the envelope function,
$\int_{\mathcal{F}_{n}\times\mathcal{F}_{n}}\mathcal{C}_{P_{n}}^{2}(f,g)dP_{n}(f)dP_{n}(g)<\infty$.
Also, the existence of an envelop function implies that
$(\mathcal{F}_{n},e_{P_{n}})$ is a compact metric space. See also Section 4.

The next result establishes consistency of the Gaussian process bootstrap
based on the truncated Karhunen-Loève expansion in (9). It is a simple
corollary of Theorem 6. Note that the intrinsic standard deviation metric
associated with (the pushforward probability measure induced by)
$\widehat{Z}_{n}^{m}$ can be expressed in terms of its covariance function as
$d_{\widehat{\mathcal{C}}_{n}^{m}}^{2}(f,g)=\widehat{\mathcal{C}}_{n}^{m}(f,f)+\widehat{\mathcal{C}}_{n}^{m}(g,g)-2\widehat{\mathcal{C}}_{n}^{m}(f,g)$.

######  Corollary 8 (Consistency of the Gaussian process bootstrap).

Let $\widehat{\mathcal{C}}_{n}$ be an admissible estimate of $\mathcal{C}_{P}$
and $\widehat{\mathcal{C}}_{n}^{m}$ its best rank-$m$ approximation. Let
$\mathcal{F}_{n}\subset L_{2}(S,\mathcal{S},P)$ have envelope $F_{n}\in
L_{3}(S,\mathcal{S},P)$. If the Gaussian processes
$\\{G_{P}(f):f\in\mathcal{F}_{n}\\}$ and
$\\{\widehat{Z}_{n}^{m}(f):f\in\mathcal{F}_{n}\\}$ have almost surely
uniformly continuous sample paths w.r.t. their respective standard deviation
metrics $d_{P}$ and $d_{\widehat{\mathcal{C}}_{n}^{m}}$, and
$\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}\vee\mathrm{E}\big{[}\|\widehat{Z}_{n}^{m}\|_{\mathcal{F}_{n}}\mid
X_{1},\ldots,X_{n}\big{]}<\infty$, then, for each $M\geq 0$,

$\displaystyle\sup_{s\geq
0}\Big{|}\mathbb{P}\left\\{\|\mathbb{G}_{n}\|_{\mathcal{F}_{n}}\leq
s\right\\}-\mathbb{P}\left\\{\|\widehat{Z}_{n}^{m}\|_{\mathcal{F}_{n}}\leq
s\mid X_{1},\ldots,X_{n}\right\\}\Big{|}$
$\displaystyle\quad{}\quad{}\quad{}\lesssim\frac{\|F_{n}\|_{P,3}}{\sqrt{n^{1/3}\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}+\frac{\|F_{n}\mathbf{1}\\{F_{n}>M\\}\|_{P,3}^{3}}{\|F_{n}\|_{P,3}^{3}}+\frac{\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}+M}{\sqrt{n\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}$
$\displaystyle\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}+\left(\frac{\sup_{f,g\in\mathcal{F}_{n}}\big{|}\mathcal{C}_{P}(f,g)-\widehat{\mathcal{C}}_{n}(f,g)\big{|}}{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}\right)^{1/3}+\left(\frac{\sup_{f,g\in\mathcal{F}_{n}}\big{|}\widehat{\mathcal{C}}_{n}(f,g)-\widehat{\mathcal{C}}_{n}^{m}(f,g)\big{|}}{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}\right)^{1/3},$

where $\lesssim$ hides an absolute constant independent of
$n,m,M,\mathcal{F}_{n},F_{n},P_{n}$, and $P$.

Moreover, if $\mathcal{F}_{n}$ is compact w.r.t.
$d_{\widehat{\mathcal{C}}_{n}}$, then the last term on the right hand side in
above display vanishes as $m\rightarrow\infty$.

In general, the strong variance $\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})$
depends on the dimension of the function class $\mathcal{F}_{n}$. Hence, the
truncation level $m$ has to be chosen (inverse) proportionate to
$\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})$ to ensure that in the upper bound
of Corollary 8 the deterministic approximation error
$\sup_{f,g\in\mathcal{F}_{n}}\big{|}\widehat{\mathcal{C}}_{n}(f,g)-\widehat{\mathcal{C}}_{n}^{m}(f,g)\big{|}$
is negligible compared to the stochastic estimation errors.

We conclude this section with a consistency result on the bootstrap
approximation of quantiles of suprema of empirical processes. This result is
the empirical process analogue to Proposition 5. It is relevant in the context
of hypothesis testing and construction of confidence regions. For
$\alpha\in(0,1)$ we denote the conditional $\alpha$-quantile of the supremum
of $\\{\widehat{Z}_{n}^{m}(f):f\in\mathcal{F}_{n}\\}$ by

$\displaystyle c_{n}(\alpha;\widehat{\mathcal{C}}_{n}^{m})$
$\displaystyle:=\inf\left\\{s\geq
0:\mathbb{P}\left\\{\|\widehat{Z}_{n}^{m}\|_{\mathcal{F}_{n}}\leq s\mid
X_{1},\ldots,X_{n}\right\\}\geq\alpha\right\\}.$

We have the following result:

######  Theorem 9 (Quantiles of the Gaussian process bootstrap).

Consider the setup of Corollary 8. Let $(\Theta_{n})_{n\geq 1}\in\mathbb{R}$
be a sequence of arbitrary random variables, not necessarily independent of
$X_{1},\ldots,X_{n}$. Then, for each $M\geq 0$,

$\displaystyle\sup_{\alpha\in(0,1)}\Big{|}\mathbb{P}\left\\{\|\mathbb{G}_{n}\|_{\mathcal{F}_{n}}+\Theta_{n}\leq
c_{n}(\alpha;\widehat{\mathcal{C}}_{n}^{m})\right\\}-\alpha\Big{|}$
$\displaystyle\quad{}\quad{}\quad{}\lesssim\frac{\|F_{n}\|_{P,3}}{\sqrt{n^{1/3}\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}+\frac{\|F_{n}\mathbf{1}\\{F_{n}>M\\}\|_{P,3}^{3}}{\|F_{n}\|_{P,3}^{3}}+\frac{\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n}}+M}{\sqrt{n\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}$
$\displaystyle\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}+\inf_{\delta>0}\left\\{\left(\frac{\delta}{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}\right)^{1/3}+\mathbb{P}\left\\{\sup_{f,g\in\mathcal{F}_{n}}\big{|}\mathcal{C}_{P}(f,g)-\widehat{\mathcal{C}}_{n}^{m}(f,g)\big{|}>\delta\right\\}\right\\}$
$\displaystyle\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}+\inf_{\eta>0}\left\\{\frac{\eta}{\sqrt{\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})}}+\mathbb{P}\left\\{|\Theta_{n}|>\eta\right\\}\right\\},$

where $\lesssim$ hides an absolute constant independent of
$n,m,M,\mathcal{F}_{n},F_{n},P_{n}$, and $P$.

In statistical applications, the statistic of interest is rarely a simple
empirical process. Instead, the empirical process usually arises as the
leading term of a (functional) Taylor expansion. The random sequence
$(\Theta_{n})_{n\geq 1}$ in above theorem can be used to capture the higher-
order approximation errors of such a expansion.

######  Remark 7 (Additional practical considerations).

It is often infeasible to draw Monte Carlo samples directly from
$\|Z_{n}^{m}\|_{\mathcal{F}_{n}}$. In practice, we suggest approximating
$\mathcal{F}_{n}$ via a finite $\delta\|F\|_{P,2}$-net
$\mathcal{G}_{n}\subseteq\mathcal{F}_{n}$ with respect to $d_{P}$. With this
additional approximation step and by Corollary 5, the conclusion of the
Corollary 8 holds with $\|\widehat{Z}_{n}^{m}\|_{\mathcal{G}_{n}}$ substituted
for $\|\widehat{Z}_{n}^{m}\|_{\mathcal{F}_{n}}$ and the additional quantity
$\big{(}\delta\|F\|_{P,2}\sup_{f\in\mathcal{F}_{n}}\sqrt{Pf^{2}}\big{)}^{1/3}\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})^{-1/3}$
in the upper bound on the Kolmogorov distance. This shows that the level of
discretization $\delta>0$ should be chosen proportional to
$\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})$.

###  3.6 Relation to previous work

The papers by Chernozhukov et al. (2014a, b, 2016) are currently the only
other existing works on Gaussian and bootstrap approximations of suprema of
empirical processes indexed by (potentially) non-Donsker function classes. In
this section, we compare their results to our Theorems 1–9. To keep the
comparison we focus on the key aspects that motivated us to write this paper.

It is important to note that the results presented in Chernozhukov et al.
(2014a, b, 2016) differ slightly in nature from ours. Instead of establishing
bounds on Kolmogorov distances, Chernozhukov and his co-authors derive
coupling inequalities. However, through standard arguments (Strassen’s theorem
and anti-concentration) these coupling inequalities indirectly yield bounds on
Kolmogorov distances. These implied bounds on the Kolmogorov distances are at
most as sharp as the ones in the coupling inequalities, up to multiplicative
constants. Since we do not care about absolute constants, but only the
dependence of the upper bounds on characteristics of the function class
$\mathcal{F}_{n}$ and the law $P$, we can meaningfully compare their findings
and ours.

  * •

_Unbounded function classes._ In classical empirical process theory the
function class $\mathcal{F}_{n}$ is typically assumed to be uniformly bounded,
i.e. $\sup_{x\in S}|f(x)|<\infty$ for all $f\in\mathcal{F}_{n}$ (e.g. van der
Vaart and Wellner, 1996). A key feature of Theorems 1–9 as well as Theorems
A.2, 2.1, and 2.1-2.3 in Chernozhukov et al. (2014a, b, 2016) is that they
hold for unbounded function classes and only require the envelope function
$F_{n}$ to have finite moments. Relaxing the uniform boundedness of the
function class is useful, among other things, for inference on high-
dimensional statistical models, functional data analysis, nonparametric
regression, and series/ sieve estimation (Chernozhukov et al., 2014a, b;
Giessing and Fan, 2023).

  * •

_Entropy conditions._ The upper bounds provided in Theorems A.2, 2.1, and
2.1-2.3 in Chernozhukov et al. (2014a, b, 2016) depend on a combination of
(truncated) second and $q$th ($q\geq 3$) moments of $\max_{1\leq i\leq
n}F_{n}$ $(X_{i})$, “local quantities” of order
$\sup_{f\in\mathcal{F}_{n}}P|f|^{3}$ and
$\sup_{f\in\mathcal{F}_{n}}\sqrt{Pf^{4}}$ (disregarding $(\log n)$-factors),
$\mathrm{E}\|G_{P}\|_{\mathcal{F}_{n,\delta}^{\prime}}\vee\big{\|}\|\mathbb{G}_{n}\|_{\mathcal{F}_{n,\delta}^{\prime}}\big{\|}_{P,1}$,
and the entropy number $\log N(\mathcal{F},e_{P},\\\ \delta\|F_{n}\|_{P,2})$,
$\delta>0$ arbitrary. These upper bounds are not only more complex than ours
but also weaker in terms of their implications: Since metric entropy with
respect to intrinsic standard deviation metric $e_{P}$ typically scales
linearly in the (VC-)dimension of the statistical model (see Appendix A for
two (counter-)examples), the upper bounds in Chernozhukov et al. (2014a, b,
2016) are vacuous in high-dimensional situations without sparsity or when the
(VC-)dimension exceeds the sample size. In contrast, the upper bounds in our
Theorems 1–9 depend only on the expected value and standard deviation of the
supremum of the Gaussian $P$-bridge process and the (truncated) third moments
of the envelope. Under mild assumptions, these quantities can be upper bounded
independently of the (VC-)dimension, thus offering useful bounds even in high-
dimensional problems (see Section 4).

The entropy term in Chernozhukov et al. (2014a, b, 2016) is due to their
proofs relying on a dimension-dependent version of Proposition 1.

  * •

_Lower bounds on weak variances._ Lemmas 2.3 and 2.4 in Chernozhukov et al.
(2014b) and all of the results in Chernozhukov et al. (2014a, 2016) require a
strictly positive lower bound on the weak variance of the Gaussian $P$-bridge
process, i.e.
$\inf_{f\in\mathcal{F}_{n}}\mathrm{Var}\big{(}G_{P}(f)\big{)}\geq\underline{\sigma}^{2}>0$.
This assumption automatically limits the applicability of these lemmas to
studentized statistics/ standardized function classes (Chernozhukov et al.,
2014a, b) and excludes relevant scenarios with variance decay (Lopes et al.,
2020; Lopes, 2022b; Lopes et al., 2023). In contrast, our Theorems 1–9 apply
to all function classes for which the strong variance of the Gaussian
$P$-bridge process, $\mathrm{Var}(\|G_{P}\|_{\mathcal{F}_{n}})$, does not
vanish “too fast”. A slowly vanishing strong variance is a weaker requirement
than a strictly positive lower bound on the weak variance (see Giessing (2023)
and also Lemmas 6 and 7).

The lower bound on the weak variance is an artifact of Chernozhukov et al.’s
(2014a; 2014b; 2016) proof technique which requires the joint distribution of
the finite-dimensional marginals of the Gaussian $P$-bridge process to be non-
degenerate.

##  4 Applications

###  4.1 Confidence ellipsoids for high-dimensional parameter vectors

Confidence regions are fundamental to uncertainty quantification in
multivariate statistics. Typically, their asymptotic validity relies on
multivariate CLTs. However, in high dimensions, when the number of parameters
exceeds the sample size, validity of confidence regions needs to be justified
differently. In this section, we show how the Gaussian process bootstrap
offers a practical solution to this problem. Existing work on bootstrap
confidence regions in high dimensions has focused exclusively on conservative,
rectangular confidence regions (e.g. Chernozhukov et al., 2013, 2014b, 2023).
For the first time, our results allow construction of tighter, elliptical
confidence regions and without sparsity assumptions.

Let $\theta_{0}\in\mathbb{R}^{d}$ be an unknown parameter and
$\hat{\theta}_{n}\equiv\hat{\theta}_{n}(X_{1},\ldots,X_{n})$ an estimator for
$\theta_{0}$ based on the simple random sample $X_{1},\ldots,X_{n}$. Then, an
asymptotic $1-\alpha$ confidence ellipsoid for $\theta_{0}$ can be constructed
as

$\displaystyle\mathcal{E}_{n}(c):=\left\\{\theta\in\mathbb{R}^{d}:\sqrt{n}\|\hat{\theta}_{n}-\theta\|_{2}\leq
c\right\\},$ (10)

where $\|\cdot\|_{2}$ denotes the Euclidean norm and $c>0$ solves

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}\left\\{\sqrt{n}\|\hat{\theta}_{n}-\theta_{0}\|_{2}\leq
c\right\\}=1-\alpha.$

In general, it is difficult to determine $c$ from above equation because the
sampling distribution of $\sqrt{n}\|\hat{\theta}_{n}-\theta_{0}\|_{2}$ often
does not have an explicit form. However, in many cases, $\hat{\theta}_{n}$
admits an asymptotically linear expansion of the form

$\displaystyle\hat{\theta}_{n}-\theta_{0}=\frac{1}{n}\sum_{i=1}^{n}\psi_{i}+R_{n},$

where $\psi_{1},\ldots,\psi_{n}$ are centered i.i.d. random vectors (influence
functions) and $R_{n}$ is a remainder term. Thus, in these cases, we have

$\displaystyle\sqrt{n}\|\hat{\theta}_{n}-\theta_{0}\|_{2}=\sup_{u\in
S^{d-1}}\left|\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi_{i}^{\prime}u\right|+\Theta_{n},\quad\quad\text{where}\quad\quad|\Theta_{n}|\leq\sqrt{n}\|R_{n}\|_{2}.$

Using this formulation, we can now apply Theorem 9 to approximate the quantile
$c$ and complete the construction of the asymptotic confidence ellipsoid. By
Algorithm 1 we need to construct a Gaussian process with index set $S^{d-1}$
and whose covariance function approximates the bi-linear map $(u,v)\mapsto
u^{\prime}\Omega_{\psi}v$, where
$\Omega_{\psi}=\mathrm{E}[\psi_{1}\psi_{1}^{\prime}]$. Clearly, the following
process does the job:

$\displaystyle\\{\widehat{Z}_{\psi}^{\prime}u:u\in
S^{d-1}\\},\quad\quad\text{where}\quad\quad\widehat{Z}_{\psi}\mid
X_{1},\ldots,X_{n}\sim N(0,\widehat{\Omega}_{\psi}),$

and
$\widehat{\Omega}_{\psi}=n^{-1}\sum_{i=1}^{n}(\psi_{i}-\overline{\psi}_{n})(\psi_{i}-\overline{\psi}_{n})^{\prime}$
and $\overline{\psi}_{n}=n^{-1}\sum_{i=1}^{n}\psi_{i}$. We denote the
$\alpha$-quantile of the supremum of this process by

$\displaystyle c_{n}(\alpha):=\inf\left\\{c\geq
0:\mathbb{P}\left\\{\|\widehat{Z}_{\psi}\|_{2}\leq c\mid
X_{1},\ldots,X_{n}\right\\}\geq\alpha\right\\}.$ (11)

To show that these quantiles uniformly approximate the quantiles of the
(asymptotic) distribution of $\sqrt{n}\|\hat{\theta}_{n}-\theta_{0}\|_{2}$
introduce the following assumptions:

######  Definition 2 (Sub-Gaussian random vector).

We call a centered random vector $X\in\mathbb{R}^{d}$ sub-Gaussian if
$\|X^{\prime}u\|_{\psi_{2}}^{2}\lesssim\mathrm{E}[(X^{\prime}u)^{2}]$ for all
$u\in\mathbb{R}^{d}$.

######  Assumption 1 (Sub-Gaussian influence functions).

The influence functions $\psi,\psi_{1},\ldots,\psi_{n}\in\mathbb{R}^{d}$ are
i.i.d. sub-Gaussian random vectors with covariance matrix $\Omega_{\psi}$ and
(i) $r(\Omega_{\psi})/r(\Omega_{\psi}^{2})=O(1)$, (ii)
$r(\Omega_{\psi})=o(n^{1/6})$, and (iii) $\sqrt{n}\|R_{n}\|_{2}=o_{p}(1)$,
where $r(A)=\mathrm{tr}(A)/\|A\|_{op}$ is the effective rank matrix $A$.

######  Assumption 2 (Heavy-tailed influence functions).

The influence functions $\psi,\psi_{1},\ldots,\psi_{n}\in\mathbb{R}^{d}$ are
centered i.i.d. random vectors with covariance matrix $\Omega_{\psi}$ such
that $M_{n}^{2}:=\mathrm{E}[\max_{1\leq i\leq
n}\|\psi_{i}\|_{2}^{2}]/\|\Omega_{\psi}\|_{op}$ $<\infty$ and
$m_{2,s}^{s}:=\mathrm{E}[\|\psi_{1}\|_{2}^{s}]/\|\Omega_{\psi}\|_{op}^{s/2}<\infty$
for some $s>3$. Furthermore, (i) $r(\Omega_{\psi})/r(\Omega_{\psi}^{2})=O(1)$,
(ii) $m_{2,s}/m_{2,3}=o(n^{1/3-1/s})$, (iii) $m_{2,3}=o(n^{1/6})$, (iv)
$M_{n}^{2}(\log n)=o(n)$, and (v) $\sqrt{n}\|R_{n}\|_{2}=o_{p}(1)$, where
$r(A)=\mathrm{tr}(A)/\|A\|_{op}$ is the effective rank matrix $A$.

Under either assumption, we have the following result:

######  Proposition 1 (Bootstrap confidence ellipsoid).

If either Assumption 1 or Assumption 2 holds, then

$\displaystyle\lim_{n\rightarrow\infty}\sup_{\alpha\in(0,1)}\left|\mathbb{P}\Big{\\{}\theta_{0}\in\mathcal{E}_{n}\big{(}c_{n}(1-\alpha)\big{)}\Big{\\}}-(1-\alpha)\right|=0,$

with the confidence ellipsoid $\mathcal{E}_{n}$ as defined in (10) and the
conditional quantile $c_{n}(1-\alpha)$ as in (11).

It is worth noticing that neither Assumption 1 nor 2 require sparsity of the
parameter $\theta_{0}$, the estimator $\hat{\theta}_{n}$, or the influence
functions $\psi_{1},\ldots,\psi_{n}$. Instead, Assumption 1 and 2 already hold
if the covariance matrix $\Omega_{\psi}$ has bounded effective rank (see
Giessing and Fan, 2023, Section 2.1 and Appendices A.1 and A.2). In this
context it is important to notice that, unlike Chernozhukov et al. (2014b), we
do not require a strictly positive lower bound on
$\inf_{\|u\|_{2}=1}\mathrm{E}[(Z_{\psi}^{\prime}u)^{2}]\equiv\lambda_{\min}(\Omega_{\psi})$,
the smallest eigenvalue of the covariance matrix $\Omega_{\psi}$. If such a
lower bound was needed, the effective rank $r(\Omega_{\psi})$ would grow
linearly in the dimension $d$ and Assumptions 1 and 2 would be violated if
$d\gg n$ (see Appendix A). Moreover, Proposition 1 cannot be deduced from
Theorem 2.1 in Chernozhukov et al. (2014b) because the upper bound in their
coupling inequality would feature the term
$\sqrt{d\log(1+1/\epsilon)}/n^{1/4}$ (and other terms as well), which is a
remnant of the entropy number of the $\epsilon$-net discretization of the
$d$-dimensional Euclidean unit-ball.

###  4.2 Inference on the spectral norm of high-dimensional covariance
matrices

Spectral statistics of random matrices play an important role in multivariate
statistical analysis. The asymptotic distributions of spectral statistics are
well established in low dimensions (e.g. Anderson, 1963; Waternaux, 1976;
Fujikoshi, 1980) and when the dimension is comparable to the sample size (e.g.
Johnstone, 2001; El Karoui, 2007; Péché, 2009; Bao et al., 2015). In the high-
dimensional case, when asymptotic arguments do not apply, bootstrap procedures
have proved to be effective in approximating distributions of certain maximum-
type spectral statistics (e.g. Han et al., 2018; Naumov et al., 2019; Lopes et
al., 2019, 2020; Lopes, 2022b). Here, we demonstrate that the Gaussian process
bootstrap is a viable alternative to these bootstrap procedures in
approximating the distribution of the spectral norm of a high-dimensional
sample covariance matrix.

Let $X_{1},\ldots,X_{n}\in\mathbb{R}^{d}$ be i.i.d. random vectors with law
$P$, mean zero, and covariance matrix $\Sigma\in\mathbb{R}^{d\times d}$.
Consider the spectral statistic

$\displaystyle
T_{n}:=\sqrt{n}\|\widehat{\Sigma}_{n}-\Sigma\|_{op},\quad\quad{}\widehat{\Sigma}_{n}=n^{-1}\sum_{i=1}^{n}X_{i}X_{i}^{\prime}.$

Since there does not exist a closed form expression of the sampling
distribution of $T_{n}$ when $d\gg n$, we apply Algorithm 1 to obtain an
approximation based on a Gaussian proxy process. We introduce the following
notation: Let $x,u,v\in\mathbb{R}^{d}$ and note that
$(x^{\prime}u)(x^{\prime}v)=\mathrm{vech}^{\prime}(xx^{\prime})H_{d}^{\prime}(v\otimes
u)$, where $\otimes$ denotes the Kronecker product, $\mathrm{vech}(\cdot)$ the
half-vectorization operator which turns a symmetric matrix
$A\in\mathbb{R}^{d\times d}$ into a $d(d+1)/2$ column vector (of unique
entries), and $H_{d}$ the duplication matrix such that
$H_{d}\>\mathrm{vech}(A)=\mathrm{vec}(A)$ where $\mathrm{vec}(\cdot)$ is the
ordinary vectorization operator. Whence,
$T_{n}\equiv\sqrt{n}\|Q_{n}-Q\|_{{\mathcal{F}_{n}}}$ where $Q_{n}$ is the
empirical measure of the collection of
$Y_{i}=\mathrm{vech}^{\prime}(X_{i}X_{i}^{\prime})H_{d}^{\prime}$, $1\leq
i\leq n$, $Q$ the pushforward of $P$ under the map $X\mapsto
Y=\mathrm{vech}^{\prime}(XX^{\prime})H_{d}^{\prime}$, and
$\mathcal{F}_{n}=\\{y\mapsto y(v\otimes u):u,v\in S^{d-1}\\}$. Since each
$f\in\mathcal{F}_{n}$ has a (not necessarily unique) representation in terms
of $u,v\in S^{d-1}$, in the following we identify $f\in\mathcal{F}_{n}$ with
pairs $(u,v)\in S^{d-1}\times S^{d-1}$. The covariance function associated
with the empirical process $\\{\sqrt{n}(Q_{n}-Q)(f):f\in\mathcal{F}_{n}\\}$ is
thus given by

$\displaystyle\big{(}(u_{1},v_{1}),(u_{2},v_{2})\big{)}$
$\displaystyle\mapsto(v_{1}\otimes u_{1})^{\prime}H_{d}\Omega
H_{d}^{\prime}(v_{2}\otimes u_{2}),$

where
$\Omega=\mathrm{E}\left[\mathrm{vech}(XX^{\prime}-\Sigma)\mathrm{vech}^{\prime}(XX^{\prime}-\Sigma)\right]\in\mathbb{R}^{d(d+1)/2\times
d(2+1)/2}$. Thus, in the light of Algorithm 1 the natural choice for the
Gaussian bootstrap process is

$\displaystyle\big{\\{}\widehat{Z}_{n}^{\prime}H_{d}^{\prime}(v\otimes
u):u,v\in S^{d-1}\big{\\}},\quad\quad\text{where}\quad\quad\widehat{Z}_{n}\mid
X_{1},\ldots,X_{n}\sim N(0,\widehat{\Omega}_{n}),$ (12)

and
$\widehat{\Omega}_{n}=n^{-1}\sum_{i=1}^{n}\mathrm{vech}(X_{i}X_{i}^{\prime}-\widehat{\Sigma}_{n})\mathrm{vech}^{\prime}(X_{i}X_{i}^{\prime}-\widehat{\Sigma}_{n})$
is the sample analogue of $\Omega$.

We make the following assumption:

######  Assumption 3 (Sub-Gaussian data).

The data $X,X_{1},\ldots,X_{n}\in\mathbb{R}^{d}$ are i.i.d. sub-Gaussian
random vectors with covariance matrix $\Sigma$ and
$\inf_{\|u\|_{2}=1}\mathrm{Var}((X^{\prime}u)^{2})\geq\kappa>0$.

######  Remark 1 (On the lower bound on the variances).

The strictly positive lower bound $\kappa>0$ on the variance of the quadratic
from $u^{\prime}X^{\prime}Xu$ is mild. The existence of the lower bound
$\kappa>0$ is equivalent to
$\mathrm{E}[(X^{\prime}u)^{4}]^{1/4}>\mathrm{E}[(X^{\prime}u)^{2}]^{1/2}$ for
all $u\in S^{d-1}$. The latter inequality holds if the law of $X$ does not
concentrate on a lower dimensional subspace of $\mathbb{R}^{d}$. Since the
bounds in Proposition 2 are explicit in $\kappa>0$, Proposition 2 also applies
to scenarios in which $\kappa\equiv\kappa(n,d)\rightarrow 0$ as
$n,d\rightarrow\infty$. Similar lower bounds on the variance of
$u^{\prime}X^{\prime}Xu$ appear in Lopes (2022b); Lopes et al. (2023).

We have the following result:

######  Proposition 2 (Bootstrap approximation of the distribution of spectral
norms of covariance matrices).

Suppose that Assumption 3 holds. Let
$\widehat{S}_{n}=\mathrm{vec}^{-1}(H_{d}\widehat{Z}_{n})\in\mathbb{R}^{d\times
d}$ with $\widehat{Z}_{n}$ as given in (12). Then,

$\displaystyle\sup_{s\geq
0}\Big{|}\mathbb{P}\left\\{\sqrt{n}\|\widehat{\Sigma}_{n}-\Sigma\|_{op}\leq
s\right\\}-\mathbb{P}\left\\{\|\widehat{S}_{n}\|_{op}\leq s\mid
X_{1},\ldots,X_{n}\right\\}\Big{|}$ $\displaystyle\quad\quad\quad\lesssim
O_{p}\left(\left(\frac{\|\Sigma\|_{op}^{2/3}}{\kappa^{2/3}}\vee\frac{\|\Sigma\|_{op}^{4/3}}{\kappa^{4/3}}\right)\left(\sqrt{\frac{(\log
en)^{2}\mathrm{r}^{7}(\Sigma)}{n}}\vee\frac{(\log
en)^{2}\mathrm{r}^{4}(\Sigma)}{n}\right)^{1/3}\right)$
$\displaystyle\quad{}\quad{}\quad{}\quad{}\quad{}\quad{}+\left(\frac{\|\Sigma\|_{op}}{\kappa}\vee\frac{\|\Sigma\|_{op}^{2}}{\kappa^{2}}\right)\frac{r^{2}(\Sigma)}{n^{1/6}},$

where $r(\Sigma)=\mathrm{tr}(\Sigma)/\|\Sigma\|_{op}$ is the effective rank of
$\Sigma$ and $\lesssim$ hides an absolute constant independent of
$n,d,\Sigma$, and $\kappa$.

######  Remark 2.

Note that the matrix $\widehat{S}_{n}$ is symmetric just as the target matrix
$\widehat{\Sigma}_{n}-\Sigma$.

We conclude this section with a comparison of Proposition 2 to existing
results in the literature: First, if
$(\|\Sigma\|_{op}/\kappa\vee\|\Sigma\|_{op}^{2}/\kappa^{2})r^{2}(\Sigma)=o(n^{1/6})$,
then the upper bound in above theorem is asymptotically negligible. In this
case, the bootstrapped distribution of the Gaussian proxy statistic
$\|\widehat{S}_{n}\|_{op}$ consistently approximates the distribution of
$\|\widehat{\Sigma}_{n}-\Sigma\|_{op}$. Since this rate depends only on the
effective rank, it is dimension-free and cannot be derived through the results
in Chernozhukov et al. (2014b). Second, unlike the results in Lopes (2022b)
and Lopes et al. (2023), Proposition 2 does not rely on specific assumptions
about the decay of the eigenvalues of $\Sigma$. And yet, under certain
circumstances, the consistency rate provided by Proposition 2 can be faster
than theirs. Specifically, the bootstrap procedure described in Lopes (2022b)
achieves consistency at a rate of $n^{-\frac{\beta-1/2}{2\beta+4+\epsilon}}$
for $\epsilon>0$ and $\beta>1/2$. In our context, the parameter $\beta$
determines the rate at which the eigenvalues of $\Sigma$ decrease, i.e.
$\lambda_{k}(\Sigma)\asymp k^{-2\beta}$, $k=1,\ldots,d$. To achieve a rate
faster than $n^{-1/6}$, $\beta$ must be greater than $(7+\epsilon)/4$ which
requires an extremely fast decay of the eigenvalues. Third, Lopes et al.
(2023) conduct extensive numerical experiments and observe that the bootstrap
approximation exhibits a sharp phase transition from accurate to inaccurate
when $\beta$ switches from greater than $1/2$ to less than $1/2$. This
observation aligns not only with their own theoretical findings but also with
the upper bound presented in Proposition 2, since, under their modeling
assumptions, the effective rank $r(\Sigma)$ remains bounded if $\beta>1/2$ but
diverges if $\beta\leq 1/2$.

###  4.3 Simultaneous confidence bands for functions in reproducing kernel
Hilbert spaces

Reproducing kernel Hilbert spaces (RKHS) are an integral part of statistics,
with applications in classical non-parametric statistics (Wahba, 1990),
machine learning (Schölkopf and Smola, 2002; Steinwart and Christmann, 2008)
and, most recently, (deep) neural nets (Belkin et al., 2018; Jacot et al.,
2018; Bohn et al., 2019; Unser, 2019; Chen and Xu, 2020). In this section, we
consider constructing simultaneous confidence bands for functions in RKHS by
bootstrapping the distribution of a bias-corrected kernel ridge regression
estimator. Recently, this problem has been addressed by Singh and Vijaykumar
(2023) with a symmetrizied multiplier bootstrap. Here, we propose an
alternative based on the Gaussian process bootstrap using a truncated
Karhunen-Loève decomposition. We point out several commonalities and
differences between the two procedures.

In the following, $\mathcal{H}$ denotes the RKHS of continuous functions from
$S$ to $\mathbb{R}$ associated to the symmetric and positive definite kernel
$k:S\times S\rightarrow\mathbb{R}$. We allow $\mathcal{H}$, $S$, and $k$ to
change with the sample size $n$, but we do not make this dependence explicit.
We write $\|\cdot\|_{\mathcal{H}}$ for the norm induced by the inner product
$\langle\cdot,\cdot\rangle_{\mathcal{H}}$ and $\|\cdot\|_{\infty}$ for the
supremum norm. For $x\in S$, we let $k_{x}:S\rightarrow\mathbb{R}$ be the
function $y\mapsto k(x,y)$. Then, $k_{x}\in\mathcal{H}$ for all $x\in S$ and,
by the reproducing property, $f(x)=\langle f,k_{x}\rangle_{\mathcal{H}}$ for
all $f\in\mathcal{H}$ and $x\in S$. The kernel induces the so-called kernel
metric $d_{k}(x,y):=\|k_{x}-k_{y}\|_{\mathcal{H}}$, $x,y\in S$. Given
$f\in\mathcal{H}$ ($z\in S$) we denote its dual by $f^{*}\in\mathcal{H}^{*}$
($z^{*}\in S^{*}\subset\mathcal{H}$). For $f,g\in\mathcal{H}$ we define the
tensor product $f\otimes g^{*}:\mathcal{H}\rightarrow\mathcal{H}$ by
$h\mapsto(f\otimes g^{*})(h):=\langle g,h\rangle_{\mathcal{H}}f$. For
operators on $\mathcal{H}$ we use $\|\cdot\|_{op}$, $\|\cdot\|_{HS}$, and
$\mathrm{tr}(\cdot)$ to denote operator norm, Hilbert-Schmidt norm, and trace,
respectively. For further details on RKHSs we refer to Berlinet and Thomas-
Agnan (2004).

Let $Y\in\mathbb{R}$ and $X\in S$ have joint law $P$. Given a simple random
sample $(Y_{1},X_{1}),\ldots,(Y_{n},$ $X_{n})$ our goal is to construct
uniform confidence bands for the conditional mean function
$x\mapsto\mathrm{E}[Y\mid X=x]$ or, rather, its best approximation in
hypothesis space $\mathcal{H}$, i.e.

$\displaystyle
f_{0}\in\operatorname*{\arg\min}_{f\in\mathcal{H}}\mathrm{E}\left[\big{(}Y-f(X)\big{)}^{2}\right].$

To this end, consider the classical kernel ridge regression estimator

$\displaystyle\widehat{f}_{n}\in\operatorname*{\arg\min}_{f\in\mathcal{H}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\big{(}Y_{i}-f(X_{i})\big{)}^{2}+\lambda\|f\|_{\mathcal{H}}^{2}\right\\},\quad\lambda>0,$

and define its bias-corrected version as

$\displaystyle\widehat{f}^{\mathrm{bc}}_{n}:=\widehat{f}_{n}+\lambda(\widehat{T}_{n}+\lambda)^{-1}\widehat{f}_{n},\quad\quad\widehat{T}_{n}=\frac{1}{n}\sum_{i=1}^{n}\left(k_{X_{i}}\otimes
k_{X_{i}}^{*}\right).$ (13)

We propose to construct simultaneous $1-\alpha$ confidence bands for $f_{0}$
based on $\widehat{f}^{\mathrm{bc}}_{n}$ via the rectangle

$\displaystyle\mathcal{R}_{n}(c):=\left\\{f\in\mathcal{F}:\sqrt{n}\|\widehat{f}^{\mathrm{bc}}_{n}-f\|_{\infty}\leq
c\right\\},$ (14)

where $c>0$ approximates the (asymptotic) $1-\alpha$ quantile of the law of
$\sqrt{n}\|\widehat{f}^{\mathrm{bc}}_{n}-f_{0}\|_{\infty}$. To compute $c>0$
we proceed in two steps: First, we show that
$\sqrt{n}\|\widehat{f}^{\mathrm{bc}}_{n}-f_{0}\|_{\infty}$ can be written as
the sum of the supremum of an empirical process and a negligible remainder
term. Then, we apply the strategy developed in Section 3.5 to bootstrap the
supremum of the empirical process.

By Lemma 20 in Appendix B.5,

$\displaystyle\widehat{f}^{\mathrm{bc}}_{n}-f_{0}=(T+\lambda)^{-2}T\left(\frac{1}{n}\sum_{i=1}^{n}\big{(}Y_{i}-f_{0}(X_{i})\big{)}k_{X_{i}}\right)+R_{n},$

where $T=\mathrm{E}[k_{X}\otimes k_{X}^{*}]$ and $R_{n}$ is a higher-order
remainder term. Since $k$ is a reproducing kernel, we have $k_{x}(z)=\langle
k_{x},k_{z}\rangle_{\mathcal{H}}=\langle k_{x},z^{*}\rangle_{\mathcal{H}}$ for
all $x,z\in S$. Hence, above expansion implies

$\displaystyle\sqrt{n}\|\widehat{f}^{\mathrm{bc}}_{n}-f_{0}\|_{\infty}=\sup_{u\in
S^{*}}\left|\left\langle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}V_{i},u\right\rangle_{\mathcal{H}}\right|+\Theta_{n},\quad\quad\text{where}\quad\quad|\Theta_{n}|\leq\sqrt{n}\|R_{n}\|_{\infty},$

and $V_{i}=(T+\lambda)^{-2}T\big{(}Y_{i}-f_{0}(X_{i})\big{)}k_{X_{i}}$, $1\leq
i\leq n$. By Lemma 21 in Appendix B.5, $\sqrt{n}\|R_{n}\|_{\infty}$ is
negligible with high probability. Moreover, since $f_{0}$ is the best
approximation in square loss, the random elements $V_{i}$’s have mean zero.
Thus, $\sup_{u\in S^{*}}|\langle
n^{-1/2}\sum_{i=1}^{n}V_{i},u\rangle_{\mathcal{H}}|\equiv\sqrt{n}\|Q_{n}-Q\|_{\mathcal{F}_{n}}$,
where $Q_{n}$ is the empirical measure of the $V_{i}$’s, $Q$ the pushforward
of $P$ under the map $(Y,X)\mapsto
V=(T+\lambda)^{-2}T\big{(}Y-f_{0}(X)\big{)}k_{X}$, and
$\mathcal{F}_{n}=\\{v\mapsto\langle v,u\rangle_{\mathcal{H}}:u\in S^{*}\\}$.
Since the functions $f\in\mathcal{F}_{n}$ are just the evaluation functionals
of $u\in S^{*}$, in the following we identify $f\in\mathcal{F}_{n}$ with its
corresponding $u\in S^{*}$. The covariance function $\mathcal{C}:S^{*}\times
S^{*}\rightarrow\mathbb{R}$ associated with the empirical process
$\\{\sqrt{n}(Q_{n}-Q)(f):f\in\mathcal{F}_{n}\\}$ is thus given by

$\displaystyle(u_{1},u_{2})\mapsto\mathrm{E}[\langle
V,u_{1}\rangle_{\mathcal{H}}\langle
V,u_{2}\rangle_{\mathcal{H}}]=\mathrm{E}\big{[}\big{\langle}(V\otimes
V^{*})u_{1},u_{2}\big{\rangle}_{\mathcal{H}}\big{]}=\big{\langle}\Omega
u_{1},u_{2}\big{\rangle}_{\mathcal{H}},$

with covariance operator

$\displaystyle\Omega:=\sigma_{0}^{2}T(T+\lambda)^{-2}T(T+\lambda)^{-2}T,\quad\quad\text{where}\quad\quad\sigma_{0}^{2}:=\mathrm{E}\big{[}\big{(}Y-f_{0}(X)\big{)}^{2}\big{]},$

where we have used that the operators $T$ and $(T+\lambda)^{-1}$ commute.

We proceed to construct a Gaussian proxy process as outlined in Section 3.5:
Let
$\widehat{\Omega}_{n}=\widehat{\sigma}_{n}^{2}\widehat{T}_{n}(\widehat{T}_{n}+\lambda)^{-2}\widehat{T}_{n}(\widehat{T}_{n}+\lambda)^{-2}\widehat{T}_{n}$
and
$\widehat{\sigma}_{n}^{2}=n^{-1}\sum_{i=1}^{n}\big{(}Y_{i}-\widehat{f}_{n}(X_{i})\big{)}^{2}$
be the plug-in estimates of the covariance operator $\Omega$ and the variance
$\sigma_{0}^{2}$, respectively. Define $\widehat{\mathcal{C}}_{n}:S^{*}\times
S^{*}\rightarrow\mathbb{R}$ by
$(u_{1},u_{2})\mapsto\big{\langle}\widehat{\Omega}_{n}u_{1},u_{2}\big{\rangle}_{\mathcal{H}}$.
Recall definition (8) of the integral operator $T_{K}$. In the present setup,
$T_{\widehat{\mathcal{C}}_{n}}=\widehat{\Omega}_{n}$ by Fubini’s theorem.
Denote by
$\big{\\{}(\widehat{\lambda}_{k},\widehat{\varphi}_{k})\big{\\}}_{k=1}^{\infty}$
the eigenvalue and eigenfunction pairs of $\widehat{\Omega}_{n}$. Further, let
$\\{\xi_{k}\\}_{k=1}^{\infty}$ be a sequence of i.i.d. standard normal random
variables. Then, for $m\geq 1$ and $u,u_{1},u_{2}\in S^{*}$ define

$\displaystyle\widehat{Z}_{n}^{m}(u):=\sum_{k=1}^{m}\xi_{k}\sqrt{\widehat{\lambda}_{k}}\widehat{\varphi}_{k}(u),\quad\quad\widehat{\mathcal{C}}_{n}^{m}(u_{1},u_{2}):=\sum_{k=1}^{m}\widehat{\lambda}_{k}\widehat{\varphi}_{k}(u_{1})\widehat{\varphi}_{k}(u_{2})=\big{\langle}\widehat{\Omega}_{n}^{m}u_{1},u_{2}\big{\rangle}_{\mathcal{H}}.$
(15)

where $\widehat{\Omega}_{n}^{m}$ is the best rank-$m$ approximation of
$\widehat{\Omega}_{n}$.

Given Proposition 7 we postulate that the process $\widehat{Z}_{n}^{\infty}$
is an almost sure version of a Gaussian process on $S^{*}$ with covariance
function $\widehat{\mathcal{C}}_{n}$ (or, equivalently, with covariance
operator $\widehat{\Omega}_{n}$. Consequently, Theorem 9 guarantees validity
of the Gaussian process bootstrap based on $\widehat{Z}_{n}^{m}$. To make all
these claims rigorous, consider the following assumptions:

######  Assumption 4 (On the kernel).

The kernel $k:S\times S\rightarrow\mathbb{R}$ is symmetric, positive semi-
definite, continuous, and bounded, i.e. $\sup_{x\in
S}\sqrt{|k(x,x)|}=:\kappa<\infty$.

######  Remark 3.

The assumptions on the kernel are standard and important (Berlinet and Thomas-
Agnan, 2004). The continuity of $k$ guarantees that $k_{X}$ is a random
element on $\mathcal{H}$ whenever $X\in S$ is a random variable. It also
implies that the RKHS $\mathcal{H}$ is separable whenever $S$ is separable.
The boundedness and the reproducing property of $k$ imply that
$\|f\|_{\infty}\leq\kappa\|f\|_{\mathcal{H}}$ for all $f\in\mathcal{H}$.

######  Assumption 5 (On the data).

The data $(Y_{1},X_{1}),\ldots,(Y_{n},X_{n})\in\mathbb{R}\times S$ are i.i.d.
random elements defined on an abstract product probability space
$(\Omega,\mathcal{A},\mathbb{P})$. The $Y_{i}$’s are almost surely bounded,
i.e. there exists an absolute constant $B>0$ such that $\max_{1\leq i\leq
n}|Y_{i}|\leq B$ almost surely.

######  Remark 4.

The almost sure boundedness of the $Y_{i}$’s is a strong assumption. We
introduce this assumption to keep technical arguments at a minimum. Singh and
Vijaykumar (2023) impose an equivalent boundedness condition on the pseudo-
residuals $\varepsilon_{i}=Y_{i}-f_{0}(X_{i})$, $1\leq i\leq n$.

######  Assumption 6 (On the population and sample covariance operators).

For all $n\geq 1$ (i) there exists $\omega_{S}>0$ such that $\inf_{u\in
S^{*}}\langle\Omega u,u\rangle_{\mathcal{H}}\geq\omega_{S}$ and (ii)
$\mathrm{tr}(\Omega)\vee\mathrm{tr}(\widehat{\Omega}_{n})<\infty$.

######  Remark 5.

Condition (i) is the Hilbert space equivalent to the lower bound on the
variance in Assumption 3. While Singh and Vijaykumar (2023) do not explicitly
impose a lower bound on the covariance operator, such a lower bound is implied
by their Assumption 5.2 (see Giessing, 2023, for details). In eq.
(LABEL:eq:theorem:Uniform-CI-Bands-RKHS-8) in Appendix 4 we provide a general
non-asymptotic complement to below Proposition 3 which applies even if
$\omega_{S}\equiv\omega_{S}(n)\rightarrow 0$ as $n\rightarrow\infty$.
Condition (ii) is a classical assumption in learning theory on RKHS
(Mendelson, 2002). Together with Conditions (i) it implies that $\Omega$ and
$\widehat{\Omega}_{n}$ are finite rank and trace class operators,
respectively.

Denote the $\alpha$-quantile of the supremum of the Gaussian proxy process
$\widehat{Z}_{n}^{m}$ in (15) by

$\displaystyle c_{n}^{m}(\alpha)$ $\displaystyle:=\inf\left\\{s\geq
0:\mathbb{P}\left\\{\sup_{u\in
S^{*}}\left|\big{\langle}\widehat{Z}_{n}^{m},u\big{\rangle}_{\mathcal{H}}\right|\leq
s\mid(Y_{1},X_{1}),\ldots,(Y_{n},X_{n})\right\\}\geq\alpha\right\\}.$ (16)

An application of Theorem 9 yields:

######  Proposition 3 (Bootstrap quantiles).

Let $(S,d_{k})$ be a compact metric space such that
$\int_{0}^{\infty}\sqrt{N(S,d_{k},\varepsilon)}\\\ d\varepsilon<\infty$. If
Assumptions 4, 5, and 6, and the rates in eq. (77) in Appendix E hold, then

$\displaystyle\sup_{\alpha\in(0,1)}\left|\mathbb{P}\left\\{\sqrt{n}\|\widehat{f}^{\mathrm{bc}}_{n}-f_{0}\|_{\infty}\leq
c_{n}^{m}(\alpha)\right\\}-\alpha\right|$
$\displaystyle\quad=o(1)+\mathbb{P}\left\\{\big{\|}\widehat{\Omega}_{n}-\widehat{\Omega}_{n}^{m}\big{\|}_{op}>\big{(}\kappa\mathfrak{n}_{1}(\lambda)+\bar{\sigma}^{2}\big{)}\|T^{3}(T+\lambda)^{-4}\|_{op}\sqrt{\frac{\log
n}{n\lambda^{2}}}\right\\},$

where
$\bar{\sigma}^{2}\geq\sigma_{0}^{2}\vee\kappa^{2}(B+\kappa\|f_{0}\|_{\mathcal{H}})^{2}\vee
1$ and
$\mathfrak{n}_{1}^{2}(\lambda)=\mathrm{tr}\left((T+\lambda)^{-2}T\right)$.

The finite metric entropy condition on the set $S$ ensures that the Gaussian
bootstrap process (15) is almost surely bounded and uniformly continuous on
$S^{*}\subset\mathcal{H}$ (as required by Theorem 9). This condition is not
merely technical but also intuitive: Since the RKHS $\mathcal{H}$ is the
completion of $\mathrm{span}(\\{k_{x}:x\in S\\})$ and $k$ is bounded and
continuous, conditions that guarantee the continuity of Gaussian bootstrap
processes on (a subset of) $\mathcal{H}$ should indeed be attributable to
properties of $S$. Importantly, the metric entropy condition on $S$ does not
impose restrictions on the dimension of $\mathcal{H}$. Only Assumption 6 (ii)
implicitly imposes restrictions on the dimension of $\mathcal{H}$.

Since under the conditions of Proposition 3,
$\lim_{m\rightarrow\infty}\|\widehat{\Omega}_{n}-\widehat{\Omega}_{n}^{m}\|_{op}=0$
almost surely for all $n\geq 1$, it follows that the bootstrap confidence band
proposed in (14) is asymptotically valid:

######  Corollary 4 (Simultaneous bootstrap confidence bands).

Under the setup of Proposition 3,

$\displaystyle\lim_{n\rightarrow\infty}\lim_{m\rightarrow\infty}\sup_{\alpha\in(0,1)}\left|\mathbb{P}\Big{\\{}f_{0}\in\mathcal{R}_{n}\big{(}c_{n}^{m}(1-\alpha)\big{)}\Big{\\}}-(1-\alpha)\right|=0,$

with the uniform confidence band $\mathcal{R}_{n}$ as defined in (14) and
quantile $c_{n}(1-\alpha)$ as in (16).

A thorough comparison of the Gaussian process bootstrap and Singh and
Vijaykumar’s (2023) symmetrized multiplier bootstrap is beyond the scope of
this paper. In practice, both methods yield biased confidence bands for
$f_{0}$, albeit for different reasons: Singh and Vijaykumar’s (2023) bias
stems from constructing a confidence band for the pseudo-true regression
function $f_{\lambda}=(T+\lambda)^{-1}Tf_{0}$ without correcting the
regularization bias induced by $\lambda>0$, ours is due to using an
$m$-truncated Karhunen-Loève decomposition based on a finite number of
eigenfunctions. In future work we will explore ways to mitigate these biases
by judiciously choosing $\lambda\equiv\lambda(n)\rightarrow 0$ and $m\equiv
m(d)\rightarrow\infty$.

##  5 Conclusion

In this paper we have developed a new approach to theory and practice of
Gaussian and bootstrap approximations of the sampling distribution of suprema
of empirical processes. We have put special emphasize on non-asymptotic
approximations that are entropy- and weak variance-free, and have allowed the
function class $\mathcal{F}_{n}$ to vary with the sample size $n$ and to be
non-Donsker. We have shown that such general approximation results are useful,
among other things, for inference on high-dimensional statistical models and
reproducing kernel Hilbert spaces. However, theory and methodology in this
paper have three limitations that need to be addressed in future work:

  * •

_Reliance on independent and identically distributed data._ All statistically
relevant results in this paper depend on Proposition 1, which heavily relies
on the assumption of independent and identically distributed data. Expanding
Proposition 1 to accommodate non-identical distributed data would be a first
step towards solving simultaneous and large-scale two-sample testing problems
and conducting inference in high-dimensional fixed design settings. Currently,
the results in this paper are exclusively applicable to one-sample testing and
unconditional inference.

  * •

_Lack of tight lower bounds on the strong variances of most Gaussian
processes._ One the most notable features of the results in this paper is the
fact that all upper bounds depend on the inverse of the strong variance of
some Gaussian proxy process. Unfortunately, in statistical applications this
poses a formidable challenge since up until now there exist only few
techniques to derive tight lower bounds on these strong variances (Giessing
and Fan, 2023; Giessing, 2023). We either need new tools or we need to develop
Gaussian and bootstrap approximations for statistics other than maxima/
suprema. The latter will require new anti-concentration inequalities.

  * •

_Biased quantile estimates due to bootstrapping non-pivotal statistic._ The
Gaussian process bootstrap is based on a non-pivotal statistic, i.e. the
sampling distribution of the supremum depends on the unknown population
covariance function. In practice, when bootstrapping non-pivotal statistics
the estimated quantiles often differ substantially from the true quantiles.
Several bias correction schemes have been proposed in the classical setting
(Davison et al., 1986; Beran, 1987; Hall and Martin, 1988; Shi, 1992). Since
the Gaussian process bootstrap is not a re-sampling procedure in the classical
sense, these techniques do not apply. In Giessing and Fan (2023) we therefore
develop the spherical bootstrap to improve accuracy and efficiency when
bootstrapping $\ell_{p}$-statistics. However, this approach does not
generalize to arbitrary empirical processes as the one in Section 4.3. Thus,
there is an urgent need for new bias correction schemes.

## Acknowledgement

Alexander Giessing is supported by NSF grant DMS-2310578.

## References

  * Anderson (2003) T. Anderson.  _An Introduction to Multivariate Statistical Analysis_.  Wiley Series in Probability and Statistics. Wiley, 2003.  ISBN 9780471360919.
  * Anderson (1963) T. W. Anderson.  Asymptotic Theory for Principal Component Analysis.  _The Annals of Mathematical Statistics_ , 34(1):122 – 148, 1963.
  * Bao et al. (2015) Z. Bao, G. Pan, and W. Zhou.  Universality for the largest eigenvalue of sample covariance matrices with general population.  _The Annals of Statistics_ , 43(1):382–421, 2015\.
  * Belkin et al. (2018) M. Belkin, S. Ma, and S. Mandal.  To understand deep learning we need to understand kernel learning.  In J. Dy and A. Krause, editors, _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pages 541–549. PMLR, 2018.
  * Bentkus (2003) V. Bentkus.  On the dependence of the Berry–Esseen bound on dimension.  _Journal of Statistical Planning and Inference_ , 113(2):385 – 402, 2003.  ISSN 0378-3758.
  * Beran (1987) R. Beran.  Prepivoting to reduce level error of confidence sets.  _Biometrika_ , 74(3):457–468, 1987.
  * Berlinet and Thomas-Agnan (2004) A. Berlinet and C. Thomas-Agnan.  _Reproducing Kernel Hilbert Spaces in Probability and Statistics_.  Springer, 2004.
  * Bhattacharya and Holmes (2010) R. Bhattacharya and S. Holmes.  An exposition of Götze’s estimation of the rate of convergence in the multivariate central limit theorem, 2010.
  * Bohn et al. (2019) B. Bohn, C. Rieger, and M. Griebel.  A representer theorem for deep kernel learning.  _The Journal of Machine Learning Research_ , 20(1):2302–2333, 2019.
  * Bolthausen (1984) E. Bolthausen.  An estimate of the remainder in a combinatorial central limit theorem.  _Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete_ , 66(3):379–386, 1984.
  * Bong et al. (2023) H. Bong, A. K. Kuchibhotla, and A. Rinaldo.  Dual induction clt for high-dimensional m-dependent data, 2023.
  * Cattaneo et al. (2022) M. D. Cattaneo, R. P. Masini, and W. G. Underwood.  Yurinskii’s coupling for martingales, 2022.
  * Chen and Xu (2020) L. Chen and S. Xu.  Deep neural tangent kernel and laplace kernel have the same rkhs.  In _International Conference on Learning Representations_ , 2020.
  * Chen (2018) X. Chen.  Gaussian and bootstrap approximations for high-dimensional u-statistics and their applications.  _Ann. Statist._ , 46(2):642–678, 04 2018.
  * Chen et al. (2015) Y.-C. Chen, C. R. Genovese, and L. Wasserman.  Asymptotic theory for density ridges.  _The Annals of Statistics_ , 43(5):1896–1928, 2015.
  * Chen et al. (2016) Y.-C. Chen, C. R. Genovese, R. J. Tibshirani, and L. Wasserman.  Nonparametric modal regression.  _The Annals of Statistics_ , 44(2):489 – 514, 2016.
  * Chernozhukov et al. (2013) V. Chernozhukov, D. Chetverikov, and K. Kato.  Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors.  _The Annals of Statistics_ , 41(6):2786–2819, 12 2013.
  * Chernozhukov et al. (2014a) V. Chernozhukov, D. Chetverikov, and K. Kato.  Anti-concentration and honest, adaptive confidence bands.  _The Annals of Statistics_ , 42(5):1787–1818, 2014a.
  * Chernozhukov et al. (2014b) V. Chernozhukov, D. Chetverikov, and K. Kato.  Gaussian approximation of suprema of empirical processes.  _The Annals of Statistics_ , 42(4):1564–1597, 08 2014b.
  * Chernozhukov et al. (2015) V. Chernozhukov, D. Chetverikov, and K. Kato.  Comparison and anti-concentration bounds for maxima of gaussian random vectors.  _Probability Theory and Related Fields_ , 162(1):47–70, Jun 2015.
  * Chernozhukov et al. (2016) V. Chernozhukov, D. Chetverikov, and K. Kato.  Empirical and multiplier bootstraps for suprema of empirical processes of increasing complexity, and related gaussian couplings.  _Stochastic Processes and their Applications_ , 126(12):3632–3651, 2016.  ISSN 0304-4149.  In Memoriam: Evarist Giné.
  * Chernozhukov et al. (2017) V. Chernozhukov, D. Chetverikov, and K. Kato.  Central limit theorems and bootstrap in high dimensions.  _The Annals of Probability_ , 45(4):2309–2352, 07 2017.
  * Chernozhukov et al. (2019) V. Chernozhukov, D. Chetverikov, K. Kato, and Y. Koike.  Improved central limit theorem and bootstrap approximations in high dimensions.  _arXiv preprint, arXiv:1912.10529_ , 12 2019.
  * Chernozhukov et al. (2020) V. Chernozhukov, D. Chetverikov, and Y. Koike.  Nearly optimal central limit theorem and bootstrap approximations in high dimensions.  2020\.
  * Chernozhukov et al. (2023) V. Chernozhukov, D. Chetverikov, K. Kato, and Y. Koike.  High-dimensional data bootstrap.  _Annual Review of Statistics and Its Application_ , 10(1):427–449, 2023.
  * Chetverikov (2019) D. Chetverikov.  Testing regression monotonicity in econometric models.  _Econometric Theory_ , 35(4):729–776, 2019.
  * Davison et al. (1986) A. C. Davison, D. V. Hinkley, and E. Schechtman.  Efficient bootstrap simulation.  _Biometrika_ , 73(3):555–566, 1986.
  * Deng and Zhang (2020) H. Deng and C.-H. Zhang.  Beyond gaussian approximation: Bootstrap for maxima of sums of independent random vectors.  _arXiv preprint, arXiv:1705.09528_ , 2020.
  * Dezeure et al. (2017) R. Dezeure, P. Bühlmann, and C.-H. Zhang.  High-dimensional simultaneous inference with the bootstrap.  _Test_ , 26(4):685–719, 2017.
  * Dudley (2002) R. Dudley.  _Real Analysis and Probability_.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2002.
  * Dudley (2014) R. Dudley.  _Uniform Central Limit Theorems_.  Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2014.
  * El Karoui (2007) N. El Karoui.  Tracy–Widom limit for the largest eigenvalue of a large class of complex sample covariance matrices.  _The Annals of Probability_ , 35(2):663 – 714, 2007.
  * Fan et al. (2018) J. Fan, Q.-M. Shao, and W.-X. Zhou.  Are discoveries spurious? Distributions of maximum spurious correlations and their applications.  _Ann. Statist._ , 46(3):989–1017, 06 2018.
  * Fang and Koike (2021) X. Fang and Y. Koike.  High-dimensional central limit theorems by Stein’s method.  _The Annals of Applied Probability_ , 31(4):1660 – 1686, 2021.
  * Folland (1999) G. Folland.  _Real Analysis: Modern Techniques and Their Applications_.  Pure and Applied Mathematics: A Wiley Series of Texts, Monographs and Tracts. Wiley, 1999.
  * Fujikoshi (1980) Y. Fujikoshi.  Asymptotic expansions for the distributions of the sample roots under nonnormality.  _Biometrika_ , 67(1):45–51, 1980.
  * Ghanem and Spanos (2003) R. Ghanem and P. Spanos.  _Stochastic Finite Elements: A Spectral Approach_.  Dover Publications, 2003.
  * Giessing (2023) A. Giessing.  Anti-Concentration of Suprema of Gaussian Processes and Gaussian Order Statistics.  _Working Paper_ , 2023.
  * Giessing and Fan (2023) A. Giessing and J. Fan.  A bootstrap hypothesis test for high-dimensional mean vectors.  _Working Paper_ , 2023.
  * Giessing and Wang (2021) A. Giessing and J. Wang.  Debiased inference on heterogeneous quantile treatment effects with regression rank-scores, 2021.
  * Giné and Nickl (2016) E. Giné and R. Nickl.  _Mathematical Foundations of Infinite-Dimensional Statistical Models_.  Cambridge University Press, 2016.
  * Götze (1991) F. Götze.  On the rate of convergence in the multivariate clt.  _The Annals of Probability_ , 19(2):724–739, 1991\.
  * Hall and Martin (1988) P. Hall and M. A. Martin.  On bootstrap resampling and iteration.  _Biometrika_ , 75(4):661–671, 1988.
  * Han et al. (2018) F. Han, S. Xu, and W.-X. Zhou.  On Gaussian comparison inequality and its application to spectral analysis of large random matrices.  _Bernoulli_ , 24(3):1787 – 1833, 2018.
  * Hardy (2006) M. Hardy.  Combinatorics of partial derivatives.  _Electronic Journal of Combinatorics_ , 2006.
  * Hsing and Eubank (2015) T. Hsing and R. Eubank.  _Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators_.  Wiley, 2015.
  * Jacot et al. (2018) A. Jacot, F. Gabriel, and C. Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  NIPS’18, page 8580–8589, 2018.
  * Jain and Kallianpur (1970) N. C. Jain and G. Kallianpur.  A Note on Uniform Convergence of Stochastic Processes.  _The Annals of Mathematical Statistics_ , 41(4):1360 – 1362, 1970.
  * Janková et al. (2020) J. Janková, R. D. Shah, P. Bühlmann, and R. J. Samworth.  Goodness-of-fit Testing in High Dimensional Generalized Linear Models.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 82(3):773–795, 05 2020.
  * Johnstone (2001) I. M. Johnstone.  On the distribution of the largest eigenvalue in principal components analysis.  _The Annals of Statistics_ , 29(2):295 – 327, 2001.
  * Koltchinskii and Lounici (2017) V. Koltchinskii and K. Lounici.  Concentration inequalities and moment bounds for sample covariance operators.  _Bernoulli_ , 23(1):110 – 133, 2017.
  * Kuchibhotla et al. (2021) A. K. Kuchibhotla, S. Mukherjee, and D. Banerjee.  High-dimensional CLT: Improvements, non-uniform extensions and large deviations.  _Bernoulli_ , 27(1):192 – 217, 2021.
  * Le Cam (1986) L. Le Cam.  _Asymptotic Methods in Statistical Decision Theory_.  Springer Series in Statistics. Springer New York, 1986.
  * Lopes (2022a) M. E. Lopes.  Central limit theorem and bootstrap approximation in high dimensions: Near $1/\sqrt{n}$ rates via implicit smoothing.  _The Annals of Statistics_ , 50(5):2492 – 2513, 2022a.
  * Lopes (2022b) M. E. Lopes.  Improved rates of bootstrap approximation for the operator norm: A coordinate-free approach, 2022b.
  * Lopes et al. (2019) M. E. Lopes, A. Blandino, and A. Aue.  Bootstrapping spectral statistics in high dimensions.  _Biometrika_ , 106(4):781–801, 09 2019.
  * Lopes et al. (2020) M. E. Lopes, Z. Lin, and H.-G. Müller.  Bootstrapping max statistics in high dimensions: Near-parametric rates under weak variance decay and application to functional and multinomial data.  _Annals of Statistics_ , 48(2):1214–1229, 04 2020\.
  * Lopes et al. (2023) M. E. Lopes, N. B. Erichson, and M. W. Mahoney.  Bootstrapping the operator norm in high dimensions: Error estimation for covariance matrices and sketching.  _Bernoulli_ , 29(1):428 – 450, 2023.
  * Mendelson (2002) S. Mendelson.  Geometric parameters of kernel machines.  In _International conference on computational learning theory_ , pages 29–43. Springer, 2002.
  * Naumov et al. (2019) A. Naumov, V. Spokoiny, and V. Ulyanov.  Bootstrap confidence sets for spectral projectors of sample covariance.  _Probability Theory and Related Fields_ , 174(3):1091–1132, 2019.
  * Nourdin and Peccati (2012) I. Nourdin and G. Peccati.  _Normal approximations with Malliavin calculus: from Stein’s method to universality_.  Number 192. Cambridge University Press, 2012.
  * Péché (2009) S. Péché.  Universality results for the largest eigenvalues of some sample covariance matrix ensembles.  _Probability Theory and Related Fields_ , 143(3):481–516, 2009.
  * Raič (2019) M. Raič.  A multivariate Berry–Esseen theorem with explicit constants.  _Bernoulli_ , 25(4A):2824–2853, 11 2019.
  * Rockafellar (1999) R. T. Rockafellar.  Second-order convex analysis.  _J. Nonlinear Convex Anal_ , 1(1-16):84, 1999\.
  * Schölkopf and Smola (2002) B. Schölkopf and A. Smola.  _Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond_.  Adaptive computation and machine learning. MIT Press, 2002.
  * Sedaghat (2003) H. Sedaghat.  _Nonlinear Difference Equations: Theory with Applications to Social Science Models_.  Springer, 2003.
  * Shi (1992) S. G. Shi.  Accurate and efficient double-bootstrap confidence limit method.  _Computational Statistics and Data Analysis_ , 13(1):21–32, 1992.
  * Singh and Vijaykumar (2023) R. Singh and S. Vijaykumar.  Kernel ridge regression inference, 2023.
  * Steinwart and Christmann (2008) I. Steinwart and A. Christmann.  _Support Vector Machines_.  Springer, 2008.
  * Tanguy (2017) K. Tanguy.  _Quelques inégalités de superconcentration: théorie et applications_.  PhD thesis, Université Paul Sabatier-Toulouse III, 2017.
  * Unser (2019) M. Unser.  A representer theorem for deep neural networks.  _Journal of Machine Learning Research_ , 20, 2019.
  * van der Vaart and Wellner (1996) A. W. van der Vaart and J. A. Wellner.  _Weak Convergence and Empirical Processes_.  Springer, 1996.
  * Vershynin (2018) R. Vershynin.  _High-Dimensional Probability: An Introduction with Applications in Data Science_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2018.
  * Wahba (1990) G. Wahba.  _Spline Models for Observational Data_.  SIAM, 1990.
  * Waternaux (1976) C. M. Waternaux.  Asymptotic distribution of the sample roots for a nonnormal population.  _Biometrika_ , 63(3):639–645, 1976.
  * Yurinsky (2006) V. Yurinsky.  _Sums and Gaussian Vectors_.  Springer, 2006.
  * Zhang and Cheng (2017) X. Zhang and G. Cheng.  Simultaneous inference for high-dimensional linear models.  _Journal of the American Statistical Association_ , 112(518):757–768, 2017.

Supplementary Materials for “Gaussian and Bootstrap Approximations for
Empirical Processes” Alexander Giessing11footnotemark: 1

Contents

[sections] [sections]l1

##  Appendix A Two toy examples

In this section we present two examples to illustrate the limitations of the
existing and the advantage of the new Gaussian approximation results. The
presentation is deliberately expository; we omit proofs as much as possible.
In both examples, we take $\mathcal{F}_{n}=\\{x\mapsto f(x)=x^{\prime}u:u\in
S^{d-1}\\}$ with $S^{d-1}=\\{u\in\mathbb{R}^{d}:\|u\|_{2}=1\\}$, which plays a
role in the construction of high-dimensional confidence regions and multiple
testing problems (see Section 4.1).

Consider a simple random sample $X_{1},\ldots,X_{n}\in\mathbb{R}^{d}$ from the
law of $X=a\cdot\xi$, where $a\in\mathbb{R}^{d}$ is a fixed vector and
$\xi\in\mathbb{R}$ a centered random variable with finite third moment. Then,
$\sup_{f\in\mathcal{F}_{n}}|\mathbb{G}_{n}(f)|=\|n^{-1/2}\sum_{i=1}^{n}X_{i}\|_{2}=\|a\|_{2}|n^{-1/2}\sum_{i=1}^{n}\xi_{i}|$.
Whence, the supremum of the empirical process reduces to the average of i.i.d.
scalar-valued random variables and the classical univariate Berry-Essen
theorem yields the following bound on the Kolmogorov distance:

$\displaystyle\varrho_{n}\lesssim\frac{\mathrm{E}[|\xi|^{3}]}{\sqrt{n}\>\mathrm{E}[|\xi|^{2}]^{3/2}}\rightarrow
0\quad\mathrm{as}\quad n\rightarrow\infty.$ (17)

This non-asymptotic bound is independent of the dimension $d$ and the weak
variance of the Gaussian proxy process. If we ignore the low-rank structure of
the data and instead use the Gaussian approximation results in Chernozhukov et
al. (2014b) (Theorem 2.1 combined with Lemma 2.3) we obtain (qualitatively)

$\displaystyle\varrho_{n}\lesssim_{\star}\frac{\mathrm{polylog}\left(N\left(\mathcal{F}_{n},e_{P},\varepsilon\right)\vee
n\right)}{n^{1/6}\sqrt{\inf_{u\in S^{d-1}}\mathrm{Var}(Z^{\prime}u)}},$ (18)

where $Z\sim N\left(0,\mathrm{E}[\xi^{2}]\cdot aa^{\prime}\right)$ and
$\lesssim_{\star}$ hides a complicated multiplicative factor which, among
other things, depends on the sample size $n$, third and higher moments of
$\|Z\|_{2}$, and the inverse of the discretization level $\varepsilon>0$ (see
Section 3.6). Since the covariance matrix $\mathrm{E}[\xi^{2}]\cdot
aa^{\prime}$ has rank one, the unit ball in $\mathcal{F}_{n}$ with respect to
the intrinsic standard deviation metric $e_{P}$ is isometrically isomorphic to
the interval $[-1,1]\subset\mathbb{R}$. Hence, the metric entropy can be upper
bounded independently of the dimension $d$; in particular, we have $\log
N\left(\mathcal{F}_{n},e_{P},\varepsilon\right)\leq\log(1+2/\varepsilon)$.
However, the low-rank structure of the data also implies that the weak
variance of the associated Gaussian proxy process vanishes; indeed,
$\inf_{u\in
S^{d-1}}\mathrm{Var}(Z^{\prime}u)=\mathrm{E}[\xi^{2}]\cdot\inf_{u\in
S^{d-1}}(a^{\prime}u)^{2}=0$. Thus, the upper bound in (18) is in fact invalid
(or trivial if we interpret $1/0=\infty$) and fails to replicate the
univariate Berry-Esseen bound in (17).

In contrast, the new Gaussian approximation inequality (Theorem 1 in Section
3.2; Theorem A.1 in Giessing and Fan (2023)) is agnostic to the covariance
structure of the data and yields

$\displaystyle\varrho_{n}\lesssim\frac{(\mathrm{E}[\|X\|_{2}^{3}])^{1/3}}{n^{1/6}\sqrt{\mathrm{Var}(\|Z\|_{2})}}+\frac{\mathrm{E}\left[\|X\|_{2}^{3}\mathbf{1}\\{\|X\|_{2}^{3}>n\>\mathrm{E}[\|X\|_{2}^{3}]\\}\right]}{\mathrm{E}\left[\|X\|_{2}^{3}\right]}+\frac{\mathrm{E}[\|Z\|_{2}]}{\sqrt{n\mathrm{Var}(\|Z\|_{2})}},$
(19)

where $\lesssim$ hides an absolute constant independent of $n,d$, and the
distribution of the $X_{i}$’s and $Z\sim N\left(0,\mathrm{E}[\xi^{2}]\cdot
aa^{\prime}\right)$. The upper bound in this inequality depends only on the
third moment of the envelop function $x\mapsto\|x\|_{2}$ and the strong
variance of the Gaussian proxy process $\mathrm{Var}(\|Z\|_{2})$. If we use
that $\|X\|_{2}=\|a\|_{2}\cdot|\xi|$ and $\|Z\|_{2}=\|a\|_{2}\cdot|g|$ for
$g\sim N(0,1)$, we obtain

$\displaystyle\varrho_{n}\lesssim\frac{\mathrm{E}[|\xi|^{3}]^{1/3}}{n^{1/6}\>\mathrm{E}[|\xi|^{2}]^{1/2}}+\frac{\mathrm{E}\left[|\xi|^{3}\mathbf{1}\\{|\xi|^{3}>n\>\mathrm{E}[|\xi|^{3}]\\}\right]}{\mathrm{E}\left[|\xi|^{3}\right]}\rightarrow
0\quad\mathrm{as}\quad n\rightarrow\infty.$ (20)

This inequality is obviously dimension- and weak variance-free. In this sense,
it recovers the essential feature of the univariate Berry-Esseen bound in (17)
and improves over the bound in (18). The dependence on the sample size $n$ and
the moments of $\xi$ is still sub-optimal; but refinements in this direction
are beyond the scope of this paper. Obviously, in this example, we have chosen
a rank one covariance matrix only to be able to compare the Gaussian
approximation result with the Berry-Esseen theorem. Any low-rank structure
implies a vanishing weak variance and, hence, a breakdown of the results in
Chernozhukov et al. (2014b).

Next, suppose that the data $X_{1},\ldots,X_{n}\in\mathbb{R}^{d}$ are a simple
random sample drawn from the law of a random vector
$X=(X^{(1)},\ldots,X^{(d)})^{\prime}$ with mean zero, element-wise bounded
entries $\max_{1\leq k\leq d}|X^{(k)}|<B$ almost surely for some $B>0$, and
equi-correlated covariance matrix
$\Sigma=(1-\rho)I_{d}+\rho\mathbf{1}_{d}\mathbf{1}_{d}^{\prime}$ for some
$\rho\in(-1/(d-1),1)$. The constraints on $\rho$ guarantee that $\Sigma$ has
full rank. Therefore, the multivariate Berry-Esseen theorem by Bentkus (2003)
(Theorem 1.1) implies

$\displaystyle\varrho_{n}\lesssim\frac{d^{1/4}\>\mathrm{E}[\|\Sigma^{-1/2}X\|_{2}^{3}]}{\sqrt{n}}.$
(21)

This upper bound is not useful in high-dimensional settings because the
expected value is polynomial in the dimension $d$, i.e.
$\mathrm{E}[\|\Sigma^{-1/2}X\|_{2}^{3}]\geq\mathrm{E}[\|\Sigma^{-1/2}X\|_{2}^{2}]^{3/2}=d^{3/2}$.

The Gaussian approximation results by Chernozhukov et al. (2014b) yield
(again) inequality (18). Since the covariance matrix $\Sigma$ has full rank,
the weak variance is now strictly positive and equal to the smallest
eigenvalue of $\Sigma$, i.e. $\inf_{u\in
S^{d-1}}\mathrm{Var}(Z^{\prime}u)=\inf_{u\in S^{d-1}}u^{\prime}\Sigma
u=1-\rho>0$. However, in this example, the metric entropy with respect to the
intrinsic standard deviation metric $e_{P}$ poses a problem. Indeed, let
$\lambda_{1}=1+(d-1)\rho$ and $\lambda_{2}=\ldots=\lambda_{d}=1-\rho$ be the
eigenvalues of $\Sigma$. Then, the unit ball in $\mathcal{F}_{n}$ with respect
to $e_{P}$ can be identified with the weighted Euclidean ball
$B_{\lambda}^{d}(0,\varepsilon):=\\{u\in\mathbb{R}^{d}:\sum_{k=1}^{d}\lambda_{k}u_{k}\leq\varepsilon^{2}\\}$.
Let $B^{d}(0,1)$ be the standard unit ball in $\mathbb{R}^{d}$ with respect to
the Euclidean distance. Then, $B^{d}(0,\varepsilon)\subseteq
B_{\lambda}^{d}(0,\varepsilon)$ whenever $\varepsilon<1-\rho$. Therefore,
standard arguments
$N\left(\mathcal{F}_{n},e_{P},\varepsilon\right)\geq\mathrm{vol}\left(B^{d}(0,1)\right)/\mathrm{vol}\left(B_{\lambda}^{d}(0,\varepsilon)\right)\geq\mathrm{vol}\left(B^{d}(0,1)\right)/\mathrm{vol}\left(B^{d}(0,\varepsilon)\right)\geq\left((1-\rho)/\varepsilon\right)^{d}$.
Thus, the metric entropy grows linear in the dimension $d$, i.e. $\log
N\left(\mathcal{F}_{n},e_{P},\varepsilon\right)\geq
d\log\left((1-\rho)/\varepsilon\right)$ for $\varepsilon\downarrow 0$. We
conclude that the results in Chernozhukov et al. (2014b) are again not useful
in high dimensions with $d\gg n$.

The new Gaussian approximation inequality implies (again) (19). Since
$\max_{1\leq k\leq d}|X^{(k)}|<B$ almost surely and $\mathrm{tr}(\Sigma)=d$ it
follows that
$\sqrt{d}=(\mathrm{E}[\|Z\|_{2}^{2}])^{1/2}=(\mathrm{E}[\|X\|_{2}^{2}])^{1/2}\leq(\mathrm{E}[\|X\|_{2}^{3}])^{1/3}\leq
B\sqrt{d}$. Moreover, by Theorem A.6 in Giessing and Fan (2023) and since
$\mathrm{tr}(\Sigma^{2})=d(1-\rho)^{2}+\rho^{2}d^{2}$, we have
$\mathrm{Var}(\|Z\|_{2})\geq\mathrm{tr}(\Sigma^{2})/\mathrm{tr}(\Sigma)=(1-\rho)^{2}+\rho^{2}d$.
Therefore, inequality (19) simplifies to

$\displaystyle\varrho_{n}\lesssim\frac{B}{n^{1/6}\rho}+B^{3}\mathbf{1}\\{B>n\\}\rightarrow
0\quad\mathrm{as}\quad n\rightarrow\infty.$ (22)

This inequality is not only dimension- and weak variance-free but also
improves qualitatively over, both, the results by Chernozhukov et al. (2014b)
and Bentkus (2003). Intuitively, the reason why we are able to shed the
$d^{1/4}$-factor from the upper bound compared to the results in Bentkus
(2003) is that we only take the supremum over all Euclidean balls with center
at the origin whereas he takes the supremum over all convex sets in
$\mathbb{R}^{d}$. (Note that to apply his bound in the context of this
example, we need to take the supremum over at least all weighted Euclidean
balls
$B_{\lambda}^{d}(0,\varepsilon):=\\{u\in\mathbb{R}^{d}:\sum_{k=1}^{d}\lambda_{k}u_{k}\leq\varepsilon^{2}\\}$.)
This second example is related to the first one in so far as the covariance
matrix has “approximately” rank one. Indeed, as the dimension increases the
law of the random vector $X$ concentrates in the neighborhood of the one-
dimensional subspace spanned by eigenvector associated to the largest
eigenvalue of $\Sigma$. This becomes even more obvious if we consider the
standardized covariance $\Sigma/\|\Sigma\|_{op}$ with eigenvalues
$\lambda_{1}=\rho+(1-\rho)/d\rightarrow\rho$ and
$\lambda_{2}=\ldots\lambda_{d}=(1-\rho)/d\rightarrow 0$ as
$d\rightarrow\infty$.

In both examples, the exact and approximate low-rank structures of the data
are crucial in order to go from the abstract bound (19) to the dimension-/
entropy-free bounds (20) and (22), respectively. This is not coincidental and
is in fact representative for the entire theory that we develop in this paper:
While our Gaussian and bootstrap approximation inequalities hold without
assumptions on the metric entropy, in concrete examples they often only yield
entropy-free upper bounds if the trace of the covariance operator is bounded
or grows at a much slower rate than the sample size $n$ (e.g. low-rank,
bounded effective rank, or variance decay, see Section 4). This is certainly a
limitation of our theory, but a low-rank covariance (function) is an
empirically well-documented property of many data sets and therefore a common
assumption in multivariate and high-dimensional statistics (e.g. Anderson,
2003; Vershynin, 2018).

##  Appendix B Auxiliary results

###  B.1 Smoothing inequalities and partial derivatives

######  Lemma 1.

Let $X,Z\in\mathbb{R}$ be arbitrary random variables. There exists a map
$h_{s,\lambda}\in C^{\infty}_{c}(\mathbb{R})$ such that

  * (i)

for all $s,t\in\mathbb{R}$ and $\lambda>0$,

$\displaystyle|D^{k}h_{s,\lambda}(t)|\leq C_{k}\lambda^{-k}\mathbf{1}\\{s\leq
t\leq s+3\lambda\\},$

where $C_{k}>0$ is a constant depending only on $k\in\mathbb{N}_{0}$; and

  * (ii)

for all $\lambda>0$,

$\displaystyle\left|\sup_{s\in\mathbb{R}}\left|\mathbb{P}\left\\{X\leq
s\right\\}-\mathbb{P}\left\\{Z\leq
s\right\\}\right|-\sup_{s\in\mathbb{R}}\big{|}\mathrm{E}[h_{s,\lambda}(X)-h_{s,\lambda}(Z)]\big{|}\right|\leq\zeta_{3\lambda}(X)\wedge\zeta_{3\lambda}(Z),$

where $\zeta_{\lambda}(V):=\sup_{s\in\mathbb{R}}\mathbb{P}\\{s\leq V\leq
s+\lambda\\}$ for real-valued $V\in\mathbb{R}$.

######  Remark 1.

We can take $C_{0}=C_{1}=1$ (see proof).

######  Lemma 2.

Let $h_{s,\lambda}\in C^{\infty}_{0}(\mathbb{R})$ be the map from Lemma 1 and
define $x\mapsto h(x):=h_{s,\lambda}(\|x\|_{\infty})$ for
$x\in\mathbb{R}^{d}$. Let $(P_{t})_{t\geq 0}$ be the Ornstein-Uhlenbeck semi-
group with stationary measure $N(0,\Sigma)$ and positive definite covariance
matrix $\Sigma$.

  * (i)

For arbitrary indices $1\leq i_{1},\ldots,i_{k}\leq d$, $k\geq 1$, and all
$x_{0}\in\mathbb{R}^{d}$,

$\displaystyle\frac{\partial^{k}}{\partial x_{i_{1}}\cdots\partial
x_{i_{k}}}\left(\int_{0}^{\infty}P_{t}h(x)dt\right)\Big{|}_{x=x_{0}}=\int_{0}^{\infty}e^{-kt}P_{t}\left(\frac{\partial^{k}h}{\partial
x_{i_{1}}\cdots\partial x_{i_{k}}}\right)(x_{0})dt;$

  * (ii)

for almost every $x_{0}\in\mathbb{R}^{d}$ the absolute value of the derivative
in (i) can be upper bounded by

$\displaystyle
C_{k}\lambda^{-k}\int_{0}^{\infty}e^{-kt}\mathrm{E}\left[\mathbf{1}\left\\{s\leq\|V_{0}^{t}\|_{\infty}\leq
s+3\lambda\right\\}\mathbf{1}\left\\{|V_{0{i_{1}}}^{t}|\geq|V_{0\ell}^{t}|,\>\ell\neq{i_{1}}\right\\}\right]dt\>\mathbf{1}\\{i_{1}=\ldots=i_{k}\\},$

where $V_{0}^{t}=e^{-t}x_{0}+\sqrt{1-e^{-2t}}Z$, $Z\sim N(0,\Sigma)$, and
$C_{k}>0$ is the absolute constant from Lemma 1.

######  Remark 2.

While claim (i) looks a lot like a “differentiating under the integral sign”
type result, it is more accurate to think of it as a specific smoothing
property of the Ornstein-Uhlenbeck semigroup with stationary measure
$N(0,\Sigma)$ when applied to (compositions of Lipschitz continuous functions
with) the map $x\mapsto\|x\|_{\infty}$.

######  Lemma 3.

Let $\varrho\in L^{1}(\mathbb{R})$ with $\int\varrho(r)dr=1$. For $\eta>0$ and
a map $h$ on $\mathbb{R}^{d}$ set

$\displaystyle(\varrho_{\eta}\ast_{(i)}h)(x):=\int\varrho(r)h(x-r\eta
e_{i})dr,$

where $e_{i}$ denotes the $i$th standard unit vector in $\mathbb{R}^{d}$. For
$k\in\mathbb{N}$, $f_{1},\ldots,f_{k}\in B(\mathbb{R}^{d})$, $g\in
L^{1}(\mathbb{R}^{d})$, and arbitrary indices $1\leq i_{1},\ldots,i_{k}\leq
d$,

$\displaystyle\left\|\prod_{j=1}^{k}(\varrho_{\eta}\ast_{(i_{j})}f_{j})g-\prod_{j=1}^{k}f_{j}g\right\|_{1}\rightarrow
0\quad{}\mathrm{as}\quad{}\eta\rightarrow 0.$

######  Remark 3.

The proof of this result on “partially regularized” functions is similar to
the one on “fully regularized” functions (Folland, 1999, Theorem 8.14 (a)).
The conditions on the functions $f_{1},\ldots,f_{k},g,\varrho$ can probably be
relaxed, but they are sufficiently general to apply to the situations that we
encounter.

######  Lemma 4 (Partial derivatives of compositions of almost everywhere
diff’able functions).

Let $g\in C^{k}(\mathbb{R})$ and $f\in
C^{k}(\mathbb{R}^{d}\setminus\mathcal{N})$ and $\mathcal{N}$ is a null set
with respect to the Lebesgue measure on $\mathbb{R}^{d}$. Then, $g\circ f\in
C^{k}(\mathbb{R}^{d}\setminus\mathcal{N})$ and, for arbitrary indices $1\leq
i_{1},\ldots,i_{k}\leq d$, $k\geq 1$, and all
$x\in\mathbb{R}^{d}\setminus\mathcal{N}$,

$\displaystyle\frac{\partial^{k}(g\circ f)}{\partial x_{i_{1}}\cdots\partial