Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors

Gaussian Approximations and Multiplier Bootstrap for Maxima of Sums of High-Dimensional Random Vectors

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:

Our goal is to obtain a distributional approximation for the statistic $T_0$ defined as the maximum coordinate of vector $X$:

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. 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 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:

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$:

We show that, under suitable moment assumptions, as $n \to \infty$ and possibly $p=p_{n} \to \infty$,

where constants $c > 0$ and $C > 0$ are independent of $n$.

Importantly, in ((ref)), $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) illustrates the result ((ref)) in a non-sub-exponential example, which is motivated by the analysis of the Dantzig selector of CandesTao2007 in non-Gaussian settings (see Section (ref)).

The proof of the Gaussian approximation result ((ref)) 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) of the Supplementary Material (SM; 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, 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 ChernozhukovChetverikovKato2012c and restated in this paper as Lemmas (ref) and (ref).

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)) 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 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) of the SM 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 Dudley1999. Third, the quality of approximation in ((ref)) 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 Leadbetter1983).

Note that the result ((ref)) 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

The 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}$:

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:

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 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 ArlotBlanchardRoquain2010a,ArlotBlanchardRoquain2010b, albeit the approach and results are quite different from those given here. In particular, in 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, CandesTao2007 and 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 RomanoWolf05. In the SM 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 (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 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., GineNickl2010, Spokoiny, and Chetverikov2012).

Organization of the paper

In Section (ref), we give the results on Gaussian approximation, and in Section (ref) on the multiplier bootstrap. In Sections (ref) and (ref), we develop applications to the Dantzig selector and multiple testing. Appendices (ref)-(ref) contain proofs for each of these sections, with Appendix (ref) stating auxiliary tools and lemmas. Due to the space limitation, we put additional results and proofs into the SM CCK13. In particular, Appendix (ref) of the SM provides additional application to adaptive specification testing. Results of Monte Carlo simulations are presented in Appendix (ref) of the SM.

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)$.

Gaussian Approximations for Maxima of Non-Gaussian Sums

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:

where vector $X$ is defined in equation ((ref)) and $Y$ in equation ((ref)). 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) below.)

The following bounds on moments will be used in stating the bounds in Gaussian approximations:

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:

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}$,

This smooth max function arises in the definition of “free energy" in spin glasses; see, for example, 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

Let $\varphi_x(u)$ be the infimum, which is attained, over all numbers $\varphi \geqslant 0$ such that

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

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:

Also, in what follows, for a smooth function $g: \mathbb{R} \to \mathbb{R}$, write

The following theorem is the main building block toward deriving a result of the form ((ref)).

We will also invoke the following lemma, which is proved in ChernozhukovChetverikovKato2012c.

By Theorem (ref) and Lemma (ref), 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.

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 VW96, Chapter 2.2, for further details.

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.

Examples

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$:

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$:

Gaussian Multiplier Bootstrap

A Gaussian-to-Gaussian Comparison Lemma

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

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

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

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,

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

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 |$.

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.

Theorem (ref) 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)), can be approximated by $T_0$ of the form ((ref)), 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)) but still can be approximated by $W_0$ of the form ((ref)).

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

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}$:

as an estimate of the $\alpha$-quantile of $T$.

The next result provides a bound on the bootstrap estimation error.

Examples Revisited

Here we revisit the examples in Section (ref) 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). The next corollary shows that the multiplier bootstrap is valid with a polynomial rate of accuracy for the significance level under weak conditions.

Application: Dantzig Selector in the Non-Gaussian Model

The purpose of this section is to demonstrate the case with which the GAR and the multiplier bootstrap theorem given in Corollaries (ref) and (ref) 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 CandesTao2007, BickelRitovTsybakov2009, YeZhang2010 in the Gaussian setting and of 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.

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

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:

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

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

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

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

(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

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

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 BelloniChernozhukov2011, but its extension to non-Gaussian cases was not available up to now).

An alternative choice of the penalty level is given by

which is the canonical choice; see CandesTao2007 and 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

Our first result below shows that the 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:

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$.

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).$

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 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$.

Heteroscedastic case

We consider the same model as above, except now the assumption on the error becomes

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

(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

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

Then we estimate $c_{Z_0} (1-\alpha)$ by

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)) with ${\lambda} = \lambda^{(1)}= c_W(1-\alpha)$ to obtain $\widehat\beta^{(1)}$.

Application: Multiple Hypothesis Testing via the Stepdown Method

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, Anderson08. We combine a general stepdown procedure described in RomanoWolf05 with the multiplier bootstrap developed in this paper. In contrast with 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:

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, HeShao00 for explicit conditions giving rise to linearizaton ((ref)). The detailed exposition of how the case of multiple linear regressions fits into this framework can be found in 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

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 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 RomanoWolf05 proved the following result. Suppose that $c_{1-\alpha,w}$ satisfy

then inequality ((ref)) 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,

Combining this chain of inequalities with ((ref)) yields ((ref)).

To obtain suitable $c_{1-\alpha,w}$ that satisfy inequalities ((ref)) and ((ref)) 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)) and ((ref)), 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,