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.
73,903 characters
Gaussian Approximations and Multiplier Bootstrap for Maxima of Sums of High-Dimensional Random Vectors
\begin{frontmatter}
\title{Gaussian Approximations and Multiplier Bootstrap for Maxima of Sums of High-Dimensional Random Vectors\thanksref{T1}}
\runtitle{Gaussian Approximations and multiplier bootstrap}
\thankstext{T1}{Date: June, 2012. Revised June, 2013. V. Chernozhukov and D. Chetverikov are supported by a National Science Foundation grant. K. Kato is supported by the Grant-in-Aid for Young Scientists (B) (25780152), the Japan Society for the Promotion of Science.}
\begin{aug}
\author{\fnms{Victor} \snm{Chernozhukov}\thanksref{m1}\ead[label=e1]{[email removed]}},
\author{\fnms{Denis} \snm{Chetverikov}\thanksref{m2}\ead[label=e2]{[email removed]}}
\and
\author{\fnms{Kengo} \snm{Kato}\thanksref{m3}
\ead[label=e3]{[email removed]}}
\runauthor{Chernozhukov Chetverikov Kato}
\affiliation{MIT\thanksmark{m1}, UCLA\thanksmark{m2}, and University of Tokyo\thanksmark{m3}}
\address{Department of Economics and\\
Operations Research Center, MIT \\
50 Memorial Drive \\
Cambridge, MA 02142, USA.\\
\printead{e1}}
\address{Department of Economics, UCLA\\
Bunche Hall, 8283 \\
315 Portola Plaza \\
Los Angeles, CA 90095, USA.\\
\printead{e2}}
\address{Graduate School of Economics \\
University of Tokyo \\
7-3-1 Hongo, Bunkyo-ku\\
Tokyo 113-0033, Japan. \\
\printead{e3}}
\end{aug}
\begin{abstract}
{We derive a Gaussian approximation result for the maximum of a sum of high dimensional random vectors.
Specifically, we establish conditions under which the distribution of the maximum is approximated by that of the maximum of a sum of the Gaussian random vectors with the same covariance matrices as the original vectors. This result applies when the dimension of random vectors ($p$) is large compared to the sample size ($n$); in fact, $p$ can be much larger than $n$, without restricting correlations of the coordinates of these vectors. We also show that the distribution of the maximum of a sum of the random vectors with unknown covariance matrices can be consistently estimated by the distribution of the maximum of a sum of the conditional Gaussian random vectors obtained by multiplying the original vectors with i.i.d. Gaussian multipliers. This is the Gaussian multiplier (or wild) bootstrap procedure. Here too, $p$ can be large or even much larger than $n$.
These distributional approximations, either Gaussian or conditional Gaussian, yield a high-quality approximation to the distribution of the original maximum,
often with approximation error decreasing polynomially in the sample size, and hence are of interest in many applications.
We demonstrate how our Gaussian approximations and the multiplier bootstrap can be used for modern high dimensional estimation, multiple hypothesis testing, and adaptive specification testing. All these results contain non-asymptotic bounds on approximation errors.}
\end{abstract}
\begin{keyword}[class=AMS]
\kwd{62E17}
\kwd{62F40}
\end{keyword}
\begin{keyword}
\kwd{Dantzig selector}
\kwd{Slepian}
\kwd{Stein method}
\kwd{maximum of vector sums}
\kwd{high dimensionality}
\kwd{anti-concentration}
\end{keyword}
\end{frontmatter}
\section{Introduction}
Let $x_{1},\dots,x_{n}$ be independent random vectors in $\mathbb{R}^{p}$, with each $x_i$ having coordinates denoted by $x_{ij}$, that is, $x_{i}= (x_{i1},\dots,x_{ip})'$.
Suppose that each $x_{i}$ is centered, namely ${\mathrm{E}}[x_i]=0$, and has a finite covariance matrix ${\mathrm{E}}[x_i x_i']$.
Consider the rescaled sum:
\begin{equation}\label{eq: define X}
X := (X_{1},\dots,X_{p})' := \frac{1}{\sqrt{n}} \sum_{i=1}^n x_{i}.
\end{equation}
Our goal is to obtain a distributional approximation for the statistic $T_0$ defined as the maximum coordinate of vector $X$:
\begin{equation*}
T_0:= \max_{1\leqslant j \leqslant p} X_{j}.
\end{equation*}
The distribution of $T_0$ is of interest in many applications.
When $p$ is fixed, this distribution can be approximated by
the classical Central Limit Theorem (CLT) applied to $X$. However, in modern applications (cf. \cite{BV11}),
$p$ is often comparable or even larger than $n$, and the classical CLT does not apply in such cases.
This paper provides a tractable approximation to the distribution of $T_0$ when $p$ can be large and possibly much larger than $n$.
The \textit{first} main result of the paper is the Gaussian approximation result (GAR), which bounds the Kolmogorov distance between
the distributions of $T_0$ and its Gaussian analog $Z_0$.
Specifically, let $y_{1},\dots,y_{n}$ be independent centered Gaussian random vectors in $\mathbb{R}^{p}$ such that each $y_{i}$ has the same covariance matrix as $x_{i}$: $y_{i} \sim N(0,{\mathrm{E}} [ x_{i} x_{i}' ])$. Consider the rescaled sum of these vectors:
\begin{equation}\label{eq: define Y}
Y := (Y_{1},\dots,Y_{p})' := \frac{1}{\sqrt{n}} \sum_{i=1}^n y_{i}.
\end{equation}
Vector $Y$ is the Gaussian analog of $X$ in the sense of sharing the same mean and covariance matrix, namely
${\mathrm{E}}[X] = {\mathrm{E}}[Y] = 0$ and ${\mathrm{E}}[XX'] = {\mathrm{E}}[YY'] = n^{-1}\sum_{i=1}^n {\mathrm{E}}[x_i x_i'].$
We then define the Gaussian analog $Z_0$ of $T_0$ as the maximum coordinate of vector $Y$:
\begin{equation}\label{eq: define Z}
Z_{0} := \max_{1 \leqslant j \leqslant p} Y_{j}.
\end{equation}
We show that, under suitable moment assumptions, as $n \to \infty$ and possibly $p=p_{n} \to \infty$,
\begin{equation}
\rho:= \sup_{t \in \mathbb{R}} \left| {\mathrm{P}}( T_0 \leqslant t ) - {\mathrm{P}} ( Z_0 \leqslant t ) \right| \leqslant C n^{- c} \to 0, \label{eq: main result}
\end{equation}
where constants $c > 0$ and $C > 0$ are independent of $n$.
Importantly, in (\ref{eq: main result}), $p$ can be large in comparison to $n$ and be as large as $e^{o(n^{c})}$ for some $c>0$. For example, if $x_{ij}$ are uniformly bounded (namely, $|x_{ij}| \leqslant C_{1}$ for some constant $C_{1} > 0$ for all $i$ and $j$) the Kolmogorov distance $\rho$ converges to zero at a polynomial rate whenever $(\log p)^7/n \to 0$ at a polynomial rate. We obtain similar results when $x_{ij}$ are sub-exponential and even non-sub-exponential under suitable moment assumptions. Figure \ref{fig: two deviations} illustrates the result (\ref{eq: main result}) in a non-sub-exponential example, which is motivated by the analysis of the Dantzig selector of \cite{CandesTao2007} in non-Gaussian settings (see Section \ref{sec: Dantzig}).
\begin{figure}
\label{fig: two deviations}
\includegraphics[width=4in]{graph11.eps}
\caption{\footnotesize
P-P plots comparing distributions of $T_0$ and $Z_0$ in the example motivated by the problem of selecting the penalty level of the Dantzig selector. Here $x_{ij}$ are generated as
$x_{ij}=z_{ij}\varepsilon_{i}$ with $\varepsilon_{i} \sim t(4),$ (a $t$-distribution with four degrees of freedom), and $z_{ij}$ are non-stochastic (simulated once using $U[0,1]$ distribution independently across $i$ and $j$). The dashed line is 45$^\circ$. The distributions of $T_0$ and $Z_0$ are close, as (qualitatively) predicted by the GAR derived in the paper. The quality of the Gaussian approximation is particularly good
for the tail probabilities, which is most relevant for practical applications.}
\end{figure}
The proof of the Gaussian approximation result (\ref{eq: main result}) builds on a number of technical tools such as Slepian's smart path interpolation (which is
related to the solution of Stein's partial differential equation; see Appendix \ref{sec: Slepian-Stein note} of the Supplementary Material (SM; \cite{CCK13})), Stein's leave-one-out method, approximation of maxima by the smooth potentials (related to ``free energy" in spin glasses) and using some fine or subtle properties of such approximation, and exponential inequalities for self-normalized sums. See, for example,
\cite{Slepian1962, Stein1981, Dudley1999, ChenGoldsteinShao2011, Talagrand2003, Chatterjee2005b, Rollin2011, Shao2009, Panchenko2013} for introduction and prior uses of some of these tools.
The proof also critically relies on the anti-concentration and comparison bounds of maxima of Gaussian vectors derived in \cite{ChernozhukovChetverikovKato2012c} and restated in this paper as Lemmas \ref{lem: anticoncentration} and \ref{lemma: distances Gaussian to Gaussian}.
Our new Gaussian approximation theorem has the following innovative features.
First, we provide a general result that establishes that maxima of sums of random vectors can be approximated in distribution by the maxima of sums of Gaussian random vectors when $p \gg n$ and especially when $p$ is of order $e^{o(n^{c})}$ for some $c > 0$. The existing techniques can also lead to results of the form (\ref{eq: main result}) when $p=p_{n} \to \infty$, but under much stronger conditions on $p$ requiring $p^{c}/n \to 0$; see Example 17 (Section 10) in \cite{Pollard2002}. Some high-dimensional cases where $p$ can be of order $e^{o(n^{c})}$ can also be handled via Hungarian couplings, extreme value theory or other methods, though special structure is required (for a detailed review, see Section \ref{sec: literature review} of the SM \cite{CCK13}). Second, our Gaussian approximation theorem covers cases where $T_{0}$ does not have a limit distribution as $n \to \infty$ and $p=p_{n} \to \infty$. In some cases, after a suitable normalization, $T_{0}$ could have an extreme value distribution as a limit distribution, but the approximation to an extreme value distribution requires some restrictions on the dependency structure among the coordinates in $x_{i}$. Our result does not limit the dependency structure. We also emphasize that our theorem specifically covers cases where the process $\{\sum_{i=1}^nx_{ij}/\sqrt{n},1\leqslant j\leqslant p\}$ is not asymptotically Donsker (i.e., can't be embedded into a path of an empirical process that is Donsker).
Otherwise, our result would follow from the classical functional central limit theorems for empirical processes, as in \cite{Dudley1999}. Third, the quality of approximation in (\ref{eq: main result}) is of polynomial order in $n$, which is better than the logarithmic in $n$ quality that we could obtain in
some (though not all) applications using the approximation of the distribution of $T_{0}$ by an extreme value distribution (see \cite{Leadbetter1983}).
Note that the result (\ref{eq: main result}) is immediately useful for inference with statistic $T_0$,
even though ${\mathrm{P}} (Z_{0} \leqslant t )$ needs not converge itself to a well-behaved distribution function. Indeed, if the covariance matrix
$n^{-1} \sum_{i=1}^n {\mathrm{E}}[x_{i} x_{i}']$ is known, then $c_{Z_0}(1-\alpha): = (1-\alpha)$-quantile of $Z_0$, can be computed numerically, and we have
\begin{equation}
\label{eq: inference}
|{\mathrm{P}} ( T_0 \leqslant c_{Z_0}(1-\alpha) ) - (1-\alpha)| \leqslant C n^{-c} \to 0.
\end{equation}
The \textit{second} main result of the paper establishes
validity of the multiplier (or Wild) bootstrap for estimating
quantiles of $Z_0$ when the covariance matrix $n^{-1} \sum_{i=1}^n {\mathrm{E}}[x_{i} x_{i}']$
is unknown. Specifically, we define
the Gaussian-symmetrized version $W_0$ of $T_0$ by multiplying $x_{i}$ with i.i.d. standard Gaussian
random variables $e_{1},\dots,e_{n}$:
\begin{equation}
W_0:=\max_{1\leqslant j\leqslant p}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}x_{ij}e_{i}. \label{eq: main quantity}
\end{equation}
We show that the conditional quantiles of $W_0$ given data $( x_{i} )_{i=1}^{n}$ are able to
consistently estimate the quantiles of $Z_0$ and hence those of $T_0$ (where the notion
of consistency used is the one that guarantees asymptotically valid inference). Here the primary factor driving the bootstrap estimation error is the maximum difference between the empirical and population covariance matrices:
\begin{equation*}
\Delta:=\max_{1 \leqslant j, k \leqslant p} \left | \frac{1}{n} \sum_{i=1}^n (x_{ij} x_{ik} - {\mathrm{E}}[x_{ij} x_{ik}] ) \right |,
\end{equation*}
which can converge to zero even when $p$ is much larger than $n$. For example, when $x_{ij}$ are uniformly bounded, the multiplier
bootstrap is valid for inference if $(\log p)^7/n \to 0$. Earlier related results on bootstrap in the ``$p \to \infty$ but $p/n \to 0$'' regime
were obtained in \cite{Mammen1993}; interesting results on inference on the mean vector of high-dimensional random vectors when $p \gg n$ based on concentration inequalities and symmetrization are obtained in \cite{ArlotBlanchardRoquain2010a,ArlotBlanchardRoquain2010b},
albeit the approach and results are quite different from those given here. In particular, in \cite{ArlotBlanchardRoquain2010a}, either Gaussianity or symmetry in distribution is imposed on the data.
The key motivating example of our analysis is the analysis of construction of one-sided or two-sided uniform
confidence band for high-dimensional means under non-Gaussian assumptions. This requires estimation of a high quantile of the maximum of sample means. We give two concrete applications. One application deals with high-dimensional sparse regression model. In this model, \cite{CandesTao2007} and \cite{BickelRitovTsybakov2009} assume Gaussian errors to analyze the Dantzig selector, where the high-dimensional means enter the constraint in the problem. Our results show that Gaussianity is not necessary and the sharp, Gaussian-like, conclusions hold approximately, with just the fourth moment of the
regression errors being bounded. Moreover, our
approximation allows to take into account correlations
among the regressors. This leads to a better choice of the penalty level
and tighter bounds on performance than those that had been available previously.
In another example we apply our results in the multiple hypothesis testing via
the step-down method of \cite{RomanoWolf05}. In the SM \cite{CCK13} we also provide an application to adaptive specification testing. In either case the number of hypotheses to be tested or the number of moment restrictions to be tested can be much larger than the sample size. Lastly, in a companion work (\cite{ChernozhukovChetverikovKato2012b}), we derive the strong coupling for suprema of general empirical processes based on the methods developed here and maximal inequalities. These results represent a useful complement to the results based on the Hungarian coupling developed by \cite{KomlosMajorTusnady1975,BretagnolleMassart1989, Koltchinskii1994,Rio1994} for the entire empirical process and have applications to inference in nonparametric problems such as construction of uniform confidence bands and testing qualitative hypotheses (see, e.g., \cite{GineNickl2010}, \cite{Spokoiny}, and \cite{Chetverikov2012}).
\subsection{Organization of the paper}In Section \ref{sec: Gaus vs NonGaus}, we give the results on Gaussian approximation, and in Section \ref{sec: multiplier bootstrap} on the multiplier bootstrap. In Sections \ref{sec: Dantzig} and \ref{sub: MHT}, we develop applications to the Dantzig selector and multiple testing. Appendices \ref{sec: auxiliary lemmas}-\ref{sec: multiplier bootstrap proofs} contain proofs for each of these sections, with Appendix \ref{sec: auxiliary lemmas} stating auxiliary tools and lemmas. Due to the space limitation, we put additional results and proofs into the SM \cite{CCK13}. In particular, Appendix \ref{sub: AST} of the SM provides additional application to adaptive specification testing. Results of Monte Carlo simulations are presented in Appendix \ref{sec: monte carlo} of the SM.
\subsection{Notation}
In what follows, unless otherwise stated, we will assume that $p \geqslant 3$. In making asymptotic statements, we assume that
$n \to \infty$ with understanding that $p$ depends on $n$ and possibly $p \to \infty$ as $n \to \infty$.
Constants $c,C,c_{1},C_{1},c_{2},C_{2},\dots$ are understood to be independent of $n$.
Throughout the paper, ${\mathbb{E}_n}[\cdot]$ denotes the
average over index $1 \leqslant i \leqslant n$, that is, it simply abbreviates the notation $n^{-1} \sum_{i=1}^n[\cdot]$. For example, ${\mathbb{E}_n}[x_{ij}^{2}]$ $=$ $n^{-1} \sum_{i=1}^{n}x_{ij}^{2}$. In addition, $\bar {\mathrm{E}}[\cdot]={\mathbb{E}_n}[{\mathrm{E}}[\cdot]]$. For example, $\bar {\mathrm{E}}[x_{ij}^2]$ $=$ $n^{-1}\sum_{i=1}^n{\mathrm{E}}[x_{ij}^2]$. For $z\in\mathbb{R}^p$, $z^\prime$ denotes the transpose of $z$.
For a function $f:\mathbb{R} \to \mathbb{R}$, we write $\partial^{k} f (x) = \partial^{k} f(x) / \partial x^{k}$ for nonnegative integer $k$; for a function $f:\mathbb{R}^p \to \mathbb{R}$, we write $\partial_{j} f (x) = \partial f(x)/ \partial x_{j}$ for $j=1,\dots,p$, where $x = (x_{1},\dots,x_{p})'$.
We denote by $C^k(\mathbb{R})$ the class of $k$ times continuously differentiable functions from $\mathbb{R}$ to itself, and denote by $C_{b}^{k}(\mathbb{R})$ the class of all functions $f \in C^{k}(\mathbb{R})$ such that $\sup_{z \in \mathbb{R}} | \partial^{j} f(z) | < \infty$ for $j=0,\dots,k$.
We write $a \lesssim b$ if $a$ is smaller than or equal to $b$ up to a universal positive constant. For $a,b \in \mathbb{R}$, we write $a \vee b = \max \{ a,b \}$.
For two sets $A$ and $B$, $A \ominus B$ denotes their symmetric difference, that is, $A \ominus B = (A \backslash B) \cup (B \backslash A)$.
\section{Gaussian Approximations for Maxima of Non-Gaussian Sums}
\label{sec: Gaus vs NonGaus}
The purpose of this section is to compare and bound the difference
between the expectations and distribution functions of the non-Gaussian to Gaussian maxima:
\begin{equation*}
T_{0} := \max_{1\leqslant j \leqslant p} X_{j} \ \text{ and } \ Z_{0}:=\max_{1\leqslant j \leqslant p} Y_{j},
\end{equation*}
where vector $X$ is defined in equation (\ref{eq: define X}) and
$Y$ in equation (\ref{eq: define Y}). Here and in what follows, without loss of generality, we will assume that $(x_{i})_{i=1}^{n}$ and $(y_{i})_{i=1}^{n}$ are independent.
In order to derive the main result of this section, we shall employ Slepian interpolation, Stein's leave-one-out method, a truncation method combined with self-normalization, as well as some fine properties of the smooth max function (such as ``stability").
(The relative complexity of the approach is justified in Comment \ref{comment: warmup} below.)
The following bounds on moments will be used in stating the bounds in Gaussian approximations:
\begin{equation}\label{define M and S}
\quad M_k:=\max_{1\leqslant j\leqslant p}(\bar {\mathrm{E}}[|x_{ij}|^k])^{1/k}.
\end{equation}
The problem of comparing distributions of maxima is of intrinsic difficulty since the maximum function $z=(z_{1},\dots,z_{p})' \mapsto \max_{1 \leqslant j \leqslant p} z_{j}$ is non-differentiable.
To circumvent the problem, we use a smooth approximation of the maximum function. For $z=(z_{1},\dots,z_{p})' \in \mathbb{R}^{p}$, consider the function:
\begin{equation*}
F_{\beta}(z):=\beta^{-1}\log\left(\sum_{j=1}^{p}\exp(\beta z_{j})\right),
\end{equation*}
where $\beta > 0$ is the smoothing parameter that controls the level of approximation (we call this function the ``smooth max function'').
An elementary calculation shows that for all $z \in \mathbb{R}^{p}$,
\begin{equation}\label{eq: smooth max property}
0 \leqslant F_{\beta}(z)- \max_{1 \leqslant j \leqslant p} z_{j} \leqslant \beta^{-1} \log p.
\end{equation}
This smooth max function arises in the definition of ``free energy" in spin glasses; see, for example, \cite{Talagrand2003}.
Some important properties of this function, such as stability, are derived in the Appendix.
Given a threshold
level $u>0$, we define a truncated version of $x_{ij}$ by
\begin{equation}
\tilde x_{ij} = x_{ij} 1 \left\{ |x_{ij}| \leqslant u (\bar {\mathrm{E}}[x_{ij}^2])^{1/2} \right \} - {\mathrm{E}}\left [x_{ij} 1 \left \{|x_{ij}| \leqslant u (\bar {\mathrm{E}}[x_{ij}^2])^{1/2} \right \} \right ]. \label{def: truncated x}
\end{equation}
Let $\varphi_x(u)$ be the infimum, which is attained, over all numbers $\varphi \geqslant 0$ such that
\begin{equation}
\bar {\mathrm{E}}\left[x_{ij}^21\left\{|x_{ij}|> u(\bar {\mathrm{E}}[x_{ij}^2])^{1/2}\right\}\right] \leqslant \varphi^{2} \bar {\mathrm{E}}[x_{ij}^{2}]. \label{def: truncation varphi}
\end{equation}
Note that the function $\varphi_{x}(u)$ is right-continuous; it measures the impact of truncation
on second moments. Define $u_x(\gamma)$ as the infimum over all numbers $u \geqslant 0$ such that
\begin{equation*}
{\mathrm{P}}\left ( |x_{ij}| \leqslant u (\bar {\mathrm{E}}[x_{ij}^2])^{1/2}, 1 \leqslant i \leqslant n, 1 \leqslant j \leqslant p \right) \geqslant 1-\gamma.
\end{equation*}
Also define $\varphi_y(u)$ and $u_y(\gamma)$ by the corresponding quantities
for the analogue Gaussian case, namely with $(x_i)_{i=1}^n$ replaced by $(y_i)_{i=1}^n$ in the above definitions.
Throughout the paper we use the following quantities:
\begin{equation*}
\varphi(u): = \varphi_x(u) \vee \varphi_y(u), \ \ u(\gamma) := u_x(\gamma) \vee u_y(\gamma).
\end{equation*}
Also, in what follows, for a smooth function $g: \mathbb{R} \to \mathbb{R}$, write
\begin{equation*}
G_{k} : = \sup_{z \in \mathbb{R}} |\partial^{k} g(z)|, \ k \geqslant 0.
\end{equation*}
The following theorem is the main building block toward deriving a result of the form (\ref{eq: main result}).
\begin{theorem}[Comparison of Gaussian to Non-Gaussian Maxima] \label{theorem:comparison non-Gaussian}
Let $\beta>0, u > 0$ and $\gamma \in (0,1)$ be such that $2\sqrt{2}u M_2 \beta/\sqrt{n} \leqslant 1$ and $u \geqslant u(\gamma)$.
Then for every $g \in C_{b}^3(\mathbb{R})$, $|{\mathrm{E}}[g(F_{\beta}(X))- g(F_{\beta}(Y))]| \lesssim D_{n}(g,\beta,u,\gamma)$, so that
\begin{align*}
&|{\mathrm{E}}[g(T_{0})- g(Z_{0})]| \lesssim D_{n}(g,\beta,u,\gamma) + \beta^{-1} G_1 \log p,
\end{align*}
where
\begin{align*}
D_{n}(g,\beta,u,\gamma) &:= n^{-1/2} (G_3 + G_2 \beta + G_1 \beta^2) M_3^3 + (G_2 + \beta G_1) M_2^2 \varphi(u)\\
&\quad + G_1 M_2 \varphi(u) \sqrt{\log (p/\gamma)} + G_0 \gamma.
\end{align*}
\end{theorem}
We will also invoke the following lemma, which is proved in \cite{ChernozhukovChetverikovKato2012c}.
\begin{lemma}[Anti-Concentration]\label{lem: anticoncentration}
(a) Let $Y_{1},\dots,Y_{p}$ be jointly Gaussian random variables with ${\mathrm{E}}[Y_{j}]=0$ and $\sigma_{j}^{2} := {\mathrm{E}} [Y_{j}^{2} ] > 0$ for all $1 \leqslant j \leqslant p$, and let $a_{p} := {\mathrm{E}} [ \max_{1 \leqslant j \leqslant p} (Y_{j}/\sigma_{j}) ]$.
Let $\underline{\sigma} = \min_{1 \leqslant j \leqslant p} \sigma_{j}$ and $\bar{\sigma} = \max_{1 \leqslant j \leqslant p} \sigma_{j}$. Then for every $\varsigma > 0$,
\begin{equation*}
\sup_{z \in \mathbb{R}} {\mathrm{P}} \left( | \max_{1 \leqslant j \leqslant p} Y_{j} - z| \leqslant \varsigma \right) \leqslant C \varsigma \{ a_p + \sqrt{1 \vee \log (\underline{\sigma}/\varsigma)} \},
\end{equation*}
where $C>0$ is a constant depending only on $\underline{\sigma}$ and $\bar{\sigma}$.
When $\sigma_{j}$ are all equal, $\log (\underline{\sigma}/\varsigma)$ on the right side can be replaced by $1$. (b) Furthermore,
the worst case bound is obtained by bounding $a_p$ by $\sqrt{2 \log p}$.
\end{lemma}
By Theorem \ref{theorem:comparison non-Gaussian} and Lemma \ref{lem: anticoncentration}, we can obtain a bound on the Kolmogorov distance, $\rho$, between the distribution functions of $T_{0}$ and $Z_{0}$, which is the main theorem of this section.
\begin{theorem}[\textbf{Main Result 1: Gaussian Approximation}]
\label{cor: Gaussian to nonGaussian KS 2}
Suppose that there are some constants $0 < c_{1} < C_{1}$ such that $c_{1} \leqslant \bar {\mathrm{E}}[x_{ij}^2] \leqslant C_{1}$ for all $1\leqslant j \leqslant p$. Then for every $\gamma \in (0,1)$,
\begin{equation*}
\rho \leqslant
C \left \{ n^{-1/8} (M_{3}^{3/4} \vee M_{4}^{1/2} ) (\log (pn/\gamma))^{7/8} + n^{-1/2} (\log (pn/\gamma))^{3/2} u(\gamma) + \gamma \right \},
\end{equation*}
where $C>0$ is a constant that depends on $c_1$ and $C_1$ only.
\end{theorem}
\begin{remark}[Removing lower bounds on the variance]\label{comment:relaxed}
The condition that $\bar {\mathrm{E}}[x_{ij}^2]\geqslant c_1$ for {\em all} $1\leqslant j\leqslant p$ can not be removed in general.
However, this condition becomes redundant, if there is at least a nontrivial fraction of components $x_{ij}$'s of vector $x_{i}$ with variance bounded away from zero and all pairwise correlations bounded away from 1: for some $J\subset\{1,\dots,p\}$,
$$
|J| \geqslant \nu p, \ \ \bar {\mathrm{E}}[x_{ij}^2]\geqslant c_1, \ \ \frac{|\bar {\mathrm{E}}[x_{ij}x_{ik}] |}{ \sqrt{\bar {\mathrm{E}}[x_{ij}^2] }\sqrt{\bar {\mathrm{E}}[x_{ik}^2]}} \leqslant 1- \nu', \ \ \forall (k, j) \in J \times J: k \neq j,
$$
where $\nu>0$ and $\nu'>0$ are some constants independent of $n$ or $p$. Section \ref{sec: low variance} of the SM \cite{CCK13} contains formal results under this condition. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
In applications, it is useful to have explicit bounds on the upper function $u(\gamma)$. To this end,
let $h: [0,\infty) \to [0,\infty)$ be a {\em Young-Orlicz modulus}, that is, a convex and strictly increasing function with $h(0) = 0$. Denote by $h^{-1}$ the inverse function of $h$. Standard examples include the power function
$h(v) =v^q$ with inverse $h^{-1}(\gamma) = \gamma^{1/q}$ and the exponential function $h(v) = \exp(v) -1$ with inverse
$h^{-1}(\gamma) = \log(\gamma+1)$. These functions describe how many moments the random variables have; for example,
a random variable $\xi$ has finite $q$th moment if ${\mathrm{E}}[|\xi|^q] < \infty$, and is sub-exponential if ${\mathrm{E}}[\exp(|\xi|/C)] < \infty$ for some $C > 0$. We refer to \cite{VW96}, Chapter 2.2, for further details.
\begin{lemma}[Bounds on the upper function $u(\gamma)$]\label{lem: bound on u}
Let $h: [0,\infty) \to [0,\infty)$ be a Young-Orlicz modulus, and let $B > 0$ and $D > 0$ be constants such that $({\mathrm{E}} [ x_{ij}^{2} ])^{1/2} \leqslant B$ for all $1 \leqslant i \leqslant n, 1 \leqslant j \leqslant p$, and $\bar {\mathrm{E}} [ h( \max_{1 \leqslant j \leqslant p} | x_{ij} |/ D )] \leqslant 1$. Then under the condition of Theorem \ref{cor: Gaussian to nonGaussian KS 2},
\begin{equation*}
u (\gamma) \leqslant C \max \{ D h^{-1}(n/\gamma), B \sqrt{\log (pn/\gamma)} \},
\end{equation*}
where $C>0$ is a constant that depends on $c_1$ and $C_1$ only.
\end{lemma}
In applications, parameters $B$ and $D$ (with $M_3$ and $M_4$ as well) are allowed to increase with $n$.
The size of these parameters and the choice of the Young-Orlicz modulus are case-specific.
\subsection{Examples}\label{sub: examples of applications GAR}
The purpose of this subsection is to obtain bounds on $\rho$ for various leading examples frequently encountered in applications.
We are concerned with simple conditions under which $\rho$ decays polynomially in $n$.
Let $c_{1} > 0$ and $C_{1} > 0$ be some constants, and let $B_{n} \geqslant 1$ be a sequence of constants.
We allow for the case where $B_{n} \to \infty$ as $n \to \infty$.
We shall first consider applications where one of the following conditions is satisfied {\em uniformly in} $1 \leqslant i \leqslant n$ and $1 \leqslant j \leqslant p$:
\begin{itemize}
\item[(E.1)] $c_{1} \leqslant \bar {\mathrm{E}}[x^2_{ij}] \leqslant C_{1}$ and $\displaystyle \max_{k =1,2}\bar {\mathrm{E}}[|x_{ij}|^{2+k}/B^k_n] + {\mathrm{E}}[\exp(|x_{ij}|/B_n)] \leqslant 4$;
\item[(E.2)] $c_{1} \leqslant \bar {\mathrm{E}}[x^2_{ij}] \leqslant C_{1}$ and $\displaystyle \max_{k =1,2}\bar {\mathrm{E}}[|x_{ij}|^{2+k}/B^k_n] +{\mathrm{E}}[ (\max_{1\leqslant j \leqslant p} |x_{ij}| / B_{n})^4] \leqslant 4$.
\end{itemize}
\begin{remark}
As a rather special case, Condition (E.1) covers vectors $x_i$ made up from sub-exponential random variables, that is, $$\bar {\mathrm{E}} [ x_{ij}^{2} ] \geqslant c_{1} \text{ and }{\mathrm{E}}[\exp ( |x_{ij}|/C_{1}) ] \leqslant 2$$
(set $B_n = C_1$), which in turn includes, as a special case, vectors $x_i$ made up from sub-Gaussian random variables. Condition (E.1) also covers the case when $|x_{ij}| \leqslant B_n$ for all $i$ and $j$, where $B_n$ may increase with $n$. Condition (E.2) is weaker than (E.1) in that it restricts only the growth of the fourth moments but stronger than (E.1) in that it restricts the growth of $\max_{1\leqslant j\leqslant p}|x_{ij}|$. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
We shall also consider regression applications where one of the following conditions is satisfied {\em uniformly in} $1 \leqslant i \leqslant n$ and $1 \leqslant j \leqslant p$:
\begin{itemize}
\item[(E.3)] \ $x_{ij}= z_{ij} \varepsilon_{ij}$, where $z_{ij}$ are non-stochastic with $|z_{ij}| \leqslant B_n$, ${\mathbb{E}_n}[z_{ij}^2]=1$, and ${\mathrm{E}}[\varepsilon_{ij}]=0$, ${\mathrm{E}}[\varepsilon^2_{ij}] \geqslant c_{1}$, and ${\mathrm{E}}[ \exp (|\varepsilon_{ij}|/C_{1}) ] \leqslant 2$; or
\item[(E.4)] \ $x_{ij}= z_{ij} \varepsilon_{ij}$, where $z_{ij}$ are non-stochastic with $|z_{ij}| \leqslant B_n$, ${\mathbb{E}_n}[z_{ij}^2]=1$, and ${\mathrm{E}}[\varepsilon_{ij}]=0$, ${\mathrm{E}}[\varepsilon^2_{ij}] \geqslant c_{1}$, and
${\mathrm{E}}[ \max_{1 \leqslant j \leqslant p} \varepsilon_{ij}^4] \leqslant C_1$.
\end{itemize}
\begin{remark} Conditions (E.3) and (E.4) cover examples that arise in high-dimensional regression, for example, \cite{CandesTao2007}, which we shall revisit later in the paper.
Typically, $\varepsilon_{ij}$'s are independent of $j$ (i.e., $\varepsilon_{ij} = \varepsilon_{i})$ and hence ${\mathrm{E}}[ \max_{1 \leqslant j \leqslant p} \varepsilon_{ij}^4] \leqslant C_1$ in condition (E.4) reduces to ${\mathrm{E}}[ \varepsilon_{i}^{4} ] \leqslant C_1$.
Interestingly, these examples are also connected to spin glasses, see, for example, \cite{Talagrand2003} and \cite{Panchenko2013} ($z_{ij}$ can be interpreted
as generalized products of ``spins" and $\varepsilon_{i}$ as their random ``interactions"). Note that conditions (E.3) and (E.4) are special cases of conditions (E.1) and (E.2) but we state (E.3) and (E.4) explicitly because these conditions are useful in applications.
\hfill{\tiny \ensuremath{\blacksquare} } \\
\end{remark}
\begin{corollary}[\textbf{Gaussian Approximation in Leading Examples}]\label{cor: central limit theorem}
Suppose that there exist constants $c_2>0$ and $C_2>0$ such that one of the following conditions is satisfied:
(i) (E.1) or (E.3) holds and $B_n^2 (\log (pn))^7/n\leqslant C_2 n^{-c_2}$ or
(ii) (E.2) or (E.4) holds and $B_n^4 (\log (pn))^7/n\leqslant C_2 n^{-c_2}$.
Then there exist constants $c > 0$ and $C>0$ depending only on $c_{1}, C_1, c_{2}$, and $C_{2}$ such that $$\rho \leqslant Cn^{-c}.$$
\end{corollary}
\begin{remark}
This corollary follows relatively directly from Theorem \ref{cor: Gaussian to nonGaussian KS 2} with help of Lemma \ref{lem: bound on u}.
Moreover, from Lemma \ref{lem: bound on u}, it is routine to find other conditions that lead to the conclusion of Corollary \ref{cor: central limit theorem}. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
\begin{remark}[The benefits from the overall proof strategy]
\label{comment: warmup}
We note in Section \ref{sec: elementary GAR} of the SM \cite{CCK13}, that it is possible to derive the following
result by a much simpler proof:
\begin{lemma}[A Simple GAR] Suppose that there are some constants $c_{1} > 0$ and $C_{1} > 0$ such that $c_{1} \leqslant \bar {\mathrm{E}}[x_{ij}^2] \leqslant C_{1}$ for all $1\leqslant j \leqslant p$. Then there exists a constant $C>0$ depending only on $c_{1}$ and $C_{1}$ such that
\begin{equation}\label{eq: warmup}
\sup_{t\in\mathbb{R}}\left|{\mathrm{P}} ( T_{0} \leqslant t ) - {\mathrm{P}} (Z_{0} \leqslant t) \right| \leqslant C (n^{-1} (\log (pn))^7)^{1/8} (\bar {\mathrm{E}} [S^3_i])^{1/4},
\end{equation}
where $S_i := \max_{1\leqslant j\leqslant p} ( |x_{ij}| +|y_{ij}|)$.
\end{lemma} This simple (though apparently new, at this level of generality) result follows from the classical Lindeberg's argument previously given in Chatterjee \cite{Chatterjee2005a} (in the special context of a spin-glass setting like (E.4) with $\epsilon_{ij} = \epsilon_i$) in combination with Lemma \ref{lem: anticoncentration} and standard kernel smoothing of indicator functions. In the SM \cite{CCK13}, we provide the proof using Slepian-Stein methods, which a reader wishing to see a simple exposition (before reading a much more involved proof of the main results) may find helpful. The bound here is only useful in some limited cases, for example, in (E.3) or (E.4) when $B_n^6 (\log (pn))^7/n \to 0$. When $B_n^6 (\log (pn))^7/n \to \infty$, the simple methods fail, requiring a more delicate argument. Note
that in applications $B_n$ typically grows at a fractional power of $n$, see, for example, \cite{ChernozhukovChetverikovKato2012b} and \cite{Chetverikov2011}, and so the limitation is rather major, and was the principal motivation for our whole paper. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
\section{Gaussian Multiplier Bootstrap}
\label{sec: multiplier bootstrap}
\subsection{A Gaussian-to-Gaussian Comparison Lemma}\label{sec: Gaussian vs Gaussian}
The proofs of the main results in this section rely on the following lemma. Let $V$ and $Y$ be centered Gaussian random vectors in $\mathbb{R}^{p}$
with covariance matrices $\Sigma^{V}$ and
$\Sigma^{Y}$, respectively. The following lemma compares the distribution functions
of $\max_{1 \leqslant j \leqslant p}V_{j} \text{and} \max_{1 \leqslant j \leqslant p} Y_{j}$
in terms of $p$ and
\begin{equation*}
\Delta_0:=\max_{1\leqslant j,k\leqslant p}\left|\Sigma^V_{jk}-\Sigma^Y_{jk}\right|.
\end{equation*}
\begin{lemma}[Comparison of Distributions of Gaussian Maxima]\label{lemma: distances Gaussian to Gaussian}
Suppose that there are some constants $0 < c_{1} < C_{1}$ such that $c_{1} \leqslant \Sigma^Y_{jj}\leqslant C_{1}$ for all $1\leqslant j \leqslant p$.
Then there exists a constant $C>0$ depending only on $c_{1}$ and $C_{1}$ such that
\begin{equation*}
\sup_{t\in\mathbb{R}}\left|{\mathrm{P}}\left(\max_{1 \leqslant j \leqslant p} V_{j}\leqslant t\right)-{\mathrm{P}}\left(\max_{1 \leqslant j \leqslant p} Y_{j}\leqslant t\right)\right| \leqslant C\Delta_{0}^{1/3}(1 \vee \log (p/\Delta_0))^{2/3}.
\end{equation*}
\end{lemma}
\begin{remark} The result is derived in \cite{ChernozhukovChetverikovKato2012c}, and extends that of \cite{Chatterjee2005b} who gave an explicit error in Sudakov-Fernique comparison of expectations of maxima of Gaussian random vectors. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
\subsection{Results on Gaussian Multiplier Bootstrap}
Suppose that we have a dataset $(x_{i})_{i=1}^{n}$ consisting of $n$ independent centered random vectors $x_{i}$ in $\mathbb{R}^{p}$.
In this section, we are interested in approximating quantiles of
\begin{equation}
T_0= \max_{1\leqslant j\leqslant p}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}x_{ij} \label{average}
\end{equation}
using the multiplier bootstrap method. Specifically, let $(e_{i})_{i=1}^{n}$ be a sequence of i.i.d. $N(0,1)$ variables independent of $(x_{i})_{i=1}^{n}$, and let
\begin{equation}
W_0= \max_{1\leqslant j\leqslant p}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}x_{ij}e_{i}. \label{average-multiplier}
\end{equation}
Then we define the multiplier bootstrap estimator of the $\alpha$-quantile of $T_0$ as the conditional $\alpha$-quantile of $W_0$ given $(x_{i})_{i=1}^{n}$, that is,
\begin{equation*}
c_{W_0}(\alpha):=\inf\{ t \in \mathbb{R}: {\mathrm{P}}_e(W_0\leqslant t) \geqslant \alpha \},
\end{equation*}
where ${\mathrm{P}}_e$ is the probability measure induced by the multiplier variables $(e_{i})_{i=1}^n$ holding $(x_{i})_{i=1}^{n}$ fixed (i.e., ${\mathrm{P}}_e(W_0\leqslant t) = {\mathrm{P}} (W_0\leqslant t \mid (x_{i})_{i=1}^{n})$).
The multiplier bootstrap theorem below provides a non-asymptotic bound on the bootstrap estimation error.
Before presenting the theorem, we first give a simple useful lemma that is helpful in the proof of the theorem and in power analysis in applications. Define
\begin{equation*}
c_{Z_0}(\alpha):=\inf\{t\in\mathbb{R}:{\mathrm{P}}(Z_0\leqslant t)\geqslant \alpha\},
\end{equation*}
where $Z_0=\max_{1\leqslant j\leqslant p}\sum_{i=1}^ny_{ij}/\sqrt{n}$ and $(y_{i})_{i=1}^n$ is a sequence of independent $N(0,{\mathrm{E}}[x_{i}x_{i}^{\prime}])$ vectors. Recall that $\Delta=\max_{1 \leqslant j, k \leqslant p} \left | {\mathbb{E}_n}[x_{ij}x_{ik}]-\bar {\mathrm{E}}[x_{ij} x_{ik}] \right |$.
\begin{lemma}[Comparison of Quantiles, I]\label{lem: quantile conditional to unconditional}
Suppose that there are some constants $0 < c_{1} < C_{1}$ such that $c_{1} \leqslant \bar {\mathrm{E}}[x_{ij}^2]\leqslant C_{1}$ for all $1\leqslant j\leqslant p$.
Then for every $\alpha \in (0,1)$,
\begin{align*}
&{\mathrm{P}}\big(c_{W_0}(\alpha)\leqslant c_{Z_0}(\alpha+\pi(\vartheta))\big) \geqslant 1- {\mathrm{P}}(\Delta> \vartheta), \\
&{\mathrm{P}}\big(c_{Z_0}(\alpha)\leqslant c_{W_0}(\alpha+\pi(\vartheta))\big) \geqslant 1- {\mathrm{P}}(\Delta> \vartheta),
\end{align*}
where, for $C_{2} > 0$ denoting a constant depending only on $c_{1}$ and $C_{1}$,
$$\pi(\vartheta):=C_2\vartheta^{1/3}(1 \vee \log(p/\vartheta))^{2/3}.$$
\end{lemma}
Recall that $\rho:=\sup_{t\in\mathbb{R}}\left|{\mathrm{P}}(T_0\leqslant t)-{\mathrm{P}}(Z_0\leqslant t)\right|.$
We are now in position to state the first main theorem of this section.
\begin{theorem}[\textbf{Main Result 2: Validity of Multiplier Bootstrap for High-Dimensional Means}]\label{thm: multiplier bootrstrap I}
Suppose that for some constants $0 < c_{1} < C_{1}$, we have $c_{1} \leqslant \bar {\mathrm{E}}[x_{ij}^2]\leqslant C_{1}$ for all $1\leqslant j\leqslant p$.
Then for every $\vartheta>0$,
\begin{equation*}
\rho_{\ominus}:=\sup_{\alpha\in(0,1)}{\mathrm{P}}(\{T_0\leqslant c_{W_0}(\alpha)\}\ominus\{T_0\leqslant c_{Z_0}(\alpha)\})\leqslant 2(\rho+\pi(\vartheta)+{\mathrm{P}}(\Delta>\vartheta)),
\end{equation*}
where $\pi(\cdot)$ is defined in Lemma \ref{lem: quantile conditional to unconditional}. In addition,
\[
\sup_{\alpha\in(0,1)} \left | {\mathrm{P}}(T_0\leqslant c_{W_0}(\alpha))-\alpha \right|\leqslant \rho_{\ominus}+\rho.
\]
\end{theorem}
Theorem \ref{thm: multiplier bootrstrap I} provides a useful result for the case where the statistics
are maxima of exact averages. There are many applications, however, where the relevant statistics
arise as maxima of approximate averages. The following result shows that
the theorem continues to apply if the approximation error of the relevant statistic
by a maximum of an exact average can be suitably controlled. Specifically, suppose
that a statistic of interest, say $T=T(x_{1}\dots,x_{n})$ which may not be of the form (\ref{average}), can be approximated by $T_0$ of the form (\ref{average}), and that
the multiplier bootstrap is performed on a statistic $W=W(x_{1},\dots,x_{n},e_{1},\dots,e_{n})$, which may be different from (\ref{average-multiplier}) but still can be approximated by $W_0$ of the form (\ref{average-multiplier}).
We require the approximation to hold in the following sense: there exist $\zeta_1 \geqslant 0$ and $\zeta_2 \geqslant 0$, depending on $n$ (and typically $\zeta_{1} \to 0, \zeta_{2} \to 0$ as $n \to \infty$),
such that
\begin{align}
&{\mathrm{P}}(|T-T_0|>\zeta_1)<\zeta_2, \label{eq: statistic approximation} \\
&{\mathrm{P}}({\mathrm{P}}_e(|W-W_0|>\zeta_1)> \zeta_2)< \zeta_2. \label{eq: conditional quantiles}
\end{align}
We use the $\alpha$-quantile of $W=W(x_{1},\dots,x_{n},e_{1},\dots,e_{n})$, computed conditional on $(x_{i})_{i=1}^{n}$:
\begin{equation*}
c_{W}(\alpha):=\inf\{t\in\mathbb{R}:{\mathrm{P}}_e(W\leqslant t)\geqslant \alpha\},
\end{equation*}
as an estimate of the $\alpha$-quantile of $T$.
\begin{lemma}[Comparison of Quantiles, II]\label{lem: quantile approximated to exact}
Suppose that condition (\ref{eq: conditional quantiles}) is satisfied. Then for every $\alpha\in(0,1)$,
\begin{align*}
&{\mathrm{P}}(c_{W}(\alpha)\leqslant c_{W_0}(\alpha+\zeta_2)+\zeta_1) \geqslant 1-\zeta_2,\\
&{\mathrm{P}}(c_{W_0}(\alpha)\leqslant c_{W}(\alpha+\zeta_2)+\zeta_1) \geqslant 1-\zeta_2.
\end{align*}
\end{lemma}
The next result provides a bound on the bootstrap estimation error.
\begin{theorem}[\textbf{Main Result 3: Validity of Multiplier Bootstrap for Approximate High-Dimensional Means}]
\label{thm: multiplier bootrstrap II}
Suppose that, for some constants $0 < c_{1} < C_{1}$, we have $c_{1} \leqslant \bar {\mathrm{E}}[x_{ij}^2]\leqslant C_{1}$ for all $1\leqslant j\leqslant p$.
Moreover, suppose that (\ref{eq: statistic approximation}) and (\ref{eq: conditional quantiles}) hold.
Then for every $\vartheta>0$,
\begin{align*}
\rho_\ominus&:=\sup_{\alpha\in(0,1)}{\mathrm{P}}(\{T\leqslant c_{W}(\alpha)\}\ominus\{T_0\leqslant c_{Z_0}(\alpha)\})\\
&\leqslant 2(\rho + \pi(\vartheta) + {\mathrm{P}}(\Delta>\vartheta)) + C_{3}\zeta_1\sqrt{1 \vee \log(p/\zeta_1)} +5\zeta_2,
\end{align*}
where $\pi(\cdot)$ is defined in Lemma \ref{lem: quantile conditional to unconditional}, and $C_{3}>0$ depends only on $c_{1}$ and $C_{1}$. In addition,
$
\sup_{\alpha\in(0,1)}\left|{\mathrm{P}}(T\leqslant c_W(\alpha))-\alpha\right|\leqslant \rho_{\ominus}+\rho.$
\end{theorem}
\begin{remark}[On Empirical and other bootstraps]
In this paper, we focus on the Gaussian multiplier bootstrap (which is a form of wild bootstrap).
This is because other exchangeable bootstrap methods
are asymptotically equivalent to this bootstrap. For example, consider the empirical (or Efron's) bootstrap which approximates the distribution of $T_{0}$ by the conditional distribution of $T_{0}^{*} = \max_{1 \leqslant j \leqslant p} \sum_{i=1}^{n}(x_{ij}^{*}-{\mathbb{E}_n}[x_{ij}])/\sqrt{n}$ where $x_{1},\dots,x_{n}^{*}$ are i.i.d. draws from the empirical distribution of $x_{1},\dots,x_{n}$.
We show in Section \ref{sec: nonparametric bootstrap} of the SM \cite{CCK13}, that the empirical bootstrap is asymptotically equivalent
to the Gaussian multiplier bootstrap, by virtue of Theorem \ref{cor: Gaussian to nonGaussian KS 2} (applied conditionally on the data). The validity of the empirical bootstrap then follows from the validity of the Gaussian multiplier method. The result is demonstrated
under a simplified condition. A detailed analysis of more sophisticated conditions, and the validity of more general exchangeably weighted bootstraps (see \cite{PW93}) in the current setting, will be pursued in future work. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
\subsection{Examples Revisited}
Here we revisit the examples in Section \ref{sub: examples of applications GAR} and see how the multiplier bootstrap works for these leading examples.
Let, as before, $c_{2} > 0$ and $C_{2} > 0$ be some constants, and let $B_{n} \geqslant 1$ be a sequence of constants.
Recall conditions (E.1)-(E.4) in Section \ref{sub: examples of applications GAR}.
The next corollary shows that the multiplier bootstrap is valid with a polynomial rate of accuracy for the significance level under weak conditions.
\begin{corollary}[\textbf{Multiplier Bootstrap in Leading Examples}]
\label{cor: multiplier bootstrap examples}
Suppose that conditions (\ref{eq: statistic approximation}) and (\ref{eq: conditional quantiles}) hold with $\zeta_1 \sqrt{\log p} + \zeta_2 \leqslant C_{2} n^{-c_{2}}$.
Moreover, suppose that one of the following conditions is satisfied:
(i) (E.1) or (E.3) holds and $B_n^2 (\log (pn))^7/n\leqslant C_2 n^{-c_2}$ or
(ii) (E.2) or (E.4) holds and $B_n^4 (\log (pn))^7/n\leqslant C_2 n^{-c_2}$.
Then there exist constants $c > 0$ and $C > 0$ depending only on $c_{1},C_1,c_{2}$, and $C_{2}$ such that
\begin{equation*}
\rho_{\ominus}=\sup_{\alpha\in(0,1)}{\mathrm{P}}(\{T\leqslant c_{W}(\alpha)\}\ominus\{T_0\leqslant c_{Z_0}(\alpha)\}) \leqslant Cn^{-c}.
\end{equation*}
In addition, $\sup_{\alpha \in (0,1)} | {\mathrm{P}}(T\leqslant c_W(\alpha)) - \alpha |\leqslant \rho_{\ominus}+\rho\leqslant Cn^{-c}$.
\end{corollary}
\section{Application: Dantzig Selector in the Non-Gaussian Model}
\label{sec: Dantzig}
The purpose of this section is to demonstrate the case with which the GAR and the multiplier bootstrap theorem given in Corollaries \ref{cor: central limit theorem} and \ref{cor: multiplier bootstrap examples} can be applied in important problems,
dealing with a high-dimensional inference and estimation. We consider
the Dantzig selector previously studied in the path-breaking works of \cite{CandesTao2007}, \cite{BickelRitovTsybakov2009}, \cite{YeZhang2010} in the Gaussian setting
and of \cite{Koltchinskii2009} in a sub-exponential setting. Here we consider the non-Gaussian case, where the errors have only four bounded moments, and derive the performance bounds that are approximately as sharp as in the Gaussian model. We consider both homoscedastic and heteroscedastic models.
\subsection{Homoscedastic case}
Let $(z_{i},y_{i})_{i=1}^{n}$ be a sample of independent observations where $z_{i} \in \mathbb{R}^{p}$ is a non-stochastic vector of regressors. We consider the model
\begin{equation*}
y_{i}=z_{i}^{\prime}\beta+\varepsilon_{i}, \ \ {\mathrm{E}}[\varepsilon_{i}]=0, \ i=1,\dots,n, \ {\mathbb{E}_n}[z_{ij}^2]=1, \ j=1,\dots,p,
\end{equation*}
where $y_{i}$ is a random scalar dependent variable, and the regressors are normalized in such a way that ${\mathbb{E}_n}[z_{ij}^2]=1$.
Here we consider the
homoscedastic case:
\begin{equation*}
{\mathrm{E}}[\varepsilon_{i}^2] = \sigma^2, \ i=1,\dots,n,
\end{equation*}
where $\sigma^{2}$ is assumed to be known (for simplicity). We allow $p$ to be substantially larger than $n$. It is well known that a condition that gives a good performance for the Dantzig selector is that $\beta$ is sparse, namely
$\|\beta\|_0 \leqslant s \ll n$ (although this assumption will not be invoked below explicitly).
The aim is to estimate the vector $\beta$ in some semi-norms of interest:
$\|\cdot\|_I$, where the label $I$ is the name of a norm of interest. For example, given an estimator $\widehat \beta$ the prediction semi-norm for $\delta = \widehat \beta - \beta$ is
\begin{equation*}
\| \delta \|_{\operatorname{pr}} := \sqrt{{\mathbb{E}_n}[ (z_{i}'\delta)^2 ]},
\end{equation*}
or the $j$th component seminorm for $\delta$ is
$\| \delta \|_{\text{jc}} := |\delta_{j}|,$ and so on.
The Dantzig selector is the estimator defined by
\begin{equation}
\label{eq: dantzig estimator}
\widehat\beta \in \arg\min_{b\in\mathbb{R}^p}\Vert b\|_{\ell_1} \ \text{subject to} \ \sqrt{ n} \max_{1\leqslant j\leqslant p}|{\mathbb{E}_n}[z_{ij}(y_{i}-z_{i}^\prime b)]|\leqslant{\lambda},
\end{equation}
where $\| \beta \|_{\ell_{1}} = \sum_{j=1}^{p} | \beta_{j} |$ is the $\ell_{1}$-norm.
An ideal choice of the penalty level $\lambda$ is meant to ensure that
\begin{equation*}
T_0 := \sqrt{n} \max_{1\leqslant j\leqslant p}|{\mathbb{E}_n}[z_{ij}\varepsilon_{i}]|\leqslant{\lambda}
\end{equation*}
with a prescribed confidence level $1-\alpha$ (where $\alpha$ is a number close
to zero.) Hence we would like to set penalty level ${\lambda}$ equal to
\begin{equation*}
c_{T_0} (1-\alpha) := \text{$(1-\alpha)$-quantile of $T_{0}$},
\end{equation*}
(note that $z_{i}$ are treated as fixed). Indeed, this penalty
would take into account the correlation amongst the regressors, thereby adapting the performance
of the estimator to the design condition.
We can approximate this quantity using the
Gaussian approximations derived in Section 2. Specifically, let
\begin{equation*}
Z_0: = \sigma \sqrt{n} \max_{1\leqslant j\leqslant p}|{\mathbb{E}_n}[z_{ij} e_{i}]|,
\end{equation*}
where $e_{i}$ are i.i.d. $N(0,1)$ random variables independent of the data. We then estimate $c_{T_0} (1-\alpha)$ by
\begin{equation*}
c_{Z_0}(1-\alpha):= \text{ $(1-\alpha)$-quantile of $Z_{0}$}.
\end{equation*}
Note that we can calculate $c_{Z_0}(1-\alpha)$ numerically with any specified precision by the simulation.
(In a Gaussian model, design-adaptive penalty level $c_{Z_0}(1-\alpha)$ was proposed in \cite{BelloniChernozhukov2011}, but its extension to non-Gaussian cases was not available up to now).
An alternative choice of the penalty level is given by
\begin{equation*}
c_{0}(1-\alpha) := \sigma \Phi^{-1}(1- \alpha/(2p)),
\end{equation*}
which is the canonical choice; see \cite{CandesTao2007} and \cite{BickelRitovTsybakov2009}.
Note that canonical choice $c_{0}(1-\alpha)$ disregards the correlation amongst the regressors, and is
therefore more conservative than $c_{Z_0}(1-\alpha)$. Indeed, by the union bound, we see that
\begin{equation*}
c_{Z_0}(1-\alpha) \leqslant c_{0}(1-\alpha).
\end{equation*}
Our first result below shows that the \textit{either} of the two penalty choices, $\lambda= c_{Z_0}(1-\alpha)$ or $\lambda = c_0(1-\alpha)$, are
approximately valid under non-Gaussian noise--under the mild moment assumption ${\mathrm{E}}[\varepsilon_{i}^4] \leqslant \text{const}.$ replacing the canonical Gaussian noise assumption. To derive this result we apply our GAR to $T_0$ to establish that the difference between distribution functions of $T_0$ and $Z_0$ approaches zero at polynomial speed.
Indeed $T_0$ can be represented as a maximum of averages, $
T_0 = \max_{1\leqslant k\leqslant 2p} n^{-1/2} \sum_{i=1}^ n \tilde z_{ik}\varepsilon_{i}$, for
$\tilde z_{i} = (z_{i}', -z_{i}')'$ where $z_i^\prime$ denotes the transpose of $z_i$.
To derive the bound on estimation error $\| \delta \|_I$ in a seminorm of interest, we employ the following identifiability factor:
\begin{equation*}
\kappa_I(\beta):=\inf_{\delta \in \mathbb{R}^p} \left \{ \max_{1\leqslant j\leqslant p}\frac{|{\mathbb{E}_n}[z_{ij}(z_{i}^\prime\delta)]|}{\| \delta \|_I }:
\delta\in \mathcal{R}(\beta), \| \delta\|_I \neq 0 \right \},
\end{equation*}
where $\mathcal{R}(\beta):= \{ \delta \in \mathbb{R}^p: \Vert\beta+\delta\|_{\ell_1}\leqslant\Vert\beta\|_{\ell_1}\}$ is the restricted
set; $\kappa_I(\beta)$ is defined as $\infty$ if $\mathcal{R}(\beta) = \{0\}$ (this happens if $\beta =0$). The factors summarize the impact of sparsity of true parameter value $\beta$ and the design on the identifiability of $\beta$ with respect to the norm $\| \cdot \|_I$.
\begin{remark}[A comment on the identifiability factor $\kappa_I(\beta)$]
The identifiability factors $\kappa_I(\beta)$ depend on the true parameter value $\beta$. These factors represent a modest generalization of the cone invertibility factors and sensitivity characteristics defined in \cite{YeZhang2010} and \cite{GautierTsybakov2011}, which are known to be quite general.
The difference is the use of a norm of interest $\|\cdot\|_I$ instead of the $\ell_q$ norms and the use of smaller (non-conic) restricted set $\mathcal{R}(\beta)$ in the definition.
It is useful to note for later comparisons that in the case of prediction norm $\|\cdot \|_{I} = \| \cdot \|_{\operatorname{pr}}$ and under the exact sparsity assumption $\|\beta\|_0 \leqslant s$, we have
\begin{equation}
\kappa_{\operatorname{pr}}(\beta) \geqslant 2^{-1} s^{-1/2} \kappa(s,1), \label{relate to RE}
\end{equation}
where $\kappa(s,1)$ is the restricted eigenvalue defined in \cite{BickelRitovTsybakov2009}. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
Next we state bounds on the estimation error for the Dantzig selector $\widehat \beta^{(0)}$ with
canonical penalty level $\lambda = \lambda^{(0)} := c_0(1-\alpha)$ and the Dantzig selector $\widehat \beta^{(1)}$
with design-adaptive penalty level $\lambda= \lambda^{(1)} := c_{Z_0}(1-\alpha).$
\begin{theorem}[Performance of Dantzig Selector in Non-Gaussian Model]\label{thm: dantzig estimator}
Suppose that there are some constants $c_{1} > 0, C_{1} > 0$ and $\sigma^2 > 0$, and a sequence $B_{n} \geqslant 1$ of constants such that for all $1 \leqslant i \leqslant n$ and $1 \leqslant j \leqslant p$:
(i) $|z_{ij}|\leqslant B_n$; (ii) ${\mathbb{E}_n}[z_{ij}^2]= 1$; (iii) $ {\mathrm{E}}[\varepsilon_{i}^2] =\sigma^2$; (iv) ${\mathrm{E}}[\varepsilon_{i}^{4}]\leqslant C_{1}$; and (v) $B_n^4(\log (pn))^7/n\leqslant C_{1} n^{-c_{1}}$.
Then there exist constants $c > 0$ and $C > 0$ depending only on $c_{1},C_{1}$ and $\sigma^{2}$ such that, with probability at least $ 1- \alpha - C n^{-c}$, for either $k=0$ or $1$,
\begin{equation*}
\Vert\widehat\beta^{(k)} -\beta \Vert_I \leqslant \frac{2 \lambda^{(k)}}{\sqrt{n} \kappa_I(\beta)}.
\end{equation*}
\end{theorem}
The most important feature of this result is that it provides Gaussian-like conclusions (as
explained below) in a model with non-Gaussian noise, having only four bounded moments. However, the probabilistic guarantee
is not $1-\alpha$ as, for example, in \cite{BickelRitovTsybakov2009}, but rather $1- \alpha - C n^{-c}$, which reflects the cost of non-Gaussianity
(along with more stringent side conditions). In what follows
we discuss details of this result. Note that the bound above holds for any semi-norm of interest $\| \cdot \|_I$.
\begin{remark}[Improved Performance from Design-Adaptive Penalty Level] The use of the design-adaptive penalty level
implies a better performance guarantee for $\widehat \beta^{(1)}$ over $\widehat \beta^{(0)}$. Indeed,
we have
\begin{equation*}
\frac{2 c_{Z_0}(1-\alpha)}{\sqrt{n} \kappa_I(\beta)} \leqslant \frac{2 c_0(1-\alpha)}{\sqrt{n} \kappa_I(\beta)}.
\end{equation*}
For example, in some designs, we can have $\sqrt{n} \max_{1 \leqslant j \leqslant p}| {\mathbb{E}_n}[z_{ij} e_{i}]| = O_{{\mathrm{P}}}(1)$,
so that $c_{Z_0}(1-\alpha) = O(1)$, whereas $c_0(1-\alpha) \propto \sqrt{\log p}$. Thus,
the performance guarantee provided by $\widehat \beta^{(1)}$ can be much better than
that of $\widehat \beta^{(0)}$. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
\begin{remark}[Relation to the previous results under Gaussianity] To compare to the previous results obtained for
the Gaussian settings, let us focus on the prediction norm and on estimator $\widehat \beta^{(1)}$ with penalty level $\lambda= c_{Z_0}(1-\alpha)$.
Suppose that the true value $\beta$ is sparse, namely $\|\beta\|_0 \leqslant s$. In this case, with probability at least $1- \alpha - C n^{-c}$,
\begin{equation}\label{eq: boundD}
\Vert\widehat\beta^{(1)} -\beta \Vert_{\operatorname{pr}} \leqslant \frac{2 c_{Z_0}(1-\alpha)}{\sqrt{n} \kappa_{\operatorname{pr}}(\beta)} \leqslant \frac{4 \sqrt{s} c_0(1-\alpha)}{\sqrt{n} \kappa(s,1)} \leqslant \frac{4 \sqrt{s} \sqrt{2 \log ( \alpha/(2p))}}{\sqrt{n} \kappa(s,1)},
\end{equation}
where the last bound is the same as in \cite{BickelRitovTsybakov2009}, Theorem 7.1, obtained for the Gaussian case.
We recover the same (or tighter) upper bound without making the Gaussianity assumption on the errors. However, the probabilistic guarantee is
not $1-\alpha$ as in \cite{BickelRitovTsybakov2009}, but rather $1- \alpha - Cn^{-c}$, which together with side conditions is the cost of non-Gaussianity. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
\begin{remark}[Other refinements] Unrelated to the main theme of this paper, we can see from (\ref{eq: boundD}) that
there is some tightening of the performance bound due to the use of the identifiability factor $\kappa_{\operatorname{pr}}(\beta)$ in place of the restricted eigenvalue $\kappa(s,1)$; for example, if $p=2$ and $s=1$ and the two regressors are identical, then $\kappa_{\operatorname{pr}}(\beta)>0$, whereas $\kappa(1,1) = 0$. There is also some tightening due to the use of $c_{Z_0}(1-\alpha)$ instead of $c_{0}(1-\alpha)$ as penalty level, as mentioned above. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
\subsection{Heteroscedastic case} We consider the same model as above, except
now the assumption on the error becomes
\begin{equation*}
\sigma_{i}^2 := {\mathrm{E}}[\varepsilon_{i}^2] \leqslant \sigma^2, \ \ i=1,\dots,n,
\end{equation*}
that is, $\sigma^{2}$ is the upper bound on the conditional variance, and we assume that this bound is known (for simplicity). As before, ideally we would like to set penalty level ${\lambda}$ equal to
\begin{equation*}
c_{T_0} (1-\alpha) := \text{$(1-\alpha)$-quantile of $T_{0}$},
\end{equation*}
(where $T_0$ is defined above, and we note that $z_{i}$ are treated as fixed).
The GAR applies as before, namely the difference of the distribution functions of $T_0$ and its Gaussian analogue $Z_0$
converges to zero. In this case, the Gaussian analogue can be represented as
\begin{equation*}
Z_0 := \sqrt{n} \max_{1\leqslant j\leqslant p}|{\mathbb{E}_n}[z_{ij} \sigma_{i} e_{i}]|.
\end{equation*}
Unlike in the homoscedastic case, the covariance structure is no longer known, since
$\sigma_{i}$ are unknown and we can no longer calculate the quantiles of $Z_0$. However,
we can estimate them using the following multiplier bootstrap procedure.
First, we estimate the residuals $\widehat \varepsilon_{i} = y_{i} - z_{i}'\widehat \beta^{(0)}$
obtained from a preliminary Dantzig selector $\widehat \beta^{(0)}$ with the conservative penalty level ${\lambda} = \lambda^{(0)} := c_0(1-1/n) := \sigma \Phi^{-1}(1- 1/(2pn))$,
where $\sigma^{2}$ is the upper bound on the error variance assumed to be known.
Let $(e_{i})_{i=1}^{n}$ be a sequence of i.i.d. standard Gaussian random variables, and let
\begin{equation*}
W: = \sqrt{n} \max_{1\leqslant j\leqslant p}|{\mathbb{E}_n}[z_{ij}\widehat\varepsilon_{i}e_{i}]|.
\end{equation*}
Then we estimate $c_{Z_0} (1-\alpha)$ by
\begin{equation*}
c_{W}(1-\alpha):= \text{$(1-\alpha)$-quantile of $W$},
\end{equation*}
defined conditional on data $(z_{i},y_{i})_{i=1}^{n}$. Note that $c_W(1-\alpha)$ can be calculated numerically with any specified precision by the simulation. Then we apply program (\ref{eq: dantzig estimator}) with ${\lambda} = \lambda^{(1)}= c_W(1-\alpha)$ to obtain $\widehat\beta^{(1)}$.
\begin{theorem}[Performance of Dantzig in Non-Gaussian Model with Bootstrap Penalty Level]\label{thm: dantzig estimator 2}
Suppose that there are some constants $c_{1} > 0, C_{1} > 0,\underline{\sigma}^2 > 0$ and $\sigma^2 > 0$, and a sequence $B_{n} \geqslant 1$ of constants such that for all $1 \leqslant i \leqslant n$ and $1 \leqslant j \leqslant p$:
(i) $|z_{ij}|\leqslant B_n$;
(ii) ${\mathbb{E}_n}[z_{ij}^2]= 1$;
(iii) $\underline{\sigma}^2 \leqslant {\mathrm{E}}[\varepsilon_{i}^2] \leqslant \sigma^2$;
(iv) ${\mathrm{E}}[\varepsilon_{i}^{4}]\leqslant C_{1}$;
(v) $B_n^4(\log (pn))^7/n\leqslant C_{1} n^{-c_{1}}$;
and (vi) $ (\log p) B_n c_0(1-1/n)/(\sqrt{n} \kappa_{\operatorname{pr}}(\beta) ) \leqslant C_{1} n^{-c_{1}}$.
Then there exist constants $c > 0$ and $C > 0$ depending only on $c_{1}, C_{1},\underline{\sigma}^2$ and $\sigma^2$ such that,
with probability at least $ 1- \alpha - \nu_{n}$ where $\nu_{n} = Cn^{-c}$, we have
\begin{equation}
\label{eq: bound1}
\Vert\widehat\beta^{(1)} -\beta \Vert_I \leqslant \frac{2 \lambda^{(1)}}{\sqrt{n} \kappa_{I}(\beta)}.
\end{equation}
Moreover, with probability at least $1- \nu_{n}$,
\begin{equation*}
\lambda^{(1)} = c_W(1-\alpha) \leqslant c_{Z_0}(1-\alpha + \nu_n),
\end{equation*}
where $c_{Z_0}(1-a) := \text{ $(1-a)$-quantile of $Z_{0}$}$; where $c_{Z_0}(1-a) \leqslant c_0(1-a)$.
\end{theorem}
\begin{remark}[A Portmanteu Signicance Test] The result above contains a practical test of joint significance of all regressors, that is, a test of the hypothesis that $\beta_0 =0$, with the exact asymptotic size $\alpha$.
\begin{corollary} Under conditions of the either of preceding two theorems, the test, that rejects the null hypothesis $\beta_0= 0$ if $\widehat \beta^{(1)} \neq 0$, has size equal to $\alpha + Cn^{-c}$.
\end{corollary}
To see this note that under the null hypothesis of $\beta_0 =0$, $\beta_0$ satisfies the constraint in (\ref{eq: dantzig estimator}) with probability $(1-\alpha - Cn^{-c})$, by construction of $\lambda$; hence $\|\widehat \beta^{(1)}\| \leqslant \| \beta_0\| = 0$ with exactly this probability. Appendix \ref{sub: AST} of the SM \cite{CCK13} generalizes this to a more general test, which tests $\beta_0 =0$ in the regression model $y_i= d_i'\gamma_0 + x_i'\beta_0 + \varepsilon_i$, where $d_i$'s are a small set of variables, whose coefficients are not known and need to be estimated. The test orthogonalizes each $x_{ij}$ with respect to $d_i$ by partialling out linearly the effect of $d_i$ on $x_{ij}$. The result similar to that in the corollary continues to hold. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
\begin{remark}[Confidence Bands]
Following Gautier and Tsybakov \cite{GautierTsybakov2011}, the bounds given in the preceding theorems can be used for Scheffe-type (simultaneous) inference on all components of $\beta_0$.
\begin{corollary} Under the conditions of either of the two preceding theorems, a $(1-\alpha - Cn^{-c})$-confidence rectangle for $\beta_0$ is given by the region $\times_{j=1}^p I_j$, where
$I_j = [ \widehat \beta^{(1)}_{j} \pm 2 \lambda^{(1)}/(\sqrt{n} \kappa_{\text{jc}}(\beta)].
$\end{corollary}
We note that $\kappa_{\text{jc}}(\beta) = 1$ if ${\mathbb{E}_n}[z_{ij}z_{ik}] =0$ for all $k \neq j$. Therefore,
in the orthogonal model of Donoho and Johnstone, where ${\mathbb{E}_n}[z_{ij}z_{ik}] =0$ for all pairs $j \neq k$, we have that
$\kappa_{\text{jc}}(\beta) = 1$ for all $1 \leqslant j \leqslant p$, so that $I_j = [ \widehat \beta^{(1)}_{j} \pm 2 \lambda^{(1)}/\sqrt{n}]$, which gives a practical simultaneous $(1-\alpha - Cn^{-c})$ confidence rectangle for $\beta$. In non-orthogonal designs,
we can rely on \cite{GautierTsybakov2011}'s tractable linear programming algorithms for computing lower bounds on $\kappa_I(\beta)$ for various norms $I$ of interest; see also \cite{JuditskyNemirovski2011}.
\hfill{\tiny \ensuremath{\blacksquare} } \end{remark}
\begin{remark}[Generalization of Dantzig Selector] There are many interesting applications where the results
given above apply. There are, for example, interesting works by \cite{AlquierHebiri2011} and \cite{FrickMarnitzMunk2012} that consider
related estimators that minimize a convex penalty subject to the multiresolution screening constraints.
In the context of the regression problem studied above, such estimators may be defined as:
\begin{equation*}
\widehat\beta \in \arg\min_{b\in\mathbb{R}^p} J(b) \text{ subject to } \sqrt{ n} \max_{1\leqslant j\leqslant p}|{\mathbb{E}_n}[z_{ij}(y_{i}-z_{i}^\prime b)]|\leqslant{\lambda},
\end{equation*}
where $J$ is a convex penalty, and the constraint is used for multiresolution screening. For example, the Lasso estimator is nested by the above formulation by using $J(b) = \|b \|_{\operatorname{pr}}$, and the previous Dantzig selector by using $J(b) = \|b \|_{\ell_1}$; the estimators can be interpreted as a point in confidence set for $\beta$, which lies closest to zero under $J$-discrepancy (see references cited above for both of these points). Our results on choosing $\lambda$ apply to this class of estimators, and the previous analysis also applies by redefining the identifiability factor $\kappa_I(\beta)$ relative to the new restricted set $\mathcal{R}(\beta):= \{ \delta \in \mathbb{R}^p: J (\beta+\delta)\leqslant J(\beta)\}$; where
$\kappa_I(\beta)$ is defined as $\infty$ if $\mathcal{R}(\beta) = \{0\}$. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}
\section{Application: Multiple Hypothesis Testing via the Stepdown Method}
\label{sub: MHT}
In this section, we study the problem of multiple hypothesis testing
in the framework of multiple means or, more generally, approximate means. The latter possibility allows us
to cover the case of testing multiple coefficients in multiple regressions, which is often required in empirical studies; see, for example, \cite{Anderson08}. We combine a general stepdown
procedure described in \cite{RomanoWolf05} with the multiplier bootstrap developed in this paper.
In contrast with \cite{RomanoWolf05}, our results do not require
weak convergence arguments, and, thus, can be applied to models with an increasing
number of means. Notably, the number of means can be large in comparison with the sample size.
Let $\beta:=(\beta_1,\dots,\beta_p)^\prime\in\mathbb{R}^p$ be a vector of parameters of interest. We are interested in simultaneously testing the set of null hypotheses $H_j:\beta_j\leqslant\beta_{0j}$ against the alternatives $H_j^\prime:\beta_j> \beta_{0j}$ for $j=1,\dots,p$ where $\beta_0:=(\beta_{01},\dots,\beta_{0p})^\prime\in\mathbb{R}^p$. Suppose that the estimator $\widehat{\beta}:=(\widehat{\beta}_1,\dots,\widehat{\beta}_p)^\prime\in\mathbb{R}^p$ is available that has an approximately linear form:
\begin{equation}\label{linearize}
\sqrt{n}(\widehat{\beta}-\beta) = \frac{1}{\sqrt{n}}\sum_{i=1}^nx_{i}+r_n,
\end{equation}
where $x_1,\dots,x_n$ are independent zero-mean random vectors in $\mathbb{R}^p$, the influence functions, and $r_n:=(r_{n1},\dots,r_{np})^\prime\in\mathbb{R}^p$ are linearization errors that are small in the sense required by condition (M) below. Vectors $x_1,\dots,x_n$ need not be directly observable. Instead, some estimators $\widehat{x}_1,\dots,\widehat{x}_n$ of influence functions $x_1,\dots,x_n$ are available, which will be used in the bootstrap simulations.
We refer to this framework as testing multiple approximate means. This framework covers the case of testing multiple means with $r_{n}=0$. More generally, this framework also covers the case of multiple linear and non-linear m-regressions; see, for example, \cite{HeShao00} for explicit conditions giving rise to linearizaton (\ref{linearize}). The detailed exposition of how the case of multiple linear regressions fits into this framework can be found in \cite{ChernozhukovChetverikovKato2012a}. Note also that this framework implicitly covers the case of testing equalities ($H_j:\beta_j=\beta_{0j}$) because equalities can be rewritten as pairs of inequalities.
We are interested in a procedure with the strong control of the family-wise error rate.
In other words, we seek a procedure that would reject at least one true null hypothesis
with probability not greater than $\alpha+o(1)$ uniformly over a large class of data-generating processes and, in particular, uniformly over the
set of true null hypotheses. More formally, let $\Omega$ be a set
of all data generating processes, and $\omega$ be the true process.
Each null hypothesis $H_{j}$ is equivalent to $\omega\in\Omega_{j}$
for some subset $\Omega_{j}$ of $\Omega$.
Let $\mathcal{W}:=\{1,\dots,p\}$ and for $w\subset \mathcal{W}$
denote $\Omega^{w}:=(\cap_{j\in w}\Omega_{j})\cap(\cap_{j\notin w}\Omega_{j}^c)$ where $\Omega_{j}^c:=\Omega\backslash\Omega_{j}$. The strong control of the family-wise error rate means
\begin{equation}
\sup_{w \subset \mathcal{W}}\sup_{\omega\in\Omega^{w}}
{\mathrm{P}}_{\omega}\{\text{reject at least one hypothesis among \ensuremath{H_{j}}, \ensuremath{j\in w}}\}
\leqslant \alpha+o(1) \label{eq: strong control}
\end{equation}
where ${\mathrm{P}}_{\omega}$ denotes the probability distribution under the data-generating process $\omega$.
This setting is clearly of interest in many empirical studies.
For $j=1,\dots,p$, denote $t_{j}:=\sqrt{n}(\widehat{\beta}_j-\beta_{0j})$. The stepdown procedure of \cite{RomanoWolf05} is described as follows. For
a subset $w\subset \mathcal{W}$, let $c_{1-\alpha,w}$ be some estimator of
the $(1-\alpha)$-quantile of $\max_{j\in w}t_{j}$. On the first step, let $w(1)=\mathcal{W}$. Reject
all hypotheses $H_{j}$ satisfying $t_{j}>c_{1-\alpha,w(1)}$.
If no null hypothesis is rejected, then stop. If some $H_{j}$ are
rejected, let $w(2)$ be the set of all null hypotheses that
were not rejected on the first step. On step $l\geqslant2$, let $w(l)\subset \mathcal{W}$
be the subset of null hypotheses that were not rejected up to step
$l$. Reject all hypotheses $H_{j}$, $j\in w(l)$, satisfying $t_{j}>c_{1-\alpha,w(l)}$.
If no null hypothesis is rejected, then stop. If some $H_{j}$ are
rejected, let $w(l+1)$ be the subset of all null hypotheses
among $j\in w(l)$ that were not rejected. Proceed in this way until
the algorithm stops.
Romano and Wolf \cite{RomanoWolf05} proved the following result. Suppose that $c_{1-\alpha,w}$
satisfy
\begin{align}
&c_{1-\alpha,w^{\prime}}\leqslant c_{1-\alpha,w^{\prime\prime}} \quad \text{whenever $w^{\prime}\subset w^{\prime\prime}$},\label{eq: critical value property 1}\\
&\sup_{w\subset \mathcal{W}}\sup_{\omega\in\Omega^{w}}{\mathrm{P}}_\omega \left ( \max_{j\in w}t_{j}>c_{1-\alpha,w} \right ) \leqslant \alpha+o(1),\label{eq: conditions MHT}
\end{align}
then inequality (\ref{eq: strong control}) holds if the stepdown procedure is used. Indeed, let $w$
be the set of true null hypotheses. Suppose that the procedure rejects
at least one of these hypotheses. Let $l$ be the step when the procedure
rejected a true null hypothesis for the first time, and let $H_{j_0}$ be
this hypothesis. Clearly, we have $w(l)\supset w$. So,
\begin{equation*}
\max_{j\in w}t_{j}\geqslant t_{j_0}>c_{1-\alpha,w(l)}\geqslant c_{1-\alpha,w}.
\end{equation*}
Combining this chain of inequalities with (\ref{eq: conditions MHT})
yields (\ref{eq: strong control}).
To obtain suitable $c_{1-\alpha,w}$ that satisfy inequalities (\ref{eq: critical value property 1}) and (\ref{eq: conditions MHT}) above, we can use the multiplier bootstrap method. Let $(e_{i})_{i=1}^n$ be an i.i.d. sequence of $N(0,1)$ random variables that are independent of the data. Let $c_{1-\alpha,w}$ be the conditional
$(1-\alpha)$-quantile of $\max_{j \in w}\sum_{i=1}^n\widehat{x}_{ij}e_i/\sqrt{n}$
given $(\widehat{x}_{i})_{i=1}^n$.
To prove that so defined critical values $c_{1-\alpha,w}$ satisfy inequalities (\ref{eq: critical value property 1}) and (\ref{eq: conditions MHT}), the following two quantities play a key role:
$$
\Delta_1:=\max_{1\leqslant j\leqslant p}|r_{nj}|\text{ and } \Delta_2:=\max_{1\leqslant j\leqslant p}{\mathbb{E}_n}[(\widehat{x}_{ij}-x_{ij})^2].
$$
We will assume the following regularity condition,
\begin{itemize}
\item[(M)] There are positive constants $c_2 $ and $C_2$: (i) ${\mathrm{P}}\left(\sqrt{\log p}\Delta_1>C_2n^{-c_2}\right)$ $<$ $C_2n^{-c_2}$ and (ii) ${\mathrm{P}}\left((\log (pn))^2\Delta_2>C_2n^{-c_2}\right)< C_2n^{-c_2}$. In addition, one of the following conditions is satisfied:
(iii) (E.1) or (E.3) holds and $B_n^2 (\log (pn))^7/n\leqslant C_2 n^{-c_2}$ or
(iv) (E.2) or (E.4) holds and $B_n^4 (\log (pn))^7/n\leqslant C_2 n^{-c_2}$.
\end{itemize}
\begin{theorem}[Strong Control of Family-Wise Error Rate]\label{thm: MHT}
Suppose that (M) is satisfied uniformly over a class of data-generating processes $\Omega$. Then the stepdown procedure with the multiplier bootstrap critical values $c_{1-\alpha,w}$ given above satisfy (\ref{eq: strong control}) for this $\Omega$ with $o(1)$ strengthened to $Cn^{-c}$ for some constants $c>0$ and $C>0$ depending only on $c_1,C_1,c_2$, and $C_2$.
\end{theorem}
\begin{remark}[The case of sample means] Let us consider the simple case of testing multiple means.
In this case, $\beta_j = {\mathrm{E}}[z_{ij}]$ and $\widehat \beta_j = {\mathbb{E}_n}[z_{ij}]$, where $z_i =(z_{ij})_{j=1}^p$ are i.i.d. vectors, so that the influence functions are $x_{ij} = z_{ij} - {\mathrm{E}}[z_{ij}]$, and the remainder is zero, $r_{n}=0$. The influence functions $x_i$ are not directly observable, though easily estimable by demeaning, $\widehat x_{ij} = z_{ij} - {\mathbb{E}_n}[z_{ij}]$ for all $i$ and $j$. It is instructive to see the implications of Theorem \ref{thm: MHT} in this simple setting. Condition (i) of assumption (M) holds trivially in this case. Condition (ii) of assumption (M) follows from Lemma \ref{lem: symmetrization inequality 1} under conditions (iii) or (iv) of assumption (M). Therefore, Theorem \ref{thm: MHT} applies provided that $\underline \sigma^2 \leqslant {\mathrm{E}}[{x}_{ij}^2] \leqslant \bar \sigma^2$, $(\log p)^7 \leqslant C_2n^{1-c_2}$ for arbitrarily small $c_2$ and, for example, either (a) ${\mathrm{E}} [ \exp (|x_{ij}|/C_{1}) ] \leqslant 2$ (condition (E.1)) or (b) ${\mathrm{E}}[\max_{1\leqslant j\leqslant p} x_{ij}^4] \leqslant C_{1}$ (condition (E.2)). Hence, the theorem implies that the Gaussian multiplier bootstrap as described above leads to a testing procedure with the strong control of the family-wise error rate for the multiple hypothesis testing problem of which the {\em logarithm} of the number of hypotheses is nearly of order $n^{1/7}$. Note here that no assumption that limits the dependence between $x_{i1},\dots,x_{ip}$ or the distribution of $x_i$ is made. Previously, \cite{ArlotBlanchardRoquain2010b} proved strong control of the family-wise error rate for the Rademacher multiplier bootstrap with some adjustment factors assuming that $x_i$'s are Gaussian with unknown covariance structure. \hfill{\tiny \ensuremath{\blacksquare} } \end{remark}
\begin{remark}[Relation to Simultaneous Testing]The question on how large $p$ can be was studied in \cite{FanHallYao07} but from a conservative perspective. The motivation there is to know how fast $p$ can grow to maintain the size of the simultaneous test when we calculate critical values (conservatively) ignoring the dependency among $t$-statistics $t_{j}$ and assuming that $t_{j}$ were distributed as, say, $N(0,1)$. This framework is conservative in
that correlation amongst statistics is dealt away by independence, namely by \v{S}id\'{a}k procedures. In contrast, our approach takes into account the correlation amongst statistics and hence is asymptotically exact, that is, asymptotically non-conservative. \hfill{\tiny \ensuremath{\blacksquare} }
\end{remark}