# Distributed Bootstrap for Simultaneous Inference Under High Dimensionality

Yang Yu<EMAIL_ADDRESS>
Department of Statistics
Purdue University
West Lafayette, IN 47907, USA Shih-Kang Chao<EMAIL_ADDRESS>
Department of Statistics
University of Missouri
Columbia, MO 65211, USA Guang Cheng<EMAIL_ADDRESS>
Department of Statistics
University of California, Los Angeles
Los Angeles, CA 90095, USA Part of this manuscript was completed while Cheng
was at Purdue.

###### Abstract

We propose a distributed bootstrap method for simultaneous inference on high-
dimensional massive data that are stored and processed with many machines. The
method produces an $\ell_{\infty}$-norm confidence region based on a
communication-efficient de-biased lasso, and we propose an efficient cross-
validation approach to tune the method at every iteration. We theoretically
prove a lower bound on the number of communication rounds $\tau_{\min}$ that
warrants the statistical accuracy and efficiency. Furthermore, $\tau_{\min}$
only increases logarithmically with the number of workers and the intrinsic
dimensionality, while nearly invariant to the nominal dimensionality. We test
our theory by extensive simulation studies, and a variable screening task on a
semi-synthetic dataset based on the US Airline On-Time Performance dataset.
The code to reproduce the numerical results is available at GitHub:
https://github.com/skchao74/Distributed-bootstrap.

Keywords: Distributed Learning, High-dimensional Inference, Multiplier
Bootstrap, Simultaneous Inference, De-biased Lasso

##  1 Introduction

Modern massive datasets with enormous sample size and tremendous
dimensionality are usually impossible to be processed with a single machine.
For remedy, a master-worker architecture is often adopted, e.g., Hadoop (Singh
and Kaur, 2014), which operates on a cluster of nodes for data storage and
processing, where the master node also contains a portion of the data; see
Figure 1. An inherent problem of this architecture is that inter-node
communication can be over a thousand times slower than intra-node computation
due to the inter-node communication protocol, which unfortunately always comes
with significant overhead (Lan et al., 2018; Fan et al., 2019a). Hence,
communication efficiency is usually a top concern for algorithm development in
distributed learning.

Figure 1: Master-worker architecture for storing and processing distributed
data.

Classical statistical methods are usually not communication-efficient as some
of them require hundreds or even thousands passes over the entire dataset. In
the last few years, active research has greatly advanced our ability to
perform distributed statistical optimization and inference in, e.g., maximum
likelihood estimation (Zhang et al., 2012; Li et al., 2013; Chen and Xie,
2014; Battey et al., 2018; Jordan et al., 2019; Huang and Huo, 2019; Chen et
al., 2018; Zhu et al., 2020), Lasso (Lee et al., 2017; Wang et al., 2017; Wang
and Zhang, 2017), partially linear models (Zhao et al., 2016), nonstandard
regression (Shi et al., 2018; Banerjee et al., 2019), quantile regression
(Volgushev et al., 2019; Chen et al., 2019), principal component analysis (Fan
et al., 2019b; Chen et al., 2020), just to name a few. However, solutions for
many other problems in the distributed framework, for example the statistical
inference for high-dimensional models, are still elusive.

Simultaneous inference for high-dimensional statistical models has been widely
considered in many applications where datasets can be handled with a
standalone computer (Cai and Sun, 2017), and many recent papers focus on
bootstrap as an effective way to implement simultaneous inference (Dezeure et
al., 2017; Zhang and Cheng, 2017; Belloni et al., 2018, 2019; Yu et al.,
2020a). These existing methods typically use the well-celebrated de-biased
Lasso (van de Geer et al., 2014; Zhang and Zhang, 2014; Javanmard and
Montanari, 2014a, b), where the de-biased score results from the KKT condition
of the Lasso optimization problem. However, de-biased Lasso is not directly
applicable in a distributed computational framework. For one thing, the
implementation of de-biased Lasso requires expensive subroutines such as
nodewise Lasso (van de Geer et al., 2014), which has to be replaced by a more
communication-efficient method. For another, the quality of the de-biased
score, which is essential to the validity of the bootstrap, is generally worse
in a distributed computational framework than that in a centralized
computational framework. In particular, it is heavily biased so the asymptotic
normality fails. However, it can possibly be improved with sufficient rounds
of communication between the master and worker nodes. The bootstrap validity
therefore critically hinges on the interplay between the dimensionality of the
model and the sparsity level, as well as the rounds of communication, the
number of worker nodes and the size of local sample that are specific to the
distributed computational framework.

In this paper, we tackle the challenges discussed above and propose a
communication-efficient simultaneous inference method for high-dimensional
models. The main component at the core of our method is a novel way to improve
the quality of the de-biased score with a carefully selected number of rounds
of communication while relaxing the constraint on the number of machines. Our
method is motivated by Wang et al. (2017), who proposed an iterative procedure
for computing the estimator but no statistical inference was provided. Note
that the de-biased Lasso has been applied by Lee et al. (2017) to obtain a
communication-efficient $\sqrt{N}$-consistent estimator, but their method
restricts the number of worker nodes to be less than the local sample size.
Next, we apply communicate-efficient multiplier bootstrap methods k-grad and
n+k-1-grad, which are originally proposed in Yu et al. (2020b) for low
dimensional models. These bootstrap methods prevent repeatedly refitting the
models and relax the constraint on the number of machines that plague the
methods proposed earlier (Kleiner et al., 2014; Sengupta et al., 2016). A key
challenge in implementation is that cross-validation, which is a popular
method for selecting tuning parameters, usually requires multiple passes of
the entire dataset and is typically inefficient in the distributed
computational framework. We propose a new cross-validation that only requires
the master node for implementation without needing to communicate with the
worker nodes.

Our theoretical study focuses on the explicit lower bounds on the rounds of
communication that warrant the validity of the bootstrap method for high-
dimensional generalized linear models; see Section 3.1 for an overview. In
short, the greater the number of worker nodes and/or the intrinsic
dimensionality, the greater the rounds of communication required for the
bootstrap validity. The bootstrap validity and efficiency are corroborated by
an extensive simulation study.

We further demonstrate the merit of our method on variable screening with a
semi-synthetic dataset, based on the large-scale US Airline On-Time
Performance dataset. By performing a pilot study on an independently sampled
subset of data, we take four key explanatory variables for flight delay, which
correspond to the dummy variables of the four years after the September 11
attacks. On another independently sampled subset of data, we combine the dummy
variables of the four years with artificial high-dimensional spurious
variables to create a design matrix. We perform our method on this artificial
dataset, and find that the relevant variables are correctly identified as the
number of iteration increases. In particular, we visualize the effect of these
four years by confidence intervals.

We go beyond our previous publication Yu et al. (2020b) in two major aspects:
(1) In this paper we focus on high-dimensional models. In particular, the
dimensionality of the model can exceed the sample size in each computing node.
We handle high dimensionality using $\ell_{1}$ penalization, and consider de-
biased Lasso under the distributed computational framework. (2) We tune the
$\ell_{1}$ penalized problem with a carefully designed cross-validation
method, which can be applied under distributed computational framework.

The rest of the paper is organized as follows. In Section 2, we introduce the
problem formulation of distributed high-dimensional simultaneous inference and
present the main bootstrap algorithm as well as the cross-validation algorithm
for hyperparameter tuning. Theoretical guarantees of bootstrap validity for
high-dimensional (generalized) linear models are provided in Section 3.
Section 4 presents simulation results that corroborate our theoretical
findings. Section 5 showcases an application on variable screening for high-
dimensional logistic regression with a big real dataset using our new method.
Finally, Section 6 concludes the paper. Technical details are in Appendices.
The proofs of the theoretical results are in Supplementary Material. The code
to reproduce the numerical results is in GitHub:
https://github.com/skchao74/Distributed-bootstrap.

Notations. We denote the $\ell_{p}$-norm ($p\geq 1$) of any vector
$v=(v_{1},\dots,v_{n})$ by $\|v\|_{p}=(\sum_{i=1}^{n}|v_{i}|^{p})^{1/p}$ and
$\|v\|_{\infty}=\max_{1\leq i\leq n}|v_{i}|$. The induced $p$-norm and the
max-norm of any matrix $M\in\mathbb{R}^{m\times n}$ (with element $M_{ij}$ at
$i$-th row and $j$-th column) are denoted by
$\left|\\!\left|\\!\left|{M}\right|\\!\right|\\!\right|_{p}=\sup_{x\in\mathbb{R}^{n};\|x\|_{p}=1}\|Mx\|_{p}$
and $\left|\\!\left|\\!\left|{M}\right|\\!\right|\\!\right|_{\max}=\max_{1\leq
i\leq m;1\leq j\leq n}|M_{i,j}|$. We write $a\lesssim b$ if $a=O(b)$, and
$a\ll b$ if $a=o(b)$.

##  2 Distributed Bootstrap for High-Dimensional Simultaneous Inference

In this section, we introduce the distributed computational framework and
present a novel bootstrap algorithm for high-dimensional simultaneous
inference under this framework. A communication-efficient cross-validation
method is proposed for tuning.

###  2.1 Distributed Computation Framework

Suppose data $\\{Z_{i}\\}_{i=1}^{N}$ are i.i.d., and $\mathcal{L}(\theta;Z)$
is a twice-differentiable convex loss function arising from a statistical
model, where $\theta=(\theta_{1},\dots,\theta_{d})\in\mathbb{R}^{d}$. Suppose
that the parameter of interest $\theta^{\ast}$ is the minimizer of an expected
loss:

$\theta^{\ast}=\operatorname*{\arg\min}_{\theta\in\mathbb{R}^{d}}\mathcal{L}^{\ast}(\theta),\mbox{
where $\mathcal{L}^{\ast}(\theta):\,=\mathbb{E}_{Z}[\mathcal{L}(\theta;Z)]$}.$

We consider a high-dimensional setting where $d>N$ is possible, and
$\theta^{\ast}$ is sparse, i.e., the support of $\theta^{\ast}$ is small.

We consider a distributed computation framework, in which the entire data are
stored distributedly in $k$ machines, and each machine has data size $n$.
Denote by $\\{Z_{ij}\\}_{i=1,\dots,n;j=1,\dots,k}$ the entire data, where
$Z_{ij}$ is $i$-th datum on the $j$-th machine $\mathcal{M}_{j}$, and $N=nk$.
Without loss of generality, assume that the first machine $\mathcal{M}_{1}$ is
the master node; see Figure 1. Define the local and global loss functions as

$\displaystyle\begin{split}\mbox{global loss:
}\mathcal{L}_{N}(\theta)&=\frac{1}{k}\sum_{j=1}^{k}\mathcal{L}_{j}(\theta),\quad\mbox{where}\\\
\mbox{local loss:
}\mathcal{L}_{j}(\theta)&=\frac{1}{n}\sum_{i=1}^{n}\mathcal{L}(\theta;Z_{ij}),\quad
j=1,\dots,k.\end{split}$ (1)

A great computational overhead occurs when the master and worker nodes
communicate. In order to circumvent the overhead, the rounds of communications
between the master and worker nodes should be minimized, and the algorithms
with reduced communication overheads are “communication-efficient”.

###  2.2 High-Dimensional Simultaneous Inference

In this paper, we focus on the simultaneous confidence region for
$\theta^{\ast}$ in a high-dimensional model, which is one of the effective
ways for variable selection and inference that are immune to the well-known
multiple testing problem. In particular, given an estimator $\widehat{\theta}$
that is $\sqrt{N}$-consistent, simultaneous confidence intervals can be found
with confidence $\alpha$, for large $\alpha\in(0,1)$, by finding the quantile

$\displaystyle c(\alpha)$
$\displaystyle:\,=\inf\\{t\in\mathbb{R}:P(\widehat{T}\leq
t)\geq\alpha\\}\quad\text{where}$ (2) $\displaystyle\widehat{T}$
$\displaystyle:\,=\big{\|}\sqrt{N}\big{(}\widehat{\theta}-\theta^{\ast}\big{)}\big{\|}_{\infty}.$
(3)

where $\widehat{\theta}$ may be computed through the de-biased Lasso (van de
Geer et al., 2014; Zhang and Zhang, 2014; Javanmard and Montanari, 2014a, b):

$\displaystyle\widehat{\theta}=\widehat{\theta}_{Lasso}-\widehat{\Theta}\nabla\mathcal{L}_{N}(\widehat{\theta}_{Lasso}),$
(4)

where

$\widehat{\theta}_{Lasso}=\operatorname*{\arg\min}_{\theta\in\mathbb{R}^{d}}\mathcal{L}_{N}(\theta)+\lambda\|\theta\|_{1}$

is the Lasso estimator with some hyperparameter $\lambda>0$,
$\widehat{\Theta}$ is a surrogate inverse Hessian matrix and
$\mathcal{L}_{N}(\theta)=N^{-1}\sum_{i=1}^{N}\mathcal{L}(\theta;Z_{i})$ is the
empirical loss.

Implementing the simultaneous inference based on $\widehat{\theta}$ and
$\widehat{T}$ in distributed computational framework inevitably faces some
computational challenges. Firstly, computing $\widehat{\theta}$ usually
involves some iterative optimization routines that can accumulate a large
communication overhead without a careful engineering. Next, some bootstrap
methods have been proposed for estimating $c(\alpha)$, e.g., the multiplier
bootstrap (Zhang and Cheng, 2017), but they cannot be straightforwardly
implemented within a distributed computational framework due to excessive
resampling and communication. Even though some communication-efficient
bootstrap methods have been proposed, e.g., Kleiner et al. (2014); Sengupta et
al. (2016); Yu et al. (2020b), they either require a large number of machines
or are inapplicable to high-dimensional models.

Because of the above-mentioned difficulties, inference based on $\widehat{T}$
is inapplicable in the distributed computational framework and is regarded as
an “oracle” in this paper. Our goal is to provide a method that is
communication-efficient while entertaining the same statistical accuracy as
that based on the oracle $\widehat{T}$.

###  2.3 High-Dimensional Distributed Bootstrap

In order to adapt (4) to the distributed computational setting, we first need
to find a good substitute $\widetilde{\theta}$ for $\widehat{\theta}_{Lasso}$
that is communication-efficient, while noting that standard algorithms for
Lasso are not communication-efficient. Fortunately, $\widetilde{\theta}$ can
be computed by the communication-efficient surrogate likelihood (CSL)
algorithm with the $\ell_{1}$-norm regularization (Wang et al., 2017; Jordan
et al., 2019), which iteratively generates a sequence of estimators
$\widetilde{\theta}^{(t)}$ with regularization parameters $\lambda^{(t)}$ at
each iteration $t=0,\dots,\tau-1$. See Remark 1 for model tuning and Lines
2-17 of Algorithm 1 for the exact implementation. Under regularity conditions,
if $t$ is sufficiently large, it is warranted that $\widetilde{\theta}$ is
close to $\widehat{\theta}_{Lasso}$.

Typical algorithms for computing $\widehat{\Theta}$, e.g., the nodewise Lasso
(van de Geer et al., 2014), cannot be extended straightforwardly to the
distributed computational framework due to the same issue of communication
inefficiency. We overcome this by performing the nodewise Lasso using only
$\mathcal{M}_{1}$ without accessing the entire dataset. This simple approach
does not sacrifice accuracy as long as a sufficient amount of communication
brings $\widetilde{\theta}$ sufficiently close to $\theta^{*}$.

Lastly, given the surrogate estimators $\widetilde{\theta}$ for
$\widehat{\theta}_{Lasso}$ and $\widetilde{\Theta}$ for $\widehat{\Theta}$, we
estimate the asymptotic quantile $c(\alpha)$ of $\widehat{T}$ by bootstrapping
$\|\widetilde{\Theta}\sqrt{N}\nabla\mathcal{L}_{N}(\widetilde{\theta})\|_{\infty}$
using the k-grad or n+k-1-grad bootstrap originally proposed by Yu et al.
(2020b) for low-dimensional models. However, the number of communication
rounds between master and worker nodes has to be carefully fine-tuned for
high-dimensional models. In particular, the k-grad algorithm computes

$\displaystyle\overline{W}^{(b)}:\,=\bigg{\|}\underbrace{-\widetilde{\Theta}\frac{1}{\sqrt{k}}\sum_{j=1}^{k}\epsilon_{j}^{(b)}\sqrt{n}(\mathbf{g}_{j}-\bar{\mathbf{g}})}_{=:\overline{A}}\bigg{\|}_{\infty},$
(5)

where $\epsilon_{j}^{(b)}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1)$
independent from the data,
$\mathbf{g}_{j}=\nabla\mathcal{L}_{j}(\widetilde{\theta})$ and
$\bar{\mathbf{g}}=k^{-1}\sum_{j=1}^{k}\mathbf{g}_{j}$. However, it is known
that k-grad does not perform well when $k$ is small (Yu et al., 2020b). The
improved algorithm n+k-1-grad computes

$\displaystyle\begin{split}\widetilde{W}^{(b)}:\,=\bigg{\|}&\underbrace{-\widetilde{\Theta}\frac{1}{\sqrt{n+k-1}}\bigg{(}\sum_{i=1}^{n}\epsilon_{i1}^{(b)}(\mathbf{g}_{i1}-\bar{\mathbf{g}})+\sum_{j=2}^{k}\epsilon_{j}^{(b)}\sqrt{n}(\mathbf{g}_{j}-\bar{\mathbf{g}})\bigg{)}}_{=:\widetilde{A}}\bigg{\|}_{\infty},\end{split}$
(6)

where $\epsilon_{i1}^{(b)}$ and $\epsilon_{j}^{(b)}$ are i.i.d.
$\mathcal{N}(0,1)$ multipliers, and
$\mathbf{g}_{i1}=\nabla\mathcal{L}(\widetilde{\theta};Z_{i1})$ is based on a
single datum $Z_{i1}$ in the master. The key advantage of k-grad or n+k-1-grad
is that once the master has the gradients $\mathbf{g}_{j}$ from the worker
nodes, the quantile of $\\{{\overline{W}}^{(b)}\\}_{b=1}^{B}$ can be computed
in the master node only, without needing to communicate with worker nodes. See
Algorithm 3 in the Appendix for the pseudocode of k-grad and n+k-1-grad.

Algorithm 1 k-grad/n+k-1-grad with de-biased $\ell_{1}$-CSL estimator

1:Require: $\tau\geq 1$ rounds of communication; hyperparameters
$\\{\lambda^{(t)}\\}_{t=0}^{\tau-1}$ , nodewise Lasso procedure
Node$(\cdot,\cdot)$ with hyperparameters $\\{\lambda_{l}\\}_{l=1}^{d}$ (see
Section B)

2:$\widetilde{\theta}^{(0)}\leftarrow\operatorname*{\arg\min}_{\theta}\mathcal{L}_{1}(\theta)+\lambda^{(0)}\|\theta\|_{1}$
at $\mathcal{M}_{1}$

3:Compute $\widetilde{\Theta}$ by running
Node$(\nabla^{2}\mathcal{L}_{1}(\widetilde{\theta}^{(0)}),\\{\lambda_{l}\\}_{l=1}^{d})$
at $\mathcal{M}_{1}$

4:for $t=1,\ldots,\tau$ do

5: Transmit $\widetilde{\theta}^{(t-1)}$ to $\\{\mathcal{M}_{j}\\}_{j=2}^{k}$

6: Compute $\nabla\mathcal{L}_{1}(\widetilde{\theta}^{(t-1)})$ at
$\mathcal{M}_{1}$

7: for $j=2,\ldots,k$ do

8: Compute $\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(t-1)})$ at
$\mathcal{M}_{j}$

9: Transmit $\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(t-1)})$ to
$\mathcal{M}_{1}$

10: end for

11: $\nabla\mathcal{L}_{N}(\widetilde{\theta}^{(t-1)})\leftarrow
k^{-1}\sum_{j=1}^{k}\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(t-1)})$ at
$\mathcal{M}_{1}$

12: if $t<\tau$ then

13:
$\widetilde{\theta}^{(t)}\leftarrow\operatorname*{\arg\min}_{\theta}\mathcal{L}_{1}(\theta)-\theta^{\top}\left(\nabla\mathcal{L}_{1}(\widetilde{\theta}^{(t-1)})-\nabla\mathcal{L}_{N}(\widetilde{\theta}^{(t-1)})\right)+\lambda^{(t)}\|\theta\|_{1}$
at $\mathcal{M}_{1}$

14: else

15:
$\widetilde{\theta}^{(\tau)}\leftarrow\widetilde{\theta}^{(\tau-1)}-\widetilde{\Theta}\nabla\mathcal{L}_{N}(\widetilde{\theta}^{(\tau-1)})$
at $\mathcal{M}_{1}$

16: end if

17:end for

18:Run DistBoots$(\text{`{k-grad}' or
`{n+k-1-grad}'},\widetilde{\theta}=\widetilde{\theta}^{(\tau)},\\{\mathbf{g}_{j}=\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(\tau-1)})\\}_{j=1}^{k},$

19: $\widetilde{\Theta}=\widetilde{\Theta})$ at $\mathcal{M}_{1}$

Algorithm 1 presents the complete statistical inference procedure. There are
two key innovative steps in Algorithm 1 that facilitate the statistical
inference for high dimensional model with a big dataset. First, we introduce
de-biased Lasso in distributed inference, which goes beyond high dimensional
model estimation considered in Jordan et al. (2019); Wang et al. (2017).
Second, we use nodewise Lasso to provide a sparse estimation of the high-
dimensional inverse Hessian matrix instead of the empirical Hessian used in Yu
et al. (2020b).

Algorithm 1 can achieve high computational efficiency due to two reasons.
First, we initialize Algorithm 1 with a warm start. Namely, we warm start with
the Lasso estimator estimated with dataset in the master node, which provides
a good initializer. Second, because the nodewise Lasso is computationally
expensive, we perform it only once at the very beginning and freeze it through
the iterations of the algorithm without updating it.

The number of iterations $\tau$ in Algorithm 1 steers the trade-off between
statistical accuracy and communication efficiency. In particular, a larger
$\tau$ leads to a more accurate coverage of the simultaneous confidence
interval, but it also induces a higher communication cost. Therefore, studying
the minimal $\tau$ that warrants the bootstrap accuracy is crucial, which is
done in Section 3.

###### Remark 1

Two groups of hyperparameters need to be chosen in Algorithm 1:
$\\{\lambda^{(t)}\\}_{t=0}^{\tau-1}$ for regularization in CSL estimation, and
$\\{\lambda_{l}\\}_{l=1}^{d}$ for regularization in nodewise Lasso (see
Algorithm 4). In Section 2.4, we propose a cross-validation method for tuning
$\\{\lambda^{(t)}\\}_{t=0}^{\tau-1}$. As to $\\{\lambda_{l}\\}_{l=1}^{d}$,
while van de Geer et al. (2014) suggests to choose the same value for all
$\lambda_{l}$ by cross-validation, a potentially better way may be to allow
$\lambda_{l}$ to be different across $l$ and select each $\lambda_{l}$ via
cross-validation for the corresponding nodewise Lasso, which is the approach
we take for a distributed variable screening task in Section 5.

###### Remark 2

There exist other options than CSL for $\widetilde{\theta}$ such as the
averaging de-biased estimator (Lee et al., 2017), but an additional round of
communication may be needed to compute the local gradients. More importantly,
their method may be inaccurate when $n<k$.

###  2.4 Communication-Efficient Cross-Validation

We propose a communication-efficient cross-validation method for tuning the
hyperparameters $\\{\lambda^{(t)}\\}_{t=0}^{\tau-1}$ in Algorithm 1. Wang et
al. (2017) proposes to hold out a validation set on each node for selecting
$\lambda^{(t)}$. However, this method requires fitting the model for each
candidate value of $\lambda^{(t)}$, which uses the same communication cost as
the complete CSL estimation procedure.

We propose a communication-efficient $K$-fold cross-validation method that
chooses $\lambda^{(t)}$ for the CSL estimation at every iteration $t$. At
iteration $t$, the master uses the gradients already communicated from the
worker nodes at iteration $t-1$. Hence, the cross-validation needs only the
master node, which circumvents costly communication between the master and the
worker nodes.

Specifically, notice that the surrogate loss (see Line 13 in Algorithm 1) is
constructed using $n$ observations $\mathcal{Z}=\\{Z_{i1}\\}_{i=1}^{n}$ in the
master node and $k-1$ gradients
$\mathcal{G}=\\{\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(t-1)})\\}_{j=2}^{k}$
from the worker nodes. We then create $K$ (approximately) equal-size
partitions to both $\mathcal{Z}$ and $\mathcal{G}$. The objective function for
training is formed using $K-1$ partitions of $\mathcal{Z}$ and $\mathcal{G}$.
In terms of the measure of fit, instead of computing the original likelihood
or loss, we calculate the unregularized surrogate loss using the last
partition of $\mathcal{Z}$ and $\mathcal{G}$, still in the master node. See
Algorithm 2 for the pseudocode.

Algorithm 2 Distributed $K$-fold cross-validation for $t$-step CSL

1:Require: $(t-1)$-step CSL estimate $\widetilde{\theta}^{(t-1)}$, set
$\Lambda$ of candidate values for $\lambda^{(t)}$, partition of master data
$\mathcal{Z}=\bigcup_{q=1}^{K}\mathcal{Z}_{q}$, partition of worker gradients
$\mathcal{G}=\bigcup_{q=1}^{K}\mathcal{G}_{q}$

2:for $q=1,\dots,K$ do

3: $\mathcal{Z}_{train}\leftarrow\bigcup_{r\neq q}\mathcal{Z}_{r}$;
$\mathcal{Z}_{test}\leftarrow\mathcal{Z}_{q}$

4: $\mathcal{G}_{train}\leftarrow\bigcup_{r\neq q}\mathcal{G}_{r}$;
$\mathcal{G}_{test}\leftarrow\mathcal{G}_{q}$

5:
$g_{1,train}\leftarrow\text{Avg}_{Z\in\mathcal{Z}_{train}}\Big{(}\nabla\mathcal{L}(\widetilde{\theta}^{(t-1)};Z)\Big{)}$;
$g_{1,test}\leftarrow\text{Avg}_{Z\in\mathcal{Z}_{test}}\Big{(}\nabla\mathcal{L}(\widetilde{\theta}^{(t-1)};Z)\Big{)}$

6:
$\bar{g}_{train}\leftarrow\text{Avg}_{g\in\\{g_{1,train}\\}\cup\mathcal{G}_{train}}(g)$;
$\bar{g}_{test}\leftarrow\text{Avg}_{g\in\\{g_{1,test}\\}\cup\mathcal{G}_{test}}(g)$

7: for $\lambda\in\Lambda_{t}$ do

8:
$\beta\leftarrow\operatorname*{\arg\min}_{\theta}\text{Avg}_{Z\in\mathcal{Z}_{train}}\big{(}\mathcal{L}(\theta;Z)\big{)}-\theta^{\top}\left(g_{1,train}-\bar{g}_{train}\right)+\lambda\|\theta\|_{1}$

9:
$Loss(\lambda,q)\leftarrow\text{Avg}_{Z\in\mathcal{Z}_{test}}\big{(}\mathcal{L}(\beta;Z)\big{)}-\beta^{\top}\left(g_{1,test}-\bar{g}_{test}\right)$

10: end for

11:end for

12:Return
$\lambda^{(t)}=\operatorname*{\arg\min}_{\lambda\in\Lambda}K^{-1}\sum_{q=1}^{K}Loss(\lambda,q)$

##  3 Theoretical Analysis

Section 3.1 provides an overview of the theoretical results. Sections 3.2 and
3.3 presents the rigorous statements for linear models and generalized linear
models (GLMs) respectively.

###  3.1 An Overview of Theoretical Results

As discussed in Section 2.3, $\tau$ has to be large enough to ensure the
bootstrap accuracy, yet it also induces a great communication cost. Hence, our
main goal is to pin down the minimal number of iterations $\tau_{\min}$
(communication rounds) sufficient for the bootstrap validity in Algorithm 1.
An overview of the theoretical results is provided in Figure 2.

As an overall trend in Figure 2, $\tau_{\min}$ is increasing logarithmically
in $k$ and decreasing in $n$ for both k-grad and n+k-1-grad in (generalized)
linear models; in addition, $\tau_{\min}$ is increasing in $\overline{s}$
logarithmically, where $\overline{s}$ is the maximum of the sparsity of the
true coefficient vector and the inverse population Hessian matrix to be
formally defined later.

By comparing the left and right panels of Figure 2 under a fixed tuple
$(n,k,\overline{s})$, the $\tau_{\min}$ for k-grad is always greater or equal
to that for n+k-1-grad, which indicates a greater communication efficiency of
n+k-1-grad. For very small $k$, n+k-1-grad can still provably work, while
k-grad cannot. Particularly, $\tau_{\min}=1$ can work for certain instances of
n+k-1-grad but is always too small for k-grad.

Regarding the comparison between high-dimensional sparse linear models (top
panels) and GLMs (bottom panels), GLMs typically require a greater $n$ than
sparse linear models, which ensures that the error between
$\widetilde{\theta}^{(t)}$ and $\theta^{\ast}$ decreases in a short transient
phase; see Section C in the Appendix for details.

Figure 2: Illustration of Theorems 3-11. Gray region are where the bootstrap
validity are not warranted by our theory, and the other area is colored blue
with varying lightness according to the lower bound of iteration $\tau$.
$\gamma_{n}=\log_{d}n$, $\gamma_{k}=\log_{d}k$ and
$\gamma_{\bar{s}}=\log_{d}\bar{s}$ are the orders of the local sample size
$n$, number of machines $k$ and the sparsity $\bar{s}$.

###  3.2 Linear Model

Suppose that $N$ i.i.d. observations are generated by a linear model
$y=x^{\top}\theta^{\ast}+e$ with an unknown coefficient vector
$\theta^{\ast}\in\mathbb{R}^{d}$, covariate random vector
$x\in\mathbb{R}^{d}$, and noise $e\in\mathbb{R}$ independent of $x$ with zero
mean and variance of $\sigma^{2}$. We consider the least-squares loss
$\mathcal{L}(\theta;z)=\mathcal{L}(\theta;x,y)=(y-x^{\top}\theta)^{2}/2$.

We impose the following assumptions on the linear model.

  * •

$x$ is sub-Gaussian, i.e.,

$\sup_{\|w\|_{2}\leq
1}\mathbb{E}\big{[}\exp((w^{\top}x)^{2}/L^{2})\big{]}=O(1),$

for some absolute constant $L>0$. Moreover,
$1/\lambda_{\tiny{\min}}(\Sigma)\leq\mu$ for some absolute constant $\mu>0$,
where $\Sigma=\mathbb{E}[xx^{\top}]$.

  * •

$e$ is sub-Gaussian, i.e.,

$\mathbb{E}\big{[}\exp(e^{2}/L^{\prime 2})\big{]}=O(1),$

for some absolute constant $L^{\prime}>0$. Moreover, $\sigma>0$ is an absolute
constant.

  * •

$\theta^{\ast}$ and $\Theta_{l,\cdot}$ are sparse for $l=1,\cdots,d$, where
$\Theta:\,=\Sigma^{-1}=\mathbb{E}[xx^{\top}]^{-1}$. Specifically, we denote by
$S:\,=\\{l:\theta^{\ast}_{l}\neq 0\\}$ the active set of covariates and its
cardinality by $s_{0}:\,=|S|$. Also, we define $s_{l}:\,=|\\{l^{\prime}\neq
l:\Theta_{l,l^{\prime}}\neq 0\\}|$, $s^{*}:\,=\max_{l}s_{l}$, and
$\overline{s}=s_{0}\vee s^{*}$.

Assumption • ‣ 3.2 ensures a restricted eigenvalue condition when
$n\gtrsim\bar{s}\log d$ by Rudelson and Zhou (2013). Under the assumptions, we
first investigate the theoretical property of Algorithm 1, where we apply
k-grad with the de-biased $\ell_{1}$-CSL estimator with $\tau$ communications.
Define

$\displaystyle T$
$\displaystyle:\,=\big{\|}\sqrt{N}\big{(}\widetilde{\theta}^{(\tau)}-\theta^{\ast}\big{)}\big{\|}_{\infty},$
(7)

where $\widetilde{\theta}^{(\tau)}$ is an output of Algorithm 1.

######  Theorem 3 (k-grad, sparse linear model)

Suppose • ‣ 3.2-• ‣ 3.2 hold, and that we run Algorithm 1 with k-grad method
in linear models. Let

$\displaystyle\lambda_{l}\asymp\sqrt{\frac{\log
d}{n}}\quad\text{and}\quad\lambda^{(t)}\asymp\sqrt{\frac{\log
d}{nk}}+\sqrt{\frac{\log d}{n}}\bigg{(}s_{0}\sqrt{\frac{\log
d}{n}}\bigg{)}^{t},$ (8)

for $l=1,\dots,d$ and $t=0,\dots,\tau-1$. Assume $n=d^{\gamma_{n}}$,
$k=d^{\gamma_{k}}$, $\overline{s}=d^{\gamma_{s}}$ for some constants
$\gamma_{n},\gamma_{k},\gamma_{s}>0$. If $\gamma_{n}>3\gamma_{s}$,
$\gamma_{k}>3\gamma_{s}$, and $\tau\geq\tau_{\min}$, where

$\displaystyle\tau_{\min}=1+\left\lfloor\max\left\\{\frac{\gamma_{k}+\gamma_{s}}{\gamma_{n}-2\gamma_{s}},1+\frac{3\gamma_{s}}{\gamma_{n}-2\gamma_{s}}\right\\}\right\rfloor,$

then for $T$ defined in (7), we have

$\displaystyle\sup_{\alpha\in(0,1)}|P(T\leq
c_{\overline{W}}(\alpha))-\alpha|=o(1).$ (9)

where
$c_{\overline{W}}(\alpha):\,=\inf\\{t\in\mathbb{R}:P_{\epsilon}(\overline{W}\leq
t)\geq\alpha\\}$, in which $\overline{W}$ is the k-grad bootstrap statistics
with the same distribution as $\overline{W}^{(b)}$ in (5) and $P_{\epsilon}$
denotes the probability with respect to the randomness from the multipliers.

In addition, (9) also holds if $T$ is replaced by $\widehat{T}$ defined in
(3).

Theorem 3 warrants the bootstrap validity for the simultaneous confidence
intervals produced by Algorithm 1 with the k-grad. Furthermore, it also
suggests that the bootstrap quantile can approximates the quantile of the
oracle statistics $T$; that is, our distributed bootstrap procedure is as
statistically efficient as the oracle centralized method.

Next, we show that the same distributed bootstrap validity and the efficiency
of the k-grad also hold for the n+k-1-grad in Algorithm 1.

######  Theorem 4 (n+k-1-grad, sparse linear model)

Suppose • ‣ 3.2-• ‣ 3.2 hold, and that we run Algorithm 1 with n+k-1-grad
method. Let $\lambda_{l}$ and $\lambda^{(t)}$ be as in (8) for $l=1,\dots,d$
and $t=0,\dots,\tau-1$. Assume $n=d^{\gamma_{n}}$, $k=d^{\gamma_{k}}$,
$\overline{s}=d^{\gamma_{s}}$ for some constants
$\gamma_{n},\gamma_{k},\gamma_{s}>0$. If $\gamma_{n}>3\gamma_{s}$,
$\gamma_{n}+\gamma_{k}>4\gamma_{s}$, and $\tau\geq\tau_{\min}$, where

$\displaystyle\tau_{\min}=1+\left\lfloor\frac{(\gamma_{k}\vee\gamma_{s})+\gamma_{s}}{\gamma_{n}-2\gamma_{s}}\right\rfloor,$

then for $T$ defined in (7), we have

$\displaystyle\sup_{\alpha\in(0,1)}|P(T\leq
c_{\widetilde{W}}(\alpha))-\alpha|=o(1).$ (10)

where

$c_{\widetilde{W}}(\alpha):\,=\inf\\{t\in\mathbb{R}:P_{\epsilon}(\widetilde{W}\leq
t)\geq\alpha\\},$

in which $\widetilde{W}$ is the n+k-1-grad bootstrap statistics with the same
distribution as $\widetilde{W}^{(b)}$ in (6) and $P_{\epsilon}$ denotes the
probability with respect to the randomness from the multipliers.

In addition, (10) also holds if $T$ is replaced by $\widehat{T}$ defined in
(3).

Note by Theorem 2.4 of van de Geer et al. (2014) that $\widehat{T}$ is well
approximated by
$\|\widehat{\Theta}\sqrt{N}\nabla\mathcal{L}_{N}(\theta^{\ast})\|_{\infty}$ ,
which is further approximated by the $\ell_{\infty}$-norm of the oracle score

$A=-\Theta\frac{1}{\sqrt{N}}\sum_{i=1}^{n}\sum_{j=1}^{k}\nabla\mathcal{L}(\theta^{\ast};Z_{ij}),$

given that $\widehat{\Theta}$ only deviates from $\Theta$ up to order
$O_{P}(s^{*}(\log d)^{1/2}N^{-1/2})$ in $\ell_{\infty}$-norm. To gain a deeper
look into the efficiency of k-grad and n+k-1-grad, we compare the difference
between the covariance of $A$ and the conditional covariance of $\overline{A}$
(for k-grad, defined in (5)), and $\widetilde{A}$ (for n+k-1-grad, defined in
(6)). In particular, conditioning on the data $Z_{ij}$, we have

$\displaystyle\left|\\!\left|\\!\left|{\operatorname{cov}_{\epsilon}(\overline{A})-\operatorname{cov}(A)}\right|\\!\right|\\!\right|_{\max}\leq$
$\displaystyle
s^{*}\|\widetilde{\theta}^{(\tau-1)}-\theta^{\ast}\|_{1}+ns^{*}\|\widetilde{\theta}^{(\tau-1)}-\theta^{\ast}\|_{1}^{2}$
$\displaystyle+O_{P}\bigg{(}\sqrt{\frac{{s^{*}}^{2}}{k}}+\sqrt{\frac{s^{*}}{n}}\bigg{)},$
(11)
$\displaystyle\left|\\!\left|\\!\left|{\operatorname{cov}_{\epsilon}(\widetilde{A})-\operatorname{cov}(A)}\right|\\!\right|\\!\right|_{\max}\leq$
$\displaystyle
s^{*}\|\widetilde{\theta}^{(\tau-1)}-\theta^{\ast}\|_{1}+(n\wedge
k)s^{*}\|\widetilde{\theta}^{(\tau-1)}-\theta^{\ast}\|_{1}^{2}$
$\displaystyle+O_{P}\bigg{(}\sqrt{\frac{{s^{*}}^{2}}{n+k}}+\sqrt{\frac{s^{*}}{n}}\bigg{)},$
(12)

up to some logarithmic terms in $d$, $n$ or $k$. Overall, n+k-1-grad in (12)
has a smaller error term than that of k-grad in (11). In particular, k-grad
requires both $n$ and $k$ to be large, while n+k-1-grad requires a large $n$
but not necessarily a large $k$. In addition, $\tau=1$ could be enough for
n+k-1-grad, but not for k-grad. To see it, if
$\|\widetilde{\theta}^{(0)}-\theta^{\ast}\|_{1}$ is of order
$O_{P}(s^{*}/\sqrt{n})$, the right-hand side of (11) can grow with $s^{*}$,
while the error in (12) still shrinks to zero as long as $k\ll n$.

###### Remark 5

Note in both Theorems 3 and 4 that the expression of $\tau_{\min}$ does not
depend on $d$, because the direct effect of $d$ only enters through an
iterative logarithmic term $\log\log d$ which is dominated by
$\log\overline{s}\asymp\log d$.

###### Remark 6

The rates of $\\{\lambda^{(t)}\\}_{t=0}^{\tau-1}$ and
$\\{\lambda_{l}\\}_{l=1}^{d}$ in Theorems 3 and 4 are motivated by those in
Wang et al. (2017) and van de Geer et al. (2014), which, unfortunately, are
not useful in practice. We therefore provide a practically useful cross-
validation method in Section 2.4.

###### Remark 7

The main result (Theorem 2.2) in Zhang and Cheng (2017) can be seen as a
justification of multiplier bootstrap for high-dimensional linear models with
data being processed in a centralized manner. Theorem 4 compliments it by
justifying a distributed multiplier bootstrap with at least one round of
communication ($\tau\geq 1$).

###### Remark 8

A rate of $\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{\overline{W}}(\alpha))-\alpha\right|$ may be shown to be polynomial in $n$
and $k$ with a more careful analysis, which is faster than the order obtained
by the extreme value distribution approach (Chernozhukov et al., 2013; Zhang
and Cheng, 2017) that is at best logarithmic.

###### Remark 9

We have not addressed the question of whether the conditions for $\tau_{\min}$
in Theorem 3 and 4 can be improved in a minimax sense. This is left for future
research. On the other hand, we remark that the total communication cost in
our algorithm is of order $\Omega(\tau_{\min}kd)$, because in each iteration
we communicate $d$-dimensional vectors between the master node and $k-1$
worker nodes, and $\tau_{\min}$ only grows logarithmically with $k$. Our order
matches those in the existing communication-efficient statistical inference
literature e.g., Jordan et al. (2019); Wang et al. (2017).

###  3.3 Generalized Linear Model

In this section, we consider GLMs, which generate i.i.d. observations
$(x,y)\in\mathbb{R}^{d}\times\mathbb{R}$. We assume that the loss function
$\mathcal{L}$ is of the form $\mathcal{L}(\theta;z)=g(y,x^{\top}\theta)$ for
$\theta,x\in\mathbb{R}^{d}$ and $y\in\mathbb{R}$ with
$g:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$, and $g(a,b)$ is three times
differentiable with respect to $b$, and denote $\frac{\partial}{\partial
b}g(a,b)$, $\left(\frac{\partial}{\partial b}\right)^{2}g(a,b)$,
$\left(\frac{\partial}{\partial b}\right)^{3}g(a,b)$ by $g^{\prime}(a,b)$,
$g^{\prime\prime}(a,b)$, $g^{\prime\prime\prime}(a,b)$ respectively. We let
$\theta^{\ast}$ be the unique minimizer of the expected loss
$\mathcal{L}^{\ast}(\theta)$.

We let $X_{1}\in\mathbb{R}^{n\times d}$ be the design matrix in the master
node $\mathcal{M}_{1}$ and $X_{1}^{*}:\,=P^{*}X_{1}$ be the weighted design
matrix with a diagonal $P^{*}\in\mathbb{R}^{n\times n}$ with elements
$\\{g^{\prime\prime}(y_{i1},x_{i1}^{\top}\theta^{\ast})^{1/2}\\}_{i=1,\dots,n}$.
We further let $(X_{1}^{*})_{-l}\varphi^{*}_{l}$ be the $L_{2}$ projection of
$(X_{1}^{*})_{l}$ on $(X_{1}^{*})_{-l}$, for $l=1,\dots,d$. Equivalently, for
$l=1,\dots,d$, we define
$\varphi^{*}_{l}:\,=\operatorname*{\arg\min}_{\varphi\in\mathbb{R}^{d-1}}\mathbb{E}[\|(X_{1}^{*})_{l}-(X_{1}^{*})_{-l}\varphi\|_{2}^{2}]$.

We impose the following assumptions on the GLM.

  * •

For some $\Delta>0$, and $\Delta^{\prime}>0$ such that
$|x^{\top}\theta^{\ast}|\leq\Delta^{\prime}$,

$\displaystyle\sup_{|b|\vee|b^{\prime}|\leq\Delta+\Delta^{\prime}}$
$\displaystyle\sup_{a}\frac{|g^{\prime\prime}(a,b)-g^{\prime\prime}(a,b^{\prime})|}{|b-b^{\prime}|}\leq
1,$ $\displaystyle\max_{|b_{0}|\leq\Delta}$
$\displaystyle\sup_{a}|g^{\prime}(a,b_{0})|=O(1),\quad\text{and}\quad\max_{|b|\leq\Delta+\Delta^{\prime}}\sup_{a}|g^{\prime\prime}(a,b)|=O(1).$

  * •

$\|x\|_{\infty}=O(1)$. Moreover, $x^{\top}\theta^{\ast}=O(1)$ and
$\max_{l}\big{|}g^{\prime\prime}(y,x^{\top}\theta^{\ast})^{1/2}x_{-l}^{\top}\varphi^{*}_{l}\big{|}=O(1)$,
where $x_{-l}$ consists of all but the $l$-th coordinate of $x$.

  * •

The least and the greatest eigenvalues of
$\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})$ and
$\mathbb{E}\left[\nabla\mathcal{L}(\theta^{\ast};Z)\nabla\mathcal{L}(\theta^{\ast};Z)^{\top}\right]$
are bounded away from zero and infinity respectively.

  * •

For some constant $L>0$,

$\max_{l}\max_{q=1,2}\mathbb{E}[|\mathbf{h}_{l}^{2+q}|/L^{q}]+\mathbb{E}[\exp(|\mathbf{h}_{l}|/L)]=O(1),\quad\text{or}$
$\max_{l}\max_{q=1,2}\mathbb{E}[|\mathbf{h}_{l}^{2+q}|/L^{q}]+\mathbb{E}[(\max_{l}|\mathbf{h}_{l}|/L)^{4}]=O(1),$

where
$\mathbf{h}=\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)$
and $\mathbf{h}_{l}$ is the $l$-th coordinate.

  * •

$\theta^{\ast}$ and $\Theta_{l,\cdot}$ are sparse, where the inverse
population Hessian matrix
$\Theta:\,=\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}$, i.e.,
$S:\,=\\{l:\theta^{\ast}_{l}\neq 0\\}$, $s_{0}:\,=|S|$,
$s_{l}:\,=|\\{l^{\prime}\neq l:\Theta_{l,l^{\prime}}\neq 0\\}|$,
$s^{*}:\,=\max_{l}s_{l}$, and $\overline{s}=s_{0}\vee s^{*}$.

Assumption • ‣ 3.3 imposes smoothness conditions on the loss function, which
is satisfied by, for example, the logistic regression. In particular, logistic
regression has $g(a,b)=-ab+\log(1+\exp(b))$, and it can be easily seen that
$|g^{\prime}(a,b)|\leq 2$, $|g^{\prime\prime}(a,b)|\leq 1$,
$|g^{\prime\prime\prime}(a,b)|\leq 1$. Assumption • ‣ 3.3 imposes some
boundedness conditions required for the validity of the nodewise Lasso
(Algorithm 4; van de Geer et al. (2014)) in the master node. Assumption • ‣
3.3 is a standard assumption in the GLM literature. Assumption • ‣ 3.3 is
required for proving the validity of multiplier bootstrap (Chernozhukov et
al., 2013).

Analogously to Theorem 3 and 4 that focus on the distributed bootstrap
validity and the efficiency of Algorithm 1 using k-grad/ n+k-1-grad for linear
models, here we extend them to the high-dimensional de-biased GLMs. See Figure
2 for a comparison between the results of high-dimensional linear models and
GLMs.

######  Theorem 10 (k-grad, sparse GLM)

Suppose • ‣ 3.3-• ‣ 3.3 hold, and that we run Algorithm 1 with k-grad method
in GLMs. Let $\lambda_{l}\asymp\sqrt{\log d/n}$ for $l=1,\dots,d$, and
$\lambda^{(t)}$ be as

$\displaystyle\lambda^{(t)}\asymp\begin{cases}\sqrt{\frac{\log
d}{nk}}+\frac{1}{s_{0}^{2}}\Big{(}s_{0}^{2}\sqrt{\frac{\log
d}{n}}\Big{)}^{2^{t}},&t\leq\tau_{0},\\\ \sqrt{\frac{\log
d}{nk}}+\frac{1}{s_{0}^{2}}\Big{(}s_{0}^{2}\sqrt{\frac{\log
d}{n}}\Big{)}^{2^{\tau_{0}}}\Big{(}s_{0}\sqrt{\frac{\log
d}{n}}\Big{)}^{t-\tau_{0}},&t>\tau_{0}+1,\end{cases}$ (13)

for $t=0,\dots,\tau-1$, where

$\displaystyle\tau_{0}=1+\left\lfloor\log_{2}\frac{\gamma_{n}-2\gamma_{s}}{\gamma_{n}-4\gamma_{s}}\right\rfloor.$
(14)

Assume $n=d^{\gamma_{n}}$, $k=d^{\gamma_{k}}$, $\overline{s}=d^{\gamma_{s}}$
for some constants $\gamma_{n},\gamma_{k},\gamma_{s}>0$. If
$\gamma_{n}>5\gamma_{s}$, $\gamma_{k}>3\gamma_{s}$, and $\tau\geq\tau_{\min}$,
where

$\displaystyle\tau_{\min}=\max\left\\{\tau_{0}+\left\lfloor\frac{\gamma_{k}+\gamma_{s}}{\gamma_{n}-2\gamma_{s}}+\nu_{0}\right\rfloor,2+\left\lfloor\log_{2}\frac{\gamma_{n}-\gamma_{s}}{\gamma_{n}-4\gamma_{s}}\right\rfloor\right\\},$
$\displaystyle\nu_{0}=2-\frac{2^{\tau_{0}}(\gamma_{n}-4\gamma_{s})}{\gamma_{n}-2\gamma_{s}}\in(0,1],$
(15)

then we have (9). In addition, (9) also holds if $T$ is replaced by
$\widehat{T}$ defined in (3).

The $\tau_{0}$ in (14) is the preliminary communication rounds needed for the
CSL estimator to go through the regions which are far from $\theta^{\ast}$. As
$\overline{s}$ grows, the time spent in these regions can increase. However,
when $n$ is large, e.g., $n\gg\overline{s}^{6}$, the loss function is more
well-behaved so that the preliminary communication round can reduce to
$\tau_{0}=1$. See Section C in the Appendix for more details.

######  Theorem 11 (n+k-1-grad, sparse GLM)

Suppose • ‣ 3.3-• ‣ 3.3 hold, and that we run Algorithm 1 with n+k-1-grad
method in GLMs. Let $\lambda_{l}\asymp\sqrt{\log d/n}$ for $l=1,\dots,d$, and
$\lambda^{(t)}$ be as in (13) for $t=0,\dots,\tau-1$. Assume
$n=d^{\gamma_{n}}$, $k=d^{\gamma_{k}}$, $\overline{s}=d^{\gamma_{s}}$ for some
constants $\gamma_{n},\gamma_{k},\gamma_{s}>0$. If $\gamma_{n}>5\gamma_{s}$
and $\tau\geq\tau_{\min}$, where

$\displaystyle\tau_{\min}=\begin{cases}\max\left\\{2+\left\lfloor\log_{2}\frac{\gamma_{k}+\gamma_{s}}{\gamma_{n}-4\gamma_{s}}\right\rfloor,1\right\\},&\text{if}\quad\gamma_{k}\leq\gamma_{n}-3\gamma_{s},\\\
\tau_{0}+\left\lfloor\frac{\gamma_{k}+\gamma_{s}}{\gamma_{n}-2\gamma_{s}}+\nu_{0}\right\rfloor,&\text{otherwise},\end{cases}$

$\tau_{0}$ and $\nu_{0}$ defined as in (14) and (15) respectively, then we
have (10). In addition, (10) also holds if $T$ is replaced by $\widehat{T}$
defined in (3).

###### Remark 12

The selection of $\\{\lambda_{l}\\}_{l=1}^{d}$ in Theorems 10 and 11 are
motivated by those in van de Geer et al. (2014),
$\\{\lambda^{(t)}\\}_{t=0}^{\tau-1}$ are motivated by Wang et al. (2017) and
Jordan et al. (2019). We perform a more careful analysis for the two phases of
model tuning as in (13).

##  4 Simulation Studies

We demonstrate the merits of our methods using synthetic data in this section.
The code to reproduce the simulation experiments, results, and plots is
available at GitHub: https://github.com/skchao74/Distributed-bootstrap.

We consider a Gaussian linear model and a logistic regression model. We fix
total sample size $N=2^{14}$ and the dimension $d=2^{10}$, and choose the
number of machines $k$ from $\\{2^{2},2^{3},\dots,2^{6}\\}$. The true
coefficient $\theta^{\ast}$ is a $d$-dimensional vector in which the first
$s_{0}$ coordinates are 1 and the rest is 0, where $s_{0}\in\\{2^{2},2^{4}\\}$
for the linear model and $s_{0}\in\\{2^{1},2^{3}\\}$ for the GLM. We generate
covariate vector $x$ independently from $\mathcal{N}(0,\Sigma)$, while
considering two different specifications for $\Sigma$:

  * •

Toeplitz: $\Sigma_{l,l^{\prime}}=0.9^{|l-l^{\prime}|}$;

  * •

Equi-correlation: $\Sigma_{l,l^{\prime}}=0.8$ for all $l\neq l^{\prime}$,
$\Sigma_{l,l}=1$ for all $l$.

For linear model, we generate the model noise independently from
$\mathcal{N}(0,1)$; for GLM, we obtain i.i.d. responses from
$y\sim\text{Ber}(1/(1+\exp[-x^{\top}\theta^{\ast}]))$. For each choice of
$s_{0}$ and $k$, we run Algorithm 1 with k-grad and n+k-1-grad on $1{,}000$
independently generated datasets, and compute the empirical coverage
probability and the average width based on the results from these $1{,}000$
replications. At each replication, we draw $B=500$ bootstrap samples, from
which we calculate the $95\%$ empirical quantile to further obtain the $95\%$
simultaneous confidence interval.

For the $\ell_{1}$-CSL computation, we choose the initial $\lambda^{(0)}$ by a
local $K$-fold cross-validation, where $K=10$ for linear regression and $K=5$
for logistic regression. For each iteration $t$, $\lambda^{(t)}$ is selected
by Algorithm 2 in Section 2.4 with $K^{\prime}$ folds with
$K^{\prime}=\min\\{k-1,5\\}$, which ensures that each partition of worker
gradients is non-empty when $k$ is small. For an efficient implementation of
the nodewise Lasso, we select a $\hat{\lambda}$ at every simulation repetition
and set $\lambda_{l}=\bar{\lambda}$ for all $l$. Specifically, for each
simulated dataset, we select
$\bar{\lambda}=10^{-1}\sum_{l=1}^{10}\hat{\lambda}_{l}$, where each
$\hat{\lambda}_{l}$ is obtained obtained by a cross-validation of nodewise
Lasso regression of $l$-th variable on the remaining variables. Since the
variables are homogeneous, these $\hat{\lambda}_{l}$’s only deviate by some
random variations, which can be alleviated by an average.

The computation of the oracle width starts with fixing $(N,d,s_{0})$ and
generating $500$ independent datasets. For each dataset, we compute the
centralized de-biased Lasso estimator $\widehat{\theta}$ as in (4). The oracle
width is defined as two times the $95\%$ empirical quantile of
$\|\widehat{\theta}-\theta^{\ast}\|_{\infty}$ of the 500 samples. The average
widths are compared against the oracle widths by taking the ratio of the two.

The empirical coverage probabilities and the average width ratios of k-grad
and n+k-1-grad are displayed for the linear model in Figures 3 (Toeplitz
design) and 4 (equi-correlation design), and for the logistic regression in
Figures 5 (Toeplitz design) and 6 (equi-correlation design), respectively.
Note that increase in $k$ indicates decrease in $n$, given the fixed $N$.

For small $k$, k-grad tends to over-cover, whereas n+k-1-grad has a more
accurate coverage. By contrast, the coverage of both algorithms fall when $k$
gets too large (or $n$ gets too small), since the estimator
$\widetilde{\theta}^{(\tau)}$ deviates from $\widehat{\theta}$ and the
deviation of the width from the oracle width, which reflects the discussion of
(11) and (12). Moreover, as $s_{0}=\|\theta^{\ast}\|_{0}$ increases, it
becomes harder for both algorithms to achieve the accurate $95\%$ coverage,
and both algorithms start to fail at a smaller $k$ (or larger $n$), which
stems from the fact that the bootstrap cannot accurately approximate variance
of the asymptotic distribution as shown in (11) and (12). Nevertheless,
raising the number of iterations improves the coverage, which verifies our
theory. We also observe an under-coverage of our bootstrap method in both the
linear regression and the logistic regression at the early stage of increasing
$k$. This is due to the loss of accuracy in estimating the inverse Hessian
matrices using only the data in the master node when $k$ increases (or $n$
decreases).

Figure 3: Empirical coverage probability (left axis, solid lines) and average
width (right axis, dashed lines) of simultaneous confidence intervals by
k-grad and n+k-1-grad in sparse linear regression with Toeplitz design and
varying sparsity. Black solid line represents the $95\%$ nominal level and
black dashed line represents 1 on the right $y$-axis.

Figure 4: Empirical coverage probability (left axis, solid lines) and average
width (right axis, dashed lines) of simultaneous confidence intervals by
k-grad and n+k-1-grad in sparse linear regression with equi-correlation design
and varying sparsity. Black solid line represents the $95\%$ nominal level and
black dashed line represents 1 on the right $y$-axis.

Figure 5: Empirical coverage probability (left axis, solid lines) and average
width (right axis, dashed lines) of simultaneous confidence intervals by
k-grad and n+k-1-grad in sparse logistic regression with Toeplitz design and
varying sparsity. Black solid line represents the $95\%$ nominal level and
black dashed line represents 1 on the right $y$-axis.

Figure 6: Empirical coverage probability (left axis, solid lines) and average
width (right axis, dashed lines) of simultaneous confidence intervals by
k-grad and n+k-1-grad in sparse logistic regression with equi-correlation
design and varying sparsity. Black solid line represents the $95\%$ nominal
level and black dashed line represents 1 on the right $y$-axis.

##  5 Variable Screening with Distributed Simultaneous Inference

Having demonstrated the performance of our method on purely synthetic data
using sparse models in the last section, in this section, we artificially
create spurious variables and mix them with the variables obtained from a real
big dataset. We check if our method can successfully select the relevant
variables associated with the response variable from the real dataset. The
code to retrieve data and reproduce the analyses, results, and plots is
available at GitHub: https://github.com/skchao74/Distributed-bootstrap.

###  5.1 Data

The US Airline On-Time Performance dataset (DVN, 2008), available at
http://stat-computing.org/dataexpo/2009, consists of flight arrival and
departure details for all commercial flights within the US from 1987 to 2008.
Given the high dimensionality after dummy transformation and the huge sample
size of the entire dataset, the most efficient way to process the data is
using a distributed computational system, with sample size on each worker node
likely to be smaller than the dimension. Our goal here is to uncover
statistically significant independent variables associated with flight delay.
We use variables Year, Month, DayOfWeek, CRSDepTime, CRSArrTime,
UniqueCarrier, Origin, Dest, and ArrDelay in our model; descriptions are
deferred to Appendix (Section D).

The response variable is labeled by $1$ to denote a delay if ArrDelay is
greater than zero, and by $0$ otherwise. The rest of the variables are treated
as categorical explanatory variables and are converted into dummy variables;
refer to Appendix (Section E) for the details of the dummy variable creation.
This results in a total of $203$ predictors. The total sample size is 113.9
million observations. We randomly sample a dataset $\mathcal{D}_{1}$ of
$N=500{,}000$ observations, and conceptually distribute them across
$k=1{,}000$ nodes such that each node receives $n=500$ observations. We
randomly sample another dataset $\mathcal{D}_{2}$ of $N=500{,}000$
observations for a pilot study to select relevant variables, where
$\mathcal{D}_{1}\cap\mathcal{D}_{2}=\emptyset$.

###  5.2 An Artificial Design Matrix and Variable Screening

In the first stage, we perform a preliminary study that informs us some
seemingly relevant variables to include in an artificial design matrix, which
will be used to demonstrate variable screening performance of our method in
the second stage. Note that the purpose of this stage is only to preliminarily
discover possibly relevant variables, rather than to select variables in a
fully rigorous manner. We perform a logistic regression in a centralized
manner with intercept and without regularization using the $N$ observations in
$\mathcal{D}_{2}$. Standard Wald tests reveal that $144$ out of $203$ slopes
are significantly non-zero ($p$-values less than $0.05$).

The four predictors with the least $p$-values correspond to the dummy
variables of years 2001–2004, and the coefficients are all negative, which
suggests less likelihood of flight delay in these years. This interesting
finding matches the results of previous study that the September 11 terrorist
attacks have negatively impacted the US airline demand (Ito and Lee, 2005),
which led to less flights and congestion. In addition, the Notice of Market-
based Actions to Relieve Airport Congestion and Delay, (Docket No.
OST-2001-9849) issued by Department of Transportation on August 21, 2001,
might also alleviate the US airline delay.

To construct the artificial design matrix, we group the $4$ predictors with
the least $p$-values mentioned above and the intercept, so the number of the
relevant columns is $5$. Given $d$, we artificially create $d-5$ columns of
binary and real valued variables by first sampling rows from
$\mathcal{N}(0,\mathcal{C}_{d-5})$, where $\mathcal{C}_{d-5}$ is a Toeplitz
matrix ($(\mathcal{C}_{d-5})_{l,l^{\prime}}=0.5^{|l-l^{\prime}|}$), and then
converting half of the columns to either 0 or 1 by their signs. Then, we
combine these $d-5$ spurious columns with a column of intercept and the $4$
columns in $\mathcal{D}_{1}$ that are associated with the selected relevant
variables to obtain an artificial design matrix.

In the second stage, using the artificial design matrix with the binary
response vector from the ArrDelay in $\mathcal{D}_{1}$, we test if our
distributed bootstrap n+k-1-grad (Algorithm 1) can screen the artificially
created spurious variables. Note that $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$
are disjoint, where $\mathcal{D}_{2}$ is used in the first stage for the
preliminary study. For model tuning, we select $\lambda^{(0)}$ by a local
$10$-fold cross-validation; for each $t\geq 1$, $\lambda^{(t)}$ is chosen by
running a distributed $10$-fold cross-validation in Algorithm 2. We select
each $\lambda_{l}$ by performing a $10$-fold cross-validation for the nodewise
Lasso of each variable. The same entire procedure is repeated under each
dimensionality $d\in\\{200,500,1{,}000\\}$.

The left panel of Figure 7 plots the number of significant variables against
the number of iterations $\tau$, which was broken down into the number
intersecting with the relevant variables (solid lines) and the number
intersecting with the spurious variables (dashed lines). First, all of the $4$
relevant variables are tested to be significant at all iterations. For the
spurious variables, we see that with $\tau=1$, the distributed bootstrap
falsely detects one of them. However, as the number of iterations increases,
less spurious variables are detected until none of them is detected. We also
see that $2$ iterations ($\tau=2$) for $d=500,1{,}000$ and $3$ iterations
($\tau=3$) for $d=200$ are sufficient, which empirically verifies that our
method is not very sensitive to the nominal dimension $d$.

As an illustration that is potentially useful in practice, the confidence
intervals computed with the simultaneous quantile for the $4$ important slopes
under $d=1{,}000$ and $\tau=2$ are plotted in the right panel of Figure 7. It
can be seen that the flights in years 2002 and 2003 are relatively less likely
to delay, which match the decreased air traffic in the aftermath of the
September 11 terrorist attacks.

Figure 7: The left panel shows the number of significant variables uncovered
by the simultaneous confidence intervals among the $4$ relevant variables and
among the $d-5$ spurious variables for $d=200,500,1{,}000$. The right panel
shows the simultaneous confidence intervals of the $4$ relevant variables for
$d=1{,}000$ and $\tau=2$.

##  6 Conclusion

We propose a distributed bootstrap method for high-dimensional simultaneous
inference based on the de-biased $\ell_{1}$-CSL estimator as well as a
distributed cross-validation method for hyperparameter tuning. The bootstrap
validity and oracle efficiency are rigorously studied, and the merits are
further shown via simulation study on coverage probability and efficiency, and
a practical example on variable screening.

Acknowledgments

Shih-Kang Chao would like to acknowledge the financial support from the
Research Council of the University of Missouri. Guang Cheng would like to
acknowledge support from the National Science Foundation (NSF – SCALE MoDL
(2134209)).

##  Appendix A Pseudocode for k-grad and n+k-1-grad

Algorithm 3
DistBoots$(\text{method},\widetilde{\theta},\\{\mathbf{g}_{j}\\}_{j=1}^{k},\widetilde{\Theta})$:
only need the master node $\mathcal{M}_{1}$

1:Require: local gradient $\mathbf{g}_{j}$ and estimate $\widetilde{\Theta}$
of inverse Hessian obtained at $\mathcal{M}_{1}$

2:$\bar{\mathbf{g}}\leftarrow k^{-1}\sum_{j=1}^{k}\mathbf{g}_{j}$

3:for $b=1,\ldots,B$ do

4: if method=‘k-grad’ then

5: Draw
$\epsilon_{1}^{(b)},\ldots,\epsilon_{k}^{(b)}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1)$
and compute $W^{(b)}$ by (5)

6: else if method=‘n+k-1-grad’ then

7: Draw
$\epsilon_{11}^{(b)},\ldots,\epsilon_{n1}^{(b)},\epsilon_{2}^{(b)},\ldots,\epsilon_{k}^{(b)}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1)$
and compute $W^{(b)}$ by (6)

8: end if

9:end for

10:Compute the quantile $c_{W}(\alpha)$ of $\\{W^{(1)},\dots,W^{(B)}\\}$ for
$\alpha\in(0,1)$

11:Return $\widetilde{\theta}_{l}\pm N^{-1/2}c_{W}(\alpha)$, $l=1,\dots,d$

###### Remark 13

Although in Algorithm 3 the same $\widetilde{\theta}$ is used for the center
of the confidence interval and for evaluating the gradients $\mathbf{g}_{ij}$,
allowing them to be different (such as in Algorithm 1) can save one round of
communication. For example, we can use $\widetilde{\theta}^{(\tau)}$ for the
center of the confidence interval, while the gradients are evaluated with
$\widetilde{\theta}^{(\tau-1)}$.

##  Appendix B Nodewise Lasso

In Algorithm 4, we state the nodewise Lasso method for constructing
approximate inverse Hessian matrix used in Section 3.1.1 of van de Geer et al.
(2014), which we apply in Algorithm 1. We define the components of
$\widehat{\gamma}_{l}$ as
$\widehat{\gamma}_{l}=\\{\widehat{\gamma}_{l,l^{\prime}};l^{\prime}=1,\dots,d,l^{\prime}\neq
l\\}$. We denote by $\widehat{M}_{l,-l}$ the $l$-th row of $\widehat{M}$
without the diagonal element $(l,l)$, and by $\widehat{M}_{-l,-l}$ the
submatrix without the $l$-th row and $l$-th column.

Algorithm 4 Node($\widehat{M}$)

1:Require: sample Hessian matrix $\widehat{M}\in\mathbb{R}^{d\times d}$,
hyperparameters $\\{\lambda_{l}\\}_{l=1}^{d}$

2:for $l=1,\ldots,d$ do

3: Compute
$\widehat{\gamma}_{l}=\operatorname*{\arg\min}_{\gamma\in\mathbb{R}^{d-1}}\widehat{M}_{l,l}-2\widehat{M}_{l,-l}\gamma+\gamma^{\top}\widehat{M}_{-l,-l}\gamma+2\lambda_{l}\|\gamma\|_{1}$

4: Compute
$\widehat{\tau}_{l}^{2}=\widehat{M}_{l,l}-\widehat{M}_{l,-l}\widehat{\gamma}_{l}$

5:end for

6:Construct $\widehat{M^{-1}}$ as
$\widehat{M^{-1}}=\begin{pmatrix}\widehat{\tau}_{1}^{-2}&0&\dots&0\\\
0&\widehat{\tau}_{2}^{-2}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\dots&\widehat{\tau}_{d}^{-2}\end{pmatrix}\begin{pmatrix}1&-\widehat{\gamma}_{1,2}&\dots&-\widehat{\gamma}_{1,d}\\\
-\widehat{\gamma}_{2,1}&1&\dots&-\widehat{\gamma}_{2,d}\\\
\vdots&\vdots&\ddots&\vdots\\\
-\widehat{\gamma}_{d,1}&-\widehat{\gamma}_{d,2}&\dots&1\end{pmatrix}.$

###### Remark 14

Throughout this paper, we fix the choice of nodewise Lasso in Algorithm 1 for
computing an approximate inverse Hessian matrix. In practice, various
approaches (e.g., Zhang and Zhang (2014); Javanmard and Montanari (2014a)) can
be chosen from in consideration of estimation accuracy and computational
efficiency.

##  Appendix C CSL Estimator for GLMs

For the $\ell_{1}$-penalized CSL estimator of generalized linear models,
Theorem 3.3 of Wang et al. (2017) states that

$\displaystyle\big{\|}\widetilde{\theta}^{(t+1)}-\theta^{\ast}\big{\|}_{1}\lesssim
s_{0}\sqrt{\frac{\log d}{N}}+s_{0}\sqrt{\frac{\log
d}{n}}\big{\|}\widetilde{\theta}^{(t)}-\theta^{\ast}\big{\|}_{1}+Ms_{0}\big{\|}\widetilde{\theta}^{(t)}-\theta^{\ast}\big{\|}_{1}^{2},$
(16)

where $M\geq 0$ is a Lipschitz constant of the $g^{\prime\prime}$, which
exists due to Assumptions • ‣ 3.3. In linear models, $g(a,b)=(a-b)^{2}/2$,
$g^{\prime\prime}$ is a constant, so $M=0$ and CSL estimator has linear
convergence to $\theta^{\ast}$ with rate $s_{0}(\log d)^{1/2}n^{-1/2}$ until
it reaches the upper bound given by the first term, which is also the rate of
the centralized (oracle) estimator. For GLMs, however, $M>0$ and the third
term can be dominant when $t$ is small. For example, when $t=0$, given that
$\|\widetilde{\theta}^{(0)}-\theta^{\ast}\|_{1}\lesssim s_{0}(\log
d)^{1/2}n^{-1/2}$, it is easy to see that the third term is always $s_{0}$
times larger than the second term (up to a constant), and a larger $n$ is
required to ensure third term is less than
$\big{\|}\widetilde{\theta}^{(t)}-\theta^{\ast}\big{\|}_{1}$ and the error is
shrinking. However, when $t$ is sufficiently large, this dominance reverses.
The threshold is given by the $\tau_{0}$ in (14), and this implies the three
phases of convergence: When $t\leq\tau_{0}$, the third term dominates and the
convergence is quadratic; when $t>\tau_{0}$, the second term dominates the
third and the linear convergence kicks in. Finally, when $t$ is sufficiently
large, the first term dominates. Our analysis complements that of Wang et al.
(2017), while in their Corollary 3.7 it is simply assumed that the second term
dominates the third.

##  Appendix D Variable Descriptions

We use the following variables in our model for the semi-synthetic study in
Section 5:

  * •

Year: from 1987 to 2008,

  * •

Month: from 1 to 12,

  * •

DayOfWeek: from 1 (Monday) to 7 (Sunday),

  * •

CRSDepTime: scheduled departure time (in four digits, first two representing
hour, last two representing minute),

  * •

CRSArrTime: scheduled arrival time (in the same format as above),

  * •

UniqueCarrier: unique carrier code,

  * •

Origin: origin (in IATA airport code),

  * •

Dest: destination (in IATA airport code),

  * •

ArrDelay: arrival delay (in minutes). Positive value means there is a delay.

The complete variable information can be found at http://stat-
computing.org/dataexpo/2009/the-data.html.

##  Appendix E Creation of Dummy Variables

We categorize CRSDepTime and CRSArrTime into $24$ one-hour time intervals
(e.g., 1420 is converted to 14 to represent the interval [14:00,15:00]), and
then treat Year, Month, DayOfWeek, CRSDepTime, CRSArrTime, UniqueCarrier,
Origin, and Dest as nominal predictors. The nominal predictors are encoded by
dummies with appropriate dimensions and merging all categories of lower counts
into “others”, and either “others” or the smallest ordinal value is treated as
the baseline.

To ensure that none of the columns of the design matrix on the master node is
completely zero so that the nodewise Lasso can be computed, we create the
dummy variables using only the observations in the master node on the dataset
$\mathcal{D}_{1}$. Specifically, for variables UniqueCarrier, Origin, and
Dest, we keep the top categories that make up $90\%$ of the data in the master
node on $\mathcal{D}_{1}$; the rest categories are merged into “others” and
are treated as baseline. For CRSDepTime and CRSArrTime, we merge the time
intervals 23:00-6:00 and 1:00-7:00 respectively (due to their low counts) and
use them as baseline. For Year, Month, and DayOfWeek, we treat year 1987,
January, and Monday as baseline respectively.

##  Appendix F Extension to Heteroscedastic Error Across Machines

As suggested by the associated editor, here we consider an extension to a more
challenging scenario for linear models where the data across machines have
heteroscedastic errors. In this scenario, Algorithm 2 can no longer apply as
it relies on the homogeneity in data across machines. We provide a new
Algorithm 5 by exploiting the multiplier bootstrap idea underlying the “High-
Dimensional Metrics” (HDM, Chernozhukov et al. (2016)).

Algorithm 5 Simultaneous inference for distributed data with
heteroscedasticity

1:Require: $\tau\geq 1$ rounds of communication; nodewise Lasso procedure
Node$(\cdot,\cdot)$ with hyperparameters $\\{\lambda_{l}\\}_{l=1}^{d}$,
theoretical constant $c$

2:$\widetilde{\theta}^{(0)}\leftarrow\operatorname*{\arg\min}_{\theta}\mathcal{L}_{1}(\theta)+\lambda^{(0)}\|\theta\|_{1}$
at $\mathcal{M}_{1}$, where $\lambda^{(0)}$ is chosen by cross-validation
using the data at $\mathcal{M}_{1}$

3:Compute $\widetilde{\Theta}$ by running
Node$(\nabla^{2}\mathcal{L}_{1}(\widetilde{\theta}^{(0)}),\\{\lambda_{l}\\}_{l=1}^{d})$
at $\mathcal{M}_{1}$

4:for $t=1,\ldots,\tau$ do

5: Transmit $\widetilde{\theta}^{(t-1)}$ to $\\{\mathcal{M}_{j}\\}_{j=2}^{k}$

6: Compute $\nabla\mathcal{L}_{1}(\widetilde{\theta}^{(t-1)})$ and
$\psi_{1}^{(t-1)}=n^{-1}\sum_{i=1}^{n}\nabla\mathcal{L}((x_{i1},y_{i1}),\widetilde{\theta}^{(t-1)})^{2}$
at $\mathcal{M}_{1}$

7: for $j=2,\ldots,k$ do

8: Compute $\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(t-1)})$ and
$\psi_{j}^{(t-1)}=n^{-1}\sum_{i=1}^{n}\nabla\mathcal{L}((x_{ij},y_{ij}),\widetilde{\theta}^{(t-1)})^{2}$
at $\mathcal{M}_{j}$

9: Transmit $\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(t-1)})$ and
$\psi_{j}^{(t-1)}$ to $\mathcal{M}_{1}$

10: end for

11: $\nabla\mathcal{L}_{N}(\widetilde{\theta}^{(t-1)})\leftarrow
k^{-1}\sum_{j=1}^{k}\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(t-1)})$ at
$\mathcal{M}_{1}$

12: if $t<\tau$ then

13: for $b=1,\ldots,B$ do

14: Draw
$\epsilon_{1}^{(b)},\ldots,\epsilon_{k}^{(b)}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1)$

15: $\Lambda_{b}^{(t)}\leftarrow
ck^{-1}\|\sum_{j=1}^{k}\epsilon_{j}^{(b)}\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(t-1)})\|_{\infty}$

16: end for

17: $\lambda^{(t)}\leftarrow 90\%$ quantile of
$\\{\Lambda_{1}^{(t)},\dots,\Lambda_{B}^{(t)}\\}$

18: for $l=1,\ldots,d$ do

19:
$\Psi_{l}^{(t)}\leftarrow\sqrt{k^{-1}\sum_{j=1}^{k}(\psi_{j}^{(t-1)})_{l}}$

20: end for

21: $\Psi^{(t)}\leftarrow\text{diag}(\Psi_{1}^{(t)},\dots,\Psi_{d}^{(t)})$

22:
$\widetilde{\theta}^{(t)}\leftarrow\operatorname*{\arg\min}_{\theta}\mathcal{L}_{1}(\theta)-\theta^{\top}\left(\nabla\mathcal{L}_{1}(\widetilde{\theta}^{(t-1)})-\nabla\mathcal{L}_{N}(\widetilde{\theta}^{(t-1)})\right)+\lambda^{(t)}\|\Psi^{(t)}\theta\|_{1}$
at $\mathcal{M}_{1}$

23: else

24:
$\widetilde{\theta}^{(\tau)}\leftarrow\widetilde{\theta}^{(\tau-1)}-\widetilde{\Theta}\nabla\mathcal{L}_{N}(\widetilde{\theta}^{(\tau-1)})$
at $\mathcal{M}_{1}$

25: end if

26:end for

27:Run DistBoots$(\text{`{k-grad}' or
`{n+k-1-grad}'},\widetilde{\theta}=\widetilde{\theta}^{(\tau)},\\{\mathbf{g}_{j}=\nabla\mathcal{L}_{j}(\widetilde{\theta}^{(\tau-1)})\\}_{j=1}^{k},$

28: $\widetilde{\Theta}=\widetilde{\Theta})$ at $\mathcal{M}_{1}$

In Algorithm 5, we select the regularization parameters
$\\{\lambda^{(t)}\\}_{t=1}^{\tau-1}$ in lines 13-17 by integrating the idea of
Spindler et al. (2016). In addition, we handle heteroscedasticity by data-
driven regularization loadings $\Psi^{(t)}$ in lines 18-22.

Under heteroscedasticity, we expect the k-grad in Algorithm 3 to continue
being valid because it treats each machine equally as an independent data
point. However, n+k-1-grad may no longer provide an accurate coverage because
each single data point in the first machine is treated as equally important as
the average of entire data in $j$ machine for $j=2,\cdots,k$, so the variance
in the first machine could dominate so that the n+k-1-grad bootstrap could
fail to precisely approximate the variance of the target empirical
distribution. A careful theoretical study deserves future research.

The empirical performance of Algorithm 5 is verified by a simulation study
based on a heteroscedastic Gaussian linear model. We fix total sample size
$N=2^{14}$ and the dimension $d=2^{10}$, and choose the number of machines $k$
from $\\{2^{2},2^{3},\dots,2^{6}\\}$. The true coefficient $\theta^{\ast}$ is
a $d$-dimensional vector in which the first $s_{0}$ coordinates are 1 and the
rest is 0, where $s_{0}\in\\{2^{2},2^{4}\\}$. We generate covariate vector $x$
independently from $\mathcal{N}(0,\Sigma)$, where $\Sigma$ is a Toeplitz
matrix with $\Sigma_{l,l^{\prime}}=0.9^{|l-l^{\prime}|}$. We introduce
heteroscedasticity across machines by first independently generating the model
noise from $\mathcal{N}(0,1)$ for all data in the master node
$\mathcal{M}_{1}$. Next, we generate model noise for each data point $i$ in
worker node $\mathcal{M}_{j}$ ($j=2,\dots,k$) independently from
$\mathcal{N}(0,\sigma_{j}^{2}+\omega_{ij})$, where the node level variance
$\sigma_{j}^{2}$ is generated independently from $\text{Unif}(2,3)$ and the
idiosyncratic variance $\omega_{ij}$ is generated independently from
$\text{Unif}(-0.2,0.2)$. For each choice of $s_{0}$ and $k$, we run Algorithm
5 with k-grad and n+k-1-grad on $1{,}000$ independently generated datasets,
and compute the empirical coverage probability and the average width based on
the results from these $1{,}000$ replications. At each replication, we draw
$B=500$ bootstrap samples, from which we calculate the $95\%$ empirical
quantile to further obtain the $95\%$ simultaneous confidence interval. For
tuning the nodewise Lasso, we use the same approach as in the main text. The
computation of the oracle width starts with fixing $(N,d,s_{0},k)$ and
generating $500$ independent datasets. For each dataset, we compute the
centralized de-biased Lasso estimator $\widehat{\theta}$. The oracle width is
defined as two times the $95\%$ empirical quantile of
$\|\widehat{\theta}-\theta^{\ast}\|_{\infty}$ of the 500 samples.

Figure 8 shows the coverage probability and efficiency in the form of relative
widths of Algorithm 5. As expected, the coverage of the simultaneous
confidence intervals is improved as the iteration goes using the new data-
driven parameter tuning and heteroscedasticity-adapted regularization. The
k-grad performs much better than the n+k-1-grad, which basically fails as the
coverage probability of n+k-1-grad is nearly zero in all cases. The failure of
the n+k-1-grad is due to the fact that it over-weigh the data in the master
node $\mathcal{M}_{1}$ which leads to an under-estimation of the variance in
other nodes, whereas in k-grad each node is weighed equally.

By comparing Figure 8 and Figure 9, we observe that our algorithm is generally
robust to the selection of $c$ as it performs similarly for $c=0.5$ and $c=1$.
However, we note that $c=0.5$ could be too small to stabilize the algorithm as
the optimization solver in 22 fails to converge in about $2\%$ of the
replications, and the divergent results are not included in Figure 8. This
suggests that the penalty at $c=0.5$ may be so small that leads to an ill-
conditioned objective function. After increasing $c$ from $0.5$ to $1$,
optimization solvers of all the replications stably converge.

Figure 8: Under $c=0.5$, empirical coverage probability (left axis, solid
lines) and average relative width (right axis, dashed lines) of simultaneous
confidence intervals by k-grad and n+k-1-grad in sparse linear regression with
Toeplitz design and varying sparsity. Black solid line represents the $95\%$
nominal level and black dashed line represents 1 on the right $y$-axis. Figure
9: Under $c=1$, empirical coverage probability (left axis, solid lines) and
average relative width (right axis, dashed lines) of simultaneous confidence
intervals by k-grad and n+k-1-grad in sparse linear regression with Toeplitz
design and varying sparsity. Black solid line represents the $95\%$ nominal
level and black dashed line represents 1 on the right $y$-axis.

## References

  * DVN (2008) Data Expo 2009: Airline on time data, 2008.  URL https://doi.org/10.7910/DVN/HG7NV7.
  * Banerjee et al. (2019) Moulinath Banerjee, Cecile Durot, Bodhisattva Sen, et al.  Divide and conquer in nonstandard problems and the super-efficiency phenomenon.  _The Annals of Statistics_ , 47(2):720–757, 2019\.
  * Battey et al. (2018) Heather Battey, Jianqing Fan, Han Liu, Junwei Lu, and Ziwei Zhu.  Distributed estimation and inference with statistical guarantees.  _Annals of Statistics_ , 46(3):1352–1382, 2018\.
  * Belloni et al. (2018) Alexandre Belloni, Victor Chernozhukov, Denis Chetverikov, and Ying Wei.  Uniformly valid post-regularization confidence regions for many functional parameters in z-estimation framework.  _Ann. Statist._ , 46(6B):3643–3675, 12 2018.  doi: 10.1214/17-AOS1671.
  * Belloni et al. (2019) Alexandre Belloni, Victor Chernozhukov, and Kengo Kato.  Valid post-selection inference in high-dimensional approximately sparse quantile regression models.  _Journal of the American Statistical Association_ , 114(526):749–758, 2019.  doi: 10.1080/01621459.2018.1442339.
  * Cai and Sun (2017) T Tony Cai and Wenguang Sun.  Large-scale global and simultaneous inference: Estimation and testing in very high dimensions.  _Annual Review of Economics_ , 9:411–439, 2017.
  * Chen et al. (2018) Xi Chen, Weidong Liu, and Yichen Zhang.  First-order newton-type estimator for distributed estimation and inference.  _arXiv preprint arXiv:1811.11368_ , 2018.
  * Chen et al. (2019) Xi Chen, Weidong Liu, and Yichen Zhang.  Quantile regression under memory constraint.  _Ann. Statist._ , 47(6):3244–3273, 12 2019.  doi: 10.1214/18-AOS1777.
  * Chen et al. (2020) Xi Chen, Jason D Lee, He Li, and Yun Yang.  Distributed estimation for principal component analysis: a gap-free approach.  _arXiv preprint arXiv:2004.02336_ , 2020.
  * Chen and Xie (2014) Xueying Chen and Min-ge Xie.  A split-and-conquer approach for analysis of extraordinarily large data.  _Statistica Sinica_ , pages 1655–1684, 2014.
  * Chernozhukov et al. (2013) Victor Chernozhukov, Denis Chetverikov, Kengo Kato, et al.  Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors.  _The Annals of Statistics_ , 41(6):2786–2819, 2013.
  * Chernozhukov et al. (2016) Victor Chernozhukov, Chris Hansen, and Martin Spindler.  hdm: High-dimensional metrics.  _The R Journal_ , 8(2):185–199, 2016.
  * Dezeure et al. (2017) Ruben Dezeure, Peter Bühlmann, and Cun-Hui Zhang.  High-dimensional simultaneous inference with the bootstrap.  _Test_ , 26(4):685–719, 2017.
  * Fan et al. (2019a) Jianqing Fan, Yongyi Guo, and Kaizheng Wang.  Communication-efficient accurate statistical estimation.  _arXiv preprint arXiv:1906.04870_ , 2019a.
  * Fan et al. (2019b) Jianqing Fan, Dong Wang, Kaizheng Wang, and Ziwei Zhu.  Distributed estimation of principal eigenspaces.  _Annals of statistics_ , 47(6):3009, 2019b.
  * Huang and Huo (2019) Cheng Huang and Xiaoming Huo.  A distributed one-step estimator.  _Mathematical Programming_ , 174(1-2):41–76, 2019\.
  * Ito and Lee (2005) Harumi Ito and Darin Lee.  Assessing the impact of the september 11 terrorist attacks on us airline demand.  _Journal of Economics and Business_ , 57(1):75–95, 2005.
  * Javanmard and Montanari (2014a) Adel Javanmard and Andrea Montanari.  Confidence intervals and hypothesis testing for high-dimensional regression.  _The Journal of Machine Learning Research_ , 15(1):2869–2909, 2014a.
  * Javanmard and Montanari (2014b) Adel Javanmard and Andrea Montanari.  Hypothesis testing in high-dimensional regression under the gaussian random design model: Asymptotic theory.  _IEEE Transactions on Information Theory_ , 60(10):6522–6554, 2014b.
  * Jordan et al. (2019) Michael I Jordan, Jason D Lee, and Yun Yang.  Communication-efficient distributed statistical inference.  _Journal of the American Statistical Association_ , 114(526):668–681, 2019.
  * Kleiner et al. (2014) Ariel Kleiner, Ameet Talwalkar, Purnamrita Sarkar, and Michael I Jordan.  A scalable bootstrap for massive data.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76(4):795–816, 2014.
  * Lan et al. (2018) Guanghui Lan, Soomin Lee, and Yi Zhou.  Communication-efficient algorithms for decentralized and stochastic optimization.  _Mathematical Programming_ , pages 1–48, 2018.
  * Lee et al. (2017) Jason D Lee, Qiang Liu, Yuekai Sun, and Jonathan E Taylor.  Communication-efficient sparse regression.  _The Journal of Machine Learning Research_ , 18(1):115–144, 2017.
  * Li et al. (2013) Runze Li, Dennis KJ Lin, and Bing Li.  Statistical inference in massive data sets.  _Applied Stochastic Models in Business and Industry_ , 29(5):399–409, 2013.
  * Rudelson and Zhou (2013) M. Rudelson and S. Zhou.  Reconstruction from anisotropic random measurements.  _IEEE Transactions on Information Theory_ , 59(6):3434–3447, 2013.  doi: 10.1109/TIT.2013.2243201.
  * Sengupta et al. (2016) Srijan Sengupta, Stanislav Volgushev, and Xiaofeng Shao.  A subsampled double bootstrap for massive data.  _Journal of the American Statistical Association_ , 111(515):1222–1232, 2016.
  * Shi et al. (2018) Chengchun Shi, Wenbin Lu, and Rui Song.  A massive data framework for m-estimators with cubic-rate.  _Journal of the American Statistical Association_ , 113(524):1698–1709, 2018.
  * Singh and Kaur (2014) Kamalpreet Singh and Ravinder Kaur.  Hadoop: addressing challenges of big data.  In _2014 IEEE International Advance Computing Conference (IACC)_ , pages 686–689. IEEE, 2014.
  * Spindler et al. (2016) Martin Spindler, Victor Chernozhukov, and Christian Hansen.  hdm: High-dimensional metrics.  _R Package Version 0.1. 0. Available at http://CRAN. R-project. org/package= hdm.[233]_ , 2016.
  * van de Geer et al. (2014) Sara van de Geer, Peter Bühlmann, Ya’acov Ritov, Ruben Dezeure, et al.  On asymptotically optimal confidence regions and tests for high-dimensional models.  _The Annals of Statistics_ , 42(3):1166–1202, 2014.
  * Volgushev et al. (2019) Stanislav Volgushev, Shih-Kang Chao, Guang Cheng, et al.  Distributed inference for quantile regression processes.  _The Annals of Statistics_ , 47(3):1634–1662, 2019.
  * Wang and Zhang (2017) Jialei Wang and Tong Zhang.  Improved optimization of finite sums with minibatch stochastic variance reduced proximal iterations.  _arXiv preprint arXiv:1706.07001_ , 2017.
  * Wang et al. (2017) Jialei Wang, Mladen Kolar, Nathan Srebro, and Tong Zhang.  Efficient distributed learning with sparsity.  In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_ , pages 3636–3645. JMLR. org, 2017.
  * Yu et al. (2020a) Ming Yu, Varun Gupta, and Mladen Kolar.  Simultaneous inference for pairwise graphical models with generalized score matching.  _Journal of Machine Learning Research_ , 21(91):1–51, 2020a.
  * Yu et al. (2020b) Yang Yu, Shih-Kang Chao, and Guang Cheng.  Simultaneous inference for massive data: Distributed bootstrap.  In _International Conference on Machine Learning_ , pages 10892–10901. PMLR, 2020b.
  * Zhang and Zhang (2014) Cun-Hui Zhang and Stephanie S Zhang.  Confidence intervals for low dimensional parameters in high dimensional linear models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76(1):217–242, 2014.
  * Zhang and Cheng (2017) Xianyang Zhang and Guang Cheng.  Simultaneous inference for high-dimensional linear models.  _Journal of the American Statistical Association_ , 112(518):757–768, 2017.
  * Zhang et al. (2012) Yuchen Zhang, Martin J Wainwright, and John C Duchi.  Communication-efficient algorithms for statistical optimization.  In _Advances in Neural Information Processing Systems_ , pages 1502–1510, 2012.
  * Zhao et al. (2016) Tianqi Zhao, Guang Cheng, and Han Liu.  A partially linear framework for massive heterogeneous data.  _Annals of Statistics_ , 44(4):1400, 2016.
  * Zhu et al. (2020) Xuening Zhu, Feng Li, and Hansheng Wang.  Least squares approximation for a distributed system.  _ArXiv Preprint arXiv:1908.04904_ , 2020.

SUPPLEMENTARY MATERIAL

## 1\. Proofs of Main Results

To simplify the notation, in the proof we denote
$\bar{\theta}=\widetilde{\theta}^{(\tau-1)}$, where
$\widetilde{\theta}^{(\tau-1)}$ is the $\ell_{1}$-penalized estimator at
$\tau-1$ iterator output by Algorithm 1. Denote
$\widetilde{\theta}=\widetilde{\theta}^{(\tau)}$ output by Algorithm 1.

Proof of Theorem 3.  We apply Theorem 3 of Wang et al. (2017), where their
Assumption 2 is inherited from Assumption • ‣ 3.2, and obtain that if $n\gg
s_{0}^{2}\log d$,

$\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}=\left\|\widetilde{\theta}^{(\tau-1)}-\theta^{\ast}\right\|_{1}=O_{P}\left(s_{0}\sqrt{\frac{\log
d}{N}}+\left(s_{0}\sqrt{\frac{\log d}{n}}\right)^{\tau}\right).$

Then, by Lemma 15, we have that $\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{\overline{W}}(\alpha))-\alpha\right|=o(1),$ as long as
$n\gg{s^{*}}^{2}\log^{3+\kappa}d+{s^{*}}\log^{5+\kappa}d+s_{0}^{2}\log d$,
$k\gg{s^{*}}^{2}\log^{5+\kappa}d$, and

$s_{0}\sqrt{\frac{\log d}{N}}+\left(s_{0}\sqrt{\frac{\log
d}{n}}\right)^{\tau}\ll\min\left\\{\frac{1}{\sqrt{k{s^{*}}}\log^{1+\kappa}d},\frac{1}{\sqrt{n{s^{*}}}\log^{1+\kappa}d}\right\\}.$

These conditions hold if
$n\gg({s^{*}}^{2}+{s^{*}}s_{0}^{2})\log^{3+\kappa}d+{s^{*}}\log^{5+\kappa}d$,
$k\gg{s^{*}}s_{0}^{2}\log^{3+\kappa}d+{s^{*}}^{2}\log^{5+\kappa}d$, and

$\tau>\max\left\\{\frac{\log k+\log{{s^{*}}}+\log(C\log^{2+\kappa}d)}{\log
n-\log(s_{0}^{2})-\log\log
d},1+\frac{\log{{s^{*}}}+\log(s_{0}^{2})+\log(C\log^{3+\kappa}d)}{\log
n-\log(s_{0}^{2})-\log\log d}\right\\}.$

If $n=d^{\gamma_{n}}$, $k=d^{\gamma_{k}}$, $\overline{s}=s_{0}\vee
s^{*}=d^{\gamma_{s}}$ for some constants $\gamma_{n}$, $\gamma_{k}$, and
$\gamma_{s}$, then a sufficient condition is $\gamma_{n}>3\gamma_{s}$,
$\gamma_{k}>3\gamma_{s}$, and

$\tau\geq
1+\left\lfloor\max\left\\{\frac{\gamma_{k}+\gamma_{s}}{\gamma_{n}-2\gamma_{s}},1+\frac{3\gamma_{s}}{\gamma_{n}-2\gamma_{s}}\right\\}\right\rfloor.$

$\blacksquare$

Proof of Theorem 4.  Similarly to the proof of Theorem 3, applying Theorem 3
of Wang et al. (2017) and Lemma 16, we have that
$\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{\widetilde{W}}(\alpha))-\alpha\right|=o(1),$ as long as
$n\gg{s^{*}}^{2}\log^{3+\kappa}d+{s^{*}}\log^{5+\kappa}d+s_{0}^{2}\log d$,
$n+k\gg{s^{*}}^{2}\log^{5+\kappa}d$, and

$s_{0}\sqrt{\frac{\log d}{N}}+\left(s_{0}\sqrt{\frac{\log
d}{n}}\right)^{\tau}\ll\min\left\\{\frac{1}{\sqrt{k{s^{*}}}\log^{1+\kappa}d},\frac{1}{{s^{*}}\sqrt{\log((n+k)d)}\log^{2+\kappa}d}\right\\}.$

These conditions hold if
$n\gg({s^{*}}^{2}+{s^{*}}s_{0}^{2})\log^{3+\kappa}d+{s^{*}}\log^{5+\kappa}d$,
$n+k\gg{s^{*}}^{2}\log^{5+\kappa}d$,
$nk\gg{s^{*}}^{2}s_{0}^{2}\log^{5+\kappa}d$, and

$\tau>\max\left\\{\frac{\log k+\log{{s^{*}}}+\log(C\log^{2+\kappa}d)}{\log
n-\log(s_{0}^{2})-\log\log
d},\frac{\log{{s^{*}}^{2}}+\log\log((n+k)d)+\log(C\log^{4+\kappa}d)}{\log
n-\log(s_{0}^{2})-\log\log d}\right\\}.$

If $n=d^{\gamma_{n}}$, $k=d^{\gamma_{k}}$, $\overline{s}=s_{0}\vee
s^{*}=d^{\gamma_{s}}$ for some constants $\gamma_{n}$, $\gamma_{k}$, and
$\gamma_{s}$, then a sufficient condition is $\gamma_{n}>3\gamma_{s}$,
$\gamma_{n}+\gamma_{k}>4\gamma_{s}$, and

$\tau\geq
1+\left\lfloor\frac{(\gamma_{k}\vee\gamma_{s})+\gamma_{s}}{\gamma_{n}-2\gamma_{s}}\right\rfloor.$

$\blacksquare$

Proof of Theorem 10.  We apply Theorem 6 of Wang et al. (2017), where their
Assumption 2 is inherited from Assumption • ‣ 3.3, and obtain that if $n\gg
s_{0}^{4}\log d$,

$\displaystyle\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}$
$\displaystyle=\left\|\widetilde{\theta}^{(\tau-1)}-\theta^{\ast}\right\|_{1}=\begin{cases}O_{P}\left(s_{0}\sqrt{\frac{\log
d}{N}}+\frac{1}{s_{0}}\left(s_{0}^{2}\sqrt{\frac{\log
d}{n}}\right)^{2^{\tau-1}}\right),&\tau\leq\tau_{0}+1,\\\
O_{P}\left(s_{0}\sqrt{\frac{\log
d}{N}}+\frac{1}{s_{0}}\left(s_{0}^{2}\sqrt{\frac{\log
d}{n}}\right)^{2^{\tau_{0}}}\left(s_{0}\sqrt{\frac{\log
d}{n}}\right)^{\tau-\tau_{0}-1}\right),&\tau>\tau_{0}+1,\end{cases}$

where $\tau_{0}$ is the smallest integer $t$ such that

$\left(s_{0}^{2}\sqrt{\frac{\log d}{n}}\right)^{2^{t}}\lesssim
s_{0}\sqrt{\frac{\log d}{n}},$

that is,

$\tau_{0}=\left\lceil\log_{2}\left(\frac{\log n-\log(s_{0}^{2})-\log(C\log
d)}{\log n-\log(s_{0}^{4})-\log\log d}\right)\right\rceil.$

Then, by Lemma 17, we have that $\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{\overline{W}}(\alpha))-\alpha\right|=o(1),$ as long as
$n\gg(s_{0}^{2}+{s^{*}}^{2})\log^{3+\kappa}d+(s_{0}+{s^{*}})\log^{5+\kappa}d+s_{0}^{4}\log
d$, $k\gg{s^{*}}^{2}\log^{5+\kappa}d$, and

$\displaystyle s_{0}\sqrt{\frac{\log
d}{N}}+\frac{1}{s_{0}}\left(s_{0}^{2}\sqrt{\frac{\log
d}{n}}\right)^{2^{\tau-1}}\ll\min\left\\{\frac{1}{\sqrt{k{s^{*}}}s_{0}\log^{1+\kappa}d},\frac{1}{\sqrt{n{s^{*}}}\log^{1+\kappa}d}\right\\},$

if $\tau\leq\tau_{0}+1$, and

$\displaystyle s_{0}\sqrt{\frac{\log
d}{N}}+\frac{1}{s_{0}}\left(s_{0}^{2}\sqrt{\frac{\log
d}{n}}\right)^{2^{\tau_{0}}}\left(s_{0}\sqrt{\frac{\log
d}{n}}\right)^{\tau-\tau_{0}-1}\ll\min\left\\{\frac{1}{\sqrt{k{s^{*}}}s_{0}\log^{1+\kappa}d},\frac{1}{\sqrt{n{s^{*}}}\log^{1+\kappa}d}\right\\},$

if $\tau>\tau_{0}+1$.

If $n=d^{\gamma_{n}}$, $k=d^{\gamma_{k}}$, $\overline{s}=s_{0}\vee
s^{*}=d^{\gamma_{s}}$ for some constants $\gamma_{n}$, $\gamma_{k}$, and
$\gamma_{s}$, then a sufficient condition is $\gamma_{n}>5\gamma_{s}$,
$\gamma_{k}>3\gamma_{s}$, and

$\displaystyle\tau$ $\displaystyle\geq
1+\left\lfloor\max\left\\{1+\log_{2}\frac{\gamma_{n}-\gamma_{s}}{\gamma_{n}-4\gamma_{s}},\tau_{0}+1+\frac{\gamma_{k}+(4\cdot
2^{\tau_{0}}+1)\gamma_{s}-2^{\tau_{0}}\gamma_{n}}{\gamma_{n}-2\gamma_{s}}\right\\}\right\rfloor$
$\displaystyle=\left\lfloor\max\left\\{2+\log_{2}\frac{\gamma_{n}-\gamma_{s}}{\gamma_{n}-4\gamma_{s}},\tau_{0}+2+\frac{\gamma_{k}+(4\cdot
2^{\tau_{0}}+1)\gamma_{s}-2^{\tau_{0}}\gamma_{n}}{\gamma_{n}-2}\right\\}\right\rfloor$
$\displaystyle=\left\lfloor\max\left\\{2+\log_{2}\frac{\gamma_{n}-\gamma_{s}}{\gamma_{n}-4\gamma_{s}},\tau_{0}+\frac{\gamma_{k}+\gamma_{s}}{\gamma_{n}-2\gamma_{s}}+\nu_{0}\right\\}\right\rfloor$
$\displaystyle=\max\left\\{\tau_{0}+\left\lfloor\frac{\gamma_{k}+\gamma_{s}}{\gamma_{n}-2\gamma_{s}}+\nu_{0}\right\rfloor,2+\left\lfloor\log_{2}\frac{\gamma_{n}-\gamma_{s}}{\gamma_{n}-4\gamma_{s}}\right\rfloor\right\\},$

where

$\tau_{0}=1+\left\lfloor\log_{2}\frac{\gamma_{n}-2\gamma_{s}}{\gamma_{n}-4\gamma_{s}}\right\rfloor,\quad\nu_{0}=2-\frac{2^{\tau_{0}}(\gamma_{n}-4\gamma_{s})}{\gamma_{n}-2\gamma_{s}}\in(0,1].$

$\blacksquare$

Proof of Theorem 11.  Similarly to the proof of Theorem 4, applying Theorem 3
of Wang et al. (2017) and Lemma 18, we have that
$\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{\overline{W}}(\alpha))-\alpha\right|=o(1),$ as long as
$n\gg(s_{0}+{s^{*}})\log^{5+\kappa}d+(s_{0}^{2}+{s^{*}}^{2})\log^{3+\kappa}d$,
$n+k\gg{s^{*}}^{2}\log^{5+\kappa}d$, and

$\displaystyle s_{0}\sqrt{\frac{\log
d}{N}}+\frac{1}{s_{0}}\left(s_{0}^{2}\sqrt{\frac{\log
d}{n}}\right)^{2^{\tau-1}}$
$\displaystyle\ll\min\left\\{\frac{n+k}{{s^{*}}\left(n+k\sqrt{\log
d}+k^{3/4}\log^{3/4}d\right)\log^{2+\kappa}d},\frac{1}{\sqrt{k{s^{*}}}s_{0}\log^{1+\kappa}d},\frac{1}{\left(nk{s^{*}}\log^{1+\kappa}d\right)^{1/4}}\right\\},$

if $\tau\leq\tau_{0}+1$, and

$\displaystyle s_{0}\sqrt{\frac{\log
d}{N}}+\frac{1}{s_{0}}\left(s_{0}^{2}\sqrt{\frac{\log
d}{n}}\right)^{2^{\tau_{0}}}\left(s_{0}\sqrt{\frac{\log
d}{n}}\right)^{\tau-\tau_{0}-1}$
$\displaystyle\ll\min\left\\{\frac{n+k}{{s^{*}}\left(n+k\sqrt{\log
d}+k^{3/4}\log^{3/4}d\right)\log^{2+\kappa}d},\frac{1}{\sqrt{k{s^{*}}}s_{0}\log^{1+\kappa}d},\frac{1}{\left(nk{s^{*}}\log^{1+\kappa}d\right)^{1/4}}\right\\},$

if $\tau>\tau_{0}+1$, where

$\tau_{0}=\left\lceil\log_{2}\left(\frac{\log n-\log(s_{0}^{2})-\log(C\log
d)}{\log n-\log(s_{0}^{4})-\log\log d}\right)\right\rceil.$

If $n=d^{\gamma_{n}}$, $k=d^{\gamma_{k}}$, $\overline{s}=s_{0}\vee
s^{*}=d^{\gamma_{s}}$ for some constants $\gamma_{n}$, $\gamma_{k}$, and
$\gamma_{s}$, then a sufficient condition is $\gamma_{n}>5\gamma_{s}$, and

Let $\overline{s}=s_{0}\vee s^{*}$. If $n=\overline{s}^{\gamma_{n}}$,
$k=\overline{s}^{\gamma_{k}}$, and $d=\overline{s}^{\gamma_{d}}$ for some
constants $\gamma_{n}$, $\gamma_{k}$, and $\gamma_{d}$, then a sufficient
condition is $\gamma_{n}>5$, and if $\tau\leq\tau_{0}+1$,

$\tau\geq\max\left\\{2+\left\lfloor\log_{2}\frac{\gamma_{k}+1}{\gamma_{n}-4}\right\rfloor,1\right\\},$

and if $\tau>\tau_{0}+1$

$\displaystyle\tau$ $\displaystyle\geq
1+\left\lfloor\tau_{0}+1+\frac{\gamma_{k}+4\cdot
2^{\tau_{0}}+1-2^{\tau_{0}}\gamma_{n}}{\gamma_{n}-2}\right\rfloor$
$\displaystyle=\left\lfloor\tau_{0}+2+\frac{\gamma_{k}+4\cdot
2^{\tau_{0}}+1-2^{\tau_{0}}\gamma_{n}}{\gamma_{n}-2}\right\rfloor$
$\displaystyle=\left\lfloor\tau_{0}+\frac{\gamma_{k}+1}{\gamma_{n}-2}+\nu_{0}\right\rfloor$
$\displaystyle=\tau_{0}+\left\lfloor\frac{\gamma_{k}+1}{\gamma_{n}-2}+\nu_{0}\right\rfloor,$

where

$\tau_{0}=1+\left\lfloor\log_{2}\frac{\gamma_{n}-2}{\gamma_{n}-4}\right\rfloor,\quad\nu_{0}=2-\frac{2^{\tau_{0}}(\gamma_{n}-4)}{\gamma_{n}-2}\in(0,1].$

$\blacksquare$

## 2\. Technical Lemmas

######  Lemma 15 (k-grad)

In sparse linear model, under Assumptions • ‣ 3.2 and • ‣ 3.2, if
$n\gg{s^{*}}^{2}\log^{3+\kappa}d+{s^{*}}\log^{5+\kappa}d$,
$k\gg{s^{*}}^{2}\log^{5+\kappa}d$, and

$\displaystyle\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}\ll\min\left\\{\frac{1}{\sqrt{k{s^{*}}}\log^{1+\kappa}d},\frac{1}{\sqrt{n{s^{*}}}\log^{1+\kappa}d}\right\\},$

for some $\kappa>0$, then we have that

$\displaystyle\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{\overline{W}}(\alpha))-\alpha\right|$ $\displaystyle=o(1),\quad\text{and}$
(17) $\displaystyle\sup_{\alpha\in(0,1)}\left|P(\widehat{T}\leq
c_{\overline{W}}(\alpha))-\alpha\right|$ $\displaystyle=o(1).$ (18)

Proof of Lemma 15.  As noted by Zhang and Cheng (2017), since
$\|\sqrt{N}(\widetilde{\theta}-\theta^{\ast})\|_{\infty}=\max_{l}\sqrt{N}|\widetilde{\theta}_{l}-\theta^{\ast}_{l}|=\sqrt{N}\max_{l}\big{(}(\widetilde{\theta}_{l}-\theta^{\ast}_{l})\vee(\theta^{\ast}_{l}-\widetilde{\theta}_{l})\big{)}$,
the arguments for the bootstrap consistency result with

$\displaystyle T$
$\displaystyle=\max_{l}\sqrt{N}(\widetilde{\theta}-\theta^{\ast})_{l}\quad\text{and}$
(19) $\displaystyle\widehat{T}$
$\displaystyle=\max_{l}\sqrt{N}(\widehat{\theta}-\theta^{\ast})_{l}$ (20)

imply the bootstrap consistency result for
$T=\|\sqrt{N}(\widetilde{\theta}-\theta^{\ast})\|_{\infty}$ and
$\widehat{T}=\|\sqrt{N}(\widehat{\theta}-\theta^{\ast})\|_{\infty}$. Hence,
from now on, we redefine $T$ and $\widehat{T}$ as (19) and (20). Define an
oracle multiplier bootstrap statistic as

$\displaystyle W^{*}:\,=\max_{1\leq l\leq
d}-\frac{1}{\sqrt{N}}\sum_{i=1}^{n}\sum_{j=1}^{k}\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z_{ij})\right)_{l}\epsilon_{ij}^{*},$
(21)

where $\\{\epsilon_{ij}^{*}\\}_{i=1,\dots,n;j=1,\dots,k}$ are $N$ independent
standard Gaussian variables, also independent of the entire dataset. The proof
consists of two steps; the first step is to show that $W^{*}$ achieves
bootstrap consistency, i.e., $\sup_{\alpha\in(0,1)}|P(T\leq
c_{W^{*}}(\alpha))-\alpha|$ converges to $0$, where
$c_{W^{*}}(\alpha)=\inf\\{t\in\mathbb{R}:P_{\epsilon}(W^{*}\leq
t)\geq\alpha\\},$ and the second step is to show the bootstrap consistency of
our proposed bootstrap statistic by showing the quantiles of $W$ and $W^{*}$
are close.

Note that
$\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)=\mathbb{E}[xx^{\top}]^{-1}x(x^{\top}\theta^{\ast}-y)=\Theta
xe$ and

$\mathbb{E}\left[\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)\right)\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)\right)^{\top}\right]=\Theta\mathbb{E}\left[xx^{\top}e^{2}\right]\Theta=\sigma^{2}\Theta\Sigma\Theta=\sigma^{2}\Theta.$

Then, under Assumptions • ‣ 3.2 and • ‣ 3.2,

$\displaystyle\min_{l}\mathbb{E}\left[\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)\right)_{l}^{2}\right]=\sigma^{2}\min_{l}\Theta_{l,l}\geq\sigma^{2}\lambda_{\tiny{\min}}(\Theta)=\frac{\sigma^{2}}{\lambda_{\tiny{\max}}(\Sigma)},$
(22)

is bounded away from zero. Under Assumption • ‣ 3.2, $x$ is sub-Gaussian, that
is, $w^{\top}x$ is sub-Gaussian with uniformly bounded $\psi_{2}$-norm for all
$w\in S^{d-1}$. To show $w^{\top}\Theta x$ is also sub-Gaussian with uniformly
bounded $\psi_{2}$-norm, we write it as

$w^{\top}\Theta x=(\Theta w)^{\top}x=\left\|\Theta
w\right\|_{2}\left(\frac{\Theta w}{\left\|\Theta
w\right\|_{2}}\right)^{\top}x.$

Since $\Theta w/\left\|\Theta w\right\|_{2}\in S^{d-1}$, we have that
$\left(\Theta w/\left\|\Theta w\right\|_{2}\right)x$ is sub-Gaussian with
$O(1)$ $\psi_{2}$-norm, and hence, $w^{\top}\Theta x$ is sub-Gaussian with
$O(\left\|\Theta
w\right\|_{2})=O(\lambda_{\tiny{\max}}(\Theta))=O(\lambda_{\tiny{\min}}(\Sigma)^{-1})=O(1)$
$\psi_{2}$-norm, under Assumption • ‣ 3.2. Since $e$ is also sub-Gaussian
under Assumption • ‣ 3.2 and is independent of $w^{\top}\Theta x$, we have
that $w^{\top}\Theta xe$ is sub-exponential with uniformly bounded
$\psi_{1}$-norm for all $w\in S^{d-1}$, and also, all
$\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)\right)_{l}$
are sub-exponential with uniformly bounded $\psi_{1}$-norm. Combining this
with (22), we have verified Assumption (E.1) of Chernozhukov et al. (2013) for
$\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)$.

Define

$\displaystyle T_{0}:\,=\max_{1\leq l\leq
d}-\sqrt{N}\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right)_{l},$
(23)

which is a Bahadur representation of $T$. Under the condition
$\log^{7}(dN)/N\lesssim N^{-c}$ for some constant $c>0$, which holds if
$N\gtrsim\log^{7+\kappa}d$ for some $\kappa>0$, applying Theorem 3.2 and
Corollary 2.1 of Chernozhukov et al. (2013), we obtain that for some constant
$c>0$ and for every $v,\zeta>0$,

$\displaystyle\begin{split}\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{W^{*}}(\alpha))-\alpha\right|&\lesssim
N^{-c}+v^{1/3}\left(1\vee\log\frac{d}{v}\right)^{2/3}+P\left(\left|\\!\left|\\!\left|{\widehat{\Omega}-\Omega_{0}}\right|\\!\right|\\!\right|_{\max}>v\right)\\\
&\quad+\zeta\sqrt{1\vee\log\frac{d}{\zeta}}+P\left(|T-T_{0}|>\zeta\right),\end{split}$
(24)

where

$\displaystyle\begin{split}\widehat{\Omega}&:\,=\operatorname{cov}_{\epsilon}\left(-\frac{1}{\sqrt{N}}\sum_{i=1}^{n}\sum_{j=1}^{k}\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z_{ij})\epsilon_{ij}^{*}\right)\\\
&=\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\left(\frac{1}{N}\sum_{i=1}^{n}\sum_{j=1}^{k}\nabla\mathcal{L}(\theta^{\ast};Z_{ij})\nabla\mathcal{L}(\theta^{\ast};Z_{ij})^{\top}\right)\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1},\quad\text{and}\end{split}$
(25) $\displaystyle\Omega_{0}$
$\displaystyle:\,=\operatorname{cov}\left(-\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)\right)=\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\mathbb{E}\left[\nabla\mathcal{L}(\theta^{\ast};Z)\nabla\mathcal{L}(\theta^{\ast};Z)^{\top}\right]\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}.$
(26)

To show the quantiles of $\overline{W}$ and $W^{*}$ are close, we first have
that for any $\omega$ such that $\alpha+\omega,\alpha-\omega\in(0,1)$,

$\displaystyle P(\\{T\leq c_{\overline{W}}(\alpha)\\}\ominus\\{T\leq
c_{W^{*}}(\alpha)\\})$ $\displaystyle\leq 2P(c_{W^{*}}(\alpha-\omega)<T\leq
c_{W^{*}}(\alpha+\omega))+P(c_{W^{*}}(\alpha-\omega)>c_{\overline{W}}(\alpha))+P(c_{\overline{W}}(\alpha)>c_{W^{*}}(\alpha+\omega)),$

where $\ominus$ denotes symmetric difference. Following the arguments in the
proof of Lemma 3.2 of Chernozhukov et al. (2013), we have that

$P(c_{\overline{W}}(\alpha)>c_{W^{*}}(\alpha+\pi(u)))\leq
P\left(\left|\\!\left|\\!\left|{\overline{\Omega}-\widehat{\Omega}}\right|\\!\right|\\!\right|_{\max}>u\right),\quad\text{and}$
$P(c_{W^{*}}(\alpha-\pi(u))>c_{\overline{W}}(\alpha))\leq
P\left(\left|\\!\left|\\!\left|{\overline{\Omega}-\widehat{\Omega}}\right|\\!\right|\\!\right|_{\max}>u\right),$

where $\pi(u):\,=u^{1/3}\left(1\vee\log(d/u)\right)^{2/3}$ and

$\displaystyle\begin{split}\overline{\Omega}&:\,=\operatorname{cov}_{\epsilon}\left(-\frac{1}{\sqrt{k}}\sum_{j=1}^{k}\widetilde{\Theta}\sqrt{n}\left(\nabla\mathcal{L}_{j}(\bar{\theta})-\nabla\mathcal{L}_{N}(\bar{\theta})\right)\epsilon_{j}\right)\\\
&=\widetilde{\Theta}\left(\frac{1}{k}\sum_{j=1}^{k}n\left(\nabla\mathcal{L}_{j}(\bar{\theta})-\nabla\mathcal{L}_{N}(\bar{\theta})\right)\left(\nabla\mathcal{L}_{j}(\bar{\theta})-\nabla\mathcal{L}_{N}(\bar{\theta})\right)^{\top}\right)\widetilde{\Theta}^{\top}.\end{split}$
(27)

By letting $\omega=\pi(u)$, we have that

$\displaystyle P(\\{T\leq c_{\overline{W}}(\alpha)\\}\ominus\\{T\leq
c_{W^{*}}(\alpha)\\})$ $\displaystyle\leq 2P(c_{W^{*}}(\alpha-\pi(u))<T\leq
c_{W^{*}}(\alpha+\pi(u)))+P(c_{W^{*}}(\alpha-\pi(u))>c_{\overline{W}}(\alpha))+P(c_{\overline{W}}(\alpha)>c_{W^{*}}(\alpha+\pi(u)))$
$\displaystyle\leq 2P(c_{W^{*}}(\alpha-\pi(u))<T\leq
c_{W^{*}}(\alpha+\pi(u)))+2P\left(\left|\\!\left|\\!\left|{\overline{\Omega}-\widehat{\Omega}}\right|\\!\right|\\!\right|_{\max}>u\right),$

where by (24),

$\displaystyle P(c_{W^{*}}(\alpha-\pi(u))<T\leq c_{W^{*}}(\alpha+\pi(u)))$
$\displaystyle=P(T\leq c_{W^{*}}(\alpha+\pi(u)))-P(T\leq
c_{W^{*}}(\alpha-\pi(u)))$
$\displaystyle\lesssim\pi(u)+N^{-c}+\zeta\sqrt{1\vee\log\frac{d}{\zeta}}+P\left(|T-T_{0}|>\zeta\right),$

and then,

$\displaystyle\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{\overline{W}}(\alpha))-\alpha\right|$ $\displaystyle\lesssim
N^{-c}+v^{1/3}\left(1\vee\log\frac{d}{v}\right)^{2/3}+P\left(\left|\\!\left|\\!\left|{\widehat{\Omega}-\Omega_{0}}\right|\\!\right|\\!\right|_{\max}>v\right)$
$\displaystyle\quad+\zeta\sqrt{1\vee\log\frac{d}{\zeta}}+P\left(|T-T_{0}|>\zeta\right)+u^{1/3}\left(1\vee\log\frac{d}{u}\right)^{2/3}+P\left(\left|\\!\left|\\!\left|{\overline{\Omega}-\widehat{\Omega}}\right|\\!\right|\\!\right|_{\max}>u\right).$
(28)

Applying Lemmas 19, 24, and 23, we have that there exist some $\zeta,u,v>0$
such that

$\displaystyle\zeta\sqrt{1\vee\log\frac{d}{\zeta}}$
$\displaystyle+P\left(|T-T_{0}|>\zeta\right)=o(1),\quad\text{and}$ (29)
$\displaystyle u^{1/3}\left(1\vee\log\frac{d}{u}\right)^{2/3}$
$\displaystyle+P\left(\left|\\!\left|\\!\left|{\overline{\Omega}-\widehat{\Omega}}\right|\\!\right|\\!\right|_{\max}>u\right)=o(1),\quad\text{and}$
(30) $\displaystyle v^{1/3}\left(1\vee\log\frac{d}{v}\right)^{2/3}$
$\displaystyle+P\left(\left|\\!\left|\\!\left|{\widehat{\Omega}-\Omega_{0}}\right|\\!\right|\\!\right|_{\max}>v\right)=o(1),$
(31)

and hence, after simplifying the conditions, obtain the first result in the
lemma. To obtain the second result, we use Lemma 20, which yields

$\displaystyle\xi\sqrt{1\vee\log\frac{d}{\xi}}+P\left(|\widehat{T}-T_{0}|>\xi\right)=o(1).$
(32)

$\blacksquare$

######  Lemma 16 (n+k-1-grad)

In sparse linear model, under Assumptions • ‣ 3.2 and • ‣ 3.2, if
$n\gg{s^{*}}^{2}\log^{3+\kappa}d+{s^{*}}\log^{5+\kappa}d$,
$n+k\gg{s^{*}}^{2}\log^{5+\kappa}d$, $nk\gtrsim\log^{7+\kappa}d$, and

$\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}\ll\min\left\\{\frac{1}{\sqrt{k{s^{*}}}\log^{1+\kappa}d},\frac{1}{{s^{*}}\sqrt{\log((n+k)d)}\log^{2+\kappa}d}\right\\},$

for some $\kappa>0$, then we have that

$\displaystyle\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{\widetilde{W}}(\alpha))-\alpha\right|$ $\displaystyle=o(1),\quad\text{and}$
(33) $\displaystyle\sup_{\alpha\in(0,1)}\left|P(\widehat{T}\leq
c_{\widetilde{W}}(\alpha))-\alpha\right|$ $\displaystyle=o(1).$ (34)

Proof of Lemma 16.  By the argument in the proof of Lemma 15, we have that

$\displaystyle\sup_{\alpha\in(0,1)}\left|P(T\leq
c_{\widetilde{W}}(\alpha))-\alpha\right|$ $\displaystyle\lesssim
N^{-c}+v^{1/3}\left(1\vee\log\frac{d}{v}\right)^{2/3}+P\left(\left|\\!\left|\\!\left|{\widehat{\Omega}-\Omega_{0}}\right|\\!\right|\\!\right|_{\max}>v\right)$
$\displaystyle\quad+\zeta\sqrt{1\vee\log\frac{d}{\zeta}}+P\left(|T-T_{0}|>\zeta\right)+u^{1/3}\left(1\vee\log\frac{d}{u}\right)^{2/3}+P\left(\left|\\!\left|\\!\left|{\widetilde{\Omega}-\widehat{\Omega}}\right|\\!\right|\\!\right|_{\max}>u\right),$
(35)

where

$\displaystyle\begin{split}\widetilde{\Omega}&:\,=\operatorname{cov}_{\epsilon}\left(-\frac{1}{\sqrt{n+k-1}}\left(\sum_{i=1}^{n}\widetilde{\Theta}\left(\nabla\mathcal{L}(\bar{\theta};Z_{i1})-\nabla\mathcal{L}_{N}(\bar{\theta})\right)\epsilon_{i1}+\sum_{j=2}^{k}\widetilde{\Theta}\sqrt{n}\left(\nabla\mathcal{L}_{j}(\bar{\theta})-\nabla\mathcal{L}_{N}(\bar{\theta})\right)\epsilon_{j}\right)\right)\\\
&=\widetilde{\Theta}\frac{1}{n+k-1}\Bigg{(}\sum_{i=1}^{n}\left(\nabla\mathcal{L}(\theta;Z_{i1})-\nabla\mathcal{L}_{N}(\theta)\right)\left(\nabla\mathcal{L}(\theta;Z_{i1})-\nabla\mathcal{L}_{N}(\theta)\right)^{\top}\\\
&\quad+\sum_{j=2}^{k}n\left(\nabla\mathcal{L}_{j}(\theta)-\nabla\mathcal{L}_{N}(\theta)\right)\left(\nabla\mathcal{L}_{j}(\theta)-\nabla\mathcal{L}_{N}(\theta)\right)^{\top}\Bigg{)}\widetilde{\Theta}^{\top},\end{split}$
(36)

if $N\gtrsim\log^{7+\kappa}d$ for some $\kappa>0$. Applying Lemmas 19, 24, and
25, we have that there exist some $\zeta,u,v>0$ such that (29),

$\displaystyle
u^{1/3}\left(1\vee\log\frac{d}{u}\right)^{2/3}+P\left(\left|\\!\left|\\!\left|{\widetilde{\Omega}-\widehat{\Omega}}\right|\\!\right|\\!\right|_{\max}>u\right)=o(1),$
(37)

and (31) hold, and hence, after simplifying the conditions, obtain the first
result in the lemma. To obtain the second result, we use Lemma 20, which
yields (32). $\blacksquare$

######  Lemma 17 (k-grad)

In sparse GLM, under Assumptions • ‣ 3.3–• ‣ 3.3, if
$n\gg(s_{0}^{2}+{s^{*}}^{2})\log^{3+\kappa}d+(s_{0}+{s^{*}})\log^{5+\kappa}d$,
$k\gg{s^{*}}^{2}\log^{5+\kappa}d$, and

$\displaystyle\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}\ll\min\left\\{\frac{1}{\sqrt{k{s^{*}}}s_{0}\log^{1+\kappa}d},\frac{1}{\sqrt{n{s^{*}}}\log^{1+\kappa}d}\right\\},$

for some $\kappa>0$, then we have that (17) and (18) hold.

Proof of Lemma 17.  We redefine $T$ and $\widehat{T}$ as (19) and (20). We
define an oracle multiplier bootstrap statistic as in (21). Under Assumption •
‣ 3.3,

$\displaystyle\min_{l}\mathbb{E}\left[\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)\right)_{l}^{2}\right]$
$\displaystyle=\min_{l}\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\mathbb{E}\left[\nabla\mathcal{L}(\theta^{\ast};Z)\nabla\mathcal{L}(\theta^{\ast};Z)^{\top}\right]\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\right)_{l,l}$
$\displaystyle\geq\lambda_{\tiny{\min}}\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\mathbb{E}\left[\nabla\mathcal{L}(\theta^{\ast};Z)\nabla\mathcal{L}(\theta^{\ast};Z)^{\top}\right]\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\right)$
$\displaystyle\geq\lambda_{\tiny{\min}}\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\right)^{2}\lambda_{\tiny{\min}}\left(\mathbb{E}\left[\nabla\mathcal{L}(\theta^{\ast};Z)\nabla\mathcal{L}(\theta^{\ast};Z)^{\top}\right]\right)$
$\displaystyle=\frac{\lambda_{\tiny{\min}}\left(\mathbb{E}\left[\nabla\mathcal{L}(\theta^{\ast};Z)\nabla\mathcal{L}(\theta^{\ast};Z)^{\top}\right]\right)}{\lambda_{\tiny{\max}}\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})\right)^{2}}$

is bounded away from zero. Combining this with Assumption • ‣ 3.3, we have
verified Assumption (E.1) of Chernozhukov et al. (2013) for
$\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}(\theta^{\ast};Z)$.
Then, we use the same argument as in the proof of Lemma 15, and obtain (28)
with

$\displaystyle\begin{split}\overline{\Omega}&:\,=\widetilde{\Theta}(\widetilde{\theta}^{(0)})\left(\frac{1}{k}\sum_{j=1}^{k}n\left(\nabla\mathcal{L}_{j}(\bar{\theta})-\nabla\mathcal{L}_{N}(\bar{\theta})\right)\left(\nabla\mathcal{L}_{j}(\bar{\theta})-\nabla\mathcal{L}_{N}(\bar{\theta})\right)^{\top}\right)\widetilde{\Theta}(\widetilde{\theta}^{(0)})^{\top},\end{split}$
(38)

under the condition $\log^{7}(dN)/N\lesssim N^{-c}$ for some constant $c>0$,
which holds if $N\gtrsim\log^{7+\kappa}d$ for some $\kappa>0$. Applying Lemmas
21, 27, and 26, we have that there exist some $\zeta,u,v>0$ such that (29),
(30), and (31) hold, and hence, after simplifying the conditions, obtain the
first result in the lemma. To obtain the second result, we use Lemma 22, which
yields (32). $\blacksquare$

######  Lemma 18 (n+k-1-grad)

In sparse GLM, under Assumptions • ‣ 3.3–• ‣ 3.3, if
$n\gg(s_{0}+{s^{*}})\log^{5+\kappa}d+(s_{0}^{2}+{s^{*}}^{2})\log^{3+\kappa}d$,
$n+k\gg{s^{*}}^{2}\log^{5+\kappa}d$, $nk\gtrsim\log^{7+\kappa}d$, and

$\displaystyle\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}$
$\displaystyle\ll\min\Bigg{\\{}\frac{n+k}{{s^{*}}\left(n+k\sqrt{\log
d}+k^{3/4}\log^{3/4}d\right)\log^{2+\kappa}d},\frac{1}{\sqrt{k{s^{*}}}s_{0}\log^{1+\kappa}d},\frac{1}{\left(nk{s^{*}}\log^{1+\kappa}d\right)^{1/4}}\Bigg{\\}},$

for some $\kappa>0$, then we have that (33) and (34) hold.

Proof of Lemma 18.  By the argument in the proof of Lemma 17, we obtain (35)
with

$\displaystyle\begin{split}\widetilde{\Omega}&:\,=\widetilde{\Theta}(\widetilde{\theta}^{(0)})\frac{1}{n+k-1}\Bigg{(}\sum_{i=1}^{n}\left(\nabla\mathcal{L}(\theta;Z_{i1})-\nabla\mathcal{L}_{N}(\theta)\right)\left(\nabla\mathcal{L}(\theta;Z_{i1})-\nabla\mathcal{L}_{N}(\theta)\right)^{\top}\\\
&\quad+\sum_{j=2}^{k}n\left(\nabla\mathcal{L}_{j}(\theta)-\nabla\mathcal{L}_{N}(\theta)\right)\left(\nabla\mathcal{L}_{j}(\theta)-\nabla\mathcal{L}_{N}(\theta)\right)^{\top}\Bigg{)}\widetilde{\Theta}(\widetilde{\theta}^{(0)})^{\top},\end{split}$
(39)

if $N\gtrsim\log^{7+\kappa}d$ for some $\kappa>0$. Applying Lemmas 21, 27, and
28, we have that there exist some $\zeta,u,v>0$ such that (29), (37), and (31)
hold, and hence, after simplifying the conditions, obtain the first result in
the lemma. To obtain the second result, we use Lemma 22, which yields (32).
$\blacksquare$

###### Lemma 19

$T$ and $T_{0}$ are defined as in (7) and (23) respectively. In sparse linear
model, under Assumptions • ‣ 3.2 and • ‣ 3.2, provided that
$\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}=O_{P}(r_{\bar{\theta}})$ and
$n\gg{s^{*}}\log d$, we have that

$|T-T_{0}|=O_{P}\left(r_{\bar{\theta}}\sqrt{{s^{*}}k\log d}+\frac{{s^{*}}\log
d}{\sqrt{n}}\right).$

Moreover, if $n\gg{s^{*}}^{2}\log^{3+\kappa}d$ and

$\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}\ll\frac{1}{\sqrt{k{s^{*}}}\log^{1+\kappa}d},$

for some $\kappa>0$, then there exists some $\zeta>0$ such that (29) holds.

Proof of Lemma 19.  First, we note that

$\displaystyle|T-T_{0}|$ $\displaystyle\leq\max_{1\leq l\leq
d}\left|\sqrt{N}(\widetilde{\theta}-\theta^{\ast})_{l}+\sqrt{N}\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right)_{l}\right|$
$\displaystyle=\sqrt{N}\left\|\widetilde{\theta}-\theta^{\ast}+\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty},$

where we use the fact that
$|\max_{l}a_{l}-\max_{l}b_{l}|\leq\max_{l}|a_{l}-b_{l}|$ for any two vectors
$a$ and $b$ of the same dimension. Next, we bound
$\left\|\widetilde{\theta}-\theta^{\ast}+\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$.
In linear model, we have that

$\widetilde{\theta}-\theta^{\ast}+\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})=\bar{\theta}+\widetilde{\Theta}\frac{X_{N}^{\top}(y_{N}-X_{N}\bar{\theta})}{N}-\theta^{\ast}-\Theta\frac{X_{N}^{\top}(y_{N}-X_{N}\theta^{\ast})}{N},$

and then,

$\displaystyle\left\|\widetilde{\theta}-\theta^{\ast}+\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$
$\displaystyle=\left\|\bar{\theta}+\widetilde{\Theta}\frac{X_{N}^{\top}(y_{N}-X_{N}\bar{\theta})}{N}-\theta^{\ast}-\Theta\frac{X_{N}^{\top}(y_{N}-X_{N}\theta^{\ast})}{N}\right\|_{\infty}$
$\displaystyle=\left\|\bar{\theta}+\widetilde{\Theta}\frac{X_{N}^{\top}(y_{N}-X_{N}\bar{\theta})}{N}-\theta^{\ast}-\widetilde{\Theta}\frac{X_{N}^{\top}(y_{N}-X_{N}\theta^{\ast})}{N}+\widetilde{\Theta}\frac{X_{N}^{\top}(y_{N}-X_{N}\theta^{\ast})}{N}-\Theta\frac{X_{N}^{\top}(y_{N}-X_{N}\theta^{\ast})}{N}\right\|_{\infty}$
$\displaystyle\leq\left\|\left(\widetilde{\Theta}\frac{X_{N}^{\top}X_{N}}{N}-I_{d}\right)(\bar{\theta}-\theta^{\ast})\right\|_{\infty}+\left\|\left(\widetilde{\Theta}-\Theta\right)\frac{X_{N}^{\top}e_{N}}{N}\right\|_{\infty}$
$\displaystyle\leq\left|\\!\left|\\!\left|{\widetilde{\Theta}\frac{X_{N}^{\top}X_{N}}{N}-I_{d}}\right|\\!\right|\\!\right|_{\max}\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}+\left|\\!\left|\\!\left|{\widetilde{\Theta}-\Theta}\right|\\!\right|\\!\right|_{\infty}\left\|\frac{X_{N}^{\top}e_{N}}{N}\right\|_{\infty},$

where we use the triangle inequality in the second to last inequality and the
fact that for any matrix $A$ and vector $a$ with compatible dimensions,
$\|Aa\|_{\infty}\leq\left|\\!\left|\\!\left|{A}\right|\\!\right|\\!\right|_{\max}\|a\|_{1}$
and
$\|Aa\|_{\infty}\leq\left|\\!\left|\\!\left|{A}\right|\\!\right|\\!\right|_{\infty}\|a\|_{\infty}$,
in the last inequality. Further applying the triangle inequality and the fact
that for any two matrices $A$ and $B$ with compatible dimensions,
$\left|\\!\left|\\!\left|{AB}\right|\\!\right|\\!\right|_{\max}\leq\left|\\!\left|\\!\left|{A}\right|\\!\right|\\!\right|_{\infty}\left|\\!\left|\\!\left|{B}\right|\\!\right|\\!\right|_{\max}$,
we have that

$\displaystyle\left|\\!\left|\\!\left|{\widetilde{\Theta}\frac{X_{N}^{\top}X_{N}}{N}-I_{d}}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle=\left|\\!\left|\\!\left|{\widetilde{\Theta}\frac{X_{N}^{\top}X_{N}}{N}-\widetilde{\Theta}\frac{X_{1}^{\top}X_{1}}{n}+\widetilde{\Theta}\frac{X_{1}^{\top}X_{1}}{n}-I_{d}}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle\leq\left|\\!\left|\\!\left|{\widetilde{\Theta}\left(\frac{X_{N}^{\top}X_{N}}{N}-\frac{X_{1}^{\top}X_{1}}{n}\right)}\right|\\!\right|\\!\right|_{\max}+\left|\\!\left|\\!\left|{\widetilde{\Theta}\frac{X_{1}^{\top}X_{1}}{n}-I_{d}}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle\leq\left|\\!\left|\\!\left|{\widetilde{\Theta}}\right|\\!\right|\\!\right|_{\infty}\left|\\!\left|\\!\left|{\frac{X_{N}^{\top}X_{N}}{N}-\frac{X_{1}^{\top}X_{1}}{n}}\right|\\!\right|\\!\right|_{\max}+\left|\\!\left|\\!\left|{\widetilde{\Theta}\frac{X_{1}^{\top}X_{1}}{n}-I_{d}}\right|\\!\right|\\!\right|_{\max}.$

Under Assumption • ‣ 3.2, $X_{N}$ has sub-Gaussian rows. Then, by Lemma 35, if
$n\gg{s^{*}}\log d$, we have that

$\left|\\!\left|\\!\left|{\widetilde{\Theta}}\right|\\!\right|\\!\right|_{\infty}=\max_{l}\left\|\widetilde{\Theta}_{l}\right\|_{1}=O_{P}\left(\sqrt{{s^{*}}}\right),$
$\left|\\!\left|\\!\left|{\widetilde{\Theta}\frac{X_{1}^{\top}X_{1}}{n}-I_{d}}\right|\\!\right|\\!\right|_{\max}=O_{P}\left(\sqrt{\frac{\log
d}{n}}\right),$

and

$\left|\\!\left|\\!\left|{\widetilde{\Theta}-\Theta}\right|\\!\right|\\!\right|_{\infty}=\max_{l}\left\|\widetilde{\Theta}_{l}-\Theta_{l}\right\|_{1}=O_{P}\left({s^{*}}\sqrt{\frac{\log
d}{n}}\right).$

It remains to bound
$\left|\\!\left|\\!\left|{\frac{X_{N}^{\top}X_{N}}{N}-\frac{X_{1}^{\top}X_{1}}{n}}\right|\\!\right|\\!\right|_{\max}$
and $\left\|\frac{X_{N}^{\top}e_{N}}{N}\right\|_{\infty}$.

Under Assumptions • ‣ 3.2, each $x_{ij,l}$ is sub-Gaussian, and therefore, the
product $x_{ij,l}x_{ij,l^{\prime}}$ of any two is sub-exponential. By
Bernstein’s inequality, we have that for any $t>0$,

$P\left(\left|\frac{(X_{N}^{\top}X_{N})_{l,l^{\prime}}}{N}-\Sigma_{l,l^{\prime}}\right|>t\right)\leq
2\exp\left(-cN\left(\frac{t^{2}}{\Sigma_{l,l^{\prime}}^{2}}\wedge\frac{t}{|\Sigma_{l,l^{\prime}}|}\right)\right),$

or for any $\delta\in(0,1)$,

$P\left(\left|\frac{(X_{N}^{\top}X_{N})_{l,l^{\prime}}}{N}-\Sigma_{l,l^{\prime}}\right|>|\Sigma_{l,l^{\prime}}|\left(\frac{\log\frac{2d^{2}}{\delta}}{cN}\vee\sqrt{\frac{\log\frac{2d^{2}}{\delta}}{cN}}\right)\right)\leq\frac{\delta}{d^{2}},$

for some constant $c>0$. Then, by the union bound, we have that

$\displaystyle
P\left(\left|\\!\left|\\!\left|{\frac{X_{N}^{\top}X_{N}}{N}-\Sigma}\right|\\!\right|\\!\right|_{\max}>\left|\\!\left|\\!\left|{\Sigma}\right|\\!\right|\\!\right|_{\max}\left(\frac{\log\frac{2d^{2}}{\delta}}{cN}\vee\sqrt{\frac{\log\frac{2d^{2}}{\delta}}{cN}}\right)\right)\leq\delta.$
(40)

Similarly, we have that

$\displaystyle
P\left(\left|\\!\left|\\!\left|{\frac{X_{1}^{\top}X_{1}}{n}-\Sigma}\right|\\!\right|\\!\right|_{\max}>\left|\\!\left|\\!\left|{\Sigma}\right|\\!\right|\\!\right|_{\max}\left(\frac{\log\frac{2d^{2}}{\delta}}{cn}\vee\sqrt{\frac{\log\frac{2d^{2}}{\delta}}{cn}}\right)\right)\leq\delta.$
(41)

Then, by the triangle inequality, we have that

$\displaystyle\left|\\!\left|\\!\left|{\frac{X_{N}^{\top}X_{N}}{N}-\frac{X_{1}^{\top}X_{1}}{n}}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle\leq\left|\\!\left|\\!\left|{\frac{X_{1}^{\top}X_{1}}{n}-\Sigma}\right|\\!\right|\\!\right|_{\max}+\left|\\!\left|\\!\left|{\frac{X_{N}^{\top}X_{N}}{N}-\Sigma}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle\lesssim\left|\\!\left|\\!\left|{\Sigma}\right|\\!\right|\\!\right|_{\max}\left(\frac{\log\frac{2d^{2}}{\delta}}{n}\vee\sqrt{\frac{\log\frac{2d^{2}}{\delta}}{n}}\right)$
$\displaystyle\lesssim\left(\frac{\log\frac{2d^{2}}{\delta}}{n}\vee\sqrt{\frac{\log\frac{2d^{2}}{\delta}}{n}}\right),$

with probability at least $1-\delta$, where we use
$\left|\\!\left|\\!\left|{\Sigma}\right|\\!\right|\\!\right|_{\max}\leq\left|\\!\left|\\!\left|{\Sigma}\right|\\!\right|\\!\right|_{2}=\lambda_{\tiny{\max}}(\Sigma)=O(1)$
under Assumption • ‣ 3.2. This implies that

$\left|\\!\left|\\!\left|{\frac{X_{N}^{\top}X_{N}}{N}-\frac{X_{1}^{\top}X_{1}}{n}}\right|\\!\right|\\!\right|_{\max}=O_{P}\left(\sqrt{\frac{\log
d}{n}}\right).$

Under Assumptions • ‣ 3.2 and • ‣ 3.2, each $x_{ij,l}$ and $e_{ij}$ are sub-
Gaussian, and therefore, their product $x_{ij,l}e_{ij}$ is sub-exponential.
Applying Bernstein’s inequality, we have that for any $\delta\in(0,1)$,

$P\left(\left|\frac{(X_{N}^{\top}e_{N})_{l}}{N}\right|>\sqrt{\Sigma_{l,l}}\sigma\left(\frac{\log\frac{2d}{\delta}}{cN}\vee\sqrt{\frac{\log\frac{2d}{\delta}}{cN}}\right)\right)\leq\frac{\delta}{d},$

for some constant $c>0$. Then, by the union bound, we have that

$\displaystyle
P\left(\left\|\frac{X_{N}^{\top}e_{N}}{N}\right\|_{\infty}>\max_{l}\sqrt{\Sigma_{l,l}}\sigma\left(\frac{\log\frac{2d}{\delta}}{cN}\vee\sqrt{\frac{\log\frac{2d}{\delta}}{cN}}\right)\right)\leq\delta,$
(42)

and then,

$\left\|\frac{X_{N}^{\top}e_{N}}{N}\right\|_{\infty}=O_{P}\left(\sqrt{\frac{\log
d}{N}}\right).$

Putting all the preceding bounds together, we obtain that

$\displaystyle\left\|\widetilde{\theta}-\theta^{\ast}+\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$
$\displaystyle\leq\left(\left|\\!\left|\\!\left|{\widetilde{\Theta}}\right|\\!\right|\\!\right|_{\infty}\left|\\!\left|\\!\left|{\frac{X_{N}^{\top}X_{N}}{N}-\frac{X_{1}^{\top}X_{1}}{n}}\right|\\!\right|\\!\right|_{\max}+\left|\\!\left|\\!\left|{\widetilde{\Theta}\frac{X_{1}^{\top}X_{1}}{n}-I_{d}}\right|\\!\right|\\!\right|_{\max}\right)\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}+\left|\\!\left|\\!\left|{\widetilde{\Theta}-\Theta}\right|\\!\right|\\!\right|_{\infty}\left\|\frac{X_{N}^{\top}e_{N}}{N}\right\|_{\infty}$
$\displaystyle=\left(O_{P}\left(\sqrt{{s^{*}}}\right)O_{P}\left(\sqrt{\frac{\log
d}{n}}\right)+O_{P}\left(\sqrt{\frac{\log
d}{n}}\right)\right)O_{P}(r_{\bar{\theta}})+O_{P}\left({s^{*}}\sqrt{\frac{\log
d}{n}}\right)O_{P}\left(\sqrt{\frac{\log d}{N}}\right)$
$\displaystyle=O_{P}\left(\sqrt{\frac{{s^{*}}\log
d}{n}}r_{\bar{\theta}}+\frac{{s^{*}}\log d}{n\sqrt{k}}\right),$

where we assume that
$\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}=O_{P}(r_{\bar{\theta}})$, and
hence,

$|T-T_{0}|=O_{P}\left(r_{\bar{\theta}}\sqrt{{s^{*}}k\log d}+\frac{{s^{*}}\log
d}{\sqrt{n}}\right).$

Choosing

$\zeta=\left(r_{\bar{\theta}}\sqrt{{s^{*}}k\log d}+\frac{{s^{*}}\log
d}{\sqrt{n}}\right)^{1-\kappa},$

with any $\kappa>0$, we deduce that

$P\left(|T-T_{0}|>\zeta\right)=o(1).$

We also have that

$\zeta\sqrt{1\vee\log\frac{d}{\zeta}}=o(1),$

provided that

$\left(r_{\bar{\theta}}\sqrt{{s^{*}}k\log d}+\frac{{s^{*}}\log
d}{\sqrt{n}}\right)\log^{1/2+\kappa}d=o(1),$

which holds if

$n\gg{s^{*}}^{2}\log^{3+\kappa}d,$

and

$r_{\bar{\theta}}\ll\frac{1}{\sqrt{k{s^{*}}}\log^{1+\kappa}d}.$

$\blacksquare$

###### Lemma 20

$\widehat{T}$ and $T_{0}$ are defined as in (20) and (23) respectively. In
sparse linear model, under Assumptions • ‣ 3.2 and • ‣ 3.2, provided that
$n\gg{s^{*}}\log d$, we have that

$|\widehat{T}-T_{0}|=O_{P}\left(\frac{\left(s_{0}\sqrt{s^{*}}+{s^{*}}\right)\log
d}{\sqrt{n}}\right).$

Moreover, if $n\gg\left(s_{0}^{2}{s^{*}}+{s^{*}}^{2}\right)\log^{3+\kappa}d$
and for some $\kappa>0$, then there exists some $\xi>0$ such that (32) holds.

Proof of Lemma 20.  By the proof of Lemma 19, we obtain that

$\displaystyle|\widehat{T}-T_{0}|$ $\displaystyle\leq\max_{1\leq l\leq
d}\left|\sqrt{N}(\widehat{\theta}-\theta^{\ast})_{l}+\sqrt{N}\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right)_{l}\right|$
$\displaystyle=\sqrt{N}\left\|\widehat{\theta}-\theta^{\ast}+\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$
$\displaystyle=\sqrt{N}\left\|\widehat{\theta}_{L}+\widetilde{\Theta}\frac{X_{N}^{\top}(y_{N}-X_{N}\widehat{\theta}_{L})}{N}-\theta^{\ast}-\Theta\frac{X_{N}^{\top}(y_{N}-X_{N}\theta^{\ast})}{N}\right\|_{\infty}$
$\displaystyle\leq\left|\\!\left|\\!\left|{\widetilde{\Theta}\frac{X_{N}^{\top}X_{N}}{N}-I_{d}}\right|\\!\right|\\!\right|_{\max}\left\|\widehat{\theta}_{L}-\theta^{\ast}\right\|_{1}+\left|\\!\left|\\!\left|{\widetilde{\Theta}-\Theta}\right|\\!\right|\\!\right|_{\infty}\left\|\frac{X_{N}^{\top}e_{N}}{N}\right\|_{\infty}$
$\displaystyle=O_{P}\left(\sqrt{{s^{*}}k\log
d}\right)\left\|\widehat{\theta}_{L}-\theta^{\ast}\right\|_{1}+O_{P}\left(\frac{{s^{*}}\log
d}{\sqrt{n}}\right).$

Since

$\left\|\widehat{\theta}_{L}-\theta^{\ast}\right\|_{1}=O_{P}\left(s_{0}\sqrt{\frac{\log
d}{N}}\right),$

we have that

$\displaystyle|\widehat{T}-T_{0}|=O_{P}\left(\frac{\left(s_{0}\sqrt{s^{*}}+{s^{*}}\right)\log
d}{\sqrt{n}}\right).$

Choosing

$\xi=\left(\frac{\left(s_{0}\sqrt{s^{*}}+{s^{*}}\right)\log
d}{\sqrt{n}}\right)^{1-\kappa},$

with any $\kappa>0$, we deduce that

$P\left(|\widehat{T}-T_{0}|>\xi\right)=o(1).$

We also have that

$\xi\sqrt{1\vee\log\frac{d}{\xi}}=o(1),$

provided that

$\left(\frac{\left(s_{0}\sqrt{s^{*}}+{s^{*}}\right)\log
d}{\sqrt{n}}\right)\log^{1/2+\kappa}d=o(1),$

which holds if

$n\gg\left(s_{0}^{2}{s^{*}}+{s^{*}}^{2}\right)\log^{3+\kappa}d.$

$\blacksquare$

###### Lemma 21

$T$ and $T_{0}$ are defined as in (7) and (23) respectively. In sparse GLM,
under Assumptions • ‣ 3.3 and • ‣ 3.3, provided that
$\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}=O_{P}(r_{\bar{\theta}})$ and
$n\gg s_{0}^{2}\log^{2}d+{s^{*}}^{2}\log d$, we have that

$|T-T_{0}|=O_{P}\left(r_{\bar{\theta}}\sqrt{{s^{*}}k\log d}+\frac{{s^{*}}\log
d}{\sqrt{n}}\right).$

Moreover, if $n\gg({s^{*}}^{2}+s_{0}^{2})\log^{3+\kappa}d$ and

$\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}\ll\min\left\\{\frac{1}{\sqrt{k{s^{*}}}s_{0}\log^{1+\kappa}d},\frac{1}{\left(nk{s^{*}}\log^{1+\kappa}d\right)^{1/4}}\right\\},$

for some $\kappa>0$, then there exists some $\zeta>0$ such that (29) holds.

Proof of Lemma 21.  Following the argument in the proof of Lemma 19, we have
that

$\displaystyle|T-T_{0}|$ $\displaystyle\leq\max_{1\leq l\leq
d}\left|\sqrt{N}(\widetilde{\theta}_{l}-\theta^{\ast}_{l})+\sqrt{N}\left(\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right)_{l}\right|$
$\displaystyle=\sqrt{N}\left\|\widetilde{\theta}-\theta^{\ast}+\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty},$

and

$\displaystyle\left\|\widetilde{\theta}-\theta^{\ast}+\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$
$\displaystyle=\left\|\bar{\theta}-\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla\mathcal{L}_{N}(\bar{\theta})-\theta^{\ast}+\Theta\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$
$\displaystyle=\left\|\bar{\theta}-\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla\mathcal{L}_{N}(\bar{\theta})-\theta^{\ast}+\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla\mathcal{L}_{N}(\theta^{\ast})-\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla\mathcal{L}_{N}(\theta^{\ast})+\Theta\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$
$\displaystyle\leq\left\|\bar{\theta}-\theta^{\ast}-\widetilde{\Theta}(\widetilde{\theta}^{(0)})\left(\nabla\mathcal{L}_{N}(\bar{\theta})-\nabla\mathcal{L}_{N}(\theta^{\ast})\right)\right\|_{\infty}+\left\|\left(\widetilde{\Theta}(\widetilde{\theta}^{(0)})-\Theta\right)\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}.$

By Taylor’s theorem, we have that

$\displaystyle\nabla\mathcal{L}_{N}(\bar{\theta})-\nabla\mathcal{L}_{N}(\theta^{\ast})=\int_{0}^{1}\nabla^{2}\mathcal{L}_{N}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast}))dt(\bar{\theta}-\theta^{\ast}),$
(43)

and then,

$\displaystyle\left\|\widetilde{\theta}-\theta^{\ast}+\nabla^{2}\mathcal{L}^{\ast}(\theta^{\ast})^{-1}\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$
$\displaystyle\leq\left\|\bar{\theta}-\theta^{\ast}-\widetilde{\Theta}(\widetilde{\theta}^{(0)})\int_{0}^{1}\nabla^{2}\mathcal{L}_{N}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast}))dt(\bar{\theta}-\theta^{\ast})\right\|_{\infty}+\left\|\left(\widetilde{\Theta}(\widetilde{\theta}^{(0)})-\Theta\right)\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$
$\displaystyle=\left\|\int_{0}^{1}\left(\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{N}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast}))-I_{d}\right)dt(\bar{\theta}-\theta^{\ast})\right\|_{\infty}+\left\|\left(\widetilde{\Theta}(\widetilde{\theta}^{(0)})-\Theta\right)\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}$
$\displaystyle\leq\int_{0}^{1}\left|\\!\left|\\!\left|{\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{N}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast}))-I_{d}}\right|\\!\right|\\!\right|_{\max}dt\left\|\bar{\theta}-\theta^{\ast}\right\|_{1}+\left|\\!\left|\\!\left|{\widetilde{\Theta}(\widetilde{\theta}^{(0)})-\Theta}\right|\\!\right|\\!\right|_{\infty}\left\|\nabla\mathcal{L}_{N}(\theta^{\ast})\right\|_{\infty}.$

By the triangle inequality, we have that

$\displaystyle\left|\\!\left|\\!\left|{\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{N}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast}))-I_{d}}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle=\bigg{|}\\!\bigg{|}\\!\bigg{|}\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{N}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast}))-\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{N}(\theta^{\ast})+\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{N}(\theta^{\ast})-\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{1}(\theta^{\ast})$
$\displaystyle\quad+\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{1}(\theta^{\ast})-\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{1}(\widetilde{\theta}^{(0)})+\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{1}(\widetilde{\theta}^{(0)})-I_{d}\bigg{|}\\!\bigg{|}\\!\bigg{|}_{\max}$
$\displaystyle\leq\left|\\!\left|\\!\left|{\widetilde{\Theta}(\widetilde{\theta}^{(0)})\left(\nabla^{2}\mathcal{L}_{N}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast}))-\nabla^{2}\mathcal{L}_{N}(\theta^{\ast})\right)}\right|\\!\right|\\!\right|_{\max}+\left|\\!\left|\\!\left|{\widetilde{\Theta}(\widetilde{\theta}^{(0)})\left(\nabla^{2}\mathcal{L}_{N}(\theta^{\ast})-\nabla^{2}\mathcal{L}_{1}(\theta^{\ast})\right)}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle\quad+\left|\\!\left|\\!\left|{\widetilde{\Theta}(\widetilde{\theta}^{(0)})\left(\nabla^{2}\mathcal{L}_{1}(\theta^{\ast})-\nabla^{2}\mathcal{L}_{1}(\widetilde{\theta}^{(0)})\right)}\right|\\!\right|\\!\right|_{\max}+\left|\\!\left|\\!\left|{\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{1}(\widetilde{\theta}^{(0)})-I_{d}}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle\leq\left|\\!\left|\\!\left|{\widetilde{\Theta}(\widetilde{\theta}^{(0)})}\right|\\!\right|\\!\right|_{\infty}\bigg{(}\left|\\!\left|\\!\left|{\nabla^{2}\mathcal{L}_{N}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast}))-\nabla^{2}\mathcal{L}_{N}(\theta^{\ast})}\right|\\!\right|\\!\right|_{\max}+\left|\\!\left|\\!\left|{\nabla^{2}\mathcal{L}_{N}(\theta^{\ast})-\nabla^{2}\mathcal{L}_{1}(\theta^{\ast})}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle\quad+\left|\\!\left|\\!\left|{\nabla^{2}\mathcal{L}_{1}(\theta^{\ast})-\nabla^{2}\mathcal{L}_{1}(\widetilde{\theta}^{(0)})}\right|\\!\right|\\!\right|_{\max}\bigg{)}+\left|\\!\left|\\!\left|{\widetilde{\Theta}(\widetilde{\theta}^{(0)})\nabla^{2}\mathcal{L}_{1}(\widetilde{\theta}^{(0)})-I_{d}}\right|\\!\right|\\!\right|_{\max}.$

Under Assumption • ‣ 3.3, we have by Taylor’s theorem that

$\displaystyle\left|g^{\prime\prime}(y_{ij},x_{ij}^{\top}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast})))-g^{\prime\prime}(y_{ij},x_{ij}^{\top}\theta^{\ast})\right|$
$\displaystyle=\left|\int_{0}^{1}g^{\prime\prime\prime}(y_{ij},x_{ij}^{\top}(\theta^{\ast}+st(\bar{\theta}-\theta^{\ast})))ds\cdot
tx_{ij}^{\top}(\bar{\theta}-\theta^{\ast})\right|$
$\displaystyle\lesssim\left|x_{ij}^{\top}(\bar{\theta}-\theta^{\ast})\right|,$

and then by the triangle inequality,

$\displaystyle\left|\\!\left|\\!\left|{\nabla^{2}\mathcal{L}_{N}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast}))-\nabla^{2}\mathcal{L}_{N}(\theta^{\ast})}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle=\left|\\!\left|\\!\left|{\frac{1}{N}\sum_{i=1}^{n}\sum_{j=1}^{k}x_{ij}x_{ij}^{\top}\left(g^{\prime\prime}(y_{ij},x_{ij}^{\top}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast})))-g^{\prime\prime}(y_{ij},x_{ij}^{\top}\theta^{\ast})\right)}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{n}\sum_{j=1}^{k}\left|\\!\left|\\!\left|{x_{ij}x_{ij}^{\top}\left(g^{\prime\prime}(y_{ij},x_{ij}^{\top}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast})))-g^{\prime\prime}(y_{ij},x_{ij}^{\top}\theta^{\ast})\right)}\right|\\!\right|\\!\right|_{\max}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{n}\sum_{j=1}^{k}\left|\\!\left|\\!\left|{x_{ij}x_{ij}^{\top}}\right|\\!\right|\\!\right|_{\max}\left|g^{\prime\prime}(y_{ij},x_{ij}^{\top}(\theta^{\ast}+t(\bar{\theta}-\theta^{\ast})))-g^{\prime\prime}(y_{ij},x_{ij}^{\top}\theta^{\ast})\right|$
$\displaystyle\lesssim\frac{1}{N}\sum_{i=1}^{n}\sum_{j=1}^{k}\|x_{ij}\|_{\infty}^{2}\left|x_{ij}^{\top}(\bar{\theta}-\theta^{\ast})\right|$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{n}\sum_{j=1}^{k}\|x_{ij}\|_{\infty}^{3}\|\bar{\theta}-\theta^{\ast}\|_{1}$
$\displaystyle\lesssim\left\|\bar{\theta}-\theta^{\ast}\right\|_{1},$ (44)

where we use that $\|x_{ij}\|_{\infty}=O(1)$ under Assumption • ‣ 3.3 in the
last inequality. Similarly, we have that

$\left|\\!\left|\\!\left|{\nabla^{2}\mathcal{L}_{1}(\theta^{\ast})-\nabla^{2}\mathcal{L}_{1}(\widetilde{\theta}^{(0)})}\right|\\!\right|\\!\right|_{\max}\lesssim\|\widetilde{\theta}^{(0)}-\theta^{\ast}\|_{1}=O_{P}\left(s_{0}\sqrt{\frac{\log
d}{n}}\right),$