Improved Central Limit Theorem and bootstrap approximations in high dimensions

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

99,742 characters

Improved Central Limit Theorem and Bootstrap Approximations in High Dimensions



\begin{frontmatter}
\title{Improved Central Limit Theorem and Bootstrap Approximations in High Dimensions}
\runtitle{Improved CLT and bootstrap in high dimensions}

\begin{aug}
\author[A]{\fnms{Victor} \snm{Chernozhuokov}\ead[label=e1]{[email removed]}},
\author[B]{\fnms{Denis} \snm{Chetverikov}\ead[label=e2,mark]{[email removed]}}
\author[C]{\fnms{Kengo} \snm{Kato}\ead[label=e3,mark]{[email removed]}} \\
\and
\author[D]{\fnms{Yuta} \snm{Koike}\ead[label=e4,mark]{[email removed]}}
\address[A]{
Department of Economics and Operations Research Center, MIT,
\printead{e1}}

\address[B]{
Department of Economics, UCLA,
\printead{e2}}

\address[C]{
Department of Statistics and Data Science, Cornell University,
\printead{e3}}

\address[D]{Department,
Mathematics and Informatics Center and Graduate School of Mathematical Sciences, The University of Tokyo,
\printead{e4}}

\end{aug}

\begin{abstract}
This paper deals with the Gaussian and bootstrap approximations to the distribution of the max statistic in high dimensions. This statistic takes the form of the maximum over components of the sum of independent random vectors and its distribution plays a key role in many high-dimensional estimation and testing problems. Using a novel iterative randomized Lindeberg method, the paper derives new bounds for the distributional approximation errors. These new bounds substantially improve upon existing ones and simultaneously allow for a larger class of bootstrap methods.
\end{abstract}

\begin{keyword}[class=MSC2020]
\kwd{60F05, 62E17}
\end{keyword}

\begin{keyword}
\kwd{bootstrap}
\kwd{central limit theorem}
\kwd{iterative randomized Lindeberg method}
\kwd{Stein kernel}
\end{keyword}

\end{frontmatter}










\section{Introduction}

Let $X_1,\dots,X_n$ be independent random vectors in $\mathbb R^p$ such that
${\mathrm{E}}[X_{i j}] = \mu_{j}$ for all $i=1,\dots,n$ and $j=1,\dots,p$,
where $X_{i j}$ denotes the $j$th component of the vector $X_i$. We are interested in approximating the distribution of the maximum coordinate of the centered sample mean of $X_{1},\dots,X_{n}$, i.e.,
\begin{equation}\label{eq: statistic first}
T_n = \max_{1\leq j\leq p}\frac{1}{\sqrt n}\sum_{i=1}^n(X_{i j} - \mu_{j}).
\end{equation}
The distribution of $T_n$ plays a particularly important role in many high-dimensional settings, where $p$ is potentially larger or much larger than $n$. For example, it appears in selecting the regularization parameters for the Lasso estimator and the Dantzig selector (\cite{CCK13}), in carrying out reality checks for data snooping and testing superior predictive ability (\cite{W00,H05}), in constructing model confidence sets (\cite{HLN11}), in testing conditional and/or many unconditional moment inequalities (\cite{BSS19, C18, CCK19, KB19}), in multiple testing with the family-wise error rate control (\cite{BCCHK18}), in constructing simultaneous confidence intervals for high-dimensional parameters (\cite{BCCW18}), in adaptive testing of regression and stochastic monotonicity (\cite{C19,CWK18}), in carrying out inference on generalized instrumental variable models (\cite{CR19}), and in constructing Lepski-type procedures for adaptive estimation and inference in nonparametric problems (\cite{CCK14}); more references can be found in \cite{DZ17} and especially in \cite{BCCHK18}. It is therefore of great interest to develop methods for obtaining feasible and accurate approximations to the distribution of $T_n$, allowing for the high-dimensional $p \gg n$ case.


Toward this goal, the first three authors of this paper obtained the following Gaussian approximation result in \cite{CCK13, CCK17}. Let $G = (G_1,\dots,G_p)'$ be a Gaussian random vector in $\mathbb R^p$ with mean $\mu = (\mu_1,\dots,\mu_p)'$ and covariance matrix $n^{-1}\sum_{i=1}^n{\mathrm{E}}[(X_i - \mu)(X_i - \mu)']$ and let the critical value $c_{1-\alpha}$ be the $(1-\alpha)$th quantile of $\max_{1\leq j\leq p}G_j$. Then under mild regularity conditions,
\begin{equation}\label{eq: basic approximation}
\Big|{\mathrm{P}}(T_n > c_{1-\alpha}) - \alpha \Big|\leq C\left(\frac{\log^7(p n)}{n}\right)^{1/6},
\end{equation}
where $C$ is a constant that is independent of $n$ and $p$. This result is important because the right-hand side of the bound \eqref{eq: basic approximation} depends on $p$ only via the logarithm of $p$, and hence it shows that the Gaussian approximation holds if $\log  p = o(n^{1/7})$, which allows $p$ to be much larger than $n$. Besides, building upon this result, the same authors have proved bounds similar to (\ref{eq: basic approximation}) for the critical values obtained by the Gaussian multiplier and empirical bootstraps in \cite{CCK17}.

Gaussian approximation of the form (\ref{eq: basic approximation}) allows us to develop powerful inference methods for high-dimensional data in applications discussed above and has stimulated further developments into dependent data \cite{ZW17,ZC17,CCK19}, $U$-statistics \cite{C18,CK19a,CK19b}, Malliavin calculus \cite{C19}, and homogeneous sums \cite{K19a}. Despite such rapid developments, the literature has left much to be desired on  coherent understanding of sharpness of the bound (\ref{eq: basic approximation}) for the Gaussian or bootstrap critical values since the first appearance of \cite{CCK17} in 2014 on arXiv. The problem can be decomposed into two parts: (i) sharpness of the bound in terms of dependence on $n$ and (ii) sharpness of the bound in terms of dependence on $p$.

There are two important developments toward the question of sharpness of the bound (\ref{eq: basic approximation}) that should be mentioned.
First, Deng and Zhang \cite{DZ17} considered direct bootstrap approximation without taking the root of Gaussian approximation, and proved the following bound for the critical value $c_{1-\alpha}$ obtained by the empirical or  third-order matching (or Mammen's \cite{M93}) multiplier bootstraps:
\begin{equation}\label{eq: dz bound}
\Big|{\mathrm{P}}(T_n > c_{1-\alpha}) - \alpha \Big|\leq C\left(\frac{\log^5(p n)}{n}\right)^{1/6}.
\end{equation}
Their bound improves the power of the logs in the previous bound (\ref{eq: basic approximation}), showing that the empirical and Mammen's bootstraps are consistent to approximate the distribution of $T_{n}$ if $\log p = o(n^{1/5})$ instead of $\log p = o(n^{1/7})$.
Second, the recent preprint by the fourth author \cite{K19b} shows that the same bound (\ref{eq: dz bound}) indeed holds for the Gaussian critical value as well.

In turn, in this paper, we show that in fact a much larger improvement is possible: under mild regularity conditions, we prove that
\begin{equation}\label{eq: empirical bootstrap introduction}
\Big|{\mathrm{P}}(T_n > c_{1-\alpha}) - \alpha \Big|\leq C\left(\frac{\log^5(p n)}{n}\right)^{1/4},
\end{equation}
both for the Gaussian and  bootstrap critical values $c_{1-\alpha}$. In comparison with the Gaussian approximation result \eqref{eq: basic approximation}, our new bound improves not only the power of the logs but also the power of the sample size $n$. Moreover, regarding the bootstrap types, we allow for not only the empirical and third-order matching multiplier bootstrap methods, but also for general multiplier bootstrap methods (with i.i.d weights), which match only two moments of the data, such as the multiplier bootstrap methods with Gaussian and Rademacher weights.

We remark that several authors have recently pointed out that an additional structural assumption on the covariance matrices of $X_i$'s can improve the bound \eqref{eq: empirical bootstrap introduction}. In particular, Fang and Koike \cite{FK20} showed that the right-hand side of \eqref{eq: empirical bootstrap introduction} can be improved to $C(\log^4(pn)/n)^{1/3}$ when the covariance matrices are non-degenerate and can be further improved to $C(\log^3(p)/n)^{1/2}\log n$ when we additionally assume that $X_i$'s have log-concave densities. The latter result is based on the fact that random vectors with log-concave densities admit Stein kernels with sub-Weibull entries, which is established by Fathi in \cite{Fa19}. Moreover, building on the important results by Lopes in \cite{L20} and Kuchibhotla and Rinaldo in \cite{KR20}, \cite{CCK20} showed that the bound $C(\log^3(p)/n)^{1/2}\log n$ can be achieved even without the assumption of log-concave densities (non-degenerate covariance matrices are still required; \cite{L20} and \cite{KR20} were the first to obtain dependence on $n$ via $1/\sqrt n$ in \eqref{eq: empirical bootstrap introduction} without requiring log-concave densities). In addition, Lopes, Lin and M\"uller \cite{LLM20} showed that the right-hand side of \eqref{eq: empirical bootstrap introduction} can be improved to $C n^{-1/2+\delta}$ for any $\delta>0$ when the coordinates of $X_i$'s have decaying variances. Compared to these results, our bound requires neither non-degenerate covariance matrices nor decaying variances.

In addition, we prove that if the distribution of the random vectors $X_1,\dots,X_n$ is symmetric around the mean, then even better approximation to the distribution of $T_n$ is possible:
\begin{equation}\label{eq: randomization introduction}
\Big|{\mathrm{P}}(T_n > c_{1-\alpha}) - \alpha \Big|\leq C\left(\frac{\log^3(p n)}{n}\right)^{1/2}
\end{equation}
as long as the critical value $c_{1-\alpha}$ is obtained via the multiplier bootstrap method with Rademacher weights. This new bound makes Rademacher weights particularly appealing in the high-dimensional settings, at least from a theoretical perspective.

We also consider bootstrap approximations with incremental factors, previously used by Andrews and Shi in \cite{AS13} in the context of testing conditional moment inequalities. Specifically, for a small but fixed constant $\eta > 0$, called an incremental factor, we derive the following bounds:
\begin{equation}\label{eq: or 1}
{\mathrm{P}}(T_n > c_{1-\alpha} + \eta) - \alpha \leq C\left(\frac{\log^3(p n)}{n}\right)^{1/2}
\end{equation}
if $c_{1-\alpha}$ is obtained via either the empirical or the third-order matching multiplier bootstrap methods and
\begin{equation}\label{eq: or 2}
{\mathrm{P}}(T_n > c_{1-\alpha} + \eta) - \alpha \leq C\left(\frac{\log^5(p n)}{n}\right)^{1/2}
\end{equation}
if $c_{1-\alpha}$ is obtained via general multiplier bootstrap methods, where the constant $C$ may depend on $\eta$. Even though these are one-sided bounds, they are useful because they show that in any test based on the statistic $T_n$, increasing the critical value $c_{1-\alpha}$ by an incremental factor $\eta$ may substantially reduce the sample complexity for over-rejection. Namely,  assuming $\log p \gtrsim \log n$ for simplicity, for the over-rejection probability to be less than or equal to a given level $0 < \Delta < 1-\alpha$, the empirical bootstrap or multiplier bootstrap (without incremental factor) requires $n \gtrsim \Delta^{-4} \log^5 p$, while adding a constant incremental factor reduces the sample complexity to $n \gtrsim (\Delta^{-2} \log^3p) \vee \log^5 p$ if we use the empirical or third-order matching bootstrap.
It is worth noting that, given that in high-dimensional settings, where $p$ is rapidly increasing together with $n$, $c_{1-\alpha}$ is typically also getting large as we increase $n$, adding an incremental factor $\eta$ may not have a large impact on the power properties of the test.

In fact, all our results apply to a more general version of the statistic $T_n$:
\begin{equation}\label{eq: statistic second}
T_n = \max_{1\leq j\leq p}\frac{1}{\sqrt n}\sum_{i=1}^n(X_{i j} - \mu_{j} + a_j),
\end{equation}
where $a=(a_1,\dots,a_p)'$ is a vector in $\mathbb R^p$, which reduces to \eqref{eq: statistic first} if we set $a = 0_p$. In most applications mentioned above, the former version \eqref{eq: statistic first} is sufficient but there are some applications where the more general version \eqref{eq: statistic second} is required; for example, the latter was used by Bai, Shaikh, and Santos in \cite{BSS19} to extend the method of testing moment inequalities proposed in \cite{RSW14} for the case of a small number of inequalities to the case of a large number of inequalities. For the rest of the paper, we will therefore work with the more general version \eqref{eq: statistic second} of the statistic $T_n$. In addition, we emphasize that our results can be equally applied with
$$
T_n = \max_{1\leq j\leq p}\left|\frac{1}{\sqrt n}\sum_{i=1}^n(X_{i j} - \mu_j + a_j)\right|
$$
by replacing the $p$-dimensional vectors $X_i - \mu + a$ with the $2p$-dimensional vectors whose first $p$ components are equal to $X_i - \mu + a$ and the last $p$ components are equal to $-(X_i - \mu + a)$.

To prove \eqref{eq: empirical bootstrap introduction}, we develop a novel and iterative version of the randomized Lindeberg method. A key feature of our approach is that we carry out a careful analysis of the coefficients in the Taylor expansion underlying the Lindeberg method. In particular, we apply the Lindeberg method iteratively in combination with an anti-concentration inequality for maxima of Gaussian processes to bound these coefficients, which substantially improves upon the original randomized Lindeberg method proposed in \cite{DZ17}. In addition, we sharpen the Gaussian approximation bounds for the multiplier processes developed in \cite{K19b} using Stein's kernels. In turn, to prove \eqref{eq: randomization introduction}, we establish a new connection between the Rademacher bootstrap and the randomization tests, as discussed in \cite{LR05}, using a recent result from the computer science literature on pseudo-random number generators by O'Donnell, Servedio, and Tan \cite{OST18}, which provides an anti-concentration inequality for maxima of Rademacher processes. Finally, to prove error bounds  \eqref{eq: or 1} and \eqref{eq: or 2}, we apply the original randomized Lindeberg method as developed in \cite{DZ17}.


Finally, we conduct a small scale simulation study. Our simulation study shows that (i) all bootstrap methods considered in this paper perform reasonably well in high dimensions; (ii) for asymmetric distributions, the empirical and
the third-order matching multiplier bootstrap methods outperform the multiplier bootstrap
methods with Gaussian and Rademacher weights; and (iii) for symmetric distributions, the multiplier bootstrap with Rademacher weights performs the best, which is consistent with Theorem \ref{thm: gauss and rademacher} ahead. See the Supplementary Material for details.







The rest of the paper is organized as follows. In the next section, we present our main results. In Section \ref{sec: main arguments}, we develop the iterative randomized Lindeberg method, which is the first key component in deriving our main results. In Section \ref{sec: stein kernels}, we provide new bounds for the Gaussian approximations using Stein's kernels, which is the second key component in deriving our main results. In Section \ref{sec: proofs}, we give proofs of the main results. In the Supplemental Material, we collect additional derivations and conduct a small simulation study.


\subsection{Notation} For any vectors $x,y\in\mathbb R^p$ and any scalar $c\in\mathbb R$, we write $x\leq y$ if $x_j \leq y_j$ for all $j=1,\dots,p$ and  write $x + c$ to denote the vector in $\mathbb R^p$ whose $j$th component is $x_j + c$ for all $j = 1,\dots,p$. Also, for any sequences of scalars $\{a_n\}_{n\geq 1}$ and $\{b_n\}_{n\geq 1}$ we write $a_n\lesssim b_n$ if $a_n \leq C b_n$ for all $n\geq 1$  for some constant $C$. Recall that,  for any random variable $T$ and a constant $\gamma\in(0,1)$, the $\gamma$th quantile of $T$ is defined as  $\inf\{t\in\mathbb R\colon  {\mathrm{P}}(T\leq t) \ge \gamma \}$. Finally, we use the notation $X_{1:n} = (X_1,\dots,X_n)$.

\section{Main Results}\label{sec: main results}
In this section, we present our main results. We first formally define all the critical values $c_{1-\alpha}$ to be used throughout the paper. We then discuss the required regularity conditions and present the results.

\subsection{Gaussian and Bootstrap Critical Values}

First, define the Gaussian critical value $c_{1-\alpha}^G$ as the $(1-\alpha)$th quantile of
\begin{equation}\label{eq: gaussian analog}
T_n^G = \max_{1\leq j\leq p}(G_j + a_j),
\end{equation}
where $G$ is a centered Gaussian random vector in $\mathbb R^p$ with the covariance matrix
\begin{equation}\label{eq: sigma definition}
\Sigma_n = \frac{1}{n}\sum_{i=1}^n{\mathrm{E}}[(X_i-\mu)(X_i-\mu)'],
\end{equation}
which coincides with the covariance matrix of $\sqrt{n}(\bar{X}_n-\mu)$.
Second, define the bootstrap critical value $c_{1-\alpha}^B$ as the $(1-\alpha)$th quantile of the conditional distribution of
\begin{equation}\label{eq: bootstrap statistic}
T_n^* = \max_{1\leq j\leq p}\frac{1}{\sqrt n}\sum_{i=1}^n(X_{i j}^* + a_j)
\end{equation}
given the data $X_1,\dots,X_n$, where $X_1^*,\dots,X_n^*$ is a (not necessarily empirical) bootstrap sample. We consider the following types of the bootstrap:
\begin{itemize}
\item  Empirical bootstrap: let $X_1^*,\dots,X_n^*$ be a sequence of i.i.d. random variables sampled from the uniform distribution on $\{X_1 - \bar X_n,\dots,X_n - \bar X_n\}$, where $\bar X_n = n^{-1}\sum_{i=1}^n X_i$ denotes the sample mean of the data $X_1,\dots,X_n$.
\item Multiplier bootstrap: let $e_1,\dots,e_n$ be a sequence of i.i.d. random variables with mean zero and unit variance, referred to as weights, which are independent of $X_1,\dots,X_n$. Define $X_i^* = e_i(X_i - \bar X_n)$ for all $i = 1,\dots,n$.
\end{itemize}


For the multiplier bootstrap, we will assume throughout the paper that the weights $e_1,\dots,e_n$  are such that
\begin{equation}\label{eq: multiplier bootstrap simplification}
\begin{array}{l}
\text{$e_i=e_{i,1}+e_{i,2}$, where $e_{i,1}$ and $e_{i,2}$ are independent, $e_{i,1}$ has}\\
\text{the $N(0,\sigma_e^2)$ distribution for some $\sigma_e\geq0$, and $|e_{i,2}|\leq3$.}
\end{array}
\end{equation}
Condition (\ref{eq: multiplier bootstrap simplification}) is mild and covers many commonly used weights, such as:
\begin{itemize}
\item Gaussian weights: $e_{i,1} \sim N(0,1)$ and $e_{i,2} = 0$.
\item Rademacher weights: $e_{i,1} = 0$ (i.e., $\sigma_{e} =0$) and ${\mathrm{P}} (e_{i,2} = \pm 1) = 1/2$.
\item Mammen's weights \cite{M93}: $e_{i,1} = 0$  and
\[
{\mathrm{P}}\left (e_{i,2} =  \frac{1 \pm \sqrt{5}}{2} \right ) =\frac{\sqrt{5} \mp 1}{2\sqrt{5}}.
\]
\end{itemize}
See Remark \ref{rem: weight discussion} for further discussion on Condition (\ref{eq: multiplier bootstrap simplification}).

Occasionally, we will also consider the weights with unit third moment, namely,
\begin{equation}\label{eq: third order matching multipliers}
{\mathrm{E}}[e_i^3] = 1,\quad\text{for all }i=1,\dots,n.
\end{equation}
The weights satisfying Condition (\ref{eq: third order matching multipliers}) correspond to  the third-order matching multiplier bootstrap mentioned in the Introduction. We note that Mammen's weights satisfy both Conditions (\ref{eq: multiplier bootstrap simplification}) and (\ref{eq: third order matching multipliers}), but neither Rademacher nor Gaussian weights satisfy Condition \eqref{eq: third order matching multipliers}.
See Lemma \ref{lem: multipliers} in the Supplemental Material, where we provide a more general class of distributions for the weights satisfying both  Conditions \eqref{eq: multiplier bootstrap simplification} and \eqref{eq: third order matching multipliers}.







Before proceeding to the regularity conditions, we also note that the multiplier bootstrap critical value $c_{1-\alpha}^B$ with Gaussian weights can be regarded as a feasible version of the Gaussian critical value $c_{1-\alpha}^G$. Indeed, it is easy to see that the former can be alternatively defined as the $(1-\alpha)$th quantile of the distribution of
$$
T_n^{\hat G} = \max_{1\leq j\leq p}(\hat G_j + a_j),$$
where
$
\hat G\sim N(0_p,\widehat\Sigma_n)$ and $\widehat \Sigma_{n}$ is the empirical covariance matrix
\begin{equation}\label{eq: sigma estimator definition}
\widehat\Sigma_n = n^{-1}\sum_{i=1}^n(X_i-\bar X_n)(X_i-\bar X_n)'.
\end{equation}
For brevity, we sometimes refer to both quantities as the Gaussian critical values.

\begin{remark}[On Condition (\ref{eq: multiplier bootstrap simplification})]
\label{rem: weight discussion}
Condition (\ref{eq: multiplier bootstrap simplification}) is technical and can be weakened depending on the moment conditions on $X_i$. A key step in the proof of Theorem \ref{cor: rejection probabilities} is to apply Theorem \ref{cor: max} ahead to approximate the conditional distribution of $T_n^*$ with that of the multiplier bootstrap statistic with weights following a Beta distribution that matches the moments of $e_i$ up to the third order (to be precise, we first replace the Gaussian components $e_{i,1}$ by bounded weights in the proof of Theorem \ref{cor: rejection probabilities}). Condition (\ref{eq: multiplier bootstrap simplification}) will be used to verify Conditions V, P, and B when we apply Theorem \ref{cor: max} there.
If, e.g., $X_i$ are bounded by $B_n$, then the conclusion of Theorem \ref{cor: rejection probabilities} continues to hold for sub-exponential weights. Since current Condition (\ref{eq: multiplier bootstrap simplification})  already covers many commonly used bootstrap weights, however, we do not pursue this generality of the weights to keep our presentation reasonable concise.
\end{remark}


\subsection{Regularity Conditions}
First, observe that given the construction of the statistic $T_n$ in \eqref{eq: statistic second} and its Gaussian and bootstrap analogs in \eqref{eq: gaussian analog} and \eqref{eq: bootstrap statistic}, it is without loss of generality to assume that $\mu_{j} = 0$ for all $j = 1,\dots,p$, which is what we do for the rest of the paper. Also, all our results follow immediately if $n= 2$, so we assume $n\geq 3$,  which  in particular implies $\log(p n)\geq 1$. In addition, since we are primarily interested in the case with large $p$, we assume $p\geq 2$.

Second, let $b_1$ and $b_2$ be some strictly positive constants such that $b_1\leq b_2$ and let $\{B_n\}_{n\geq 1}$ be a sequence of constants such that $B_n \geq 1$ for all $n\geq 1$. Here, the sequence $\{B_n\}_{n\geq 1}$ can diverge to infinity as the sample size $n$ increases.


\medskip

\noindent
{\bf Condition E:} {\em For all $i=1,\dots,n$ and $j = 1,\dots,p$, we have
\[
{\mathrm{E}}[\exp(|X_{i j}|/B_n)]\leq 2.
\]
}


\noindent
{\bf Condition M:} {\em For all $j=1,\dots,p$, we have
$$
b_1^2\leq \frac{1}{n}\sum_{i=1}^n {\mathrm{E}}[X_{i j}^2]\quad\text{and}\quad \frac{1}{n}\sum_{i=1}^n {\mathrm{E}}[X_{i j}^4]\leq B_n^2 b_2^2.
$$
}

\noindent
{\bf Condition S:} {\em For all $i=1,\dots,n$, the distribution of $X_i$ is symmetric in the sense that $X_i$ and $-X_i$ are identically distributed.}

\medskip


Condition E implies that the random variables $X_{i j}$ are sub-exponential with the Orlicz $\psi_1$-norm bounded by $B_n$; see \cite{V11} for details. The same sub-exponential condition was assumed in e.g. \cite{CCK17} and \cite{DZ17}; see Condition (E.1) in \cite{CCK17} and (E.1) in \cite{DZ17}.
The first part of Condition M, which we refer to as the variance lower bound condition,  requires that each component of the random vectors $X_i$ is scaled properly. The variance lower bound condition is needed to apply the anti-concentration inequalities (cf.  Lemmas \ref{lem: anticoncetration} and \ref{lem: rademacher anticoncentration} in the Supplemental Material) but can be dropped in Theorem \ref{thm: infinitesimal factors} ahead. Also, at least for Theorems \ref{thm: gaussian approximation main} and \ref{cor: rejection probabilities}, it can be relaxed by using Theorem 10 in \cite{DZ17}. However, to consistently state all the results, we work with the present assumption. Given the first part, the second part of Condition M holds if, for example, all random variables $X_{i j}$ are bounded  by $B_n$ and $n^{-1}\sum_{i=1}^n{\mathrm{E}}[X_{ij}^2]\leq b_2^2$ for all $j=1,\dots,p$.
Condition S means that the distribution of each $X_i$ is symmetric around the mean.
Importantly, none of these conditions restrict the correlation matrices of $X_i$, and so our results do not follow from the classical results in  empirical process theory.

In what follows, we will always  maintain Conditions E and M and will assume Condition S only in Theorem \ref{thm: gauss and rademacher}, which shows that imposing the symmetric distributions improves the approximation bound for the multiplier bootstrap with Rademacher weights.


\subsection{Main Results} We first present a non-asymptotic bound on the error of the Gaussian approximation to the distribution of the statistic $T_n$:
\begin{theorem}[Gaussian Approximation]\label{thm: gaussian approximation main}
Suppose that Conditions E and M are satisfied. Then
\begin{equation}\label{eq: gaussian bound main}
\left|{\mathrm{P}}\left(T_n > c^G_{1-\alpha}\right) - \alpha\right| \leq C\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4},
\end{equation}
where $C$ is a constant depending only on $b_1$ and $b_2$.
\end{theorem}
This result improves upon the bound in \cite{K19b}, who obtained a similar result with the rate $1/6$ instead of $1/4$. Since $a\in\mathbb R^p$ in the definition of $T_n$ in \eqref{eq: statistic second} is arbitrary, the bound \eqref{eq: gaussian bound main} can be equivalently stated as
$$
\sup_{A\in\mathcal A}\left|{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n X_i \in A\right) - {\mathrm{P}}(G\in A)\right|\leq C\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4},
$$
where $G\sim N(0_p,\Sigma_n)$ and $\mathcal A$ is the class of all hyper-rectangles in $\mathbb R^p$, i.e. sets of the form
$$
A = \Big\{ w = (w_1,\dots,w_p)'\in\mathbb R^p\colon a_{lj} \leq w_j \leq a_{rj}\text{ for all }j=1,\dots,p \Big\},
$$
for some constants $-\infty\leq a_{lj} \leq a_{rj} \leq \infty$ with $j=1,\dots,p$. This gives a quantitative Central Limit Theorem (CLT) over the hyper-rectangles in high dimensions.


The proof of Theorem \ref{thm: gaussian approximation main}, which is deferred to Section \ref{sec: proofs}, is fairly complicated and goes somewhat backward: (i) we first compare the conditional distribution of a third-order matching bootstrap statistic $T_n^*$ with that of the Gaussian multiplier bootstrap statistic $T_n^{\hat{G}}$, and then compare the conditional distribution of $T_n^{\hat{G}}$ with the distribution of $T_n^G$. These two comparisons rely on the Gaussian approximation via Stein kernel (Theorem \ref{thm: stein kernel}). Then, (ii) we use the preceding comparison between $T_n^*$ and $T_n^{G}$ to verify the anti-concentration for $T_n^*$ to invoke Theorem \ref{cor: max} and compare the conditional distribution of $T_n^*$ with the distribution of $T_n$. The proof of Theorem \ref{cor: max} relies on a novel technique which we call the \textit{iterative randomized Lindeberg method}.
The conclusion of Theorem \ref{thm: gaussian approximation main} follows from combining the results in Steps (i) and (ii) and the triangle inequality.


Comparison of the the Gaussian multiplier bootstrap statistic $T_n^{\hat{G}}$ with $T_n^G$ relies on  the following Gaussian-to-Gaussian comparison inequality, which can be of independent interest and whose proof is presented in Section \ref{sec: stein kernels} as a consequence of Theorem \ref{thm: stein kernel}:


\begin{proposition}[Gaussian-to-Gaussian Comparison]\label{coro:g-g-comparison}
If $Z_1$ and $Z_2$ are centered Gaussian random vectors in $\mathbb R^p$ with covariance matrices $\Sigma^1$ and $\Sigma^2$, respectively, and $\Sigma^2$ is such that $\Sigma^2_{jj}\geq c$ for all $j=1,\dots,p$ for some constant $c>0$, then
$$
\sup_{y\in\mathbb R^p}\Big|{\mathrm{P}}(Z_1\leq y) - {\mathrm{P}}(Z_2\leq y)\Big| \leq C\Big(\Delta \log^2 p\Big)^{1/2},
$$
where $C$ is a constant depending only on $c$ and $\Delta = \max_{1\leq j,k\leq p}|\Sigma^1_{jk} - \Sigma^2_{jk}|$.
\end{proposition}
\begin{remark}
Two comments on Proposition \ref{coro:g-g-comparison} are warranted. First, Proposition \ref{coro:g-g-comparison} improves upon Theorem 2 in \cite{CCK15}, which shows that
$$
\sup_{x\in\mathbb R}\left|{\mathrm{P}}\left(\max_{1\leq j\leq p}Z_{1j}\leq x\right) - {\mathrm{P}}\left(\max_{1\leq j\leq p}Z_{2j}\leq x\right)\right| \leq C\Big(\Delta \log^2 p\Big)^{1/3},
$$
under the same conditions. Second, the bound in this proposition is sharp in the sense that there exists a constant $c>0$ such that for infinitely many values of $p$, there exist centered Gaussian random vectors $Z_1$ and $Z_2$ in $\mathbb R^p$ such that the covariance matrix $\Sigma^2$ of $Z_2$ satisfies $\Sigma_{jj}^2 = 1$ for all $j=1,\dots,p$ and
$$
\sup_{y\in\mathbb R^p}\Big|{\mathrm{P}}(Z_1\leq y) - {\mathrm{P}}(Z_2\leq y)\Big| \geq c\Big(\Delta \log^2 p\Big)^{1/2}.
$$
The latter claim is proven in Appendix \ref{sec: sharpness} of the Supplemental Material.\hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}


Comparison of the conditional distribution of the third-order matching bootstrap statistic $T_n^*$ with that of $T_n$ (Theorem \ref{cor: max}) relies on the iterative randomized Lindeberg method.
An intuition behind the iterative randomized Lindeberg method goes as follows. Recall that, for any smooth function $g\colon\mathbb R^p\to\mathbb R$ and any two sequences of independent random vectors $X_1,\dots,X_n$ and $Y_1,\dots,Y_n$ in $\mathbb R^p$, in order to approximate ${\mathrm{E}}[g(X_1+\dots+X_n)]$ by ${\mathrm{E}}[g(Y_1+\dots+Y_n)]$, the original Lindeberg method constructs an interpolation path from ${\mathrm{E}}[g(X_1+\dots+X_n)]$ to ${\mathrm{E}}[g(Y_1+\dots+Y_n)]$ by replacing $X_i$'s with $Y_i$'s one-by-one in a given order and uses Taylor's expansion to show that the change in the expectation at each step is sufficiently small; see \cite{C06} for example. The randomized Lindeberg method, introduced in \cite{DZ17}, is similar to the original Lindeberg method but it replaces $X_i$'s with $Y_i$'s in a randomly selected order. It turns out that this randomization may bring substantial benefits to the final bound. In turn, to improve upon this version of the randomized Lindeberg method, we carry out a careful analysis of the coefficients in the Taylor's expansions underlying the method. In particular, given that $k$th order coefficients take the form of ${\mathrm{E}}[g^{(k)}(Z_1+\dots+Z_n)]$, up to some approximation error, where $g^{(k)}$ is a vector of the $k$th partial derivatives of $g$ and $Z_1,\dots,Z_n$ is a sequence such that some of its elements are given by $X_i$'s and others by $Y_i$, and using the fact that it is easier in our setting to bound ${\mathrm{E}}[g^{(k)}(Y_1+\dots+Y_n)]$, we apply the randomized Lindeberg method once again to approximate ${\mathrm{E}}[g^{(k)}(Z_1+\dots+Z_n)]$ by ${\mathrm{E}}[g^{(k)}(Y_1+\dots+Y_n)]$. Here, since a new application of the method will bring new Taylor's coefficients, we apply the same method over and over again until the approximation error becomes sufficiently small. We demonstrate that this iterative use of the randomized Lindeberg method gives further substantial benefits to the final bound. See also the discussion before Lemma \ref{lem: main} concerning comparisons of the iterative randomized Lindeberg method with the randomized Lindeberg method used in \cite{DZ17} and the related Slepian-Stein method used in our earlier work \cite{CCK13,CCK17}.

Our second main result gives a non-asymptotic bound on the deviation of the bootstrap rejection probabilities ${\mathrm{P}}(T_n > c_{1-\alpha}^B)$ from the nominal level $\alpha$ for the empirical and the multiplier bootstrap methods:
\begin{theorem}[Bootstrap Approximation]\label{cor: rejection probabilities}
Suppose that Conditions E and M are satisfied and that $c_{1-\alpha}^B$ is obtained via either the empirical bootstrap or the multiplier bootstrap with weights satisfying \eqref{eq: multiplier bootstrap simplification}. Then
\begin{equation}\label{eq: empirical bootstrap main}
\left|{\mathrm{P}}\left(T_n > c_{1-\alpha}^B\right) - \alpha\right| \leq C\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4},
\end{equation}
where $C$ is a constant depending only on $b_1$ and $b_2$.
\end{theorem}

This theorem improves upon the bounds in \cite{DZ17}, who obtained a similar result with the rate $1/6$ instead of $1/4$. In addition, we allow for a larger class of multiplier bootstrap methods. In particular, we do not require the weights $e_1,\dots,e_n$ to satisfy \eqref{eq: third order matching multipliers}. The proof of this theorem is given in Section \ref{sec: proofs}.







Our third main result gives a non-asymptotic bound on the deviation of the bootstrap rejection probabilities from the nominal level for the multiplier bootstrap method with Rademacher weights in the case of symmetric distributions:

\begin{theorem}[Rademacher Bootstrap Approximation in Symmetric Case]\label{thm: gauss and rademacher}
Suppose that Conditions E, M, and S are satisfied and that $c_{1-\alpha}^B$ is obtained via the multiplier bootstrap with Rademacher weights. Then
\begin{equation}\label{eq: rademacher symmetric}
\left|{\mathrm{P}}\left(T_n > c_{1-\alpha}^B\right) - \alpha\right| \leq C\left(\frac{B_n^2\log^3(p n)}{n}\right)^{1/2},
\end{equation}
where $C$ is a constant depending only on $b_1$ and $b_2$.
\end{theorem}
This theorem implies that the multiplier bootstrap with Rademacher weights is very accurate in the symmetric case. To prove it, we note that under the assumption of symmetric distributions, one can construct the randomization critical value $c_{1-\alpha}^R$ such that ${\mathrm{P}}(T_n > c_{1-\alpha}^R) = \alpha$, up to possible mass points in the distribution of $T_n$. Thus, given that the critical value based on the multiplier bootstrap with Rademacher weights turns out to be a feasible version of this randomization critical value and the two are close to each other, \eqref{eq: rademacher symmetric} follows if we can show that the distribution of $T_n$ is not too concentrated. To this end, we use an anti-concentration inequality for maxima of Rademacher processes derived in \cite{OST18}. The proof of Theorem \ref{thm: gauss and rademacher} is given in Appendix \ref{sec: proof of theorem 23} of the Supplemental Material.








Our fourth and final result shows that one-sided bounds in the bootstrap approximation can be substantially improved if we allow for incremental factors:

\begin{theorem}[Bootstrap Approximation with Incremental Factors]\label{thm: infinitesimal factors}
Suppose that Conditions E and M are satisfied and let $\eta > 0$ be a constant that may depend on $n$ and $p$. Then there exists a constant C depending only $b_1$ and $b_2$ such that the following hold.
\begin{enumerate}
\item[(i)] If $B_n^2\log^5(pn)\leq n$ and $c_{1-\alpha}^B$ is obtained via either the empirical bootstrap or the multiplier bootstrap with weights satisfying \eqref{eq: multiplier bootstrap simplification} and \eqref{eq: third order matching multipliers}, then we have
$$
{\mathrm{P}}(T_n > c_{1-\alpha}^B + \eta) \leq \alpha + C(1\vee\eta^{-4})\left(\frac{B_n^2\log^3(p n)}{n}\right)^{1/2}.
$$
\item[(ii)] If $n^{-1}\sum_{i=1}^n{\mathrm{E}}[X_{ij}^2]\leq b_2^2$ for all $j=1,\dots,p$ and $c_{1-\alpha}^B$ is obtained via the multiplier bootstrap with weights satisfying \eqref{eq: multiplier bootstrap simplification}, then
$$
{\mathrm{P}}(T_n > c_{1-\alpha}^B + \eta) \leq \alpha + C(1\vee\eta^{-4})\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/2}.
$$
\end{enumerate}
\end{theorem}

Theorem \ref{thm: infinitesimal factors} allows $\eta$ to (slowly) decrease with $n$ and/or $p$. For example, if we choose $\eta \sim (\log n)^{-1}$, then the over-rejection probability is of order $n^{-1/2}$ in $n$ up to log factors, while only requiring  $p$ to be $\log p = o\big(n^{1/3}/\text{polylog} (n)\big)$ in (i) and $\log p = o\big(n^{1/5}/\text{polylog}(n)\big)$ in (ii) provided that $B_n$ is bounded in $n$.

To prove this theorem, we use the randomized Lindeberg method but with an important simplification that the incremental factor $\eta$ now absorbs all the terms arising from smoothing the functions of the form $x\mapsto 1\{\max_{1\leq j\leq p}x_j > c\}$, which is used in the Lindeberg method. As discussed in the Introduction, Theorem \ref{thm: infinitesimal factors} is useful if one is concerned with the finite-sample over-rejection of tests based on the statistic $T_n$ as it says that adding an incremental factor $\eta$ to the critical value $c_{1-\alpha}^B$ may substantially reduce over-rejection, with a minimal effect on the power of the test. The proof of Theorem \ref{thm: infinitesimal factors} is given in Appendix \ref{sec: proof of theorem 24} of the Supplemental Material.


We conclude this section with a few remarks on cases with
 approximate sample means.
 In many applications (such as simultaneous inference for high-dimensional statistical models; cf. \cite{BCK15}), the statistic $T_n$ can only be asymptotically approximated by the maximum coordinate of the sample mean of independent random vectors. Also, those random vectors, often corresponding to the influence functions, may not be directly observable but have to be estimated. We emphasize here  that all our results can be extended to such approximate sample mean cases using the same arguments as those used in \cite{BCCHK18}; however, we have opted not to carry out the extension here for brevity of the paper.





\subsection{Gaussian and Bootstrap Approximations under Polynomial Moment Conditions}\label{sec: polynomial moment conditions}


So far we have assumed the sub-exponential condition (Condition E) for $X_i$. It turns out that combining some elements of the proof of Lemma \ref{lem: main} below and a truncation argument leads to analogs of the Gaussian and bootstrap approximation results  under polynomial moment conditions, which are given next. The proofs of Theorems \ref{cor: ga under polynomial conditions} and \ref{cor: boot app pol mom cond} can be found in Appendix \ref{sec: proof of poly moment} in the Supplementary Material.




\begin{theorem}[Gaussian Approximation under Polynomial Moment Conditions]\label{cor: ga under polynomial conditions}
Suppose that Condition M is satisfied and that for some $q>2$, we have
\begin{equation}
{\mathrm{E}}\left[\max_{1\leq j\leq p}|X_{ij}|^q\right]\leq B_n^q
\label{eq: poly}
\end{equation}
for all $i=1,\dots,n$. Then
$$
\left|{\mathrm{P}}\left(T_n > c^G_{1-\alpha}\right) - \alpha\right| \leq C\left\{\left(\frac{B_n^2\log^5 p}{n}\right)^{1/4}+\sqrt{\frac{B_n^2(\log p)^{3-2/q}}{n^{1-2/q}}}\right\},
$$
where $C$ is a constant depending only on, $q$, $b_1$, and $b_2$.
\end{theorem}
This theorem improves on the corresponding result obtained by \cite{K19b} by the fourth author. For bootstrap approximation, we focus on the Gaussian multiplier and empirical bootstraps for simplicity.



\begin{theorem}[Bootstrap Approximation under Polynomial Moment Conditions]\label{cor: boot app pol mom cond}
Suppose that Condition M is satisfied and that Condition (\ref{eq: poly}) holds for all $i=1,\dots,n$ for some $q>2$.  Let $c_{1-\alpha}^B$ be the critical value obtained via either the empirical bootstrap or the Gaussian multiplier bootstrap. Then
\begin{equation*}
\left|{\mathrm{P}}\left(T_n > c_{1-\alpha}^B\right) - \alpha\right| \leq C\left\{\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4}+\sqrt{\frac{B_n^2\log^{3-2/q}(p n)}{n^{1-2/q}}}\right\},
\end{equation*}
where $C$ is a constant depending only on $q$,  $b_1$, and $b_2$.
\end{theorem}

This theorem improves on the error bound for the empirical bootstrap given in \cite{DZ17} under the polynomial moment condition.


\section{Main Theoretical Arguments}

\subsection{Iterative Randomized Lindeberg Method}\label{sec: main arguments}
In this section, we derive a distributional approximation result, Theorem \ref{cor: max}, using a novel proof technique, which we call the iterative randomized Lindeberg method. We will use this result in Section \ref{sec: proofs} to prove our main results on the Gaussian and bootstrap approximations in high dimensions, as stated in Section \ref{sec: main results}.

Let $V_1,\dots,V_n, Z_1,\dots,Z_n$ be a sequence of independent random vectors in $\mathbb R^p$ such that ${\mathrm{E}}[V_{i j}] = {\mathrm{E}}[Z_{i j}] = 0$ for all $i = 1,\dots,n$ and $j = 1,\dots,p$, where $V_{i j}$ and $Z_{i j}$ denote the $j$th components of $V_{i j}$ and $Z_{i j}$, respectively. We will assume that these vectors obey the following conditions:

\medskip
\noindent
{\bf Condition V:} {\em There exists a constant $C_v>0$ such that for all $j = 1,\dots,p$, we have
$$
\frac{1}{n}\sum_{i=1}^n{\mathrm{E}}\left[V_{i j}^4 + Z_{i j}^4\right] \leq C_v B_n^2.
$$

}
\noindent
{\bf Condition P:} {\em There exists a constant $C_p\geq 1$ such that for all $i = 1,\dots,n$, we have
$$
{\mathrm{P}}\Big(\|V_i\|_{\infty}\vee\|Z_i\|_{\infty} > C_p B_n\log(p n)\Big) \leq 1/n^4.
$$

}
\medskip
\noindent
{\bf Condition B:} {\em There exists a constant $C_b>0$ such that for all $i = 1,\dots,n$, we have
$$
{\mathrm{E}}\Big[\|V_i\|_{\infty}^8 + \|Z_i\|_{\infty}^8\Big]\leq C_bB_n^8\log^8(p n).
$$

}
\medskip
\noindent
{\bf Condition A:} {\em
There exist constants $C_a > 0$ and $\delta\geq0$ such that for all $(y,t)\in\mathbb R^p\times (0,\infty)$, we have
$$
{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n Z_i \leq y + t \right) - {\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n Z_i \leq y\right) \leq C_a\left(t\sqrt{\log p}+\delta\right).
$$}

\noindent
Note that the constants $C_v$, $C_p$, $C_b$, and $C_a$ appearing in these conditions are not supposed to be dependent on their indices, e.g. $C_p$ here is not allowed to change with $p$; the indices are introduced with the only goal to differentiate between the constants.

The following is the main result of this section:

\begin{theorem}[Distributional Approximation via Iterative Randomized Lindeberg Method]\label{cor: max}
Suppose that Conditions V, P, B, and A are satisfied. In addition, suppose that
\begin{equation}\label{eq: bn bounds 1}
\max_{1\leq j,k\leq p}\left| \frac{1}{\sqrt n}\sum_{i=1}^n({\mathrm{E}}[V_{i j}V_{i k}] - {\mathrm{E}}[Z_{i j}Z_{i k}]) \right| \leq C_mB_n\sqrt{\log(p n)}
\end{equation}
and
\begin{equation}\label{eq: bn bounds 2}
\max_{1\leq j,k,l\leq p}\left| \frac{1}{\sqrt n}\sum_{i=1}^n ({\mathrm{E}}[V_{i j}V_{i k}V_{i l}] - {\mathrm{E}}[Z_{i j}Z_{i k}Z_{i l}]) \right| \leq C_m B_n^2\sqrt{\log^3(p n)}
\end{equation}
for some constant $C_m$. Then
\[
\sup_{y\in\mathbb R^p}\left| {\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n V_i \leq y\right) - {\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n Z_i \leq y\right) \right|
 \le C\left(\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4}+\delta\right),
\]
where $C$ is a constant depending only on $C_v$, $C_p$, $C_b$, $C_a$, and $C_m$.
\end{theorem}

\begin{remark}[On Sharpness of Theorem \ref{cor: max}]
We do not claim sharpness of Theorem \ref{cor: max} in the high-dimensional case $p\gg n$ (when $p$ is fixed, the theorem is not sharp in view of the classical Berry-Esseen bound). On one hand, classical Edgeworth expansions in the low-dimensional case suggest that conditions like \eqref{eq: bn bounds 2} should lead to better distributional approximation results than the corresponding Gaussian approximation results, which we do not observe in Theorem \ref{cor: max} since Theorem \ref{thm: gaussian approximation main} gives the same dependence on both $n$ and $p$ for the Gaussian approximation. On the other hand, to the best of our knowledge, there exist no analogs of  Edgeworth expansions in high dimensions. The question whether conditions like \eqref{eq: bn bounds 2} can be used to improve distributional approximations (relative to the Gaussian approximations) thus remains open.
\hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}


To prove this result, we will need additional notation. For all $\epsilon\in\{0,1\}^n$, define
\begin{equation}\label{eq: rho eps definition}
\varrho_{\epsilon} = \sup_{y\in\mathbb R^p}\left|{\mathrm{P}}\left(S_{n,\epsilon}^V \leq y\right) - {\mathrm{P}}\left(S_n^Z \leq y\right)\right|,
\end{equation}
where
$$
S_{n,\epsilon}^V = \frac{1}{\sqrt n}\sum_{i=1}^n (\epsilon_iV_i + (1 - \epsilon_i)Z_i)\quad\text{and}\quad S_n^Z = \frac{1}{\sqrt n}\sum_{i=1}^n Z_i.
$$
We will replace $\epsilon$ with a certain sequence of random vectors $\epsilon^0,\dots,\epsilon^D \in \{ 0,1 \}^n$,  independent of $V_1,\dots,V_n,Z_1,\dots,Z_n$, and derive  recursive bounds for $\rho_{\epsilon^d}$ for $d=0,\dots,D$, which lead to the desired bound in Theorem \ref{cor: max}. Such a sequence of random vectors $\epsilon^0,\dots,\epsilon^D \in \{ 0,1 \}^n$ is constructed as follow:
\begin{itemize}
\item Set $D = [4\log n]+1$ and initialize $\epsilon^0 = (1,\dots,1)$.
\item Let $U_1,\dots,U_D$ be a sequence of independent uniform $[0,1]$ random variables that are independent of $V_1,\dots,V_n,Z_1,\dots,Z_n$.
\item For $d=1,\dots,D$: conditionally on $\epsilon^{d-1}$ and $U_1,\dots,U_D$, set $\epsilon^d_i = 0$  if $\epsilon^{d-1}_i = 0$, and generate $\{\epsilon^d_i\}_{i\in I_{d-1}}$ with $I_{d-1} = \{ i=1,\dots, n : \epsilon^{d-1}_i = 1 \}$  as i.i.d.~Bernoulli($U_d$) random variables.
\end{itemize}
It is not difficult to see that for each $d=1,\dots,D$, the random vector $\epsilon^d$ satisfies the following properties:
\begin{enumerate}
\item[(i)]  for all $i = 1,\dots,n$, $\epsilon^d_i = 0$ if $\epsilon^{d-1}_i = 0$, and
\item[(ii)] for $I_{d-1} = \{i=1,\dots,n\colon \epsilon^{d-1}_i = 1\}$, the random variables $\{\epsilon^d_i\}_{i\in I_{d-1}}$ are exchangeable conditional on $\epsilon^{d-1}$ and satisfy
\begin{equation}\label{eq: recursive epsilon}
{\mathrm{P}}\left( \sum_{i\in I_{d-1}}\epsilon^d_i = s \mid \epsilon^{d-1} \right) = \frac{1}{|I_{d-1}| + 1},\quad\text{for all }s = 0,\dots,|I_{d-1}|.
\end{equation}
\end{enumerate}
Indeed, to see that (\ref{eq: recursive epsilon}) holds, observe that,  conditional on $\epsilon^{d-1}$ and $U_d$,  $\sum_{i \in I_{d-1}} \epsilon_i^d$ follows the binomial distribution with parameters $|I_{d-1}|$ and (success probability) $U_d$, so that
\[
\begin{split}
{\mathrm{P}}\left( \sum_{i\in I_{d-1}}\epsilon^d_i = s \mid \epsilon^{d-1} \right) &=
\binom{|I_{d-1}|}{s}  \int_0^1 u^{s} (1-u)^{|I_{d-1}|-s} du \\
&=\binom{|I_{d-1}|}{s} \frac{s! (|I_{d-1}|-s)!}{(|I_{d-1}|+1)!}  = \frac{1}{|I_{d-1}|+1}.
\end{split}
\]
Also, two properties (i) and (ii) ensure that $S_{n,\epsilon^d}^V-n^{-1/2}\sum_{i\notin I_{d-1}}Z_i$ is the randomized Lindeberg interpolant between $n^{-1/2}\sum_{i\in I_{d-1}}V_i$ and $n^{-1/2}\sum_{i\in I_{d-1}}Z_i$; see Lemma \ref{lem: rand lind} and the discussion at the beginning of Step 1 of the proof of Lemma \ref{lem: main}.



Further, for all $i=1,\dots,n$ and $j,k,l = 1,\dots,p$, define
$$
\mathcal E_{i, j k}^V = {\mathrm{E}}[V_{i j}V_{i k}], \ \mathcal E_{i, j k l}^V = {\mathrm{E}}[V_{i j} V_{i k}V_{i l}],
$$
$$
\mathcal E_{i, j k}^Z = {\mathrm{E}}[Z_{i j}Z_{i k}], \ \mathcal E_{i, j k l}^Z = {\mathrm{E}}[Z_{i j} Z_{i k}Z_{i l}].
$$
For all $n\geq 1$ and $d = 0,\dots,D$, let $\mathcal B_{n,1,d}$ and $\mathcal B_{n,2,d}$ be some strictly positive constants, and define the event $\mathcal A_d$ by
\[
\begin{split}
\mathcal A_d &= \left \{
\max_{1\leq j,k\leq p}\left| \frac{1}{\sqrt n}\sum_{i=1}^n\epsilon^d_i(\mathcal E_{i, j k}^V - \mathcal E_{i, jk}^Z) \right| \leq \mathcal B_{n,1,d} \right \} \\
&\quad \bigcap
\left \{
\max_{1\leq j,k,l\leq p}\left| \frac{1}{\sqrt n}\sum_{i=1}^n \epsilon_i^d (\mathcal E_{i,j k l}^V - \mathcal E_{i, j k l}^Z) \right| \leq \mathcal B_{n,2,d}
\right \}.
\end{split}
\]

The proof of Theorem \ref{cor: max} proceeds as follows. In Lemma \ref{lem: main} and Corollary \ref{cor: main}, we establish a recursive inequality for ${\mathrm{E}}[\varrho_{\epsilon^d}1\{\mathcal A_d\}]$, $d=0,\dots,D$. Next, we show in Lemma \ref{lem: closing} that ${\mathrm{E}}[\varrho_{\epsilon^D}1\{\mathcal A_D\}]$ is bounded by $1/n$. Then, we use an induction argument backward to derive a bound for ${\mathrm{E}}[\varrho_{\epsilon^0}1\{\mathcal A_0\}]$. Since $\epsilon^0_i=1$ for all $i$, this gives the claim of the theorem once we appropriately choose the constants $\mathcal B_{n,1,d}$ and $\mathcal B_{n,2,d}$. The proof of Lemma \ref{lem: main} is long and is given in Appendix \ref{sec: proof of main lemma} of the Supplemental Material.


The derivation of the recursive inequality is based on connecting $S_{n,\epsilon^d}^V$ with $S_n^Z$ by the randomized Lindeberg method originally developed by \cite{DZ17}.
A similar approach was used in \cite{CCK13,CCK17} to connect $S_{n,\epsilon^0}^V$ with $G$, where the Slepian--Stein method was applied instead. Unlike the latter approach, the randomized Lindeberg method allows us to match the moments of $S_{n,\epsilon^d}^V$ and $S_n^Z$ up to the third order rather than the second order. This leads to improvement on the power of $\log(pn)$ factors. In addition, we incorporate a smoothing effect induced by $Z_i$ via Condition A into our argument. This along with the higher-order moment matching lead to improvement on the power of the sample size $n$.


\begin{lemma}\label{lem: main}
Suppose that Conditions V, P, B, and A are satisfied. Then for any $d = 0,\dots,D-1$ and any constant $\phi>0$ such that
\begin{equation}\label{eq: phi restriction}
C_pB_n\phi \log^2(p n)\leq \sqrt n,
\end{equation}
we have on the event $\mathcal A_d$,
\begin{align*}
\varrho_{\epsilon^d} &\lesssim \frac{\sqrt{\log p}}{\phi} + \delta + \frac{B_n^2\phi^4\log^5(p n)}{n^2} +  \left( {\mathrm{E}}[\varrho_{\epsilon^{d+1}}\mid \epsilon^d] + \frac{\sqrt{\log p}}{\phi}+\delta \right)\\
&\quad \times \left( \frac{\mathcal B_{n,1,d}\phi^2\log p}{\sqrt n} + \frac{\mathcal B_{n,2,d}\phi^3\log^2p}{n} + \frac{B_n^2\phi^4\log^3(p n)}{n} \right)
\end{align*}
up to a constant depending only on $C_v$, $C_p$, $C_b$, and $C_a$.
\end{lemma}

\begin{remark}[Choice of $\phi$]
We will choose $\phi$ to depend on $n$ via $n^{1/4}$ when applying this lemma.\hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}



\begin{corollary}\label{cor: main}
Suppose that all assumptions of Lemma \ref{lem: main} are satisfied. Then there exists a constant $K>0$ depending only on $C_v$, $C_p$, and $C_b$ such that for all $d = 0,\dots,D-1$, if $\mathcal B_{n,1,d+1} \geq \mathcal B_{n,1,d} + KB_n\log^{1/2}(p n)$ and $\mathcal B_{n,2,d+1} \geq \mathcal B_{n,2,d} + KB_n^2\log^{3/2}(pn)$, then for any constant $\phi>0$ satisfying \eqref{eq: phi restriction}, we have
\begin{align}
{\mathrm{E}}[\varrho_{\epsilon^d}1\{\mathcal A_d\}] &\lesssim \frac{\sqrt{\log p}}{\phi} + \delta+ \frac{B_n^2\phi^{4}\log^5(p n)}{n^2} + \left( {\mathrm{E}}[\varrho_{\epsilon^{d+1}} 1\{\mathcal A_{d+1}\}] + \frac{\sqrt{\log p}}{\phi} +\delta\right)  \nonumber\\
&\quad\times \left( \frac{\mathcal B_{n,1,d}\phi^2\log p}{\sqrt n} + \frac{\mathcal B_{n,2,d}\phi^3\log^2p}{n} + \frac{B_n^2\phi^4\log^3(p n)}{n} \right)\label{eq: max induction}
\end{align}
up to a constant depending only on $C_v$, $C_p$, $C_b$, and $C_a$.
\end{corollary}

\begin{proof}
Since we assume throughout the paper that $p\geq 2$, the conclusion is trivial if $\phi<1$. We will therefore assume in the proof that $\phi\geq 1$. In turn, $\phi\geq 1$ together with \eqref{eq: phi restriction} imply that
\begin{equation}\label{eq: lazy condition 1}
C_p B_n \log^2(p n)\leq \sqrt n.
\end{equation}
This condition will be useful in the proof.


Fix $d = 0,\dots,D-1$. Then, given that $\mathcal A_d$ depends only on $\epsilon^d$, we have by Lemma \ref{lem: main} that
\begin{align*}
{\mathrm{E}}[\varrho_{\epsilon^d}1\{\mathcal A_d\}] &\lesssim \frac{\sqrt{\log p}}{\phi} + \delta + \frac{B_n^2\phi^{4}\log^5(pn)}{n^2} +  \left( {\mathrm{E}}[\varrho_{\epsilon^{d+1}} 1\{\mathcal A_{d}\}] + \frac{\sqrt{\log p}}{\phi} + \delta\right)\\
&\quad \times \left( \frac{\mathcal B_{n,1,d}\phi^2\log p}{\sqrt n} + \frac{\mathcal B_{n,2,d}\phi^3\log^2p}{n} + \frac{B_n^2\phi^4\log^3(p n)}{n} \right)
\end{align*}
up to a constant depending only on $C_v$, $C_p$, $C_b$, and $C_a$. Thus, given that \eqref{eq: phi restriction} implies that $\sqrt{\log p}/\phi\geq 1/n$, the conclusion of the corollary will follow if we can show that
\begin{equation}\label{eq: event switch}
{\mathrm{E}}[\varrho_{\epsilon^{d+1}}1\{\mathcal A_d\}] \leq {\mathrm{E}}[\varrho_{\epsilon^{d+1}}1\{\mathcal A_{d+1}\}] + 4/n.
\end{equation}
To this end, we first observe that , as $\varrho_{\epsilon^{d+1}}\in[0,1]$,
\begin{equation}
\begin{split}
{\mathrm{E}}[\varrho_{\epsilon^{d+1}}1\{\mathcal A_d\}]
&= {\mathrm{E}}[\varrho_{\epsilon^{d+1}}1\{\mathcal A_d\}1\{\mathcal A_{d+1}\}] + {\mathrm{E}}[\varrho_{\epsilon^{d+1}}1\{\mathcal A_d\}(1 - 1\{\mathcal A_{d+1}\})]\\
& \leq {\mathrm{E}}[\varrho_{\epsilon^{d+1}}1\{\mathcal A_{d+1}\}]
+ {\mathrm{E}}[1\{\mathcal A_d\}(1 - 1\{\mathcal A_{d+1}\})]\\
&={\mathrm{E}}[\varrho_{\epsilon^{d+1}}1\{\mathcal A_{d+1}\}] + \underbrace{{\mathrm{P}}(\mathcal A_d) - {\mathrm{P}}(\mathcal A_d \cap \mathcal A_{d+1})}_{={\mathrm{P}}(\mathcal A_d) \left ( 1-{\mathrm{P}}(\mathcal A_{d+1}\mid \mathcal A_d) \right)} \\
& \leq {\mathrm{E}}[\varrho_{\epsilon^{d+1}}1\{\mathcal A_{d+1}\}]
+ 1 - {\mathrm{P}}(\mathcal A_{d+1}\mid \mathcal A_d).
\end{split}
\label{eq: event switch proof}
\end{equation}
Moreover, by Lemma \ref{lem: exponential inequality} in the Supplemental Material, for all $j,k = 1,\dots,p$ and $t > 0$, we have
\begin{align*}
&{\mathrm{P}}\left(\left| \frac{1}{\sqrt n}\sum_{i=1}^n\epsilon_i^{d+1}(\mathcal E_{i,jk}^V - \mathcal E_{i,jk}^Z) \right| > \left| \frac{1}{\sqrt n}\sum_{i=1}^n\epsilon_i^{d}(\mathcal E_{i,jk}^V - \mathcal E_{i,jk}^Z) \right| + t \mid \epsilon^d \right)\\
&\quad \leq 2\exp\left(-\frac{n t^2}{32\sum_{i=1}^n(\mathcal E_{i,jk}^V - \mathcal E_{i,jk}^Z)^2}\right) \leq 2\exp\left(-\frac{t^2}{128 C_v B_n^2}\right),
\end{align*}
where the second inequality follows from Condition V. Applying this inequality with $t = 8B_n\sqrt{6C_v\log(p n)}$ and using the fact that
$$
\max_{1\leq j,k\leq p}\left|\frac{1}{\sqrt n}\sum_{i=1}^n \epsilon_i^d(\mathcal E_{i,jk}^V - \mathcal E_{i,jk}^Z)\right| \leq \mathcal B_{n,1,d} \quad \text{on $\mathcal A_d$},
$$
 we have by the union bound that for any $\mathcal B_{n,1,d+1} \geq \mathcal B_{n,1,d} + t$,
$$
{\mathrm{P}}\left(\max_{1\leq j,k\leq p}\left|\frac{1}{\sqrt n}\sum_{i=1}^n \epsilon_i^{d+1}(\mathcal E_{i,jk}^V - \mathcal E_{i,jk}^Z)\right| > \mathcal B_{n,1,d+1}\mid \mathcal A_d\right) \leq \frac{2p^2}{(pn)^3}\leq \frac{2}{n}.
$$
In addition, for all $i=1,\dots,n$ and $j,k,l=1,\dots,p$, setting $\tilde V_i = 1\{\|V_i\|_{\infty}\leq C_p B_n\log(pn)\}$, we have that
\begin{equation}
\begin{split}
|\mathcal E_{i,jkl}^V| &\leq {\mathrm{E}}[|V_{ij}V_{ik}V_{il}|] = {\mathrm{E}}\Big[\tilde V_i |V_{ij}V_{ik}V_{il}|\Big] + {\mathrm{E}}\Big[(1 - \tilde V_i) |V_{ij}V_{ik}V_{il}|\Big]\\\
&\leq C_pB_n\log(p n){\mathrm{E}}[|V_{ij}V_{ik}|] + ({\mathrm{E}}[1-\tilde V_i])^{1/2}({\mathrm{E}}[\|V_i\|_{\infty}^6])^{1/2}\\
&\leq C_pB_n\log(p n){\mathrm{E}}[|V_{ij}V_{ik}|] + C_b^{3/8}B_n^3\log^{3}(p n)/n^2
\end{split}
\label{eq: details moment calculations}
\end{equation}
and similarly
$$
|\mathcal E_{i,jkl}^Z| \leq C_pB_n\log(p n){\mathrm{E}}[|Z_{ij}Z_{ik}|] + C_b^{3/8}B_n^3\log^{3}(p n)/n^2
$$
by Conditions P and B. Hence, by Condition V and \eqref{eq: lazy condition 1}, there exists a constant $C$ depending only on $C_v$, $C_p$, and $C_b$ such that
$$
\frac{32}{n}\sum_{i=1}^n(\mathcal E_{i,jkl}^V - \mathcal E_{i,jkl}^Z)^2 \leq C B_n^4\log^2(p n).
$$
Thus, by the same argument as above, for all $j,k,l=1,\dots,p$ and $t>0$,
\begin{align*}
&{\mathrm{P}}\left(\left| \frac{1}{\sqrt n}\sum_{i=1}^n\epsilon_i^{d+1}(\mathcal E_{i,jkl}^V - \mathcal E_{i,jkl}^Z) \right| > \left| \frac{1}{\sqrt n}\sum_{i=1}^n\epsilon_i^{d}(\mathcal E_{i,jkl}^V - \mathcal E_{i,jkl}^Z) \right| + t \mid \epsilon^d \right)\\
&\quad \leq 2\exp\left(-\frac{n t^2}{32\sum_{i=1}^n(\mathcal E_{i,jkl}^V - \mathcal E_{i,jkl}^Z)^2}\right) \leq 2\exp\left(-\frac{t^2}{C B_n^4\log^2(p n)}\right).
\end{align*}
Applying this inequality with $t = \sqrt{3C}B_n^2\log^{3/2}(p n)$ shows that for any $\mathcal B_{n,2,d+1} \geq \mathcal B_{n,2,d} + t$, we have
$$
{\mathrm{P}}\left(\max_{1\leq j,k,l\leq p}\left|\frac{1}{\sqrt n}\sum_{i=1}^n \epsilon_i^{d+1}(\mathcal E_{i,jkl}^V - \mathcal E_{i,jkl}^Z)\right| > \mathcal B_{n,2,d+1} \mid \mathcal A_d\right) \leq \frac{2p^3}{(p n)^3}\leq \frac{2}{n}.
$$
Thus,
$
1 - {\mathrm{P}}(\mathcal A_{d+1}\mid \mathcal A_d) \leq 4/n,
$
which in combination with \eqref{eq: event switch proof} implies \eqref{eq: event switch} and completes the proof.
\end{proof}


\begin{lemma}\label{lem: closing}
For any constant $\phi > 0$ such that \eqref{eq: phi restriction} holds, we have ${\mathrm{E}}[\varrho_{\epsilon^D}1\{\mathcal A_D\}] \leq 1/n$.
\end{lemma}
\begin{proof}
Recall that $D = [4\log n] + 1$ and note that $\varrho_{\epsilon^D} = 0$ if $\epsilon^D = (0,\dots,0)'$. Moreover, by Markov's inequality,
\begin{align*}
&{\mathrm{P}}(\epsilon^D \neq (0,\dots,0)') =  {\mathrm{P}}\left(\sum_{i=1}^n\epsilon^D_i \geq 1\right) \leq {\mathrm{E}}\left[\sum_{i=1}^n\epsilon^D_i\right] = {\mathrm{E}}\left[{\mathrm{E}}\left[\sum_{i=1}^n\epsilon^D_i \mid \sum_{i=1}^n\epsilon^{D-1}_i\right]\right] \\
&\quad    = {\mathrm{E}}\left[\frac{1}{2}\sum_{i=1}^n\epsilon^{D-1}_i\right]
 = \dots = {\mathrm{E}}\left[\frac{1}{2^D}\sum_{i=1}^n\epsilon^0_i\right] = \frac{n}{2^D} \leq \frac{n}{2^{4\log n}} \leq \frac{1}{n},
\end{align*}
where the equalities on the second line follow from \eqref{eq: recursive epsilon}. Hence,
$$
{\mathrm{E}}[\varrho_{\epsilon^D}1\{\mathcal A_D\}] \leq {\mathrm{E}}[\varrho_{\epsilon^D}] \leq {\mathrm{P}}(\epsilon^D \neq (0,\dots,0)') \leq 1/ n,
$$
as desired.
\end{proof}
\begin{proof}[Proof of Theorem \ref{cor: max}]
Throughout the proof, we will assume that
\begin{equation}\label{eq: natural bound}
C_p^4B_n^2\log^5(p n)\leq n
\end{equation}
since otherwise the conclusion  of the theorem is trivial.

Let $K$ be the constant from Corollary \ref{cor: main} and for all $d = 0,\dots,D$, define
\begin{equation}\label{eq: fancy b bounds}
\mathcal B_{n,1,d}= C_1(d+1)B_n\log^{1/2}(p n) \quad \text{and} \quad \mathcal B_{n,2,d}=C_1(d+1)B_n^2\log^{3/2}(p n),
\end{equation}
where $C_1 = C_m + K$, so that $\mathcal A_0$ holds by \eqref{eq: bn bounds 1} and \eqref{eq: bn bounds 2} and, in addition, the requirements of Corollary \ref{cor: main} on $\mathcal B_{n,1,d}$ and $\mathcal B_{n,2,d}$ also hold.

Now, for all $d = 0,\dots,D$, define
$$
f_d = \inf\left\{x\geq 1\colon {\mathrm{E}}[\varrho_{\epsilon^d}1\{\mathcal A_d\}] \leq  x\left(\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4}+\delta\right)\right\}.
$$
Note that $f_d<\infty$ because $\varrho_{\epsilon^d}\leq1$. Then, for all $d = 0,\dots,D-1$, apply Corollary \ref{cor: main} with
$$
\phi = \phi_d = \frac{n^{1/4}}{B_n^{1/2}\log^{3/4}(p n)((d+1)f_{d+1})^{1/3}},
$$
which satisfies the required condition \eqref{eq: phi restriction} since we assume \eqref{eq: natural bound}. Since
\begin{align*}
\frac{B_n^2\phi_d^4\log^{5}(pn)}{n^2}
& \leq\frac{\log^{2}(p n)}{n}
\leq\frac{\log^{1/4}(p n)}{n^{1/4}}
\leq \frac{C_p B_n^{1/2} \log^{1/4}(p n)}{n^{1/4}}\\
& \leq\frac{C_p\sqrt{\log p}}{\phi_d} \leq C_p((d+1)f_{d+1})^{1/3}\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4}, \\
\frac{\mathcal B_{n,1,d}\phi_d^2\log p}{\sqrt n} &\leq \frac{C_1(d+1)}{((d+1)f_{d+1})^{2/3}}, \quad \text{and} \\
\frac{\mathcal B_{n,2,d}\phi^3_d\log^2p}{n} &\bigvee \frac{B_n^2\phi_d^4\log^3(p n)}{n} \leq \frac{C_1\vee 1}{f_{d+1}},
\end{align*}
we have by  Corollary \ref{cor: main}
$$
{\mathrm{E}}[\rho_{\epsilon^d}1\{\mathcal A_d\}] \leq C_2\Big(f_{d+1}^{2/3} + (d+1)^{2/3} + 1\Big)\left(\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4}+\delta\right)
$$
for some constant $C_2\geq 1$ depending only on $C_v$, $C_p$, $C_b$, $C_a$, and $C_m$. Hence,
$$
f_d \leq C_2\Big(f_{d+1}^{2/3} + (d+1)^{2/3} + 1\Big),\quad\text{for all }d = 0,\dots,D-1.
$$
Here, we have $f_D = 1$ by Lemma \ref{lem: closing} since $B_n\geq 1$ by assumption. Therefore, by a simple induction argument, we conclude that there exists a constant $C\geq 1$ depending only on $C_2$ such that
$$
f_d \leq C(d+1),\quad\text{for all }d = 0,\dots,D.
$$
In particular, it follows that
$$
\varrho_{\epsilon^0}1\{\mathcal A_0\} = {\mathrm{E}}[\varrho_{\epsilon^0}1\{\mathcal A_0\}] \leq
C\left(\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4}+\delta\right).
$$
Since $\mathcal A_0$ holds by construction, so that $1\{\mathcal A_0\} = 1$, the desired bound follows by combining this inequality and the definition of $\varrho_{\epsilon^0}$.
\end{proof}


\subsection{Stein Kernels and Gaussian Approximation}\label{sec: stein kernels}
Let $C_b^2(\mathbb R^p)$ be the class of twice continuously differentiable functions $\varphi$ on $\mathbb R^p$ such that $\varphi$ and all its partial derivatives up to the second order are bounded where $p\geq 2$.
Let $V$ be a centered random vector in $\mathbb R^p$ and assume that there exists a measurable function $\tau:\mathbb R^p\to\mathbb R^{p\times p}$ such that
$$
\sum_{j=1}^p{\mathrm{E}}[\partial_j \varphi(V)V_j] = \sum_{j,k=1}^p{\mathrm{E}}[\partial_{jk}\varphi(V)\tau_{jk}(V)]
$$
for all $\varphi\in C_b^2(\mathbb R^p)$. This function $\tau$ is called a \textit{Stein kernel} for the random vector $V$.  Also, let $Z$ be a centered Gaussian random vector in $\mathbb R^p$ with covariance matrix $\Sigma$.
\begin{theorem}[Gaussian Approximation via Stein Kernels]\label{thm: stein kernel}
If $\Sigma_{jj}\geq c$ for all $j=1,\dots,p$ and some constant $c>0$, then
$$
\sup_{y\in\mathbb R^p}\Big|{\mathrm{P}}(V\leq y) - {\mathrm{P}}(Z\leq y)\Big| \leq C\Big(\Delta \log^2 p\Big)^{1/2},
$$
where $C$ is a constant depending only on $c$ and
$
\Delta = {\mathrm{E}}\left[\max_{1\leq j,k\leq p}|\tau_{jk}(V) - \Sigma_{jk}|\right].
$
\end{theorem}
\begin{remark}
This theorem improves upon Proposition 4.1 in \cite{K19a}, which shows that
$$
\sup_{y\in\mathbb R^p}\Big|{\mathrm{P}}(V\leq y) - {\mathrm{P}}(Z\leq y)\Big| \leq C\Big(\Delta \log^2 p\Big)^{1/3}
$$
under the same conditions.\hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}

Theorem \ref{thm: stein kernel} is proven in Appendix \ref{sec: proof of kernel result} of the Supplemental Material. It has two important corollaries.
The first  is Proposition \ref{coro:g-g-comparison}, a sharp Gaussian-to-Gaussian comparison inequality stated in Section \ref{sec: main results}:

\begin{proof}[Proof of Proposition \ref{coro:g-g-comparison}]
If $V$ is a centered Gaussian random vector, then by the multivariate Stein identity, its Stein kernel coincides with its covariance matrix. Hence, Theorem \ref{thm: stein kernel} immediately implies the conclusion of Proposition \ref{coro:g-g-comparison}.
\end{proof}







Second, combining Theorem \ref{thm: stein kernel} with Lemma 4.6 in \cite{K19b} gives the following result:
\begin{corollary}[Multiplier-Bootstrap-to-Gaussian Comparison]\label{coro:beta-comparison}
Let $a_1,\dots,a_n$ be vectors  in $\mathbb R^p$ such that
$$
\min_{1\leq j\leq p}\frac{1}{n}\sum_{i=1}^n a_{ij}^2 \geq c\quad\text{and}\quad\max_{1\leq j\leq p}\frac{1}{n}\sum_{i=1}^na_{ij}^4 \leq B^2
$$
for some constants $c,B>0$.
Also, let $\varepsilon_1,\dots,\varepsilon_n$ be independent $N(0,1)$ random variables.
Moreover, for some constants $\alpha,\beta>0$, let $e_1,\dots,e_n$ be independent standardized Beta$(\alpha,\beta)$ random variables so that
\begin{equation}\label{eq: beta distribution mean variance}
{\mathrm{E}}[e_i]=0 \quad\text{and}\quad {\mathrm{E}}[e_i^2]=1,\quad\text{for all }i=1,\dots,n.
\end{equation}
Then, for the random vectors
$$
V = \frac{1}{\sqrt n}\sum_{i=1}^n e_i a_i\quad\text{and}\quad Z = \frac{1}{\sqrt n}\sum_{i=1}^n \varepsilon_i a_i
$$
we have
\begin{equation}\label{eq: beta distribution comparison}
\sup_{y\in\mathbb R^p}\Big|{\mathrm{P}}(V\leq y) - {\mathrm{P}}(Z\leq y)\Big| \leq C\left(\frac{B^2\log^5 p}{n}\right)^{1/4},
\end{equation}
where $C$ is a constant depending only on $c$, $\alpha$ and $\beta$.
\end{corollary}
\begin{proof}
Recall that $\eta\sim\text{Beta}(\alpha,\beta)$ has  density function $f_{\alpha,\beta}(x) \propto x^{\alpha - 1}(1-x)^{\beta - 1}$ for $x\in[0,1]$, mean $\mu = \alpha/(\alpha + \beta)$, and variance $\sigma^2 = \alpha\beta/((\alpha + \beta)^2(\alpha + \beta + 1))$.
By definition, the common distribution of the random variables $e_1,\dots,e_n$ equals that of $(\eta - \mu)/\sigma$.


Define
$$
\tau(x) = -\frac{\int_{-\mu/\sigma}^xsf(s)ds}{f(x)} = \frac{\int_x^{(1-\mu)/\sigma}sf(s)ds}{f(x)} \quad \text{for} \ x\in\left(-\frac{\mu}{\sigma},\frac{1-\mu}{\sigma}\right),
$$
where $f(x)=\sigma f_{\alpha,\beta}(\sigma x + \mu)$ for $x\in\big (-\frac{\mu}{\sigma},\frac{1-\mu}{\sigma}\big)$
is the density function of $(\eta - \mu)/\sigma$. From L'Hospital's rule, there exists a constant $C_1$ depending only on $\alpha$ and $\beta$ such that
$
|\tau(x)|\leq C_1$ for all $x\in\big(-\frac{\mu}{\sigma}, \frac{1-\mu}{\sigma}\big)$. Also, by integration by parts,
$
{\mathrm{E}}[e_1\varphi(e_1)]={\mathrm{E}}[\varphi'(e_1)\tau(e_1)]
$
for any continuously differentiable function $\varphi\colon\mathbb R\to\mathbb R$. Then, by Lemma 4.6 in \cite{K19b}, a Stein kernel $\tau^V$ for the random vector $V$ satisfies
$$
{\mathrm{E}}\left[ \max_{1\leq j,k\leq p}\left|\tau_{jk}^V(V) - \frac{1}{n}\sum_{i=1}^n a_{ij}a_{ik} \right| \right] \leq C_2\sqrt{\frac{\log p}{n}}\times\max_{1\leq j\leq p}\sqrt{\frac{1}{n}\sum_{i=1}^n a_{i j}^4}
$$
for some constant $C_2$ depending only on $C_1$. The desired conclusion  \eqref{eq: beta distribution comparison} follows from combining this bound with Theorem \ref{thm: stein kernel} and observing that ${\mathrm{E}}[Z_{j}Z_{k}] =n^{-1}\sum_{i=1}^na_{ij}a_{ik}$ for all $j,k=1,\dots,p$.
\end{proof}

\section{Proofs of Theorems \ref{thm: gaussian approximation main} and \ref{cor: rejection probabilities}}\label{sec: proofs}
In this section, we provide proofs of Theorems \ref{thm: gaussian approximation main} and \ref{cor: rejection probabilities}. Proofs of Theorems \ref{thm: gauss and rademacher} and \ref{thm: infinitesimal factors} will be given in Appendices \ref{sec: proof of theorem 23} and \ref{sec: proof of theorem 24} of the Supplemental Material. To simplify notation, we write
\[
\delta_n=\left(\frac{B_n^2\log^5(pn)}{n}\right)^{1/4}\quad\text{and}\quad
\upsilon_n=\sqrt{\frac{B_n^2\log^3(pn)}{n}}.
\]

Our proof strategy for Theorems \ref{thm: gaussian approximation main} and \ref{cor: rejection probabilities} is summarized as follows. First, we consider the multiplier bootstrap statistic $T_n^*$ with the weights $e_i$ constructed from the standardized Beta($\alpha$,$\beta$) distribution and parameters $\alpha$ and $\beta$ chosen so that ${\mathrm{E}}[e_i^3] = 1$. Thanks to Corollary \ref{coro:beta-comparison} and Proposition \ref{coro:g-g-comparison}, we have Gaussian approximation to this statistic with the rate $\delta_n$. This implies that
Condition A in Section \ref{sec: main arguments} is satisfied with $Z_i=e_i(X_i-\bar X_n)$ and $\delta=\delta_n$ due to the Gaussian anti-concentration inequality in Lemma \ref{lem: anticoncetration} of the Supplemental Material. In turn, the latter allows us to invoke Theorem \ref{cor: max}, which gives the approximation to $T_n$ by $T_n^*$ with the rate $\delta_n$. (Note that having ${\mathrm{E}}[e_i^3]=1$ is important here since otherwise Theorem \ref{cor: max} would give a slower approximation rate.) Combining this result with the aforementioned Gaussian approximation for $T_n^*$, we obtain the Gaussian approximation for $T_n$ with the rate $\delta_n$. This is done in Lemma \ref{coro: gaussian approximation} and gives Theorem \ref{thm: gaussian approximation main}.

Second, we consider the empirical bootstrap statistic $T_n^*$. Since we now have the Gaussian approximation for $T_n$ with the rate $\delta_n$, it follows that Condition A is satisfied with $Z_i=X_i$ and $\delta=\delta_n$. Hence, applying Theorem \ref{cor: max} with $V_i=X^*_i$ and $Z_i=X_i$, we can verify the empirical bootstrap approximation for $T_n$ with the rate $\delta_n$. This is done in Lemma \ref{thm: empirical bootstrap} and gives one part of Theorem \ref{cor: rejection probabilities}.

Third, we consider the multiplier bootstrap statistic $T_n^*$ with arbitrary weights $e_i$  satisfying \eqref{eq: multiplier bootstrap simplification}. By choosing parameters $\alpha$ and $\beta$ appropriately, we can match the first three moments of these weights by weights constructed from the standardized  Beta($\alpha$,$\beta$) distribution. Thus, yet another application of Theorem \ref{cor: max} allows us to link the distribution of any multiplier bootstrap statistic to the distribution of the multiplier bootstrap statistic with weights constructed from the standardized  Beta($\alpha$,$\beta$) distribution and further, via Corollary \ref{coro:beta-comparison} and Proposition \ref{coro:g-g-comparison}, to the Gaussian distribution. This leads to the Gaussian approximation for the multiplier bootstrap statistic $T_n^*$ with the rate $\delta_n$. This is done in Lemma \ref{eq: second order matching multiplier bootstrap ks metric} and gives the other part of Theorem \ref{cor: rejection probabilities}.


Before proceeding to the main body of the proofs, we present a few auxiliary results.
\begin{lemma}\label{lem: subgaussian growth}
Suppose that Condition E is satisfied. Then
\begin{equation}\label{eq: subgaussian bound on x}
\max_{1\leq i\leq n}\|X_i\|_{\infty} \leq 5B_n\log(p n)
\end{equation}
with probability at least $1 - 1/(2n^4)$. In addition,
$$
\max_{1\leq i\leq n}{\mathrm{E}}\left[\|X_i\|_{\infty}^8\right] \leq CB_n^8\log^8(p n),
$$
where $C$ is a universal constant.
\end{lemma}
\begin{proof}
By the union bound, Markov's inequality, and Condition E, we have for any $x>0$ that
\begin{align*}
&{\mathrm{P}}\left(\max_{1\leq i\leq n}\max_{1\leq j\leq p}|X_{i j}| > x\right)
 \leq pn \max_{1\leq i\leq n}\max_{1\leq j\leq p} {\mathrm{P}}(|X_{i j}| > x)\\
&\quad
\leq pn \max_{1\leq i\leq n}\max_{1\leq j\leq p} \frac{{\mathrm{E}}[\exp(|X_{i j}| / B_n)]}{\exp(x/B_n)} \leq 2pn \exp(-x/B_n).
\end{align*}
Substituting here $x = 5B_n\log(p n)$ gives the first asserted claim. The second asserted claim follows from combining Condition E, inequalities on page 95 in \cite{VW96}, and Lemma 2.2.2 in \cite{VW96}.
\end{proof}


\begin{lemma}\label{lem: all conditions}
Suppose that Conditions E and M are satisfied and set $\tilde X_i = X_i - \bar X_n$ for all $i = 1,\dots,n$.
Then there exist a universal constant $c\in(0,1]$ and constants $C>0$ and $n_0\in\mathbb N$ depending only on $b_1$ and $b_2$ such that for all $n\geq n_0$, if the inequality
\begin{equation}\label{eq: bn restriction once again}
B_n^2\log^5(p n)\leq c n
\end{equation}
holds, then the following events hold jointly with probability at least $1 - 1/n-3\upsilon_n$:
\begin{align}
&\frac{b_1^2}{2}\leq \frac{1}{n}\sum_{i=1}^n \tilde X_{i j}^2\quad \text{and} \quad \frac{1}{n}\sum_{i=1}^n \tilde X_{i j}^4 \leq 2B_n^2 b_2^{2}, \quad \text{for all }j=1,\dots,p,
\label{eq: upper and lower bounds for second moment} \\
&\max_{1\leq j,k\leq p}\left|\frac{1}{\sqrt n}\sum_{i=1}^n(\tilde X_{i j}\tilde X_{i k} - {\mathrm{E}}[X_{ij}X_{i k}])\right| \leq CB_n\sqrt{\log(p n)}, \label{eq: second order deviation} \\
&\max_{1\leq j,k,l\leq p}\left|\frac{1}{\sqrt n}\sum_{i=1}^n(\tilde X_{i j}\tilde X_{i k}\tilde X_{i l} - {\mathrm{E}}[X_{ij}X_{i k}X_{i l}])\right| \leq CB_n^2\sqrt{\log^3(p n)}. \label{eq: third order deviation}
\end{align}
\end{lemma}

The proof of this lemma is rather standard but long, and so is deferred to Appendix \ref{sec: proof lemma 52} of the Supplemental Material.






\begin{lemma}\label{coro: gaussian approximation}
Suppose that Conditions E and M are satisfied. Then
\begin{equation}\label{eq: gaussian approximation ks metric}
\sup_{x\in\mathbb{R}}|{\mathrm{P}}(T_n\leq x)-{\mathrm{P}}(T^G_n\leq x)|\leq C\left(\frac{B_n^2\log^5(pn)}{n}\right)^{1/4},
\end{equation}
where $C$ is a constant depending only on $b_1$ and $b_2$.
\end{lemma}

\begin{proof}
Without loss of generality, we may assume that \eqref{eq: bn restriction once again} holds and that $n$ is large enough so that $n\geq n_0$ for $n_0$ from Lemma \ref{lem: all conditions}, since otherwise the conclusion of the lemma is trivial by taking $C$ large enough. This will justify an application of Lemma \ref{lem: all conditions} when needed. In addition, by again taking $C$ large enough, we may assume that $1/n^4 + 2/n + 3v_n < 1$.

Let $\mathcal{A}_n$ be the event that \eqref{eq: subgaussian bound on x} and \eqref{eq: upper and lower bounds for second moment}--\eqref{eq: third order deviation} hold jointly. By Lemmas \ref{lem: subgaussian growth} and \ref{lem: all conditions}, ${\mathrm{P}}(\mathcal A_n)\geq 1 - 1/(2n^4) - 1/n -3v_n > 0$. Further, let $e_1,\dots,e_n$ be independent standardized Beta$(1/2,3/2)$ random variables, standardized in such a way that they have mean zero and unit variance (cf. Corollary \ref{coro:beta-comparison}), that are independent of $X_{1:n} = (X_1,\dots,X_n)$. It is not difficult to check that  ${\mathrm{E}}[e_i^3]=1$ for all $i=1,\dots,n$.

Let $T_n^*$ be the multiplier bootstrap statistic with weights $e_1,\dots,e_n$.
Condition on $X_{1:n}$ such that $\mathcal A_n$ holds.
Then, by Corollary \ref{coro:beta-comparison} and the definition of $\mathcal A_n$, we have
\begin{equation}\label{beta-to-gauss}
\sup_{y\in\mathbb{R}^p}\left|{\mathrm{P}}\left( \frac{1}{\sqrt n}\sum_{i=1}^n e_i(X_i - \bar X_n) \leq y \mid X_{1:n} \right)-{\mathrm{P}}(\hat{G}\leq y\mid X_{1:n})\right|\leq C_1\delta_n,
\end{equation}
while by Proposition \ref{coro:g-g-comparison}, we have
\begin{equation*}
\sup_{y\in\mathbb{R}^p}|{\mathrm{P}}(\hat{G}\leq y\mid X_{1:n})-{\mathrm{P}}(G\leq y)|\leq C_2\delta_n,
\end{equation*}
where  $C_1$ and $C_2$ are constants depending only on $b_1$ and $b_2$.

Next, we shall invoke Theorem \ref{cor: max} to  compare the distribution of $T_n$ with the conditional distribution of $T_n^*$. Formally, let $Y_1,\dots,Y_n$ be independent copies of $X_1,\dots,X_n$ that are independent of $X_{1:n}$, and define $T_n'$ by $T_n$ with $X_i$'s replaced by $Y_i$'s. Then, ${\mathrm{P}} (T_n \le x) = {\mathrm{P}} (T_n' \le x \mid X_{1:n})$. Condition on $X_{1:n}$ such that $\mathcal A_n$ holds and apply Theorem \ref{cor: max} with $V_i = Y_i$ and $Z_i = e_i \tilde X_i$ for all $i=1,\dots,n$. Since ${\mathrm{E}}[e_i]=0$ and ${\mathrm{E}}[e_i^2] = {\mathrm{E}}[e_i^3]=1$ for all $i=1,\dots,n$, it is not difficult to see from the definition of $\mathcal A_n$ that Conditions V, P, and B, as well as inequalities (\ref{eq: bn bounds 1}) and (\ref{eq: bn bounds 2}) of Theorem \ref{cor: max} are satisfied with appropriate constants $C_v$, $C_p$, $C_b$, and $C_m$ that depend only on $b_1,b_2$.
It remains to verify Condition A in Theorem \ref{cor: max}.
Observe that for  any $y\in\mathbb R^p$ and $t>0$,
\begin{equation}
 \label{beta-anti}
\begin{split}
&{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n e_i(X_i - \bar X_n) \leq y + t \mid X_{1:n}\right)  \leq {\mathrm{P}}\left(\hat G \leq y+t \mid X_{1:n}\right) + C_1\delta_n \quad (\text{by (\ref{beta-to-gauss})}) \\
&\quad  \leq {\mathrm{P}}\left(\hat G \leq y\mid X_{1:n}\right) + K_1t\sqrt{\log p} + C_1\delta_n \quad (\text{by Lemma \ref{lem: anticoncetration} and (\ref{eq: upper and lower bounds for second moment})}) \\
&\quad  \leq{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n e_i(X_i - \bar X_n) \leq y \mid X_{1:n}\right)
+K_1t\sqrt{\log p} + 2C_1\delta_n, \quad (\text{by (\ref{beta-to-gauss})})
\end{split}
\end{equation}
where $K_1>0$ is a constant depending only on $b_1$. Thus, applying Theorem \ref{cor: max}, we conclude that
\begin{equation*}
\sup_{x\in\mathbb{R}}|{\mathrm{P}}(T_n\leq x)-{\mathrm{P}}(T^*_n\leq x\mid X_{1:n})| = \sup_{x\in\mathbb{R}}|{\mathrm{P}}(T_n'\leq x \mid X_{1:n})-{\mathrm{P}}(T^*_n\leq x\mid X_{1:n})|\leq C_3\delta_n
\end{equation*}
 for some constant $C_3$ depending only on $b_1$ and $b_2$. The asserted claim follows from these bounds via the triangle inequality by noting that the left-hand side of \eqref{eq: gaussian approximation ks metric} is non-stochastic, so that if \eqref{eq: gaussian approximation ks metric} holds with strictly positive probability (recall that ${\mathrm{P}}(\mathcal A_n)>0$), then it holds with probability one.
\end{proof}

\begin{lemma}\label{thm: empirical anticoncentration}
Suppose that Conditions E and M are satisfied. Then for any $y\in\mathbb R^p$ and $t > 0$,
\begin{align*}
&{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n X_i \leq y+t\right) - {\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n X_i \leq y\right) \leq C\left(t\sqrt{\log p} + \left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4}\right),
\end{align*}
where $C$ is a constant depending only on $b_1$ and $b_2$.
\end{lemma}


\begin{proof}
Fix $y\in\mathbb R^p$ and $t>0$. Then for some constant $C$ depending only on $b_1$ and $b_2$,
\begin{align*}
{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n X_i \leq y + t\right) &\leq {\mathrm{P}}(G \leq y + t) + C\delta_n
\leq {\mathrm{P}}(G \leq y ) + Ct\sqrt{\log p} + C\delta_n\\
&\leq {\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n X_i \leq y\right) + Ct\sqrt{\log p} + 2C\delta_n,
\end{align*}
where the first and the third inequalities follow from Lemma \ref{coro: gaussian approximation} and the second from Lemma \ref{lem: anticoncetration} of the Supplemental Material. This gives the asserted claim.
\end{proof}




\begin{lemma}\label{thm: empirical bootstrap}
Suppose that Conditions E and M are satisfied and that the random variables $X_1^*,\dots,X_n^*$ are obtained via the empirical bootstrap. Then with probability at least $1 - 2/n-3\upsilon_n$, we have
$$
\sup_{x\in\mathbb R}\left| {\mathrm{P}}\left(T_n \leq x\right) - {\mathrm{P}}\left(T_n^* \leq x \mid X_{1:n}\right) \right| \leq C\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4},
$$
where $C$ is a constant depending only on $b_1$ and $b_2$.
\end{lemma}



\begin{proof}
As before, we may assume that \eqref{eq: bn restriction once again} holds and that $n$ is large enough so that $n\geq n_0$ for $n_0$ from Lemma \ref{lem: all conditions}, since otherwise the conclusion of the lemma is trivial by taking $C$ large enough. This will justify an application of Lemma \ref{lem: all conditions} when needed.

Let $Y_1,\dots,Y_n$ be vectors in $\mathbb R^p$ such that
\begin{equation}\label{eq: Y infty bound}
\|Y_i\|_{\infty} \leq 10B_n\log(p n)\quad\text{for all }i=1,\dots,n,
\end{equation}
\begin{equation}\label{eq: Y2 bound lower and upper}
b_1^2/2\leq \frac{1}{n}\sum_{i=1}^n Y_{i j}^2 \quad \text{and} \quad \frac{1}{n}\sum_{i=1}^n Y_{i j}^4 \leq 2B_n^2 b_2^{2},\quad\text{for all }j=1,\dots,p,
\end{equation}
\begin{equation}\label{eq: Y2 bound plus}
\max_{1\leq j,k\leq p}\left|\frac{1}{\sqrt n}\sum_{i=1}^n(Y_{i j}Y_{i k} - {\mathrm{E}}[X_{ij}X_{i k}])\right| \leq C_mB_n\sqrt{\log(p n)},
\end{equation}
and
\begin{equation}\label{eq: Y3 bound}
\max_{1\leq j,k,l\leq p}\left|\frac{1}{\sqrt n}\sum_{i=1}^n(Y_{i j}Y_{i k}Y_{i l} - {\mathrm{E}}[X_{ij}X_{i k}X_{i l}])\right| \leq C_mB_n^2\sqrt{\log^3(p n)},
\end{equation}
where $C_m$ is the constant $C$ from Lemma \ref{lem: all conditions}. Also, let $Y_1^*,\dots,Y_n^*$ be independent random vectors with each $Y_i^*$ having uniform distribution on $\{Y_1,\dots,Y_n\}$.

To prove the asserted claim, we will apply Theorem \ref{cor: max} with $V_i = Y_i^*$ and $Z_i = X_i$ for all $i=1,\dots,n$. Conditions V, P, and B with constants $C_v$, $C_p$, and $C_b$ depending only on $b_1$ and $b_2$ follow immediately from Conditions E and M, Lemma \ref{lem: subgaussian growth}, and the inequalities in \eqref{eq: Y infty bound} and \eqref{eq: Y2 bound lower and upper}. Also, Condition A with $\delta=\delta_n$ and $C_a$ depending only on $b_1$ and $b_2$ follows from Lemma \ref{thm: empirical anticoncentration}. Hence, an application of Theorem \ref{cor: max} is justified if we can verify \eqref{eq: bn bounds 1} and \eqref{eq: bn bounds 2} but these inequalities follow from \eqref{eq: Y2 bound plus} and \eqref{eq: Y3 bound} by noting that
$$
\frac{1}{\sqrt n}\sum_{i=1}^n({\mathrm{E}}[V_{i j}V_{i k}] - Y_{i j}Y_{i k}) = 0 \quad \text{and} \quad \frac{1}{\sqrt n}\sum_{i=1}^n({\mathrm{E}}[V_{i j}V_{i k}V_{i l}] - Y_{i j}Y_{i k}Y_{i l}) = 0
$$
for all $j,k,l=1,\dots,p$.
Now, applying Theorem \ref{cor: max} shows that for all $y\in\mathbb R^p$, we have
$$
\left|{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n V_i\leq y\right) -  {\mathrm{P}}\left( \frac{1}{\sqrt n}\sum_{i=1}^n X_i \leq y \right)\right|\leq K_1\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4}
$$
for some constant $K_1$ depending only on $b_1$, $b_2$, and $C_m$. The asserted claim follows from this bound by setting $Y_i = X_i - \bar X_n$ for all $i = 1,\dots,n$, and noting that in this case \eqref{eq: Y infty bound} holds with probability at least $1 - 1/(2n^4)$ by Lemma \ref{lem: subgaussian growth} and  \eqref{eq: Y2 bound lower and upper}, \eqref{eq: Y2 bound plus}, and \eqref{eq: Y3 bound} hold jointly with probability at least $1 - 1/n-3\upsilon_n$ by Lemma \ref{lem: all conditions}.
\end{proof}

\begin{lemma}\label{eq: second order matching multiplier bootstrap ks metric}
Suppose that Conditions E and M are satisfied and that the random variables $X_1^*,\dots,X_n^*$ are obtained via the multiplier bootstrap with weights $e_1,\dots,e_n$ satisfying \eqref{eq: multiplier bootstrap simplification}. Then with probability at least $1-2/n-3v_n$, we have
\[
\sup_{x\in\mathbb{R}}|{\mathrm{P}}(T_n\leq x)-{\mathrm{P}}(T^*_n\leq x\mid X_{1:n})|\leq C\left(\frac{B_n^2\log^5(pn)}{n}\right)^{1/4},
\]
where $C$ is a constant depending only on ${\mathrm{E}}[e_1^3]$, $b_1$ and $b_2$.
\end{lemma}

\begin{remark}
The constant $C$ in this result depends on ${\mathrm{E}}[e_1^3]$ continuously, and so we can take $C$ independent of ${\mathrm{E}}[e_1^3]$ under the implicitly maintained assumption that \eqref{eq: multiplier bootstrap simplification} holds.
\end{remark}

\begin{proof}
As before, we may assume that \eqref{eq: bn restriction once again} holds and that $n$ is large enough so that $n\geq n_0$ for $n_0$ from Lemma \ref{lem: all conditions}, since otherwise the conclusion of the lemma is trivial by taking $C$ large enough. This will justify an application of Lemma \ref{lem: all conditions} when needed.

Let $\mathcal{A}_n$ be the event that \eqref{eq: subgaussian bound on x} and \eqref{eq: upper and lower bounds for second moment}--\eqref{eq: third order deviation} hold jointly.
By Lemmas \ref{lem: subgaussian growth} and \ref{lem: all conditions}, we have ${\mathrm{P}}(\mathcal{A}_n)\geq1-2/n-3\upsilon_n$. Moreover, by Proposition \ref{coro:g-g-comparison},
\begin{equation}\label{eq: lem 56 beginning}
\sup_{y\in\mathbb{R}^p}|{\mathrm{P}}(\hat G\leq y\mid X_{1:n})-{\mathrm{P}}(G \leq y)|\leq C_1\delta_n
\end{equation}
on the event $\mathcal A_n$, where $C_1$ is a constant depending only on $b_1$ and $b_2$.

Next, we claim that the case with $\sigma_e > 0$ can be reduced to the case with $\sigma_e = 0$ (and the constant 3 appearing in \eqref{eq: multiplier bootstrap simplification} replaced by some other universal constant). To prove this claim, define random variables $e'_1,\dots,e'_n$ as in Corollary \ref{coro:beta-comparison} with $\alpha=\beta=1$ such that they are independent of everything else. Then on the event $\mathcal A_n$, by Corollary \ref{coro:beta-comparison}, we have that
\[
\sup_{y\in\mathbb{R}^p}\left|{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n e'_i\tilde X_i\leq y\mid X_{1:n}\right)-{\mathrm{P}}\left(\frac{1}{\sigma_e\sqrt n}\sum_{i=1}^ne_{i,1}\tilde X_i\leq y\mid X_{1:n}\right)\right|\leq C_2\delta_n,
\]
where $\tilde X_i=X_i - \bar X_n$ for all $i=1,\dots,n$ and $C_2$ is a constant depending only on $b_1$ and $b_2$. Therefore, noting that the sequences $\{e_{i,1}\}_{i=1}^n$, $\{e_{i,2}\}_{i=1}^n$, and $\{e'_{i}\}_{i=1}^n$ are independent, we have on $\mathcal A_n$ that
\begin{align*}
&\sup_{y\in\mathbb{R}^p}\left|{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^ne_{i}\tilde X_i\leq y\mid X_{1:n}\right)-{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n(\sigma_ee'_i+e_{i,2})\tilde X_i\leq y\mid X_{1:n}\right)\right|\\
&\quad\leq{\mathrm{E}}\Bigg[\sup_{y\in\mathbb{R}^p}\Bigg|{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^ne_{i,1}\tilde X_i\leq y\mid X_{1:n},\{e_{i,2}\}_{i=1}^n\right) \\
&\qquad\qquad\qquad-{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n\sigma_e e'_i\tilde X_i\leq y\mid X_{1:n},\{e_{i,2}\}_{i=1}^n\right)\Bigg|\mid X_{1:n}\Bigg] \leq C_2\delta_n.
\end{align*}
Thus, it suffices to prove the asserted claim with $e_i$'s replaced by $\sigma_ee'_i+e_{i,2}$'s, which are bounded by a universal constant (note that $\sigma_e \le 1$ since $e_i$ has unit variance).

Further, define the function $f\colon(0,1)\to\mathbb{R}$ by
\[
f(\alpha)=\frac{2\sqrt{2}(1-2\alpha)}{3\sqrt{\alpha(1-\alpha)}},\qquad\text{for all }\alpha\in(0,1).
\]
One can directly check that $f(\alpha)$ is the skewness of the Beta$(\alpha,1-\alpha)$ distribution for all $\alpha\in(0,1)$. Since $\lim_{\alpha\to0}f(\alpha)=\infty$, $\lim_{\alpha\to1}f(\alpha)=-\infty$ and $f$ is continuous, there is an $\alpha^*\in(0,1)$ satisfying $f(\alpha^*)={\mathrm{E}}[e_{1}^3]$.
We define random variables $\tilde e_1,\dots,\tilde e_n$ as in Corollary \ref{coro:beta-comparison} with $\alpha=\alpha^*$ and $\beta=1-\alpha^*$ such that they are independent of everything else. It is then easy to check that ${\mathrm{E}}[\tilde e_i]=0$, ${\mathrm{E}}[\tilde e_i^2]=1$, and ${\mathrm{E}}[\tilde e_i^3]={\mathrm{E}}[e_i^3]$ for all $i = 1,\dots,n$. Also, applying Corollary \ref{coro:beta-comparison}, we have on $\mathcal{A}_n$ that
\begin{equation}\label{eq: third order to gaussian multipliers}
\sup_{y\in\mathbb{R}^p}\left|{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n \tilde e_i \tilde X_i \leq y\mid X_{1:n}\right)-{\mathrm{P}}\left(\hat G \leq y\mid X_{1:n}\right)\right|\leq C_3\delta_n,
\end{equation}
where $C_3$ is a constant depending only on $\alpha^*$, $b_1$ and $b_2$.

We now apply Theorem \ref{cor: max} with $V_i = e_i \tilde X_i$ and $Z_i = \tilde e_i \tilde X_i$ for all $i=1,\dots,n$ conditional on $X_{1:n}$ on the event $\mathcal A_n$. Conditions V, P, and B with $C_v$, $C_p$, and $C_b$ depending only on $\alpha^*$, $b_1$ and $b_2$ follow immediately from the inequalities \eqref{eq: subgaussian bound on x} and \eqref{eq: upper and lower bounds for second moment} and the boundedness of $e_{i}$'s and $\tilde e_i$'s. Condition A with $\delta=\delta_n$ follows from \eqref{eq: third order to gaussian multipliers} and the derivation in \eqref{beta-anti}. Moreover, \eqref{eq: bn bounds 1} and \eqref{eq: bn bounds 2} are evident by construction. Thus, by Theorem \ref{cor: max}, we have on $\mathcal{A}_n$  that
\begin{equation}\label{eq: lemma 56 end}
\sup_{y\in\mathbb{R}^p}\left|{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n \tilde e_i \tilde X_i \leq y\mid X_{1:n}\right)-{\mathrm{P}}\left(\frac{1}{\sqrt n}\sum_{i=1}^n  e_i \tilde X_i  \leq y\mid X_{1:n}\right)\right|\leq C_4\delta_n,
\end{equation}
where $C_4$ is a constant depending only on $\alpha^*$, $b_1$ and $b_2$. The asserted claim now follows from combining \eqref{eq: lem 56 beginning}, \eqref{eq: third order to gaussian multipliers}, and \eqref{eq: lemma 56 end} via the triangle inequality and using Lemma \ref{coro: gaussian approximation}.
\end{proof}

We are now in the position to prove the main results from Section \ref{sec: main results}.


\begin{proof}[Proof of Theorem \ref{thm: gaussian approximation main}]
The asserted claim follows immediately from Lemma \ref{coro: gaussian approximation} by applying \eqref{eq: gaussian approximation ks metric} with $x = c_{1-\alpha}^G$.
\end{proof}


\begin{proof}[Proof of Theorem \ref{cor: rejection probabilities}]
Let $C_1$, $C_2$, and $C_3$ be the constants $C$ in Lemmas \ref{thm: empirical anticoncentration}, \ref{thm: empirical bootstrap}, and \ref{eq: second order matching multiplier bootstrap ks metric}, respectively. Set
$$
\beta_n = (1\vee C_1\vee C_2\vee C_3)\left( \frac{B_n^2\log^5(p n)}{n} \right)^{1/4}.
$$
By Lemmas \ref{thm: empirical bootstrap} and \ref{eq: second order matching multiplier bootstrap ks metric}, we have
$
\sup_{x\in\mathbb R}\left| {\mathrm{P}}(T_n \leq x) - {\mathrm{P}}(T_n^* \leq x\mid X_{1:n}) \right| \leq \beta_n
$
with probability at least $1 - 2/n-3\upsilon_n$. Hence, letting $c_{1-\gamma}$ be the $(1-\gamma)$th quantile of $T_n$ for all $\gamma\in(0,1)$, we have with the same probability that
$$
{\mathrm{P}}(T_n^* \leq c_{1-\alpha + \beta_n}\mid X_{1:n}) \geq {\mathrm{P}}(T_n \leq c_{1-\alpha + \beta_n}) - \beta_n \geq 1 - \alpha, \quad \text{and}
$$
\begin{align*}
{\mathrm{P}}(T_n^* \leq c_{1-\alpha - 3\beta_n} \mid X_{1:n}) & \leq {\mathrm{P}}(T_n \leq c_{1-\alpha - 3\beta_n}) + \beta_n\\
& \leq 1-\alpha - 2\beta_n + C_{1}\left(\frac{B_n^2\log^5(p n)}{n}\right)^{1/4} < 1 - \alpha,
\end{align*}
where the second inequality follows from Lemma \ref{thm: empirical anticoncentration}. Therefore,
$$
{\mathrm{P}}(c_{1-\alpha - 3\beta_n} < c_{1-\alpha}^B \leq c_{1-\alpha + \beta_n}) \geq 1 - 2/n-3\upsilon_n\geq1-5\upsilon_n,
$$
so that
$$
{\mathrm{P}}(T_n > c_{1-\alpha}^B) \leq {\mathrm{P}}(T_n > c_{1-\alpha - 3\beta_n}) + 5\upsilon_n \leq \alpha + 3\beta_n + 5\upsilon_n \leq \alpha + 8\beta_n \quad \text{and}
$$
\begin{align*}
{\mathrm{P}}(T_n > c_{1-\alpha}^B) &\geq {\mathrm{P}}(T_n > c_{1-\alpha + \beta_n}) - 5\upsilon_n \\
& \geq \alpha - \beta_n - C_{1}\left( \frac{B_n^2\log^5(p n)}{n} \right)^{1/4} - 5\upsilon_n \geq \alpha - 7\beta_n,
\end{align*}
where the second inequality follows from Lemma \ref{thm: empirical anticoncentration}. Combining these inequalities gives the asserted claim.
\end{proof}





\begin{acks}[Acknowledgments]
We are grateful to Tim Armstrong, Matias Cattaneo, Xiaohong Chen, and Tengyuan Liang for helpful discussions. We also thank seminar participants at  the University of Pennsylvania and Yale University.
\end{acks}

\begin{funding}
 K. Kato is suport by the NSF DMS-1952306 and DMS-2014636.
\end{funding}

\begin{supplement}
The Supplementary Material contains proofs omitted in the main text as well as several technical tools, and the simulation results.
\end{supplement}


\begin{thebibliography}{99}
\bibitem[Andrews and Shi(2013)]{AS13}
Andrews, D. and Shi, X. (2013). Inference based on conditional moment inequalities. {\em Econometrica} \textbf{81} 609-666.


\bibitem[Bai et al.(2019)]{BSS19}
Bai, Y., Santos, A., and Shaikh, A. (2019). A practical method for testing many moment inequalities. {\em Becker Friedman Institute for Economics Working Paper}.

\bibitem[Belloni et al.(2018)]{BCCHK18}
Belloni, A., Chernozhukov, V., Chetverikov, D., Hansen, C., and Kato, K. (2018). High-dimensional econometrics and regularized GMM. {\em arXiv:1806.01888}.

\bibitem[Belloni et al.(2018)]{BCCW18}
Belloni, A., Chernozhukov, V., Chetverikov, D., and Wei, Y. (2018). Uniformly valid post-regularization confidence regions for many functional parameters in Z-estimation framework. {\em Annals of Statistics} \textbf{46} 3643-3675.



\bibitem[Belloni et al.(2015)]{BCK15}
Belloni, A., Chernozhukov, V., and Kato, K. (2015). Uniform post-selection inference for least absolute deviation regression and other Z-estimation problems. {\em Biometrika} \textbf{102} 77--94.

\bibitem[Bentkus(1984)]{B85}
Bentkus, V. (1985). Lower bounds for the rate of convergence in the central limit theorem in Banach spaces. {\em Litovskii Matematicheskii Sbornik} \textbf{25} 10-21.








\bibitem[Canay et al.(2016)]{CRS16}
Canay, I., Romano, J. and Shaikh, A. (2016). Randomization tests under an approximate symmetry assumption. {\em Econometrica} \textbf{85} 1013-1030.

\bibitem[Chatterjee(2006)]{C06}
Chatterjee, S. (2006). A generalization of the Lindeberg principle. {\em Annals of Probability} \textbf{34} 2061-2076.


\bibitem[Chen(2018)]{Ch18}
Chen, X. (2018). Gaussian and bootstrap approximations for high-dimensional U-statistics and their applications. \textit{Annals of Statistics} \textbf{46} 642-678.

\bibitem[Chen and Kato(2019)]{CK19a}
Chen, X. and Kato, K. (2019). Randomized incomplete $U$-statistics in high dimensions. \textit{Annals of Statistics} \textbf{47}  3127-3156.

\bibitem[Chen and Kato(2020)]{CK19b}
Chen, X. and Kato, K. (2020). Jackknife multiplier bootstrap: finite sample approximations to the U-process supremum with applications. \textit{Probability Theory and Related Fields} \textbf{176} 1097-1163.



\bibitem[Chernozhukov et al.(2013)]{CCK13}
Chernozhukov, V., Chetverikov, D. and Kato, K. (2013). Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors. {\em Annals of Statistics} \textbf{41} 2786-2819.

\bibitem[Chernozhukov et al.(2014)]{CCK14}
Chernozhukov, V., Chetverikov, D. and Kato, K. (2014). Anti-concentration and honest, adaptive confidence bands. {\em Annals of Statistics} \textbf{42} 1787-1818.

\bibitem[Chernozhukov et al.(2015)]{CCK15}
Chernozhukov, V., Chetverikov, D. and Kato, K. (2015). Comparison and anti-concentration bounds for maxima of Gaussian random vectors. {\em Probability Theory and Related Fields} \textbf{162} 47-70.

\bibitem[Chernozhukov et al.(2017)]{CCK17}
Chernozhukov, V., Chetverikov, D. and Kato, K. (2017). Central limit theorems and bootstrap in high dimensions. {\em Annals of Probability} \textbf{45} 2309-2352.

\bibitem[Chernozhukov et al.(2019)]{CCK19}
Chernozhukov, V., Chetverikov, D. and Kato, K. (2019). Inference on causal and structural parameters using many moment inequalities. {\em Review of Economic Studies} \textbf{86} 1867-1900.

\bibitem[Chernozhukov et al.(2020)]{CCK20}
Chernozhukov, V., Chetverikov, D. and Koike, Y. (2020). Nearly optimal central limit theorem and bootstrap approximations in high dimensions. {\em arXiv:2012.09513}.

\bibitem[Chesher and Rosen(2019)]{CR19}
Chesher, A. and Rosen, A. (2019). Generalized instrumental variable models, methods, and applications. {\em Cemmap working paper CWP41/19}.

\bibitem[Chetverikov(2018)]{C18}
Chetverikov, D. (2018). Adaptive tests of conditional moment inequalities. {\em Econometric Theory} \textbf{34} 186-227.

\bibitem[Chetverikov(2019)]{C19}
Chetverikov, D. (2019). Testing regression monotonicity in econometric models. {\em Econometric Theory} \textbf{35} 729-776.

\bibitem[Chetverikov et al.(2018)]{CWK18}
Chetverikov, D., Wilhelm, D., and Kim, D. (2018). An adaptive test of stochastic monotonicity. {\em Cemmap working paper CWP24/18}.

\bibitem[Deng and Zhang(2020)]{DZ17}
Deng, H. and Zhang, C.-H. (2020). Beyond Gaussian approximation: bootstrap for maxima of sums of independent random vectors. \textit{Annals of Statistics} \textbf{48}  3643-3671.



\bibitem[Fang and Koike(2020)]{FK20}
Fang, X. and Koike, Y. (2020). High-dimensional Central Limit Theorems by Stein's Method.
To appear in \textit{Annals of Applied Probability}.

\bibitem[Fathi(2019)]{Fa19}
Fathi, M. (2020). Stein kernels and moment maps.
\textit{Annals of Probability} \textbf{47} 2172--2185.


\bibitem[Hansen(2005)]{H05}
Hansen, P. (2005). A test for superior predictive ability. {\em Journal of Business and Economic Statistics} \textbf{23} 365-380.

\bibitem[Hansen et al.(2011)]{HLN11}
Hansen, P., Lunde, A. and Nason, J. (2011). The model confidence set. {\em Econometrica} \textbf{79} 453-497.


\bibitem[Horowitz(2001)]{H01}
Horowitz, J. (2001). The bootstrap. {\em Handbook of Econometrics, Volume 5} 3159-3228.


\bibitem[Koike(2019a)]{K19}
Koike, Y. (2019a). Gaussian approximation of maxima of {W}iener functionals and its application to high-frequency data. \textit{Annals of Statistics} \textbf{47} 1663--1687.

\bibitem[Koike(2019b)]{K19a}
Koike, Y. (2019b). High-dimensional central limit theorems for homogeneous sums. {\em arXiv:1902.03809}.

\bibitem[Koike(2019c)]{K19b}
Koike, Y. (2021). Notes on the dimension dependence in high-dimensional central limit theorems for hyperrectangles. {\em Japanese Journal of Statistics and Data Science} \textbf{1}  257--297.

\bibitem[Koning and Bekker(2019)]{KB19}
Koning, N. and Bekker, P. (2019). Exact testing of many moment inequalities against multiple violations. {\em arXiv:1904.12775}.

\bibitem[Kuchibhotla and Rinaldo(2020)]{KR20}
Kuchibhotla, A. and Rinaldo, A. (2020). High-dimensional CLT for sums of non-degenerate random vectors: $n^{-1/2}$ rate. {\em arXiv:2009.13673}.


\bibitem[Lehmann and Romano(2005)]{LR05}
Lehmann, E. and Romano, J. (2005). {\em Testing Statistical Hypotheses}. Springer Texts in Statistics.

\bibitem[Lopes(2020)]{L20}
Lopes, M. (2020). Central limit theorem and bootstrap approximation in high dimensions with near $1/\sqrt n$ rates. {\em arXiv:2009.06004}.

\bibitem[Lopes, Lin and M\"uller(2020)]{LLM20}
Lopes, M., Lin, Z., and M\"uller, H.-G. (2020). Bootstrapping max statistics in high dimensions: Near-parametric rates under weak variance decay and application to functional and multinomial data.
{\em Annals of Statistics} \textbf{48} 1214--1229.

\bibitem[Mammen(1993)]{M93}
Mammen, E. (1993). Bootstrap and wild bootstrap for high dimensional linear models. {\em Annals of Statistics} \textbf{21} 255-285.

\bibitem[Meckes(2009)]{M09}
Meckes, E. (2009). On Stein's method for multivariate normal approximation. In: {\em High Dimensional Probability V: The Luminy Volume}, IMS Collections, Vol.5,  pp.159--178.
\bibitem[Nazarov(2003)]{Na03}
Nazarov, F. (2003). On the maximal perimeter of a convex set in $\mathbb R^n$ with respect to a Gaussian measure.
In {\em Geometric Aspects of Functional Analysis. Lecture Notes in Math.} \textbf{1807} 169--187. Springer.





\bibitem[O'Donnell et al.(2018)]{OST18}
O'Donnell, R., Servedio, R., and Tan, L. (2019). Fooling polytopes. {\em STOC 2019: Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing}, pp.  614--625


\bibitem[Romano et al.(2014)]{RSW14}
Romano, J., Shaikh, A. and Wolf, M. (2014). A practical two-step method for testing moment inequalities. {\em Econometrica} \textbf{82} 1979-2002.

\bibitem[Tanguy (2015)]{Tanguy15}
Tanguy, K. (2015). Some superconcentration inequalities for extrema of stationary Gaussian processes.
{\em Statistics and Probability Letters} \textbf{106} 239-246.




\bibitem[van der Vaart and Wellner(1996)]{VW96}
van der Vaart, A.W. and Wellner, J.A. (1996). {\em Weak Convergence and Empirical Processes: With Applications to Statistics}. Springer.


\bibitem[Vershynin(2018)]{V11}
Vershynin, R. (2018). {\em High-Dimensional Probability}. Cambridge University Press.


\bibitem[Wainwright(2019)]{Wa19}
Wainwright, M. J. (2019). {\em High-Dimensional Statistics}. Cambridge University Press.


\bibitem[White(2000)]{W00}
White, H. (2000). A reality check for data snooping. {\em Econometrica} \textbf{5} 1007-1126.

\bibitem[Zhang and Cheng(2017)]{ZC17}
Zhang, X. and Cheng, G. (2017). Gaussian approximation for high dimensional vector under physical dependence. \textit{Bernoulli} \textbf{24} 2640--2675.
\bibitem[Zhang and Wu(2017)]{ZW17}
Zhang, D. and Wu, W.-B. (2017). Gaussian approximation for high dimensional time series. \textit{Annals of Statistics} \textbf{45} 1895-1919.
\end{thebibliography}

\newpage
\centerline{\textbf{\Large{Supplementary Material}}}