# StarTrek: Combinatorial Variable Selection with False Discovery Rate Control

Lu Zhang and Junwei Lu Department of Statistics, Harvard University,
Cambridge, MA 02138; Email<EMAIL_ADDRESS>of
Biostatistics, Harvard School of Public Health, Boston, MA 02115; Email:
<EMAIL_ADDRESS>
###### Abstract

Variable selection on the large-scale networks has been extensively studied in
the literature. While most of the existing methods are limited to the local
functionals especially the graph edges, this paper focuses on selecting the
discrete hub structures of the networks. Specifically, we propose an
inferential method, called StarTrek filter, to select the hub nodes with
degrees larger than a certain thresholding level in the high dimensional
graphical models and control the false discovery rate (FDR). Discovering hub
nodes in the networks is challenging: there is no straightforward statistic
for testing the degree of a node due to the combinatorial structures;
complicated dependence in the multiple testing problem is hard to characterize
and control. In methodology, the StarTrek filter overcomes this by
constructing p-values based on the maximum test statistics via the Gaussian
multiplier bootstrap. In theory, we show that the StarTrek filter can control
the FDR by providing accurate bounds on the approximation errors of the
quantile estimation and addressing the dependence structures among the maximal
statistics. To this end, we establish novel Cramér-type comparison bounds for
the high dimensional Gaussian random vectors. Comparing to the Gaussian
comparison bound via the Kolmogorov distance established by Chernozhukov et
al. (2014), our Cramér-type comparison bounds establish the relative
difference between the distribution functions of two high dimensional Gaussian
random vectors, which is essential in the theoretical analysis of FDR control.
Moreover, the StarTrek filter can be applied to general statistical models for
FDR control of discovering discrete structures such as simultaneously testing
the sparsity levels of multiple high dimensional linear models. We illustrate
the validity of the StarTrek filter in a series of numerical experiments and
apply it to the genotype-tissue expression dataset to discover central
regulator genes.

Keywords. Graphical models, multiple testing, false discovery rate control,
combinatorial inference, Gaussian multiplier bootstrap, comparison bounds.

##  1 Introduction

Graphical models are widely used for real-world problems in a broad range of
fields, including social science, economics, genetics, and computational
neuroscience (Newman et al., 2002; Luscombe et al., 2004; Rubinov and Sporns,
2010). Scientists and practitioners aim to understand the underlying network
structure behind large-scale datasets. For a high-dimensional random vector
$\bm{X}=(\bm{X}_{1},\cdots,\bm{X}_{d})\in\mathbb{R}^{d}$, we let
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ be an undirected graph, which encodes
the conditional dependence structure among $\bm{X}$. Specifically, each
component of $\bm{X}$ corresponds to some vertex in
$\mathcal{V}=\\{1,2\cdots,d\\}$, and $(j,k)\notin\mathcal{E}$ if and only if
$\bm{X}_{j}$ and $\bm{X}_{k}$ are conditionally independent given the rest of
variables. Many existing works in the literature seek to learn the structure
of $\mathcal{G}$ via estimating the weight matrix $\bm{\Theta}$. For example,
Meinshausen and Bühlmann (2006); Yuan and Lin (2007); Friedman et al. (2008);
Rothman et al. (2008); Peng et al. (2009); Lam and Fan (2009); Ravikumar et
al. (2011); Cai et al. (2011); Shen et al. (2012) focus on estimating the
precision matrix in a Gaussian graphical model. Further, there is also a line
of work developing methodology and theory to assess the uncertainty of edge
estimation, i.e., constructing hypothesis tests and confidence intervals on
the network edges, see Cai and Ma (2013); Gu et al. (2015); Ren et al. (2015);
Cai and Zhang (2016); Janková and van de Geer (2017); Yang et al. (2018); Feng
and Ning (2019); Ding and Zhou (2020). Recently, simultaneously testing
multiple hypotheses on edges of the graphical models has received increasing
attention (Liu, 2013; Cai et al., 2013; Xia et al., 2015, 2018; Li and
Maathuis, 2019; Eisenach et al., 2020).

Most of the aforementioned works formulate the testing problems based on
continuous parameters and local properties. For example, Liu (2013) proposes a
method to select edges in Gaussian graphical models with asymptotic FDR
control guarantees. Testing the existence of edges concerns the local
structure of the graph. Under certain modeling assumptions, its null
hypothesis can be translated into a single point in the continuous parameter
space, for example, $\bm{\Theta}_{jk}=0$ where $\bm{\Theta}$ is the precision
matrix or the general weight matrix. However, for many scientific questions
involving network structures, we need to detect and infer discrete and
combinatorial signals in the networks, which does not follow from single edge
testing. For example, in the study of social networks, it is interesting to
discover active and impactful users, usually called “hub users,” as they are
connected to many other nodes in the social network (Ilyas et al., 2011; Lee
et al., 2019). In gene co-expression network analysis, identifying central
regulators/hub genes (Yuan et al., 2017; Liu et al., 2019b, a) is known to be
extremely useful to the study of progression and prognosis of certain cancers
and can support the treatment in the future. In neuroscience, researchers are
interested in identifying the cerebral areas which are intensively connected
to other regions (Shaw et al., 2008; van den Heuvel and Sporns, 2013; Power et
al., 2013) during certain cognitive processes. The discovery of such
central/hub areas can provide scientists better understanding of the
mechanisms of human cognition.

Motivated by these applications in various areas, in this paper, we consider
the hub node selection problem from the network models. In specific, given a
graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is the
vertex set and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the edge
set, we consider multiple hypotheses on whether the degree of some node
$j\in\mathcal{V}$ exceeds a given threshold $k_{\tau}$:

$H_{0j}:\text{degree of node }j<k_{\tau}\text{ v.s. }H_{1j}:\text{degree of
node }j\geq k_{\tau}.$

Throughout the paper, these nodes with large degrees will be called hub nodes.
For each $j\in[d]$, let $\psi_{j}=1$ if $H_{0j}$ is rejected and $\psi_{j}=0$
otherwise. When selecting hub nodes, we would like to control the false
discovery rate, as defined below:

${\rm
FDR}=\mathbb{E}\left[{\frac{\sum_{j\in\mathcal{H}_{0}}\psi_{j}}{\max\big{\\{}\sum_{j=1}^{d}\psi_{j},1\big{\\}}}}\right],$

where $\mathcal{H}_{0}=\\{j\mid\text{degree of node }j<k_{\tau}\\}$. Remark
the hypotheses $H_{0j},j\in[d]$ are not based on continuous parameters. They
instead involve the degrees of the nodes, which are intrinsically
discrete/combinatorial functionals. To the best of our knowledge, there is no
existing literature studying such combinatorial variable selection problems.
The most relevant work turns out to be Lu et al. (2017), which proposes a
general framework for inference about graph invariants/combinatorial
quantities on undirected graphical models. However, they study single
hypothesis testing and have to decide which subgraph to be tested before
running the procedure.

The combinatorial variable selection problems bring many new challenges.
First, most of the existing work focus on testing continuous parameters (Liu,
2013; Javanmard and Montanari, 2013; Javanmard and Montanari, 2014a, b;
Belloni et al., 2014; Van de Geer et al., 2014; Xia et al., 2015, 2018;
Javanmard and Javadi, 2019; Sur and Candès, 2019; Zhao et al., 2020). For
discrete functionals, it is more difficult to construct appropriate test
statistics and estimate its quantile accurately, especially in high
dimensions. Second, many multiple testing procedures rely on an independence
assumption (or certain dependence assumptions) on the null p-values (Benjamini
and Hochberg, 1995; Benjamini and Yekutieli, 2001; Benjamini, 2010). However,
the single hypothesis here is about the global property of the graph, which
means that any reasonable test statistic has to involve the whole graph.
Therefore, complicated dependence structures exist inevitably, which presents
another layer of difficulty for controlling the false discoveries. Now we
summarize the motivating question for this paper: how to develop a
combinatorial selection procedure to discover nodes with large degrees on a
graph with FDR control guarantees?

This paper introduces the StarTrek filter to select hub nodes. The filter is
based on the maximum statistics, whose quantiles are approximated by the
Gaussian multiplier bootstrap procedure. Briefly speaking, the Gaussian
multiplier bootstrap procedure estimates the distribution of a given maximum
statistic of general random vectors with unknown covariance matrices by the
distribution of the maximum of a sum of the conditional Gaussian random
vectors. The validity of high dimensional testing problems, such as family-
wise error rate (FWER) control, relies on the non-asymptotic bounds of the
Kolmogorov distance between the true distribution of the maximum statistics
and the Gaussian multiplier bootstrap approximation, which is established in
Chernozhukov et al. (2013). However, in order to control the FDR in the
context of combinatorial variable selection, a more refined characterization
of the quantile approximation errors is required. In specific, we need the so
called Cramér-type comparison bounds quantifying the accuracy of the p-values
in order to control the FDR in the simultaneous testing procedures (Chang et
al., 2016). In our context, consider two centered Gaussian random vectors
$U,V\in\mathbb{R}^{d}$ with different covariance matrices $\bm{\Sigma}^{U}$,
$\bm{\Sigma}^{V}$ and denote the $\ell_{\infty}$ norms of $U,V$ by
$||U||_{\infty},||V||_{\infty}$ respectively, then the Cramér-type comparison
bounds aim to control the relative error
$\left|\frac{\mathbb{P}(||U||_{\infty}>t)}{\mathbb{P}(||V||_{\infty}>t)}-1\right|$
for certain range of $t$. Comparing to the Kolmogorov distance
$\sup_{t\in\mathbb{R}}\left|{\mathbb{P}(||U||_{\infty}>t)}-{\mathbb{P}(||V||_{\infty}>t)}\right|$
(Chernozhukov et al., 2015), the Cramér-type comparison bound leads to the
relative error between two cumulative density functions, which is necessary to
guarantee the FDR control. In specific, we show in this paper a novel Cramér-
type Gaussian comparison bound

$\sup_{0\leq t\leq C_{0}\sqrt{\log
d}}\left|\frac{\mathbb{P}(||U||_{\infty}>t)}{\mathbb{P}(||V||_{\infty}>t)}-1\right|=O\left(\min\Big{\\{}(\log
d)^{5/2}\Delta_{\infty}^{1/2},\frac{\Delta_{0}\log
d}{\mathfrak{p}}\Big{\\}}\right),$ (1.1)

for some constant $C_{0}>0$, where
$\Delta_{\infty}:=||\bm{\Sigma}^{U}-\bm{\Sigma}^{V}||_{\max}$ is the entrywise
maximum norm difference between the two covariance matrices,
$\Delta_{0}:=||\bm{\Sigma}^{U}-\bm{\Sigma}^{V}||_{0}$ with
$\left\|\cdot\right\|_{0}$ being the entrywise $\ell_{0}$-norm of the matrix,
and $\mathfrak{p}$ is the number of connected subgraphs in the graph whose
edge set $\mathcal{E}=\\{(j,k):\bm{\Sigma}^{U}_{jk}\neq 0\text{ or
}\bm{\Sigma}^{V}_{jk}\neq 0\\}$. This comparison bound in (1.1) characterizes
the relative errors between Gaussian maxima via two types of rates: the
$\ell_{\infty}$-norm $\Delta_{\infty}$ and the $\ell_{0}$-norm $\Delta_{0}$.
This implies a new insight that the Cramér type bound between two Gaussian
maxima is small as long as either their covariance matrices are uniformly
close or only sparse entries of the two covariance matrices differ. As far as
we know, the second type of rate in (1.1) has not been developed even in
Kolmogorov distance results of high dimensional Gaussian maxima. In the study
of FDR control, we need both types of rates: the $\Delta_{\infty}$ rate is
used to show that the Gaussian multiplier bootstrap procedure is an accurate
approximation for the maximum statistic quantiles and the $\Delta_{0}$ rate is
used to quantify the complicated dependence structure of the p-values for the
single tests on the degree of graph nodes. In order to prove the Cramér-type
comparison bound in (1.1), we develop two novel theoretic techniques to prove
the two types of rates separately. For the $\Delta_{\infty}$ rate, we
reformulate the Slepian’s interpolation (Slepian, 1962) into an ordinary
differential inequality such that the relative error can be controlled via the
Grönwall’s inequality (Grönwall, 1919). To control the $\Delta_{0}$ rate, the
anti-concentration inequality of Gaussian maxima developed in Chernozhukov et
al. (2015) is no longer sufficient, we establish a new type of anti-
concentration inequality for the derivatives of the soft-max of high
dimensional Gaussian vectors. The existing works on the Cramér type comparison
bounds such as Liu and Shao (2010, 2014); Chang et al. (2016) does not cover
the high dimensional maximum statistics. Therefore, their techniques can not
be directly extended to our case. To the best of our knowledge, it is the
first time in our paper to prove the Cramér-type Gaussian comparison bounds
(1.1) for high dimensional Gaussian maxima.

In summary, our paper makes the following major contributions. First, we
develop a novel StarTrek filter to select combinatorial statistical signals:
the hub nodes with the FDR control. This procedure involves maximum statistic
and Gaussian multiplier bootstrap for quantile estimation. Second, in theory,
the proposed method is shown to be valid for many different models with the
network structures. In this paper, we provide two examples, the Gaussian
graphical model and the bipartite network in the multiple linear models.
Third, we prove a new Cramér-type Gaussian comparison bound with two types of
rates: the maximum norm difference and $\ell_{0}$ norm difference. These
results are quite generic and has its own significance in the probability
theory.

###  1.1 Related work

Canonical approaches to FDR control and multiple testing (Benjamini and
Hochberg, 1995; Benjamini and Yekutieli, 2001; Benjamini, 2010) require that
valid p-values are available, and they only allow for certain forms of
dependence between these p-values. However, obtaining asymptotic p-values with
sufficient accuracy is generally non-trivial for high dimensional hypothesis
testing problems concerning continuous parameters (Javanmard and Montanari,
2013; Javanmard and Montanari, 2014a, b; Belloni et al., 2014; Van de Geer et
al., 2014; Sur and Candès, 2019; Zhao et al., 2020), not even to mention
discrete/combinatorial functionals.

Recently, there is a line of work conducting variable selection without
needing to act on a set of valid p-values, including Barber and Candès (2015,
2019); Candès et al. (2018); Xing et al. (2019); Dai et al. (2020a, b). These
approaches take advantage of the symmetry of the null test statistics and
establish FDR control guarantee. As their single hypothesis is often
formulated as conditional independence testing, it is challenging to apply
those techniques to select discrete signals for the problem studied in this
paper.

Another line of work develops multiple testing procedures based on asymptotic
p-values for specific high dimensional models (Liu, 2013; Liu and Luo, 2014;
Javanmard and Javadi, 2019; Xia et al., 2015, 2018; Liu et al., 2020). Among
them, Liu (2013) studies the edge selection problem on Gaussian graphical
models, which turns out to be the most relevant work to our paper. However,
their single hypothesis is about the local property of the graph. Our problem
of discovering nodes with large degrees concerns the global property of the
whole network, therefore requiring far more work.

There exists some recent work inferring combinatorial functionals. For
example, the method proposed in Ke et al. (2020) provides a confidence
interval for the number of spiked eigenvalues in a covariance matrix. Jin et
al. (2020) focuses on estimating the number of communities in a network and
yields confidence lower bounds. Neykov et al. (2019); Lu et al. (2017) propose
a general framework for conducting inference on graph invariants/combinatorial
quantities, such as the maximum degree, the negative number of connected
subgraphs, and the size of the longest chain of a given graph. Shen and Lu
(2020) develops methods for testing the general community combinatorial
properties of the stochastic block model. Regarding the hypothesis testing
problem, all these works only deal with a single hypothesis and establish
asymptotic type-I error rate control. While simultaneously testing those
combinatorial hypotheses is also very interesting and naturally arises from
many practical problems.

###  1.2 Outline

In Section 2, we set up the general testing framework and introduce the
StarTrek filter for selecting hub nodes. In Section 3, we present our core
probabilistic tools: Cramér-type Gaussian comparison bounds in terms of
maximum norm difference and $\ell_{0}$ norm difference. To offer a relatively
simpler illustration of our generic theoretical results, we first consider the
hub selection problem on a bipartite network (multitask regression with linear
models). Specifically, the input of the general StarTrek filter is chosen to
be the estimators and quantile estimates described in Section 4. Applying the
probabilistic results under this model, we establish FDR control guarantees
under certain conditions. Then we move to the Gaussian graphical model in
Section 5. In Section 6, we demonstrate StarTrek’s performance through
empirical simulations and a real data application.

###  1.3 Notations

Let $\phi(x),\Phi(x)$ be the probability density function (PDF) and the
cumulative distribution function (CDF) respectively of the standard Gaussian
distribution and denote $\bar{\Phi}(x)=1-\Phi(x)$. Let $\mathbf{1}_{d}$ be the
vector of ones of dimension $d$. We use $\mathds{1}{(\cdot)}$ to denote the
indicator function of a set and $|\cdot|$ to denote the cardinality of a set.
For two sets $A$ and $B$, denote their symmetric difference by $A\ominus B$,
i.e., $A\ominus B=(A\setminus B)\cup(B\setminus A)$; let $A\times B$ be the
Cartesian product. For two positive sequences $\\{x_{n}\\}_{n=1}^{\infty}$ and
$\\{y_{n}\\}_{n=1}^{\infty}$, we say $x_{n}=O\left(y_{n}\right)$ if $x_{n}\leq
Cy_{n}$ holds for any $n$ with some large enough $C>0$. And we say
$x_{n}=o\left(y_{n}\right)$ if $x_{n}/y_{n}\rightarrow 0$ as
$n\rightarrow\infty$. For a sequence of random variables
$\\{X_{n}\\}_{n=1}^{\infty}$ and a scalar $a$, we say $X_{n}\leq
a+o_{\mathbb{P}}(1)$ if for all $\epsilon>0$,
$\lim_{n\rightarrow\infty}\mathbb{P}\left({X_{n}-a>\epsilon}\right)=0$. Let
$[d]$ denote the set $\\{1,\dots,d\\}$. The $\ell_{\infty}$ norm and the
$\ell_{1}$ norm on $\mathbb{R}^{d}$ are denoted by $||\cdot||_{\infty}$ and
$||\cdot||_{1}$ respectively. For a random vector $X$, let $||X||_{\infty}$ be
its $\ell_{\infty}$ norm. For a matrix $\mathbf{A}\in\mathbb{R}^{d_{1}\times
d_{2}}$, we denote its minimal and maximal eigenvalues by
$\lambda_{\min}(\mathbf{A}),\lambda_{\max}(\mathbf{A})$ respectively, the
elementwise max norm by
$\left\|\mathbf{A}\right\|_{\max}=\max_{i\in[d_{1}],j\in[d_{2}]}|\mathbf{A}_{ij}|$
and the elementwise $\ell_{0}$ norm by
$\left\|\mathbf{A}\right\|_{0}=\sum_{i\in[d_{1}],j\in[d_{2}]}\mathds{1}{(\mathbf{A}_{ij}\neq
0)}$. Throughout this paper,
$C,C^{\prime},C^{\prime\prime},C_{0},C_{1},C_{2},\dots$ are used as generic
constants whose values may vary across different places.

##  2 Methodology

Before introducing our method, we set up the problem with more details.
Specifically, we consider a graph
$\mathcal{G}=(\mathcal{V}_{1},\mathcal{V}_{2},\mathcal{E})$ with the node sets
$\mathcal{V}_{1},\mathcal{V}_{2}$ and the edge set $\mathcal{E}$. Let
$d_{1}=|\mathcal{V}_{1}|$, $d_{2}=|\mathcal{V}_{2}|$ and denote its weight
matrix by $\bm{\Theta}\in\mathbb{R}^{d_{1}\times d_{2}}$. In the undirected
graph where $\mathcal{V}_{1}=\mathcal{V}_{2}:=\mathcal{V}$, $\bm{\Theta}$ is a
square matrix and its element $\bm{\Theta}_{jk}$ is nonzero when there is an
edge between node $j$ and node $k$, zero when there is no edge. In a bipartite
graph where $\mathcal{V}_{1}\neq\mathcal{V}_{2}$, elements of $\bm{\Theta}$
describe the existence of an edge between node j in $\mathcal{V}_{1}$ and node
$k$ in $\mathcal{V}_{2}$. Without loss of generality, we focus on one of the
node sets and denote it by $\mathcal{V}$ with $|\mathcal{V}|:=d$. We would
like to select those nodes among $\mathcal{V}$ whose degree exceeds a certain
threshold $k_{\tau}$. And the selection problem is equivalent to
simultaneously testing $d$ hypotheses:

$H_{0j}:\text{degree of node }j<k_{\tau}\text{ v.s. }H_{1j}:\text{degree of
node }j\geq k_{\tau},$ (2.1)

for $j\in[d]$. Let $\psi_{j}=1$ if $H_{0j}$ is rejected and $\psi_{j}=0$
otherwise, then for some multiple testing procedure with output
$\\{\psi_{j}\\}_{j\in[d]}$, the false discovery proportion (FDP) and FDR can
be defined as below:

${\rm
FDP}=\frac{\sum_{j\in\mathcal{H}_{0}}^{d}\psi_{j}}{\max\left\\{1,\sum_{j=1}^{d}\psi_{j}\right\\}},\quad{\rm
FDR}:=\mathbb{E}[{\rm FDP}],$

where $\mathcal{H}_{0}=\\{j\mid\text{degree of node }j<k_{\tau}\\}$. We aim to
propose a multiple testing procedure such that the FDP or FDR can be
controlled at a given level $0<q<1$.

We illustrate the above general setup in two specific examples. In multitask
regression with linear models, we are working with the bipartite graph case,
then the weight matrix $\bm{\Theta}$ corresponds to the parameter matrix whose
row represents the linear coefficients for one given response variable. Given
a threshold $k_{\tau}$, we want to select those rows (response variables) with
$\ell_{0}$ norm being at least $k_{\tau}$. In the context of Gaussian
graphical models where $\mathcal{V}_{1}=\mathcal{V}_{2}$, $\bm{\Theta}$
represents the precision matrix, and we want to select those hub nodes i.e.,
whose degree is larger than or equal to $k_{\tau}$.

###  2.1 StarTrek filter

Letting $\bm{\Theta}_{j}$ be the $j$-th row of $\bm{\Theta}$ and
$\bm{\Theta}_{j,-j}$ be the vector $\bm{\Theta}_{j}$ excluding its $j$-th
element, we formulate the testing problem for each single node as below,

$H_{0j}:\|\bm{\Theta}_{j,-j}\|_{0}<k_{\tau}\text{ v.s.
}H_{1j}:\|\bm{\Theta}_{j,-j}\|_{0}\geq k_{\tau}.$

To test the above hypothesis, we need some estimator of the weight matrix
$\bm{\Theta}$. In Gaussian graphical model, it is natural to use the estimator
of a precision matrix. In the bipartite graph (multiple response model),
estimated parameter matrix will suffice. Denote this generic estimator by
$\widetilde{\bm{\Theta}}$ (without causing confusion in notation), the maximum
test statistic over a given subset $E$ of $\mathcal{V}\times\mathcal{V}$ will
be

$T_{E}:=\max_{(j,k)\in E}\sqrt{n}\left|\widetilde{\bm{\Theta}}_{jk}\right|$

and its quantile is defined as
${c}(\alpha,E)=\inf\left\\{t\in\mathbb{R}\;|\;\mathbb{P}\left(T_{E}\leq
t\right)\geq 1-\alpha\right\\}$, which is often unknown. Assume it can be
estimated by $\widehat{c}(\alpha,E)$ from some procedure such as the Gaussian
multiplier bootstrap, a generic method called skip-down procedure can be used,
which was originally proposed in Lu et al. (2017) for testing a family of
monotone graph invariants. When applied to the specific degree testing
problem, it leads to the following algorithm.

Algorithm 1 Skip-down Method in Lu et al. (2017) (for testing the degree of
node $j$)

Input: $\\{\widetilde{\bm{\Theta}}_{e}\\}_{e\in\mathcal{V}\times\mathcal{V}}$,
significance level $\alpha$.

Initialize $t=0,E_{0}=\\{(j,k):k\in[d],k\neq j\\}$.

repeat

$t\leftarrow t+1$;

Select the rejected edges $\mathcal{R}\leftarrow\\{(j,k)\in
E_{t-1}\mid\sqrt{n}|\widetilde{\bm{\Theta}}_{jk}|>\widehat{c}(\alpha,E_{t-1})\\}$;

$E_{t}\leftarrow E_{t-1}\backslash\mathcal{R}$;

until $|E_{t}^{c}|\geq k$ or $\mathcal{R}=\varnothing$

Output: $\psi_{j,\alpha}=1$ if $|E_{t}^{c}|\geq k$ and $\psi_{j,\alpha}=0$
otherwise.

To conduct the node selection over the whole graph, we need to determine an
appropriate threshold $\widehat{\alpha}$ then reject $H_{0j}$ if
$\psi_{j,\widehat{\alpha}}=1$. A desirable choice of $\widehat{\alpha}$ should
be able to discover as many as hub nodes with the FDR remaining controlled
under the nominal level $q$. For example, if the BHq procedure is considered,
$\widehat{\alpha}$ can be defined as follows:

$\widehat{\alpha}=\sup\left\\{\alpha\in(0,1):\frac{\alpha
d}{\max\left\\{1,\sum_{j\in[d]}{\psi_{j,\alpha}}\right\\}}\leq q\right\\}.$
(2.2)

The above range of $\alpha$ is $(0,1)$, it will be very computationally
expensive if we do an exhaustive search since for each $\alpha$, we have to
recompute the quantiles $\widehat{c}(\alpha,E)$ for a lot of sets $E$.

We overcome the computational difficulty and propose a efficient procedure
called StarTrek filter, which is presented in Algorithm 2.

Algorithm 2 StarTrek Filter

Input: $\\{\widetilde{\bm{\Theta}}_{e}\\}_{e\in\mathcal{V}\times\mathcal{V}}$,
nominal FDR level $q$.

for $j\in[d]$ do

We order the elements in $\\{|\widetilde{\bm{\Theta}}_{j\ell}|:\ell\neq j\\}$
as
$|\widetilde{\bm{\Theta}}_{j,(1)}|\geq|\widetilde{\bm{\Theta}}_{j,(2)}|\geq\ldots\geq|\widetilde{\bm{\Theta}}_{j,(d-1)}|,$
where $|\widetilde{\bm{\Theta}}_{j,(\ell)}|$ is the $\ell$th largest entry.
Compute $\alpha_{j}=\max_{1\leq s\leq
k_{\tau}}\widehat{c}^{-1}(\sqrt{n}|\widetilde{\bm{\Theta}}_{j,(s)}|,E^{(s)}_{j})$
where $E^{(s)}_{j}:=\\{(j,\ell):\ell\neq
j,|\widetilde{\bm{\Theta}}_{j\ell}|\leq|\widetilde{\bm{\Theta}}_{j,(s)}|\\}$.

end for

Order $\alpha_{j}$ as $\alpha_{(1)}\leq\alpha_{(2)}\leq\dots\leq\alpha_{(d)}$
and set $\alpha_{(0)}=0$, let $j_{\max}=\max\\{0\leq j\leq
d:\alpha_{(j)}\leq{qj}/{d}\\}$.

Output: $S=\\{j:\alpha_{j}\leq\alpha_{(j_{\max})}\\}$ if $j_{\max}>0$;
$S=\varnothing$ otherwise.

Remark it only involves estimating $k_{\tau}$ different quantiles of some
maximum statistics per node, which is more efficient than the Skip-down
procedure (Lu et al., 2017) in terms of computation.

###  2.2 Accuracy of approximate quantiles

Before diving into the theoretical results, we pause to give specific forms of
the estimator of $\bm{\Theta}$ and how to compute the estimated quantiles of
the maximum statistic. Take the Gaussian graphical model as an example,
suppose that
$\bm{X}_{1},\ldots,\bm{X}_{n}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}N_{d}(0,\bm{\Sigma})$.
Let $\bm{\Theta}=\bm{\Sigma}^{-1}$, which will have the same $\ell_{0}$
elementwise norm as the adjacency matrix $\bm{\Theta}$. Denote
$\mathbf{e}_{k}$ be the $k$th canonical basis in $\mathbb{R}^{d}$, we consider
the following one-step estimator of $\bm{\Theta}_{jk}$,

$\widehat{\bm{\Theta}}^{\text{d}}_{jk}:=\widehat{\bm{\Theta}}_{jk}-\frac{\widehat{\bm{\Theta}}_{j}^{\top}\left(\widehat{\bm{\Sigma}}\widehat{\bm{\Theta}}_{k}-\mathbf{e}_{k}\right)}{\widehat{\bm{\Theta}}_{j}^{\top}\widehat{\bm{\Sigma}}_{j}},$
(2.3)

where $\widehat{\bm{\Theta}}$ could be either the graphical Lasso (GLasso)
estimator (Friedman et al., 2008) or the CLIME estimator (Cai et al., 2011).
Let
$\widetilde{\bm{\Theta}}^{\text{d}}_{jk}:={\widehat{\bm{\Theta}}^{\text{d}}_{jk}}/{\sqrt{\widehat{\bm{\Theta}}^{\text{d}}_{jj}\widehat{\bm{\Theta}}^{\text{d}}_{kk}}}$
and the standardized version
$\\{\widetilde{\bm{\Theta}}^{\text{d}}_{e}\\}_{\mathcal{V}\times\mathcal{V}}$
will be the input
$\\{\widetilde{\bm{\Theta}}_{e}\\}_{\mathcal{V}\times\mathcal{V}}$ of
Algorithm 2. Then the maximum test statistics (over the subset $E$) is defined
as $T_{E}=\underset{(j,k)\in
E}{\max}\;\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{jk}|$. To estimate its
quantile, we construct the following Gaussian multiplier bootstrap

$T^{\mathcal{B}}_{E}:=\underset{(j,k)\in
E}{\max}\;\frac{1}{\sqrt{n~{}{\widehat{\bm{\Theta}}_{jj}\widehat{\bm{\Theta}}_{kk}}}}\bigg{|}\sum_{i=1}^{n}{\widehat{\bm{\Theta}}^{\top}_{j}\left(\bm{X}_{i}\bm{X}_{i}^{\top}\widehat{\bm{\Theta}}_{k}-\mathbf{e}_{k}\right)}\xi_{i}\bigg{|},$
(2.4)

where $\xi_{i}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}N(0,1)$, which
produces
$\widehat{c}(\alpha,E)=\inf\left\\{t\in\mathbb{R}:\mathbb{P}_{\xi}\left(T^{\mathcal{B}}_{E}\leq
t\right)\geq 1-\alpha\right\\}$ as the quantile estimate. We also denote the
standardized true precision matrix
$(\bm{\Theta}_{jk}/\sqrt{\bm{\Theta}_{jj}\bm{\Theta}_{kk}})_{j,k\in[d]}$ by
$\bm{\Theta}^{\star}$. The theoretical results for Gaussian multiplier
bootstrap developed in Chernozhukov et al. (2013) basically imply the above
quantile estimates are accurate in the following sense:

######  Lemma 2.1.

Suppose $\bm{\Theta}\in\mathcal{U}(M,s,r_{0})$ and
$(\log(dn))^{7}/n+s^{2}(\log dn)^{4}/{n}=o(1)$, for any edge set
$E\subseteq\mathcal{V}\times\mathcal{V}$, we have

$\lim_{(n,d)\rightarrow\infty}\sup_{\bm{\Theta}\in\mathcal{U}(M,s,r_{0})}\sup_{\alpha\in(0,1)}\left|\mathbb{P}\left(\max_{e\in
E}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{e}-\bm{\Theta}^{\star}_{e}|>\widehat{c}(\alpha,E)\right)-\alpha\right|=0.$
(2.5)

where $\widetilde{\bm{\Theta}}^{\text{d}}_{e}$ is the standardized version of
the one-step estimator (2.3).

Note that $\mathcal{U}(M,s,r_{0})$ denotes the parameter space of precision
matrices and is defined as below:

$\displaystyle\mathcal{U}(M,s,r_{0})$
$\displaystyle=\Big{\\{}\bm{\Theta}\in\mathbb{R}^{d\times
d}\,\big{|}\,\lambda_{\min}(\bm{\Theta})\geq
1/r_{0},\lambda_{\max}(\bm{\Theta})\leq
r_{0},\max_{j\in[d]}\|\bm{\Theta}_{j}\|_{0}\leq s,\|\bm{\Theta}\|_{1}\leq
M\Big{\\}}.$

The proof of Lemma 2.1 can be found in Appendix D.2. However, Lemma 2.1 is not
sufficient for our multiple testing problem. Generally speaking, the
probabilistic bounds in Chernozhukov et al. (2013) are in terms of Kolmogorov
distance, which only provides a uniform characterization for the deviation
behaviors. Their results can be used to establish FWER control for global
testing problems based on the maximum test statistics. However, in order to
establish FDR control, we have to show that the estimation of number of false
discoveries is sufficiently accurate enough in the following sense, i.e.,
uniformly over certain range of $\alpha$,

$\frac{\alpha d_{0}}{\sum_{j\in\mathcal{H}_{0}}\psi_{j,\alpha}}\rightarrow
1,\quad\text{ in probability}$

where $\mathcal{H}_{0}=\\{j:\|\bm{\Theta}_{j,-j}\|_{0}<k_{\tau}\\}$. The above
result is different from the one needed for FWER control:
$\mathbb{E}\left[{{\psi_{j,\alpha}}}\right]=\alpha+o(1),j\in\mathcal{H}_{0}$.
In the context of our node selection problem, it can be reduced to the
following,

$\left|\frac{\sum_{j\in\mathcal{H}_{0}}\mathds{1}{(\max_{e\in
E}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{e}-\bm{\Theta}^{\star}_{e}|\geq\widehat{c}(\alpha,E))}}{d_{0}\alpha}-1\right|\rightarrow
0~{}~{}~{}\text{in probability}$

uniformly over certain range of $\alpha$ for some subset $E$. The above ratio
is closedly related to the ratio in Cramér-type moderation deviation results
(Liu and Shao, 2010, 2014; Liu, 2013). To this end, we establish the Cramér-
type deviation bounds for the Gaussian multiplier bootstrap procedure. This
type of results is built on two types of Cramér-type Gaussian comparison
bounds, which are presented in Section 3.

##  3 Cramér-type comparison bounds for Gaussian maxima

In this section, we present the theoretic results on the Cramér-type
comparison bounds for Gaussian maxima. Let $U,V\in\mathbb{R}^{d}$ be two
centered Gaussian random vectors with different covariance matrices
$\bm{\Sigma}^{U}=(\sigma_{jk}^{U})_{1\leq j,k\leq
d},\bm{\Sigma}^{V}=(\sigma_{jk}^{V})_{1\leq j,k\leq d}$. Recall that the
maximal difference of the covariance matrices is
$\Delta_{\infty}:=||\bm{\Sigma}^{U}-\bm{\Sigma}^{V}||_{\max}$ and the
elementwise $\ell_{0}$ norm difference of the covariance matrices is denoted
by
$\Delta_{0}:=\left\|\bm{\Sigma}^{U}-\bm{\Sigma}^{V}\right\|_{0}=\sum_{j,k\in[d]}\mathds{1}{(\sigma^{U}_{jk}\neq\sigma^{U}_{jk})}$.
The Gaussian maxima of $U$ and $V$ are denoted as $||U||_{\infty}$ and
$||V||_{\infty}$. Now we present a Cramér-type comparison bound (CCB) between
Gaussian maxima in terms of the maximum norm difference $\Delta_{\infty}$.

######  Theorem 3.1 (CCB with maximum norm difference).

Suppose $(\log d)^{5}\Delta_{\infty}=O(1)$, then we have

$\sup_{0\leq t\leq C_{0}\sqrt{\log
d}}\left|\frac{\mathbb{P}(||U||_{\infty}>t)}{\mathbb{P}(||V||_{\infty}>t)}-1\right|=O\left((\log
d)^{5/2}\Delta_{\infty}^{1/2}\right),$ (3.1)

for some constant $C_{0}>0$.

######  Remark 3.1.

We can actually prove a more general form (see Theorem B.2 in the appendix) of
the upper bound on the above term, without the assumption on
$\Delta_{\infty}$. In fact, we bound the right hand side of (3.1) as
$M_{3}(\log d)^{3/2}A(\Delta_{\infty})e^{M_{3}(\log
d)^{3/2}A(\Delta_{\infty})},$ where $A(\Delta_{\infty})=M_{1}\log
d\Delta_{\infty}^{1/2}\exp{(M_{2}\log^{2}d\Delta_{\infty}^{1/2})}$ with the
constants $M_{1},M_{2}$ only depending on the variance terms $\min_{1\leq
j\leq d}\\{\sigma^{U}_{jj},\sigma^{V}_{jj}\\},\max_{1\leq j\leq
d}\\{\sigma^{U}_{jj},\sigma^{V}_{jj}\\}$ and $M_{3}$ being a universal
constant.

When applying Theorem 3.1 to Gaussian multiplier bootstrap,
$||\Delta||_{\infty}$ actually controls the maximum differences between the
true covariance matrix and the empirical covariance matrix. Based on the bound
of $||\Delta||_{\infty}$, we can show that the Cramér-type comparison bound in
(3.1) will be $O((\log d)^{3/2}n^{-1/4})$ with high probability.

The proof can be found in Appendix B.1. The above result bounds the relative
difference between the distribution functions of the two Gaussian maxima.
Compared with the bound in terms of Kolmogorov distance, it has more refined
characterization when $t$ is large, which benefits from our iterative use of
the Slepian interpolation. We denote the interpolation between $U$ and $V$ as
$W(s)=\sqrt{s}U+\sqrt{1-s}V,s\in[0,1]$ and let
$Q_{t}(s)=\mathbb{P}(||W(s)||_{\infty}>t)$. Existing results (Chernozhukov et
al., 2013, 2014) quantify the difference between $Q_{t}(1)$ and $Q_{t}(0)$
uniformly over $t\in\mathbb{R}$, which leads to a bound on the Kolmogorov
distance between Gaussian maxima. Our main innovation is to consider
$R_{t}(s)=Q_{t}(s)/Q_{t}(0)-1$ and show that for any given $t$,
$\mathcal{R}_{t}:s\in[0,1]\mapsto|R_{t}(s)|$ is a contraction mapping with $0$
being its fixed point. Specifically, we have the following upper bound on
$|R_{t}(s)|$,

$|R_{t}(s)|\leq AB\int_{0}^{s}|R_{t}(\mu)|d\mu+AB\cdot s+A,$

where $AB$ and $A$ can be controlled via the bound on the maximal difference
of the covariance matrices $\Delta_{\infty}$ and converge to $0$ under certain
conditions. By Grönwall’s inequality (Grönwall, 1919), we then derive the
bound on $R_{t}(1)$ explicitly in terms of $A$ and $B$, which finally lead to
the desired Cramér-type comparison bound in (3.1).

The above theorem is a key ingredient for deriving Cramér-type deviation
results for the Gaussian multiplier bootstrap procedure. However, in certain
situations, comparison bounds in terms of maximum norm difference may not be
appropriate. There exist cases where the covariance matrices of two Gaussian
random vectors are not uniformly closed to each other, but have lots of
identical entries. In particular, for the combinatorial variable selection
problem in this paper, there exist complicated dependence structures between
the maximum statistic for different nodes, since each time when the degree of
one single node is tested, the statistic is computed based on the whole graph.
Again, this highlights the challenge of the multiple testing problem in our
paper. To establish FDR control, we need to deal with the dependence between
the maximum statistic of pairs of non-hub nodes. By the definition of non-hub
nodes, the covariance matrix difference between each pair of the involving
Gaussian vectors actually has lots of zero entries. We would like to take
advantage of this sparsity pattern when applying the comparison bound.
However, the bound in (3.1) is not sharp when $\Delta_{\infty}$ is not
negligible but $\Delta_{0}$ is small. To this end, we develop a different
version of the Cramér-type comparison bound as below.

######  Theorem 3.2 (CCB with elementwise $\ell_{0}$-norm difference).

Assume the Gaussian random vectors $U$ and $V$ have unit variances, i.e.,
$\sigma^{U}_{jj}=\sigma^{V}_{jj}=1,j\in[d]$ and there exists some
$\sigma_{0}<1$ such that
$|\sigma^{V}_{jk}|\leq\sigma_{0},|\sigma^{U}_{jk}|\leq\sigma_{0}$ for any
$j\neq k$. Suppose there exists a disjoint $\mathfrak{p}$-partition of nodes
$\cup_{\ell=1}^{\mathfrak{p}}\mathcal{C}_{\ell}=[d]$ such that
$\sigma^{U}_{jk}=\sigma^{V}_{jk}=0$ when $j\in\mathcal{C}_{\ell}$ and
$k\in\mathcal{C}_{\ell^{\prime}}$ for some $\ell\neq\ell^{\prime}$. We have

$\sup_{0\leq t\leq C_{0}\sqrt{\log
d}}\left|\frac{\mathbb{P}(||U||_{\infty}>t)}{\mathbb{P}(||V||_{\infty}>t)}-1\right|\leq
O\left(\frac{\Delta_{0}\log d}{\mathfrak{p}}\right),$ (3.2)

for some constant $C_{0}>0$.

When applying the above result to our multiple degree testing problem,
specifically the covariance of maximum test statistics for pairs of non-hub
nodes, $\Delta_{0}$ can be controlled as $k_{\tau}^{2}$ which is in a constant
order. In Theorem 3.2, the quantity $\mathfrak{p}$ represents the number of
connected subgraphs shared by the coviarance matrix networks of $U$ and $V$.
We refer to Theorem B.4 in the appendix for a generalized definition of
$\mathfrak{p}$ to strengthen the results in (3.2). The $\mathfrak{p}$ in the
denominator of the right hand side of Cramér-type comparison bound in (3.2) is
necessary: it is possible that even if $\Delta_{0}$ is small, when
$\mathfrak{p}$ is large, the Camér-type Gaussian comparison bound is not
converging to zero. For example, consider Gaussian vectors with unit variances
$U=(X_{1},X_{2},Z,\ldots,Z)\in\mathbb{R}^{d}$,
$V=(Y_{1},Y_{2},Z,\ldots,Z)\in\mathbb{R}^{d}$, where
$\text{corr}(X_{1},X_{2})=0.9,\text{corr}(Y_{1},Y_{2})=0$ and
$(X_{1},X_{2})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}Z$, $(Y_{1},Y_{2})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}Z$. For this case, the Camér-type Gaussian
comparison bound

$\sup_{0\leq t\leq C_{0}\sqrt{\log
d}}\left|\frac{\mathbb{P}(||U||_{\infty}>t)}{\mathbb{P}(||V||_{\infty}>t)}-1\right|=\sup_{0\leq
t\leq C_{0}\sqrt{\log
d}}\left|\frac{\mathbb{P}(\max\\{|X_{1}|,|X_{2}|,|Z|\\}>t)}{\mathbb{P}(\max\\{|Y_{1}|,|Y_{2}|,|Z|\\}>t)}-1\right|$

is not converging to zero as $d$ goes to infinity even if the corresponding
$\Delta_{0}$ is 1 but $\mathfrak{p}=2$.

Compared with Theorem 3.1, the above theorem provides a sharper comparison
bound for large $\mathfrak{p}$ and small $\Delta_{0}$. The two theorems
together describe a interesting phase transition phenomenon, i.e., the
dependence on $\bm{\Sigma}^{U}-\bm{\Sigma}^{V}$ of the Cramér-type comparison
bound exhibits a difference behavior in the regime of large $\mathfrak{p}$ and
small $\Delta_{0}$ versus the regime of small $\Delta_{\infty}$.

The proof of Theorem 3.2 can be found in Appendix B.2. Our main technical
innovation is to establish a new type of anti-concentration bound for
“derivatives” of Gaussian maxima. Since both the indicator function and
maximum function are discontinuous, we follows the idea of using smoothing
approximation as in the proof of Theorem 3.1, specifically, we bound the
following term

$\mathbb{E}[|\partial_{j}\partial_{k}\varphi(W(s))|\cdot\mathds{1}{(t-\epsilon\leq||W(s)||_{\infty}\leq
t+\epsilon)}],$ (3.3)

where $\varphi$ is the same approximation function of the indicator of
$\ell_{\infty}$ norm with certain smoothing parameter $\beta$. Note that
$\mathbb{E}[\mathds{1}{(t-\epsilon\leq||W(s)||_{\infty}\leq t+\epsilon)}]$ is
the anti-concentration bound for Gaussian maxima (Chernozhukov et al., 2014).
A non-uniform version is also established in Kuchibhotla et al. (2021). (3.3)
can be viewed as the anti-concentration bound on the second order partial
derivatives of the smooth approximation function. When deriving the comparison
bound in terms of $\ell_{0}$ norm difference, we have to deal with such terms
as (3.3) when $\sigma^{U}_{jk}\neq\sigma^{V}_{jk}$. We show (3.3) can be
controlled as

$\mathbb{E}[|\partial_{j}\partial_{k}\varphi(W(s))|\cdot\mathds{1}{(t-\epsilon\leq||W(s)||_{\infty}\leq
t+\epsilon)}]\lesssim\frac{\mathbb{P}\left({||V||_{\infty}>t}\right)(\log
d)^{2}}{\epsilon\beta\mathfrak{p}}.$

The above anti-concentration bound is non-uniform and has only a logarithm
dependence on the dimension $d$. It provides a relatively sharp
characterization when $t$ is large and the Gaussian graphical model is not
highly connected (i.e., the number of connected components $\mathfrak{p}$
being large).

##  4 Discovering hub responses in multitask regression

The theoretical results presented in Section 3 will be the cornerstone for
establishing FDR control of the multiple testing problem described in Section
2. As seen previously, the testing problem (2.1) is set up in a quite general
way: $\bm{\Theta}$ is a weight matrix, and we would like to select rows whose
$\ell_{0}$ norm exceeds some threshold. This section considers the specific
application to multitask/multiple response regression, which turns out to be
less involved. We take advantage of it and demonstrate how to utilize the
probabilistic tools in Section 3. After that, the theoretical results on FDR
control for the Gaussian graphical models are presented and discussed in
Section 5.

In multitask regression problem, multiple response variables are regressed on
a common set of predictors. We can view this example as a bipartite graph
$\mathcal{G}=(\mathcal{V}_{1},\mathcal{V}_{2},\mathcal{E}),|\mathcal{V}_{1}|=d_{1},|\mathcal{V}_{2}|=d_{2}$,
where $\mathcal{V}_{1}$ contains the response variables and $\mathcal{V}_{2}$
represents the common set of predictors. Each entry of the weight matrix
$\bm{\Theta}$ indicates whether a given predictor is non-null or not for a
given response variable. In the case of parametric model,
$\bm{\Theta}\in\mathbb{R}^{d_{1}\times d_{2}}$ corresponds to the parameter
matrix. One might be interested in identifying shared sparsity patterns across
different response variables. It can be solved by selecting a set of
predictors being non-null for all response variables (Obozinski et al., 2006;
Dai and Barber, 2016). This section problem is column-wise in the sense that
we want to select columns of $\bm{\Theta}$, denoted by $\bm{\Theta}_{\cdot
j}$, such that $||\bm{\Theta}_{\cdot j}||_{0}=d_{1}$. It is also interesting
to consider a row-wise selection problem formalized in (2.1). Under the
multitask regression setup, we would like to select response variables with at
least a certain amount of non-null predictors. We will call this type of
response variables hub responses throughout the section. This has practical
applications in real-world problems such as the gene-disease network.

Consider the multitask regression problem with linear models, we have $n$
i.i.d. pairs of the response vector and the predictor vector, denoted by
$(\bm{Y}_{1},\bm{X}_{1}),(\bm{Y}_{2},\bm{X}_{2}),\dots,(\bm{Y}_{n},\bm{X}_{n})$,
where $\bm{Y}_{i}\in\mathbb{R}^{d_{1}},\bm{X}_{i}\in\mathbb{R}^{d_{2}}$
satisfy the following relationship,

$\displaystyle\bm{Y}_{i}=\bm{\Theta}\bm{X}_{i}+\bm{E}_{i},\text{ where
}\bm{E}_{i}\sim\mathcal{N}(0,\mathbf{D}_{d_{1}\times d_{1}})\text{ and
}\bm{X}_{i}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\bm{E}_{i},$ (4.1)

where $\bm{\Theta}\in\mathbb{R}^{d_{1}\times d_{2}}$ is the parameter matrix
and $\mathbf{D}$ is a $d_{1}$ by $d_{1}$ diagonal matrix whose diagonal
elements $\sigma_{j}^{2}$ is the noise variance for response variable
$\bm{Y}^{(j)}$. Let $\bm{X}$ be the design matrix with rows
$\bm{X}_{1}^{\top},\dots,\bm{X}_{n}^{\top}$, shared by different response
variables, and assume the noise variables are independent conditional on the
design matrix $\bm{X}$. Let $s=\max_{j\in[d_{1}]}||\bm{\Theta}_{j}||_{0}$ be
the sparsity level of the parameter matrix $\bm{\Theta}$, we want to select
columns of the parameter matrix which has at least $k_{\tau}$ nonzero entries,
i.e., select nodes with large degree among $[d_{1}]$ in the bipartite graph
$\mathcal{G}=(\mathcal{V}_{1},\mathcal{V}_{2},\mathcal{E})$.

As mentioned in Section 2, some estimator of the parameter matrix is needed to
conduct hypothesis testing. Debiased Lasso is widely used for parameter
estimation and statistical inference in high dimensional linear models
(Javanmard and Montanari, 2014a, b). For each response variable
$\bm{Y}^{(j)},j\in[d_{1}]$, we compute the debiased Lasso estimator, denoted
by $\widetilde{\bm{\Theta}}^{\text{d}}_{j}$ as

$\displaystyle\widetilde{\bm{\Theta}}^{\text{d}}_{j}=\widehat{\bm{\Theta}}_{j}+\frac{1}{n}\,\mathbf{M}\bm{X}^{\top}(\bm{Y}^{(j)}-\bm{X}\widehat{\bm{\Theta}}_{j}),\text{
where
}\widehat{\bm{\Theta}}_{j}=\arg\min_{\beta\in\mathbb{R}^{d_{2}}}\Big{\\{}\frac{1}{2n}\|\bm{Y}^{(j)}-\bm{X}\beta\|^{2}_{2}+\lambda\|\beta\|_{1}\Big{\\}}\,.$
(4.2)

Note the above $\mathbf{M}$ is defined as
$\mathbf{M}=(m_{1},\dots,m_{d_{2}})^{\top}$ where

$\displaystyle m_{i}=$
$\displaystyle\mathop{\mathrm{argmin}}_{m}m^{\top}\widehat{\bm{\Sigma}}m,\quad\text{s.t.
}\|\widehat{\bm{\Sigma}}m-e_{i}\|_{\infty}\leq\mu\,,$ (4.3)

and here $\widehat{\bm{\Sigma}}=(\bm{X}^{\top}\bm{X})/n$.

Then the debiased estimator of the parameter matrix, defined by
$\widetilde{\bm{\Theta}}^{\text{d}}:=(\widetilde{\bm{\Theta}}^{\text{d}}_{1},\cdots,\widetilde{\bm{\Theta}}^{\text{d}}_{d_{1}})^{\top}$,
will be used the input
$\\{\widetilde{\bm{\Theta}}_{e}\\}_{e\in\mathcal{V}_{1}\times\mathcal{V}_{2}}$
of Algorithm 2. In addition, we also need to compute the quantile of the
maximum statistics. There exist many work studying the asymptotic distribution
of the debiased Lasso estimator. Among them, the results in Javanmard and
Montanari (2014a) (when translated into our multitask regression setup) imply,
for each response variable $\bm{Y}^{(j)},j\in[d_{1}]$,

$\sqrt{n}(\widetilde{\bm{\Theta}}^{\text{d}}_{j}-\bm{\Theta}_{j})=Z+\Xi,\quad
Z|\bm{X}\sim\mathcal{N}(0,\sigma_{j}^{2}M\widehat{\bm{\Sigma}}M^{\top}),$
(4.4)

under proper assumptions. Additionally with a natural probabilistic model of
the design matrix, the bias term can be showed to be
$||\Xi||_{\infty}=O(\frac{s\log d_{2}}{\sqrt{n}})$ with high probability. As
discussed in (Javanmard and Montanari, 2014a), the asymptotic normality result
can be used for deriving confidence intervals and statistical hypothesis
tests. As the noise variance $\sigma_{j}$ is unknown, the scaled Lasso is used
for its estimation (Javanmard and Montanari, 2014a; Sun and Zhang, 2012),
given by the following joint optimization problem,

$\\{\widehat{\bm{\Theta}}_{j},\widehat{\sigma}_{j}\\}=\arg\min_{\beta\in\mathbb{R}^{d_{2}},\sigma>0}\Big{\\{}\frac{1}{2\sigma
n}\|\bm{Y}^{(j)}-\bm{X}\beta\|^{2}_{2}+\frac{\sigma}{2}+\lambda\|\beta\|_{1}\Big{\\}}.$
(4.5)

Regarding our testing problem, intuitively we can use the quantile of the
Gaussian maxima of
$\mathcal{N}(0,\widehat{\sigma}_{j}^{2}M\widehat{\bm{\Sigma}}M^{\top})$ to
approximate the quantile of maximum statistic $T_{E}=\underset{(j,k)\in
E}{\max}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{jk}|$ for some given
subset $E$. Specifically, let
$Z_{j}\mid\bm{X},\bm{Y}^{(j)}\sim\mathcal{N}(0,\widehat{\sigma}_{j}^{2}M\widehat{\bm{\Sigma}}M^{\top})$
where $Z_{j}\in\mathbb{R}^{d_{2}}$ and consider the subset
$E\subset\\{j\\}\times\mathcal{V}_{2}$, we approximate the quantile of $T_{E}$
by the following

$T^{\mathcal{N}}_{E}:=\underset{(j,k)\in
E}{\max}\;|Z_{jk}|,\quad\widehat{c}(\alpha,E)=\inf\left\\{t\in\mathbb{R}:\mathbb{P}_{Z}\left(T^{\mathcal{N}}_{E}\leq
t\right)\geq 1-\alpha\right\\}.$ (4.6)

Indeed, under proper scaling conditions, similar results as (2.5) can be
established, i.e., as ${n,d\rightarrow\infty}$,

$\sup_{\alpha\in(0,1)}\left|\mathbb{P}\left(\max_{(j,k)\in
E}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{jk}-{\bm{\Theta}}_{jk}|>\widehat{c}(\alpha,E)\right)-\alpha\right|\rightarrow
0.$ (4.7)

The above result is based on two ingredients: the asymptotic normality result
and the control of the bias term $\Xi$. Below we list the required assumptions
for those two ingredients, i.e., (4.4) and $||\Xi||_{\infty}=O(\frac{s\log
d_{2}}{\sqrt{n}})$.

######  Assumption 4.1 (Debiased Lasso with random designs).

The following assumptions are from the ones of Theorems 7 and 8 in Javanmard
and Montanari (2014a).

  * •

Let
$\bm{\Sigma}=\mathbb{E}\left[{\bm{X}_{1}\bm{X}_{1}^{\top}}\right]\in\mathbb{R}^{d_{2}\times
d_{2}}$ be such that $\sigma_{\min}(\bm{\Sigma})\geq C_{\min}>0$, and
$\sigma_{\max}(\bm{\Sigma})\leq C_{\max}<\infty$, and
$\max_{j\in[d_{2}]}\bm{\Sigma}_{jj}\leq 1$. Assume $\bm{X}\bm{\Sigma}^{-1/2}$
have independent subgaussian rows, with zero mean and subgaussian norm
$\|\bm{\Sigma}^{-1/2}\bm{X}_{i}\|_{\psi_{2}}=\kappa$, for some constant
$\kappa\in(0,\infty)$.

  * •

$\mu=a\sqrt{(\log d_{2})/n}$, and $n\geq\max(\nu_{0}s\log(d_{2}/s),\nu_{1}\log
d_{2})$, $\nu_{1}=\max(1600\kappa^{4},a/4)$, and
$\lambda=\sigma\sqrt{(c^{2}\log d_{2})/n}$.

Remark that there may exist other ways of obtaining a consistent estimator of
$\bm{\Theta}$ and sufficiently accurate quantile estimates under different
assumptions. Since it is not the main focus of this paper, we will not
elaborate on it. As mentioned before, the Kolmogorov type result in (4.7) can
be immediately applied to the global testing problem to guarantee FWER
control. However, it is not sufficient for FDR control of the multiple testing
problem in this paper. And this is when the Cramér-type comparison bound for
Gaussian maxima established in Section 3 play its role. In addition, signal
strength condition is needed. Recall that
$\mathcal{H}_{0}=\\{j\in[d_{1}]:||\bm{\Theta}_{j}||_{0}<k_{\tau}\\}$ with
$d_{0}=|\mathcal{H}_{0}|$, we consider the following rows of $\bm{\Theta}$,

$\mathcal{B}:=\\{j\in\mathcal{H}_{0}^{c}:\forall
k\in\text{supp}(\bm{\Theta}_{j}),|\bm{\Theta}_{jk}|>c\sqrt{{\log
d_{2}}/{n}}\\},$ (4.8)

and define the proportion of such rows as $\rho=|\mathcal{B}|/d_{1}$. In the
context of multitask regression, $\rho$ measures the proportion of hub
response variables whose non-null parameter coefficients all exceed certain
thresholds, thus characterizes the overall signal strength. Below we present
our result on FDP/FDR control under appropriate assumptions.

######  Theorem 4.2 (FDP/FDR control).

Under Assumption 4.1 and the scaling condition $\frac{d_{2}\log
d_{2}+d_{0}}{d_{0}d_{2}\rho}+\frac{s\log^{2}d_{2}}{n^{1/2}}+\frac{\log^{2}d_{2}}{(n\rho)^{1/5}}=o(1)$,
if we implement the StarTrek procedure in Algorithm 2 with $\bm{\Theta}$
estimated by (4.2) and the quantiles approximated by (4.6), as
$(n,d_{1},d_{2})\rightarrow\infty$, we have

$\text{FDP}\leq
q\frac{d_{0}}{d_{1}}+o_{\mathbb{P}}(1)\quad\text{and}\quad\lim_{(n,d_{1},d_{2})\rightarrow\infty}\text{FDR}\leq
q\frac{d_{0}}{d_{1}}.$ (4.9)

The proof of Theorem 4.2 can be found in Appendix A.3. Note that signal
strength conditions which require some entries of parameter matrix
$\bm{\Theta}$ have magnitudes exceeding $c\sqrt{\log d_{2}/n}$ are usually
assumed in existing work studying FDR control problem for high dimensional
models (Liu, 2013; Liu and Shao, 2014; Liu and Luo, 2014; Xia et al., 2015,
2018; Javanmard and Javadi, 2019).

##  5 Discovering hub nodes in Gaussian graphical models

This section focuses on the hub node selection problem on Gaussian graphical
models. Recall in Section 2, we first compute the one-step estimator
$\\{\widehat{\bm{\Theta}}^{\text{d}}_{e}\\}_{e\in\mathcal{V}\times\mathcal{V}}$
in (2.3) then take its standardized version
$\\{\widetilde{\bm{\Theta}}^{\text{d}}_{e}\\}_{e\in\mathcal{V}\times\mathcal{V}}$
as the input of Algorithm 2 i.e.,

$\widehat{\bm{\Theta}}^{\text{d}}_{jk}:=\widehat{\bm{\Theta}}_{jk}-\frac{\widehat{\bm{\Theta}}_{j}^{\top}\left(\widehat{\bm{\Sigma}}\widehat{\bm{\Theta}}_{k}-\mathbf{e}_{k}\right)}{\widehat{\bm{\Theta}}_{j}^{\top}\widehat{\bm{\Sigma}}_{j}},\quad\widetilde{\bm{\Theta}}^{\text{d}}_{jk}:={\widehat{\bm{\Theta}}^{\text{d}}_{jk}}/{\sqrt{\widehat{\bm{\Theta}}^{\text{d}}_{jj}\widehat{\bm{\Theta}}^{\text{d}}_{kk}}}.$
(5.1)

Our StarTrek filter selects nodes with large degrees based on the maximum
statistics $T_{E}=\underset{(j,k)\in
E}{\max}\;\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{jk}|$ over certain
subset $E$. We use the Gaussian multiplier bootstrap (2.4) to approximate the
quantiles, specifically,

$\widehat{c}(\alpha,E)=\inf\left\\{t\in\mathbb{R}:\mathbb{P}_{\xi}\left(T^{\mathcal{B}}_{E}\leq
t\right)\geq 1-\alpha\right\\}.$ (5.2)

Chernozhukov et al. (2013) shows that this quantile approximation is accurate
enough for FWER control in modern high dimensional simultaneous testing
problems. Their results are based on the control of the non-asymptotic bounds
in a Kolmogorov distance sense. Lu et al. (2017) also takes advantage of this
result to test single hypothesis of graph properties or derive confidence
bounds on graph invariants.

However, in order to conduct combinatorial variable selection with FDR control
guarantees, we need more refined studies about the accuracy of the quantile
approximation. This is due to the ratio nature of the definition of FDR, as
explained in Section 2.2. Compared with the results in Chernozhukov et al.
(2013), we provide a Cramér-type control on the approximation errors of the
Gaussian multiplier bootstrap procedure. This is built on the probabilistic
tools in Section 3, in particular, the Cramér-type Gaussian comparison bound
with max norm difference in Theorem 3.1. Due to the dependence structure
behind the hub selection problem in Graphical models, we also have to utilize
Theorem 3.2. In a bit more detail, computing the maximum test statistic for
testing node node actually involves the whole graph, resulting complicated
dependence among the test statistics. The non-differentiability of the maximum
function makes it very difficult to track this dependence. Also note that,
this type of difficulty can not be easily circumvented by alternative methods,
due to the discrete nature of the combinatorial inference problem. However, we
figure out that the Cramér-type Gaussian comparison bound with $\ell_{0}$ norm
difference plays an important role in handling this challenge.

In general, the sparsity/density of the graph is closed related to the
dependence level of multiple testing problem on graphical models. For example,
Liu (2013); Xia et al. (2015, 2018) make certain assumptions on the sparsity
level and control the dependence of test statistics when testing multiple
hypotheses on graphical models/networks. For the hub node selection problem in
this paper, a new quantity is introduced, and we will explain why it is
suitable. Recall that we define the set of non-hub response variables in
Section 4. Similarly, the set of non-hub nodes is denoted by
$\mathcal{H}_{0}=\\{j\in[d]:||\bm{\Theta}_{j}||_{0}<k_{\tau}\\}$ with
$d_{0}=|\mathcal{H}_{0}|$. Now we consider the following set,

$S=\\{(j_{1},j_{2},k_{1},k_{2}):j_{1},j_{2}\in\mathcal{H}_{0},j_{1}\neq
j_{2},k_{1}\neq
k_{2},\bm{\Theta}_{j_{1}j_{2}}=\bm{\Theta}_{j_{1}k_{1}}=\bm{\Theta}_{j_{2}k_{2}}=0,\bm{\Theta}_{j_{1}k_{2}}\neq
0,\bm{\Theta}_{j_{2}k_{1}}\neq 0\\}.$ (5.3)

Figure 1: Left panel: a graphical demonstration of the definition of $S$ via
four examples of a $4$-vertex graph; Right panel: four different graph
patterns with $6$ vertices. Calculating $|S|$ yields $10,15,24,51$ for
(a),(b),(c),(d) respectively.

Remark that in the above definition, $k_{1}$ can be the same as $j_{2}$ and
$k_{2}$ can be the same as $j_{1}$. If there exists a large number of nodes
which are neither connected to $j_{1}$ nor $j_{2}$, we then do not need to
worry much about the dependence between the test statistics for non-hub nodes.
Therefore, $|S|$ actually measures the dependence level via checking how a
pair of non-hub nodes interact through other nodes. Liu (2013); Cai et al.
(2013) also examine the connection structures in the $4$-vertex graph and
control the dependence level by carefully bounding the number of the
$4$-vertex graphs with different numbers of edges.

We provide a graphical demonstration of $S$ and show how $|S|$ looks like in
certain types of graph patterns via some simple examples. Though the
definition of $S$ does not exclude the possibility of
$(j_{1},j_{2},k_{1},k_{2})$ being a graph with $2$ or $3$ vertices, we only
draw $4$-vertex graph in Figure 1 for convenience. In the left panel of Figure
1, we consider four different cases of the $4$-vertex graph. The upper two
belong to the set $S$, while the lower three do not. In the right panel, we
consider four graphs which all have $6$ vertices. They have different graph
patterns. For example, (a) clearly has a hub structure. All of the non-hub
nodes are only connected to the hub node. While in (d), the edges are evenly
distributed and each node are connected to its two nearest neighbours. For
each graph, we count the value of $|S|$ and obtain $10,15,24,51$ respectively,
which show a increasing trend of $|S|$. This sort of matches our intuition
that it is relatively easier to discover hub nodes on graph (a) compared with
graph (d). See more evidence in the empirical results of Section 6.

In addition to $|S|$, we also characterize the dependence level via the
connectivity of the graph, specifically let $p$ be the number of connected
components. And similarly as in Section 4, we define $\rho$ to measure the
signal strength, i.e., $\rho=|\mathcal{B}|/d$, where
$\mathcal{B}:=\\{j\in\mathcal{H}_{0}^{c}:\forall
k\in\text{supp}(\bm{\Theta}_{j}),|\bm{\Theta}_{jk}|>c\sqrt{{\log d}/{n}}\\}$.
In the following, we list our assumptions needed for FDR control.

######  Assumption 5.1.

Suppose that $\bm{\Theta}\in\mathcal{U}(M,s,r_{0})$ and the following
conditions hold:

  1. (i)

Signal strength and scaling condition.

$\displaystyle\frac{\log d}{\rho}\left(\frac{(\log
d)^{19/6}}{n^{1/6}}+\frac{(\log d)^{11/6}}{\rho^{1/3}n^{1/6}}+\frac{{s(\log
d)^{3}}}{{n}^{1/2}}\right)=o(1).$ (5.4)

  2. (ii)

Dependency and connectivity condition.

$\displaystyle\frac{\log d}{\rho d_{0}}+\frac{({\log d})^{2}|S|}{\rho
d_{0}^{2}p}=o(1).$ (5.5)

In the above assumption, (5.4) places conditions on the signal strength and
scaling. The first and the second term come from the Cramér-type large
deviation bounds in the high dimensional CLT setting (Kuchibhotla et al.,
2021) and the Cramér-type Gaussian comparison bound established in Theorem
3.1. And the third term comes from the fact that the relevant test statistics
arise as maxima of approximate averages instead of the exact averages and thus
the approximation error needs to be controlled. See similar discussions about
this in (Chernozhukov et al., 2013). Remark that the signal strength condition
is mild here, due to similar reasons as the discussion in Section 4. Regarding
(5.5), there is a trade-off between the dependence level and connectivity
level of the topological structure. $|S|/d_{0}^{2}$ characterizes how the test
statistics of non-hub nodes are correlated to each other in average. $p$ by
definition describes the level of connectivity. Due to the condition (5.5),
larger signal strength generally makes the hub selection problem easier. And
when $|S|/d_{0}^{2}$ is small, the graph is allowed to be more connected. When
there exist more sub-graphs, we allow higher correlations between the non-hub
nodes. Note that the cardinality of $S$ is directly related to the $\ell_{0}$
norm covariance matrix difference term $\Delta_{0}$, and arises from the
application of Theorem 3.2. In the following, we present our core theoretical
result on FDP/FDR control for hub selection using the StarTrek filter on
Gaussian graphical models.

######  Theorem 5.2 (FDP/FDR control).

Under Assumption 5.1, the StarTrek procedure in Algorithm 2 with (5.1) as
input and the quantiles approximated by (5.2) satisfies: as
$(n,d)\rightarrow\infty$,

$\text{FDP}\leq
q\frac{d_{0}}{d}+o_{\mathbb{P}}(1)\quad\text{and}\quad\lim_{(n,d)\rightarrow\infty}\text{FDR}\leq
q\frac{d_{0}}{d}.$ (5.6)

The proof can be found in Appendix A.1. Remark that control of the FDR does
not prohibit the FDP from varying. Therefore our result on FDP provides a
stronger guarantee on controlling the false discoveries. See clear empirical
evidence in Section 6.1. To the best of our knowledge, the proposed StarTrek
filter in Section 2 and the above FDP/FDR control result are the first
Algorithm and theoretical guarantee for the problem of simultaneously
selecting hub nodes. Existing work like Liu (2013); Liu and Luo (2014); Xia et
al. (2015, 2018); Javanmard and Javadi (2019) focus on the discovery of
continuous signals and their tools are not applicable to the problem here.

##  6 Numerical results

###  6.1 Synthetic data

In this section, we apply the StarTrek filter to synthetic data and
demonstrate the performance of our method. The synthetic datasets are
generated from Gaussian graphical models. The corresponding precision matrices
are specified based on four different types of graphs. Given the number of
nodes $d$ and the number of connected components $p$, we will randomly assign
those nodes into $p$ groups. Within each group (sub-graph), the way of
assigning edges for different graph types will be explained below in detail.
After determinning the adjacency matrix of the graph, we follow Zhao et al.
(2012) to construct the precision matrix, more specifically, we set the off-
diagonal elements to be of value $v$ which control the magnitude of partial
correlations and is closely related to the signal strength. In order to ensure
positive-definiteness, we add some value $v$ together with the absolute value
of the minimal eigenvalues to the diagonal terms. In the following
simulations, $v$ and $u$ are set to be $0.4$ and $0.1$ respectively. Now we
explain how to determine the edges within each group (sub-graph) for four
different graph patterns.

  * •

Hub graph. We randomly pick one node as the hub node of the sub-graph, then
the rest of the nodes are made to connect with this hub node. There is no edge
between the non-hub nodes.

  * •

Random graph. This is the Erdös-Rényi random graph. There is an edge between
each pair of nodes with certain probability independently. In the following
simulations, we will set this probability to be $0.15$ unless stated
otherwise.

  * •

Scale-free graph. In this type of graphs, the degree distribution follows a
power law. We construct it by the Barabási-Albert algorithm: starting with two
connected nodes, then adding each new node to be connected with only one node
in the existing graph; and the probability is proportional to the degree of
the each node in the existing graph. The number of the edges will be the same
as the number of nodes.

  * •

K-nearest-neighbor (knn) graph. For a given number of $k$, we add edges such
that each node is connected to another $k$ nodes. In our simulations, $k$ is
sampled from $\\{1,2,3,4\\}$ with probability mass $\\{0.4,0.3,0.2,0.1\\}$.

See a visual demonstration of the above four different graph patterns in
Appendix E.1. Throughout the simulated examples, we fix the number of nodes
$d$ to be $300$ and vary other quantities such as sample size $n$ or the
number of connected components $p$. To estimate the precision matrix, we run
the graphical Lasso algorithm with 5-fold cross-validation. Then we obtain the
standardized debiased estimator as described in (2.3). To obtain the quantile
estimates, we use the Gaussian multiplier bootstrap with 4000 bootstrap
samples. The threshold $k_{\tau}$ for determining hub nodes is set to be $3$.
And all results (of FDR and power) are averaged over 64 independent
replicates.

As we can see from Table 1, the FDRs of StarTrek filter for different types of
graph are well controlled below the nominal levels. In hub graph, the FDRs are
relatively small but the power is still pretty good. Similar phenomenon for
multiple edge testing problem is observed (Liu, 2013). In the context of node
testing, it is also unsurprising. These empirical results actually match our
demonstration about $|S|$ in Figure 1: hub graphs have a relatively weaker
dependence structure (smaller $S$ values) and make it is easier to discover
true hub nodes without making many errors.

Table 1:  Empirical FDR $d=300~{}$ | $q=0.1$ | $q=0.2$
---|---|---
$n$ | 200 | 300 | 400 | 200 | 300 | 400
| $\quad\quad p=20$ |
hub | 0.0000 | 0.0000 | 0.0007 | 0.0000 | 0.0000 | 0.0029
random | 0.0255 | 0.0383 | 0.0467 | 0.0521 | 0.0770 | 0.0833
scale-free | 0.0093 | 0.0211 | 0.0282 | 0.0352 | 0.0486 | 0.0581
knn | 0.0101 | 0.0296 | 0.0370 | 0.0228 | 0.0620 | 0.0769
| $\quad\quad p=30$ |
hub | 0.0013 | 0.0000 | 0.0016 | 0.0027 | 0.0054 | 0.0036
random | 0.0347 | 0.0359 | 0.0568 | 0.0725 | 0.0753 | 0.0963
scale-free | 0.0215 | 0.0335 | 0.0317 | 0.0521 | 0.0624 | 0.0584
knn | 0.0297 | 0.0420 | 0.0563 | 0.0504 | 0.0857 | 0.1030

The power performance of the StarTrek filter is showed in Table 2. As the
sample size grows, we see the power is increasing for all four different types
of graphs. When $p$ is larger, there are more hub nodes in general due to the
way of constructing the graphs, and we find the power is higher. Among
different types of graphs, the power in hub graph and scale-free graph is
higher than that in random and knn graph since the latter two are relatively
denser and have more complicated topological structures.

Table 2: Power $d=300~{}$ | $q=0.1$ | $q=0.2$
---|---|---
$n$ | 200 | 300 | 400 | 200 | 300 | 400
| $\quad\quad p=20$ |
hub | 0.7109 | 0.9453 | 0.9898 | 0.7805 | 0.9648 | 0.9938
random | 0.3343 | 0.7815 | 0.9408 | 0.4520 | 0.8514 | 0.9604
scale-free | 0.4524 | 0.8145 | 0.9363 | 0.5281 | 0.8614 | 0.9568
knn | 0.0905 | 0.5306 | 0.8067 | 0.1634 | 0.6511 | 0.8630
| $\quad\quad p=30$ |
hub | 0.6848 | 0.9244 | 0.9706 | 0.7588 | 0.9459 | 0.9784
random | 0.4882 | 0.8863 | 0.9790 | 0.5770 | 0.9225 | 0.9870
scale-free | 0.6472 | 0.9047 | 0.9810 | 0.7197 | 0.9331 | 0.9870
knn | 0.2409 | 0.6841 | 0.8922 | 0.3298 | 0.7706 | 0.9241

In Figure 2 and 3, we demonstrate the performance of our method in the random
graph with different parameters. Specifically, we vary the connecting
probability changing from $0.1$ to $0.3$ in the x-axis. In those plots, we see
the FDRs are all well controlled below the nominal level $q=0.1$. As the
connecting probability of the random graph grows, the graph gets denser,
resulting more hub nodes. Thus we can see the height of the short blue solids
lines (representing $qd_{0}/d$) is decreasing. Based on our results in Theorem
5.2, the target level of FDP/FDR control is $qd_{0}/d$. This is why we find
the mean and median of each box-plot is getting smaller as the connecting
probability increases (hence $d_{0}$ decreases).

The box-plots and the jittering points show that our StarTrek procedure not
only controls the FDR but also prohibit it from varying too much, as implied
by the theoretical results on FDP control in Section 5. Regarding the power
plots, we see that the power is smaller when the graph is denser since the hub
selection problem becomes more difficult with more disturbing factors. Plots
with nominal FDR level $q=0.2$ are deferred to Appendix E.3.

###  6.2 Application to gene expression data

We also apply our method to the Genotype-Tissue Expression (GTEx) data studied
in Lonsdale et al. (2013). Beginning with a 2.5-year pilot phase, the GTEx
project establishes a great database and associated tissue bank for studying
the relationship between certain genetic variations and gene expressions in
human tissues. The original dataset involves 54 non-diseased tissue sites
across 549 research subjects. Here we only focus on analyzing the breast
mammary tissues. It is of great interest to identify hub genes over the gene
expression network.

Figure 2: FDP and power plots for the StarTrek filter in the random graph. The connecting probability is varied on the x-axis. The number of samples $n$ is chosen to be $300$ and the number of connected components $p$ equals 20. The nominal FDR level is set to be $q=0.1$; the short blue solid lines correspond to $q{d_{0}}/{d}$, calculated by averaging over the $64$ replicates. For both panels, the box plots are plotted with the black points representing the outliers. Colored points are jittered around, demonstrating how the FDP and power distribute. Figure 3: FDP and power plots for the StarTrek filter in the random graph. The other setups are the same as Figure 3 except for $p=30$. Male |  |  |  |
---|---|---|---|---
Female |  |  |  |
Figure 4: The above graphs are based the estimated precision matrices (the
left two plots). The adjacency matrices of the other six plots are based on
the standardized estimated precision matrices but thresholded at
$0.025,0.05,0.075$ respectively. Blue vertices represent the selected hub
genes.

First we calculate the variances of the gene expression data and focus on the
top $100$ genes in the following analysis. The data involves $n=291$ samples
for male individuals and $n=168$ samples for female individuals. The original
count data is log-transformed and scaled. We then obtain the estimator of the
precision matrix by the Graphical Lasso with 2-fold cross-validation. As for
the hub node criterion, we set $k_{\tau}$ as the 50% quantile of the node
degrees in the estimated precision matrix. We run StarTrek filter with $2000$
bootstrap samples and nominal FDR level $q=0.1$ to select hub genes for both
the male and female datasets.

Figure 5: Plots of the sorted p-values ($\alpha_{j},j\in[d]$) in Algorithm 2.
Those blue points correspond to selected hub genes. The blue line is the
rejection line of the BHq procedure. The coordinates of the plots are flipped.
We abbreviate the names of the 100 genes and only show selected ones with blue
colored text. The upper panel and the lower panel are based on male and female
data respectively.

Figure 4 shows that the selected hub genes by the StarTrek filter also have
large degrees on the estimated gene networks (based on the estimated precision
matrices). In Figure 5, the results for male and female dataset agree with
each other except that the number of selected hub genes using female dataset
is smaller due to a much smaller sample size. The selected hub genes are found
to play an important role in breast-related molecular processes, either as
central regulators or their abnormal expressions are considered as the causes
of breast cancer initiation and progression, see relevant literature in
genetic research such as Hellwig et al. (2016); Blein et al. (2015); Chen et
al. (2016); Li et al. (2019); Lou et al. (2020); Mohamed et al. (2014); Bai et
al. (2019); Sirois et al. (2019); Marino et al. (2020); Malvia et al. (2019).
Therefore, our proposed method for selecting hub nodes can be applied to the
hub gene identification problem. It may improve our understanding of the
mechanisms of breast cancer and provide valuable prognosis and treatment
signature.

## References

  * Bai et al. (2019) Bai, J., Zhang, X., Kang, X., Jin, L., Wang, P. and Wang, Z. (2019).  Screening of core genes and pathways in breast cancer development via comprehensive analysis of multi gene expression datasets.  Oncology Letters 18 5821–5830.
  * Barber and Candès (2015) Barber, R. F. and Candès, E. J. (2015).  Controlling the false discovery rate via knockoffs.  The Annals of Statistics 43 2055 – 2085.
  * Barber and Candès (2019) Barber, R. F. and Candès, E. J. (2019).  A knockoff filter for high-dimensional selective inference.  The Annals of Statistics 47 2504 – 2537.
  * Belloni et al. (2014) Belloni, A., Chernozhukov, V. and Hansen, C. (2014).  Inference on treatment effects after selection among high-dimensional controls.  The Review of Economic Studies 81 608–650.
  * Benjamini (2010) Benjamini, Y. (2010).  Discovering the false discovery rate.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72 405–416.
  * Benjamini and Hochberg (1995) Benjamini, Y. and Hochberg, Y. (1995).  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  Journal of the Royal statistical society: series B (Methodological) 57 289–300.
  * Benjamini and Yekutieli (2001) Benjamini, Y. and Yekutieli, D. (2001).  The control of the false discovery rate in multiple testing under dependency.  The Annals of Statistics 1165–1188.
  * Bentkus (1990) Bentkus, V. (1990).  Smooth approximations of the norm and differentiable functions with bounded support in banach space $l_{\infty}^{k}$.  Lithuanian Mathematical Journal 30 223–230.
  * Blein et al. (2015) Blein, S., Barjhoux, L., investigators, G., Damiola, F., Dondon, M.-G., Eon-Marchais, S., Marcou, M., Caron, O., Lortholary, A., Buecher, B., Vennin, P., Berthet, P., Noguès, C., Lasset, C., Gauthier-Villars, M., Mazoyer, S., Stoppa-Lyonnet, D., Andrieu, N., Thomas, G., Sinilnikova, O. M. and Cox, D. G. (2015).  Targeted sequencing of the mitochondrial genome of women at high risk of breast cancer without detectable mutations in brca1/2.  PLoS One 10 e0136192.
  * Cai et al. (2011) Cai, T., Liu, W. and Luo, X. (2011).  A constrained $\ell_{1}$ minimization approach to sparse precision matrix estimation.  Journal of the American Statistical Association 106 594–607.
  * Cai et al. (2013) Cai, T., Liu, W. and Xia, Y. (2013).  Two-sample covariance matrix testing and support recovery in high-dimensional and sparse settings.  Journal of the American Statistical Association 108 265–277.
  * Cai and Ma (2013) Cai, T. T. and Ma, Z. (2013).  Optimal hypothesis testing for high dimensional covariance matrices.  Bernoulli 19 2359 – 2388.
  * Cai and Zhang (2016) Cai, T. T. and Zhang, A. (2016).  Inference for high-dimensional differential correlation matrices.  Journal of Multivariate Analysis 143 107–126.
  * Candès et al. (2018) Candès, E. J., Fan, Y., Janson, L. and Lv, J. (2018).  Panning for gold: Model-X knockoffs for high-dimensional controlled variable selection.  Journal of the Royal Statistical Society: Series B 80 551–577.
  * Chang et al. (2016) Chang, J., Shao, Q.-M. and Zhou, W.-X. (2016).  Cramér-type moderate deviations for Studentized two-sample $U$-statistics with applications.  The Annals of Statistics 44 1931 – 1956.
  * Chen et al. (2016) Chen, W.-C., Wang, C.-Y., Hung, Y.-H., Weng, T.-Y., Yen, M.-C. and Lai, M.-D. (2016).  Systematic analysis of gene expression alterations and clinical outcomes for long-chain acyl-coenzyme a synthetase family in cancer.  PLoS One 11 e0155660.
  * Chernozhukov et al. (2013) Chernozhukov, V., Chetverikov, D. and Kato, K. (2013).  Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors.  The Annals of Statistics 41 2786–2819.
  * Chernozhukov et al. (2014) Chernozhukov, V., Chetverikov, D. and Kato, K. (2014).  Anti-concentration and honest, adaptive confidence bands.  The Annals of Statistics 42 1787–1818.
  * Chernozhukov et al. (2015) Chernozhukov, V., Chetverikov, D. and Kato, K. (2015).  Comparison and anti-concentration bounds for maxima of gaussian random vectors.  Probability Theory and Related Fields 162 47–70.
  * Dai et al. (2020a) Dai, C., Lin, B., Xing, X. and Liu, J. S. (2020a).  False discovery rate control via data splitting.  arXiv preprint arXiv:2002.08542 .
  * Dai et al. (2020b) Dai, C., Lin, B., Xing, X. and Liu, J. S. (2020b).  A scale-free approach for false discovery rate control in generalized linear models.  arXiv preprint arXiv:2007.01237 .
  * Dai and Barber (2016) Dai, R. and Barber, R. (2016).  The knockoff filter for FDR control in group-sparse and multitask regression.  In International Conference on Machine Learning. PMLR.
  * Ding and Zhou (2020) Ding, X. and Zhou, Z. (2020).  Estimation and inference for precision matrices of nonstationary time series.  The Annals of Statistics 48 2455 – 2477.
  * Eisenach et al. (2020) Eisenach, C., Bunea, F., Ning, Y. and Dinicu, C. (2020).  High-dimensional inference for cluster-based graphical models.  Journal of Machine Learning Research 21.
  * Feng and Ning (2019) Feng, H. and Ning, Y. (2019).  High-dimensional mixed graphical model with ordinal data: Parameter estimation and statistical inference.  In The 22nd International Conference on Artificial Intelligence and Statistics. PMLR.
  * Friedman et al. (2008) Friedman, J., Hastie, T. and Tibshirani, R. (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics 9 432–441.
  * Grönwall (1919) Grönwall, T. H. (1919).  Note on the derivatives with respect to a parameter of the solutions of a system of differential equations.  Annals of Mathematics 20 292–296.
  * Gu et al. (2015) Gu, Q., Cao, Y., Ning, Y. and Liu, H. (2015).  Local and global inference for high dimensional nonparanormal graphical models.  arXiv preprint arXiv:1502.02347 .
  * Hellwig et al. (2016) Hellwig, B., Madjar, K., Edlund, K., Marchan, R., Cadenas, C., Heimes, A.-S., Almstedt, K., Lebrecht, A., Sicking, I., Battista, M. J., Micke, P., Schmidt, M., Hengstler, J. G. and Rahnenführer, J. (2016).  Epsin family member 3 and ribosome-related genes are associated with late metastasis in estrogen receptor-positive breast cancer and long-term survival in non-small cell lung cancer using a genome-wide identification and validation strategy.  PLoS One 11 1–18.
  * Ilyas et al. (2011) Ilyas, M. U., Shafiq, M. Z., Liu, A. X. and Radha, H. (2011).  A distributed and privacy preserving algorithm for identifying information hubs in social networks.  In 2011 Proceedings IEEE INFOCOM. IEEE.
  * Isserlis (1918) Isserlis, L. (1918).  On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables.  Biometrika 12 134–139.
  * Janková and van de Geer (2017) Janková, J. and van de Geer, S. (2017).  Honest confidence regions and optimality in high-dimensional precision matrix estimation.  Test 26 143–162.
  * Janková and van de Geer (2018) Janková, J. and van de Geer, S. (2018).  Inference in high-dimensional graphical models.  arXiv preprint arXiv:1801.08512 .
  * Javanmard and Javadi (2019) Javanmard, A. and Javadi, H. (2019).  False discovery rate control via debiased lasso.  Electronic Journal of Statistics 13 1212 – 1253.
  * Javanmard and Montanari (2013) Javanmard, A. and Montanari, A. (2013).  Nearly optimal sample size in hypothesis testing for high-dimensional regression.  In 2013 51st Annual Allerton Conference on Communication, Control, and Computing (Allerton). IEEE.
  * Javanmard and Montanari (2014a) Javanmard, A. and Montanari, A. (2014a).  Confidence intervals and hypothesis testing for high-dimensional regression.  Journal of Machine Learning Research 15 2869–2909.
  * Javanmard and Montanari (2014b) Javanmard, A. and Montanari, A. (2014b).  Hypothesis testing in high-dimensional regression under the gaussian random design model: Asymptotic theory.  IEEE Transactions on Information Theory 60 6522–6554.
  * Jin et al. (2020) Jin, J., Ke, Z. T., Luo, S. and Wang, M. (2020).  Estimating the number of communities by stepwise goodness-of-fit.  arXiv preprint arXiv:2009.09177 .
  * Ke et al. (2020) Ke, Z. T., Ma, Y. and Lin, X. (2020).  Estimation of the number of spiked eigenvalues in a covariance matrix by bulk eigenvalue matching analysis.  arXiv preprint arXiv:2006.00436 .
  * Kuchibhotla et al. (2021) Kuchibhotla, A. K., Mukherjee, S. and Banerjee, D. (2021).  High-dimensional CLT: Improvements, non-uniform extensions and large deviations.  Bernoulli 27 192 – 217.
  * Lam and Fan (2009) Lam, C. and Fan, J. (2009).  Sparsistency and rates of convergence in large covariance matrix estimation.  The Annals of Statistics 37 4254–4278.
  * Lee et al. (2019) Lee, R. K.-W., Hoang, T.-A. and Lim, E.-P. (2019).  Discovering hidden topical hubs and authorities across multiple online social networks.  IEEE Transactions on Knowledge and Data Engineering 33 70–84.
  * Li and Maathuis (2019) Li, J. and Maathuis, M. H. (2019).  GGM knockoff filter: False discovery rate control for gaussian graphical models.  arXiv preprint arXiv:1908.11611 .
  * Li et al. (2019) Li, Y., Giorgi, E. E., Beckman, K. B., Caberto, C., Kazma, R., Lum-Jones, A., Haiman, C. A., Marchand, L. L., Stram, D. O., Saxena, R. and Cheng, I. (2019).  Association between mitochondrial genetic variation and breast cancer risk: The multiethnic cohort.  PLoS One 14 1–14.
  * Liu et al. (2020) Liu, M., Xia, Y., Cai, T. and Cho, K. (2020).  Integrative high dimensional multiple testing with heterogeneity under data sharing constraints.  arXiv preprint arXiv:2004.00816 .
  * Liu (2013) Liu, W. (2013).  Gaussian graphical model estimation with false discovery rate control.  The Annals of Statistics 41 2948–2978.
  * Liu and Luo (2014) Liu, W. and Luo, S. (2014).  Hypothesis testing for high-dimensional regression models.
  * Liu and Shao (2010) Liu, W. and Shao, Q.-M. (2010).  Cramér-type moderate deviation for the maximum of the periodogram with application to simultaneous tests in gene expression time series.  The Annals of Statistics 38 1913 – 1935.
  * Liu and Shao (2014) Liu, W. and Shao, Q.-M. (2014).  Phase transition and regularized bootstrap in large-scale $t$-tests with false discovery rate control.  The Annals of Statistics 42 2003 – 2025.
  * Liu et al. (2019a) Liu, Y., Gu, H.-Y., Zhu, J., Niu, Y.-M., Zhang, C. and Guo, G.-L. (2019a).  Identification of hub genes and key pathways associated with bipolar disorder based on weighted gene co-expression network analysis.  Frontiers in Physiology 10 1081\.
  * Liu et al. (2019b) Liu, Y., Yi, Y., Wu, W., Wu, K. and Zhang, W. (2019b).  Bioinformatics prediction and analysis of hub genes and pathways of three types of gynecological cancer.  Oncology Letters 18 617–628.
  * Lonsdale et al. (2013) Lonsdale, J., Thomas, J., Salvatore, M., Phillips, R., Lo, E., Shad, S., Hasz, R., Walters, G., Garcia, F., Young, N. et al. (2013).  The genotype-tissue expression (GTEx) project.  Nature Genetics 45 580–585.
  * Lou et al. (2020) Lou, W., Ding, B., Wang, S. and Fu, P. (2020).  Overexpression of gpx3, a potential biomarker for diagnosis and prognosis of breast cancer, inhibits progression of breast cancer cells in vitro.  Cancer Cell International 20 1–15.
  * Lu et al. (2017) Lu, J., Neykov, M. and Liu, H. (2017).  Adaptive inferential method for monotone graph invariants.  arXiv preprint arXiv:1707.09114 .
  * Luscombe et al. (2004) Luscombe, N. M., Babu, M. M., Yu, H., Snyder, M., Teichmann, S. A. and Gerstein, M. (2004).  Genomic analysis of regulatory network dynamics reveals large topological changes.  Nature 431 308–312.
  * Malvia et al. (2019) Malvia, S., Bagadi, S. A. R., Pradhan, D., Chintamani, C., Bhatnagar, A., Arora, D., Sarin, R. and Saxena, S. (2019).  Study of gene expression profiles of breast cancers in indian women.  Scientific Reports 9 1–15.
  * Marino et al. (2020) Marino, N., German, R., Rao, X., Simpson, E., Liu, S., Wan, J., Liu, Y., Sandusky, G., Jacobsen, M., Stoval, M., Cao, S. and Storniolo, A. M. V. (2020).  Upregulation of lipid metabolism genes in the breast prior to cancer diagnosis.  NPJ Breast Cancer 6 1–13.
  * Meinshausen and Bühlmann (2006) Meinshausen, N. and Bühlmann, P. (2006).  High-dimensional graphs and variable selection with the lasso.  The Annals of Statistics 34 1436–1462.
  * Mohamed et al. (2014) Mohamed, M. M., Sabet, S., Peng, D.-F., Nouh, M. A., El-Shinawi, M. and El-Rifai, W. (2014).  Promoter hypermethylation and suppression of glutathione peroxidase 3 are associated with inflammatory breast carcinogenesis.  Oxidative Medicine and Cellular Longevity 2014.
  * Newman et al. (2002) Newman, M. E., Watts, D. J. and Strogatz, S. H. (2002).  Random graph models of social networks.  Proceedings of the National Academy of Sciences 99 2566–2572.
  * Neykov et al. (2019) Neykov, M., Lu, J. and Liu, H. (2019).  Combinatorial inference for graphical models.  The Annals of Statistics 47 795–827.
  * Obozinski et al. (2006) Obozinski, G., Taskar, B. and Jordan, M. (2006).  Multi-task feature selection.  Statistics Department, UC Berkeley, Tech. Rep 2 2\.
  * Peng et al. (2009) Peng, J., Wang, P., Zhou, N. and Zhu, J. (2009).  Partial correlation estimation by joint sparse regression models.  Journal of the American Statistical Association 104 735–746.
  * Power et al. (2013) Power, J. D., Schlaggar, B. L., Lessov-Schlaggar, C. N. and Petersen, S. E. (2013).  Evidence for hubs in human functional brain networks.  Neuron 79 798–813.
  * Ravikumar et al. (2011) Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B. (2011).  High-dimensional covariance estimation by minimizing $\ell_{1}$-penalized log-determinant divergence.  Electronic Journal of Statistics 5 935–980.
  * Ren et al. (2015) Ren, Z., Sun, T., Zhang, C.-H. and Zhou, H. H. (2015).  Asymptotic normality and optimalities in estimation of large Gaussian graphical models.  The Annals of Statistics 43 991 – 1026.
  * Rothman et al. (2008) Rothman, A. J., Bickel, P. J., Levina, E. and Zhu, J. (2008).  Sparse permutation invariant covariance estimation.  Electronic Journal of Statistics 2 494–515.
  * Rubinov and Sporns (2010) Rubinov, M. and Sporns, O. (2010).  Complex network measures of brain connectivity: uses and interpretations.  Neuroimage 52 1059–1069.
  * Shaw et al. (2008) Shaw, P., Kabani, N. J., Lerch, J. P., Eckstrand, K., Lenroot, R., Gogtay, N., Greenstein, D., Clasen, L., Evans, A., Rapoport, J. L., Giedd, J. N. and Wise, S. P. (2008).  Neurodevelopmental trajectories of the human cerebral cortex.  Journal of Neuroscience 28 3586–3594.
  * Shen and Lu (2020) Shen, S. and Lu, J. (2020).  Combinatorial-probabilistic trade-off: Community properties test in the stochastic block models.  arXiv preprint arXiv:2010.15063 .
  * Shen et al. (2012) Shen, X., Pan, W. and Zhu, Y. (2012).  Likelihood-based selection and sharp parameter estimation.  Journal of the American Statistical Association 107 223–232.
  * Sirois et al. (2019) Sirois, I., Aguilar-Mahecha, A., Lafleur, J., Fowler, E., Vu, V., Scriver, M., Buchanan, M., Chabot, C., Ramanathan, A., Balachandran, B., Légaré, S., Przybytkowski, E., Lan, C., Krzemien, U., Cavallone, L., Aleynikova, O., Ferrario, C., Guilbert, M.-C., Benlimame, N., Saad, A., Alaoui-Jamali, M., Saragovi, H. U., Josephy, S., O’Flanagan, C., Hursting, S. D., Richard, V. R., Zahedi, R. P., Borchers, C. H., Bareke, E., Nabavi, S., Tonellato, P., Roy, J.-A., Robidoux, A., Marcus, E. A., Mihalcioiu, C., Majewski, J. and Basik, M. (2019).  A unique morphological phenotype in chemoresistant triple-negative breast cancer reveals metabolic reprogramming and PLIN4 expression as a molecular vulnerability.  Molecular Cancer Research 17 2492–2507.
  * Slepian (1962) Slepian, D. (1962).  The one-sided barrier problem for gaussian noise.  Bell System Technical Journal 41 463–501.
  * Sun and Zhang (2012) Sun, T. and Zhang, C.-H. (2012).  Scaled sparse linear regression.  Biometrika 99 879–898.
  * Sur and Candès (2019) Sur, P. and Candès, E. J. (2019).  A modern maximum-likelihood theory for high-dimensional logistic regression.  Proceedings of the National Academy of Sciences 116 14516–14525.
  * Van de Geer et al. (2014) Van de Geer, S., Bühlmann, P., Ritov, Y. and Dezeure, R. (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  The Annals of Statistics 42 1166–1202.
  * van den Heuvel and Sporns (2013) van den Heuvel, M. P. and Sporns, O. (2013).  Network hubs in the human brain.  Trends in cognitive sciences 17 683–696.
  * Van Der Vaart and Wellner (1996) Van Der Vaart, A. W. and Wellner, J. A. (1996).  Weak convergence.  In Weak convergence and empirical processes. Springer, 16–28.
  * van Handel (2014) van Handel, R. (2014).  Probability in high dimension.  Tech. rep., Princeton University.
  * Xia et al. (2015) Xia, Y., Cai, T. and Cai, T. T. (2015).  Testing differential networks with applications to the detection of gene-gene interactions.  Biometrika 102 247–266.
  * Xia et al. (2018) Xia, Y., Cai, T. and Cai, T. T. (2018).  Multiple testing of submatrices of a precision matrix with applications to identification of between pathway interactions.  Journal of the American Statistical Association 113 328–339.
  * Xing et al. (2019) Xing, X., Zhao, Z. and Liu, J. S. (2019).  Controlling false discovery rate using gaussian mirrors.  arXiv preprint arXiv:1911.09761 .
  * Yang et al. (2018) Yang, Z., Ning, Y. and Liu, H. (2018).  On semiparametric exponential family graphical models.  Journal of Machine Learning Research 19 2314–2372.
  * Yuan et al. (2017) Yuan, L., Chen, L., Qian, K., Qian, G., Wu, C.-L., Wang, X. and Xiao, Y. (2017).  Co-expression network analysis identified six hub genes in association with progression and prognosis in human clear cell renal cell carcinoma (ccrcc).  Genomics Data 14 132–140.
  * Yuan and Lin (2007) Yuan, M. and Lin, Y. (2007).  Model selection and estimation in the Gaussian graphical model.  Biometrika 94 19–35.
  * Zhao et al. (2020) Zhao, Q., Sur, P. and Candes, E. J. (2020).  The asymptotic distribution of the mle in high-dimensional logistic models: Arbitrary covariance.  arXiv preprint arXiv:2001.09351 .
  * Zhao et al. (2012) Zhao, T., Liu, H., Roeder, K., Lafferty, J. and Wasserman, L. (2012).  The huge package for high-dimensional undirected graph estimation in R.  Journal of Machine Learning Research 13 1059–1062.

Supplementary material to

StarTrek: Combinatorial Variable Selection with False Discovery Rate Control

This document contains the supplementary material to the paper “StarTrek:
Combinatorial Variable Selection with False Discovery Rate Control”. Appendix
A presents the proofs of the FDR control results. In Appendix B, we provide
the proofs of two types of Cramér-type comparison bounds for Gaussian maxima.
Appendix C proves the Cramér-type deviation bounds for the Gaussian multiplier
bootstrap. In Appendix D, we establish the validity and a power result of our
test on the degree of a single node. Appendix E contains some plots and tables
deferred from the main paper.

##  Appendix A Proofs for FDR control

In this section, we aim to prove Theorem 5.2. In order to prove the theorem,
we need Lemma A.1 which is about the test of single node degree. Remark that
this lemma proves the asymptotic validity of the test in Algorithm 1 and
provides a power analysis. The signal strength condition is only required for
the power analysis part. To see why Lemma A.1 is useful for establishing FDR
control for our StarTrek procedure in Algorithm 2, we notice the following
equivalence:

$\\{\psi_{j,\alpha}=1\\}=\\{\alpha_{j}\leq\alpha\\},$ (A.1)

where $\alpha$ is a given type-I error level, $\psi_{j,\alpha}$ is the test
described in Algorithm 1, and $\alpha_{j}$ is defined in Algorithm 2. First,
we show $\\{\alpha_{j}\leq\alpha\\}\subset\\{\psi_{j,\alpha}=1\\}$. Note

$\displaystyle\\{\alpha_{j}\leq\alpha\\}$ $\displaystyle=$
$\displaystyle\bigcap_{1\leq s\leq
k_{\tau}}\\{\widehat{c}^{-1}(\sqrt{n}|\widetilde{\bm{\Theta}}_{j,(s)}|,E^{(s)}_{j})\leq\alpha\\}$
(A.2) $\displaystyle=$ $\displaystyle\bigcap_{1\leq s\leq
k_{\tau}}\\{\sqrt{n}|\widetilde{\bm{\Theta}}_{j,(s)}|\geq\widehat{c}(\alpha,E^{(s)}_{j})\\},$

where $E^{(s)}_{j}:=\\{(j,\ell):\ell\neq
j,|\widetilde{\bm{\Theta}}_{j\ell}|\leq|\widetilde{\bm{\Theta}}_{j,(s)}|\\}$.
The first equality is due to the definition of $\alpha_{j}$ and the second
equality holds by the definition of $\widehat{c}^{-1}$. Examining (A.2), we
immediately know
$\sqrt{n}|\widetilde{\bm{\Theta}}_{j,(1)}|\geq\widehat{c}(\alpha,E^{(1)}_{j})$
(here $E^{(1)}_{j}=E_{0}=\\{(k,j):k\in[d],k\neq j\\}$), thus the edge
corresponding to $\widetilde{\bm{\Theta}}_{j,(1)}$ will be rejected in the
first iteration of Algorithm 1. Regarding the edge corresponding to
$\widetilde{\bm{\Theta}}_{j,(2)}$, if
$\sqrt{n}|\widetilde{\bm{\Theta}}_{j,(2)}|\geq\widehat{c}(\alpha,E^{(1)}_{j})$,
then it will be rejected in the first iteration, too. Otherwise, Algorithm 1
enters the second iteration. Since (A.2) implies
$\sqrt{n}|\widetilde{\bm{\Theta}}_{j,(2)}|\geq\widehat{c}(\alpha,E^{(2)}_{j})$,
we know the edge corresponding to $\widetilde{\bm{\Theta}}_{j,(2)}$ must be
rejected in the second iteration of Algorithm 1. Following this kind of
argument, we are able to show that (A.2) implies that all those edges
corresponding to $\\{\widetilde{\bm{\Theta}}_{j,(s)},1\leq s\leq k_{\tau}\\}$
will be rejected according to Algorithm 1. Since the number of rejected edges
is at least $k_{\tau}$, we have $\psi_{j,\alpha}=1$. Second, we show
$\\{\psi_{j,\alpha}=1\\}\subset\\{\alpha_{j}\leq\alpha\\}$. If
$\psi_{j,\alpha}=1\\}$, we know the edges corresponding to
$\\{\widetilde{\bm{\Theta}}_{j,(s)},1\leq s\leq k_{\tau}\\}$ will be rejected,
which immediately imply
$\sqrt{n}|\widetilde{\bm{\Theta}}_{j,(1)}|\geq\widehat{c}(\alpha,E^{(1)}_{j})$.
Regarding the edge corresponding to $\widetilde{\bm{\Theta}}_{j,(2)}$, it must
get rejected in the first two iterations of Algorithm 1. In either cases, we
always have
$\sqrt{n}|\widetilde{\bm{\Theta}}_{j,(2)}|\geq\widehat{c}(\alpha,E^{(2)}_{j})$
due to $E^{(2)}_{j}\subset E^{(1)}_{j}$ and the fact that
$\widehat{c}(\alpha,E)\leq\widehat{c}(\alpha,E^{\prime})$ when $E\subset
E^{\prime}$. Finally, we establish (A.1).

######  Lemma A.1.

Under the same conditions as Lemma 2.1, given some $1\leq j\leq d$, we have
the following results.

  1. (i)

Under the alternative hypothesis $H_{1j}:\|\bm{\Theta}_{j,-j}\|_{0}\geq
k_{\tau}$, we have for any $\alpha\in(0,1)$,

$\lim_{(n,d)\rightarrow\infty}\mathbb{P}(\psi_{j,\alpha}=1)=1.$

  2. (ii)

Additionally, suppose for any $|\bm{\Theta}_{jk}|>0$, we also have
$|\bm{\Theta}_{jk}|\geq c\sqrt{\log d/n}$ for some constant $c>0$. Under the
null hypothesis $H_{0j}:\|\bm{\Theta}_{j,-j}\|_{0}<k_{\tau}$, we then have for
any $u\in(0,1)$,

$\lim_{(n,d)\rightarrow\infty}\mathbb{P}\left({\psi_{j,\alpha}=1}\right)\leq\alpha.$

The proof of the above lemma is deferred to Section D.1. The maximum statistic
used in our testing procedure takes the form of $T_{E}=\max_{(j,k)\in
E}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{jk}|$. In our key proof
procedure, we deal with the case where $E=\\{(j,k):\bm{\Theta}_{jk}=0\\}$.
Since some of the results hold for general $E$, we will work with the general
notations. Specifically, through out Appendices A.1 and A.2, we introduce the
following notations: in order to approximate

$T_{E}:=\max_{(j,k)\in
E}\sqrt{n}\left|(\widehat{\bm{\Theta}}^{\text{d}}_{jk}/\sqrt{\widehat{\bm{\Theta}}^{\text{d}}_{jj}\widehat{\bm{\Theta}}^{\text{d}}_{kk}}-{\bm{\Theta}_{jk}}/{\sqrt{\bm{\Theta}_{jj}\bm{\Theta}_{jk}}})\right|$
(A.3)

by the multiplier bootstrap process

$T_{E}^{\mathcal{B}}:=\max_{(j,k)\in
E}\frac{1}{\sqrt{n~{}\widehat{\bm{\Theta}}_{jj}\widehat{\bm{\Theta}}_{kk}}}\bigg{|}\sum_{i=1}^{n}\widehat{\bm{\Theta}}_{j}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}\widehat{\bm{\Theta}}_{k}-\mathbf{e}_{k})\xi_{i}\bigg{|},$
(A.4)

we define two intermediate processes:

$\displaystyle\breve{T}_{E}$ $\displaystyle:=\max_{(j,k)\in
E}\bigg{|}\frac{1}{\sqrt{n~{}{\bm{\Theta}}_{jj}{\bm{\Theta}}_{kk}}}\sum_{i=1}^{n}{\bm{\Theta}}_{j}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k}-\mathbf{e}_{k})\bigg{|},$
(A.5) $\displaystyle\breve{T}_{E}^{\mathcal{B}}$
$\displaystyle:=\max_{(j,k)\in
E}\bigg{|}\frac{1}{\sqrt{n~{}\bm{\Theta}_{jj}\bm{\Theta}_{kk}}}\sum_{i=1}^{n}{\bm{\Theta}}_{j}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k}-\mathbf{e}_{k})\xi_{i}\bigg{|}.$
(A.6)

###  A.1 Proof of Theorem 5.2

###### Proof of Theorem 5.2.

Given some $j\in\mathcal{H}_{0}$, denote
$N_{0j}=\\{(j,k):\bm{\Theta}_{jk}=0\\}$. By the first part of Lemma A.1, we
have

$\frac{\sum_{j\in\mathcal{B}}{\psi_{j,\alpha}}}{|\mathcal{B}|}\rightarrow
1\quad\text{in probability},$ (A.7)

when $\alpha=\Omega(1/d)$, where
$\mathcal{B}:=\\{j\in\mathcal{H}_{0}^{c}:\forall
k\in\text{supp}(\bm{\Theta}_{j}),|\bm{\Theta}_{jk}|>c\sqrt{{\log d}/{n}}\\}$.
Note that we have

$\mathbb{P}\left(\frac{q|\mathcal{B}|}{d}\leq\widehat{\alpha}\leq
1\right)\geq\mathbb{P}\left({\frac{{q|\mathcal{B}|}/{d}\cdot
d}{\sum_{j\in[d]}\psi_{j,{q|\mathcal{B}|}/{d}}}\leq
q}\right)=\mathbb{P}\left({\frac{|\mathcal{B}|}{\sum_{j\in[d]}\psi_{j,{q|\mathcal{B}|}/{d}}}\leq
1}\right)\rightarrow 1$ (A.8)

in probability, where the first inequality is by (2.2) and the last
convergence in probability is due to $q\frac{|\mathcal{B}|}{d}=\Omega(1/d)$
and (A.7). Rewrite the FDP (with $\widehat{\alpha}$) as

$\displaystyle\mathrm{FDP}(\widehat{\alpha}):=~{}\frac{\sum_{j\in\mathcal{H}_{0}}{\psi_{j,\widehat{\alpha}}}}{\max\left\\{1,\sum_{j\in[d]}\psi_{j,\widehat{\alpha}}\right\\}}=~{}\frac{\widehat{\alpha}d}{\max\left\\{1,\sum_{j\in[d]}\psi_{j,\widehat{\alpha}}\right\\}}\cdot{\frac{\sum_{j\in\mathcal{H}_{0}}\psi_{j,\widehat{\alpha}}}{d_{0}\widehat{\alpha}}}\cdot\frac{d_{0}}{d},$

and notice that

$\frac{\widehat{\alpha}d}{\max\left\\{1,\sum_{j\in[d]}\psi_{j,\widehat{\alpha}}\right\\}}\cdot\frac{d_{0}}{d}\leq\frac{qd_{0}}{d}\leq
q.$

Then it suffices to control the $\mathrm{FDP}(\widehat{\alpha})$ by dealing
with
${\left(\sum_{j\in\mathcal{H}_{0}}\psi_{j,\widehat{\alpha}}\right)}/{d_{0}\widehat{\alpha}}$.
By (A.8), the FDP control problem is now reduced to showing

$\sup_{\alpha\in[\alpha_{L},1]}\frac{\sum_{j\in\mathcal{H}_{0}}\psi_{j,\alpha}}{d_{0}\alpha}\leq
1+o_{\mathbb{P}}(1),$

where $\alpha_{L}=q|\mathcal{B}|/d$, By (D.3) in the proof of the second part
of Lemma A.1, $\psi_{j,\alpha}=1$ implies that $\max_{e\in
N_{0j}}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{e}-\bm{\Theta}^{\star}_{e}|\geq\widehat{c}(\alpha,N_{0j})$,
where
$N_{0j}=\\{(j,k):\bm{\Theta}_{jk}=0\\}=\\{(j,k):\bm{\Theta}^{\star}_{jk}=0\\}$.
Therefore, we have

$\frac{\sum_{j\in\mathcal{H}_{0}}\psi_{j,\alpha}}{d_{0}\alpha}\leq\frac{\sum_{j\in\mathcal{H}_{0}}\mathds{1}{(\max_{e\in
N_{0j}}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{e}-\bm{\Theta}^{\star}_{e}|\geq\widehat{c}(\alpha,N_{0j}))}}{d_{0}\alpha}.$

Hence it suffices to prove that

$\sup_{\alpha\in[\alpha_{L},1]}\left|\frac{\sum_{j\in\mathcal{H}_{0}}\mathds{1}{(\max_{e\in
N_{0j}}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{e}-\bm{\Theta}^{\star}_{e}|\geq\widehat{c}(\alpha,N_{0j}))}}{d_{0}\alpha}-1\right|\rightarrow
0~{}~{}~{}\text{in probability}.$ (A.9)

In order to prove (A.9), we construct a discrete grid of the interval
$[\alpha_{L},1]$. The number of grid points is denoted by $\lambda_{d}$ and
will be decided later. First, we let $t_{1}:=\widehat{c}(1,N_{0j})=0$,
$t_{\lambda_{d}}:=\widehat{c}(\alpha_{L},N_{0j})$. Here
$\widehat{c}(\alpha_{L},N_{0j})=\inf\left\\{t\in\mathbb{R}:\mathbb{P}_{\xi}\left(T^{\mathcal{B}}_{N_{0j}}\leq
t\right)\geq 1-\alpha\right\\}$ is the quantile based on the Gaussian
multiplier bootstrap process and depends on the data $\bm{X}$. Note that the
involving random vectors in the Gaussian multiplier bootstrap process are
Gaussian conditioning on the data $\bm{X}$ and have bounded variances with
probability growing to $1$. Since $\alpha_{L}=\Omega(1/d)$, then by the
maximal inequalities for sub-Gaussian random variables (Lemma 5.2 in van
Handel (2014)), we have $t_{\lambda_{d}}=O(\sqrt{\log d})$ with probability
growing to $1$. Second, note there exists $h_{d}$ such that
$h_{d}t_{\lambda_{d}}=o(1)$ and $t_{\lambda_{d}}/h_{d}=O(\log d)$. Based on
such $h_{d}$, we construct equally spaced sequences
$\\{t_{m}\\}_{m=1}^{\lambda_{d}}$ over the range
$[t_{1},t_{\lambda_{d}}]=[0,t_{\lambda_{d}}]$ with $t_{m}-t_{m-1}=h_{d}$. Then
by setting $\alpha_{m}$ such that $t_{m}=\widehat{c}(\alpha_{m},N_{0j})$, we
obtain a discrete grid $\\{\alpha_{m}\\}_{m=1}^{\lambda_{d}}$ of the interval
$[\alpha_{L},1]$. For such $\alpha_{m},1\leq m\leq\lambda_{d}$, we have

$\displaystyle\max_{1\leq
m\leq\lambda_{d}}\left|\frac{\alpha_{m-1}}{\alpha_{m}}-1\right|$
$\displaystyle=\max_{1\leq
m\leq\lambda_{d}}\left|\frac{\mathbb{P}\left({T^{\mathcal{B}}_{N_{0j}}>t_{m-1}}\right)}{\mathbb{P}\left({T^{\mathcal{B}}_{N_{0j}}>t_{m}}\right)}-1\right|$
$\displaystyle\leq\max_{1\leq
m\leq\lambda_{d}}C^{\prime\prime}(t_{m}-t_{m-1})(t_{m}+1)\exp(C^{\prime}(t_{m}-t_{m-1})(t_{m}+1))=o(1)$
(A.10)

with probability growing to 1, where the first equality holds by the
definition of $\alpha_{m}$, the first inequality holds due to part $2$ and $3$
of Theorem 2.1 in Kuchibhotla et al. (2021) (by first choosing
$r-\epsilon,r+\epsilon$ in part $3$ to be $t_{m-1},t_{m}$ respectively then
letting $r-\epsilon,r$ in part $2$ to be $t_{m-1},t_{m}$ respectively). And
the right hand side of the inequality is $o(1)$ since
$(t_{m}-t_{m-1})t_{m}\leq h_{d}t_{\lambda_{d}}=o(1)$ with probability growing
to 1.

Denote $I_{j}(\alpha)=\mathds{1}{(\max_{e\in
N_{0j}}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{e}-\bm{\Theta}^{\star}_{e}|\geq\widehat{c}(\alpha,N_{0j}))}$.
Then given $\alpha_{m}\leq\alpha\leq\alpha_{m-1}$, for
$m=1,\cdots,\lambda_{d}$, we have

$\displaystyle\frac{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})}{d_{0}\alpha_{m}}\cdot\frac{\alpha_{m}}{\alpha_{m-1}}\leq\frac{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha)}{d_{0}\alpha}\leq\frac{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m-1})}{d_{0}\alpha_{m-1}}\cdot\frac{\alpha_{m-1}}{\alpha_{m}}.$
(A.11)

Hence by (A.10) and (A.11), showing (A.9) is reduced to proving

$\max_{1\leq
m\leq\lambda_{d}}\left|\frac{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})}{d_{0}\alpha_{m}}-1\right|\rightarrow
0,~{}~{}~{}\text{in probability}.$ (A.12)

Then it suffices to show that, for any $\epsilon>0$,

$\mathbb{P}\left({\max_{1\leq
m\leq\lambda_{d}}\left|\frac{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})}{d_{0}\alpha_{m}}-1\right|\geq\epsilon}\right)\rightarrow
0.$

By the union bound argument and Chebyshev’s inequality, we have

$\displaystyle~{}~{}~{}~{}~{}\mathbb{P}\left({\max_{1\leq
m\leq\lambda_{d}}\left|\frac{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})}{d_{0}\alpha_{m}}-1\right|\geq\epsilon}\right)$
$\displaystyle\leq~{}\sum_{m=1}^{\lambda_{d}}\mathbb{P}\left({\left|\frac{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})}{d_{0}\alpha_{m}}-1\right|\geq\epsilon}\right)$
$\displaystyle\leq~{}\sum_{m=1}^{\lambda_{d}}\frac{\mathbb{E}\left[{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})-d_{0}\alpha_{m}}\right]^{2}}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}}$
(A.13)
$\displaystyle\begin{split}&=~{}\underbrace{\sum_{m=1}^{\lambda_{d}}\frac{\sum_{j\in\mathcal{H}_{0}}\mathrm{Var}\left({I_{j}(\alpha_{m})-d_{0}\alpha_{m}}\right)}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}}}_{\mathrm{III}_{1}}+\underbrace{\sum_{m=1}^{\lambda_{d}}\frac{\left(\mathbb{E}\left[{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})-d_{0}\alpha_{m}}\right]\right)^{2}}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}}}_{\mathrm{III}_{2}}\\\
&~{}~{}~{}~{}~{}+\underbrace{\sum_{m=1}^{\lambda_{d}}\frac{\sum_{j_{1},j_{2}\in\mathcal{H}_{0},j_{1}\neq
j_{2}}\mathrm{Cov}\left({I_{j_{1}}(\alpha_{m})},{I_{j_{2}}(\alpha_{m})}\right)}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}}}_{\mathrm{III}_{3}}.\end{split}$
(A.14)

By Lemma A.2 and Lemma A.3, we have

$\displaystyle\mathrm{III}_{1}+\mathrm{III}_{2}+\mathrm{III}_{3}$
$\displaystyle\leq$
$\displaystyle\frac{C^{\prime}t_{\lambda_{d}}}{\epsilon^{2}h_{d}}\left(\frac{d}{d_{0}|\mathcal{B}|}+\eta^{2}(d,n)\right)+\frac{C^{\prime\prime\prime}d}{\epsilon^{2}|\mathcal{B}|d_{0}}\cdot\frac{t_{\lambda_{d}}}{h_{d}}\cdot\left(1+\eta(d,n)d_{0}+\frac{|S|\log
d}{d_{0}p}\right)$ (A.15) $\displaystyle\leq$
$\displaystyle\frac{C_{1}t_{\lambda_{d}}\eta^{2}(d,n)}{\epsilon^{2}h_{d}}+\frac{C_{2}}{\epsilon^{2}\rho
d_{0}}\cdot\frac{t_{\lambda_{d}}}{h_{d}}\cdot\left(1+\eta(d,n)d_{0}+\frac{|S|\log
d}{d_{0}p}\right),$

where we substitute $\zeta_{1}={s(\log d)^{2}}/\sqrt{n}$, $\zeta_{2}=1/d^{2}$
and $\alpha_{L}=q|\mathcal{B}|/d=\Omega(\rho)$ in
$\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})$ of Lemma A.2 and note
$|\mathcal{B}|>0$ then obtain the concise form $\eta(d,n)$ below,

$\displaystyle\eta(d,n)=\frac{(\log d)^{19/6}}{n^{1/6}}+\frac{(\log
d)^{11/6}}{\rho^{1/3}n^{1/6}}+\frac{{s(\log d)^{3}}}{n^{1/2}}+\frac{1}{d}.$

Recall that $t_{\lambda_{d}}=q(\alpha_{L};T^{B}_{N_{0j}})=O\left(\sqrt{\log
d}\right)$ with probability growing to 1 and $t_{\lambda_{d}}/h_{d}=O(\log
d)$. Under Assumption 5.1, we have

$\displaystyle\frac{\log d}{\rho}\left(\frac{(\log
d)^{19/6}}{n^{1/6}}+\frac{(\log d)^{11/6}}{\rho^{1/3}n^{1/6}}+\frac{{s(\log
d)^{3}}}{{n}^{1/2}}\right)=o(1),\quad\frac{\log d}{\rho d_{0}}+\frac{({\log
d})^{2}|S|}{\rho d_{0}^{2}p}=o(1),$

and thus $\mathrm{III}_{1}+\mathrm{III}_{2}+\mathrm{III}_{3}=o(1)$ with
probability growing to 1. Therefore, we have proved $\eqref{eq:fdp_term2}$,
and finally establish the FDP control result below,

$\text{FDP}(\widehat{\alpha})\leq q\frac{d_{0}}{d}+o_{\mathbb{P}}(1).$

In order to establish FDR control, it remains to check the uniformly
integrability of the random variable sequence in (A.12). Note for a sequence
of random variable $R_{1},R_{2},\cdots$, we have
$\sup_{n}\mathbb{E}\left[{|R_{n}|\mathds{1}{(|R_{n}|>x)}}\right]\leq
x^{-1}\sup_{n}\mathbb{E}\left[{R_{n}^{2}}\right]$ by Markov’s inequality. Then
to show the uniform integrability of the random variable sequence
$\\{R_{n}\\}_{n=1}^{\infty}$, where $R_{n}=\max_{1\leq
m\leq\lambda_{d}}\left|\frac{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})}{d_{0}\alpha_{m}}-1\right|$,
it suffices to show $\sup_{n}\mathbb{E}\left[{R^{2}_{n}}\right]<\infty$.
Indeed, we have

$\displaystyle\sup_{n}\mathbb{E}\left[{\left(\max_{1\leq
m\leq\lambda_{d}}\left|\frac{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})}{d_{0}\alpha_{m}}-1\right|\right)^{2}}\right]$
$\displaystyle\leq$
$\displaystyle\sup_{n}\sum_{m=1}^{\lambda_{d}}\frac{\mathbb{E}\left[{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})-d_{0}\alpha_{m}}\right]^{2}}{d_{0}^{2}\alpha_{m}^{2}}$
$\displaystyle=$
$\displaystyle\sup_{n}\epsilon^{2}(\mathrm{III}_{1}+\mathrm{III}_{2}+\mathrm{III}_{3}).$

Since $\mathrm{III}_{1}+\mathrm{III}_{2}+\mathrm{III}_{3}=o(1)$ with
probability growing to 1, we immediately have
$\sup_{n}\mathbb{E}\left[{R^{2}_{n}}\right]<\infty$, thus finally establish
the FDR control result:

$\lim_{(n,d)\rightarrow\infty}\text{FDR}\leq q\frac{d_{0}}{d}.$

∎

###  A.2 Ancillary lemmas for Theorem 5.2

######  Lemma A.2.

Recalling the definitions of $\mathrm{III}_{1},\mathrm{III}_{2}$ in
(LABEL:eq:seq_III_3), we have

$\mathrm{III}_{1}+\mathrm{III}_{2}\leq\frac{C^{\prime}t_{\lambda_{d}}}{\epsilon^{2}h_{d}}\left(\frac{1}{\rho
d_{0}}+\eta^{2}(d,n,\zeta_{1},\zeta_{2},\alpha_{L})\right),$

where $\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})=O\big{(}\frac{(\log
d)^{19/6}}{n^{1/6}}+\frac{(\log
d)^{11/6}}{n^{1/6}\alpha_{L}^{1/3}}+\zeta_{1}\log
d+\frac{\zeta_{2}}{\alpha_{L}}\big{)}$ with $\zeta_{1}={s(\log
d)^{2}}/\sqrt{n}$, $\zeta_{2}=1/d^{2}$.

###### Proof of Lemma A.2.

First note the definitions of $T_{E},\breve{T}_{E},T_{E}^{\mathcal{B}}$ and
$\breve{T}_{E}^{\mathcal{B}}$ in (A.3), (A.5), (A.4) and (A.6) respectively,
then we apply Proposition C.2 to
$T=T_{E},T_{\bm{Y}}=\breve{T}_{E},T^{\mathcal{B}}=T_{E}^{\mathcal{B}},T_{\mathcal{W}}=\breve{T}_{E}^{\mathcal{B}}$
with $E=N_{0j}$. And we can find the terms $\zeta_{1},\zeta_{2}$ in (C.4),
(C.5) to be ${s(\log d)^{2}}/\sqrt{n},1/d^{2}$ respectively, due to(D.25) and
(D.26) (i.e., the bound on the differences $T_{E}-T_{0}$,
$T_{E}^{\mathcal{B}}-T_{0}^{\mathcal{B}}$) in the proof of Lemma 2.1. Thus we
have

$\left|\frac{\mathbb{P}(\max_{e\in
N_{0j}}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{e}-\bm{\Theta}^{\star}_{e}|\geq\widehat{c}(\alpha,N_{0j}))}{\alpha}-1\right|=\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L}),$
(A.16)

where $\bm{\Theta}^{\star}_{e}=0,e\in N_{0j}$ and
$\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})=O\left(\frac{(\log
d)^{19/6}}{n^{1/6}}+\zeta_{1}\log d+\frac{\zeta_{2}}{\alpha_{L}}\right)$ with
$\zeta_{1}={s(\log d)^{2}}/\sqrt{n}$, $\zeta_{2}=1/d^{2}$. Recalling the
definition of $\mathrm{III}_{2}$ in (LABEL:eq:seq_III_3), we have

$\mathrm{III}_{2}=\sum_{m=1}^{\lambda_{d}}\frac{\left(\mathbb{E}\left[{\sum_{j\in\mathcal{H}_{0}}I_{j}(\alpha_{m})-d_{0}\alpha_{m}}\right]\right)^{2}}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}},$

where $I_{j}(\alpha)=\mathds{1}{(\max_{e\in
N_{0j}}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{e}-\bm{\Theta}^{\star}_{e}|\geq\widehat{c}(\alpha,N_{0j}))}$.
Note that $\alpha_{m}\in[\alpha_{L},1],$ $\forall~{}1\leq m\leq\lambda_{d}$,
then we arrive at the following bound

$\mathrm{III}_{2}\leq\frac{\lambda_{d}}{\epsilon^{2}}\cdot\eta^{2}(d,n,\zeta_{1},\zeta_{2},\alpha_{L})\leq\frac{t_{\lambda_{d}}}{\epsilon^{2}h_{d}}\cdot\eta^{2}(d,n,\zeta_{1},\zeta_{2},\alpha_{L})$
(A.17)

up to some constant, where the first inequality holds by (A.16). As for the
second inequality, we recall the construction of
$\\{t_{m}\\}_{m=1}^{\lambda_{d}}$ (over the course of derivations from (A.9)
to (A.10)) in the proof of Theorem 5.2 thus note $\alpha_{1}=1,t_{1}=0$ and
$t_{\lambda_{d}}-t_{1}=\sum_{m=2}^{\lambda_{d}}(t_{m}-t_{m-1})=(\lambda_{d}-1)h_{d}$.
Regarding the term $\mathrm{III}_{1}$, we have

$\displaystyle\mathrm{III}_{1}$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{\lambda_{d}}\frac{\sum_{j\in\mathcal{H}_{0}}\mathrm{Var}\left({I_{j}(\alpha_{m})-d_{0}\alpha_{m}}\right)}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}}$
(A.18) $\displaystyle=$
$\displaystyle\sum_{m=1}^{\lambda_{d}}\frac{\sum_{j\in\mathcal{H}_{0}}\mathbb{E}(I_{j}(\alpha_{m}))(1-\mathbb{E}(I_{j}(\alpha_{m}))}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}}\leq\frac{1}{\epsilon^{2}d_{0}}\sum_{m=1}^{\lambda_{d}}\frac{C}{\alpha_{m}}\leq\frac{C}{\epsilon^{2}d_{0}\alpha_{L}}\cdot\frac{t_{\lambda_{d}}}{h_{d}},$

where the first inequality holds due to $\eqref{eq:fdp_cgmb}$ and the second
inequality holds since $\alpha_{m}\geq\alpha_{L}$ $\forall~{}1\leq
m\leq\lambda_{d}$ and $t_{\lambda_{d}}=(\lambda_{d}-1)h_{d}$. Therefore,
combining (A.17) with (A.18), we obtain

$\displaystyle\mathrm{III}_{1}+\mathrm{III}_{2}\leq\frac{1}{\epsilon^{2}}\cdot\frac{t_{\lambda_{d}}}{h_{d}}\left(\frac{C}{d_{0}\alpha_{L}}+\eta^{2}(d,n,\zeta_{1},\zeta_{2},\alpha_{L})\right)\leq\frac{C^{\prime}t_{\lambda_{d}}}{\epsilon^{2}h_{d}}\left(\frac{1}{\rho
d_{0}}+\eta^{2}(d,n,\zeta_{1},\zeta_{2},\alpha_{L})\right)$

for some constant $C^{\prime}$, where the second inequality holds by the
definition $\alpha_{L}=q|\mathcal{B}|/d$ in the proof of Theorem 5.2 and the
definition $\rho=|\mathcal{B}|/d$ in Section 5. ∎

######  Lemma A.3.

Recalling the definition of $\mathrm{III}_{3}$ in (LABEL:eq:seq_III_3), we
have

$\mathrm{III}_{3}\leq\frac{C^{\prime\prime\prime}t_{\lambda_{d}}}{\rho\epsilon^{2}d_{0}h_{d}}\left(1+\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})d_{0}+\frac{|S|\log
d}{d_{0}p}\right),$

where $\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})=O\big{(}\frac{(\log
d)^{19/6}}{n^{1/6}}+\frac{(\log
d)^{11/6}}{n^{1/6}\alpha_{L}^{1/3}}+\zeta_{1}\log
d+\frac{\zeta_{2}}{\alpha_{L}}\big{)}$ with $\zeta_{1}={s(\log
d)^{2}}/\sqrt{n}$, $\zeta_{2}=1/d^{2}$.

###### Proof of Lemma A.3.

Note that $\mathrm{III}_{3}$ in (LABEL:eq:seq_III_3) equals

$\displaystyle\begin{split}\mathrm{III}_{3}&=\sum_{m=1}^{\lambda_{d}}\frac{\sum_{j_{1},j_{2}\in\mathcal{H}_{0},j_{1}\neq
j_{2}}\mathrm{Cov}\left({I_{j_{1}}(\alpha_{m})},{I_{j_{2}}(\alpha_{m})}\right)}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}},\\\
&\quad\text{ where $I_{j}(\alpha)=\mathds{1}{(\max_{e\in
N_{0j}}\sqrt{n}|\widetilde{\bm{\Theta}}^{\text{d}}_{e}-\bm{\Theta}^{\star}_{e}|\geq\widehat{c}(\alpha,N_{0j}))}$
}\end{split}$ (A.19)

for $j\in\\{j_{1},j_{2}\\}$. To quantify the covariance between
$I_{j_{1}}(\alpha_{m})$ and $I_{j_{2}}(\alpha_{m})$ for
$j_{1},j_{2}\in\mathcal{H}_{0},j_{1}\neq j_{2}$, we define

$W_{j}(\alpha)=\mathds{1}{(\max_{e\in N_{0j}}|Z_{e}|\geq c(\alpha,N_{0j}))},$
(A.20)

where $(Z_{e})_{e\in E}$ (with $E=N_{0j}$) is a Gaussian random vector and
shares the same mean vector and covariance matrix as the term
$(\frac{1}{\sqrt{n~{}{\bm{\Theta}}_{jj}{\bm{\Theta}}_{kk}}}\sum_{i=1}^{n}{\bm{\Theta}}_{j}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k}-\mathbf{e}_{k}))_{(j,k)\in
E}$ in $\breve{T}_{E}$. Here $\breve{T}_{E}$ (with $E=N_{0j}$) has the
explicit form below

$\breve{T}_{E}=\max_{(j,k)\in
E}\frac{1}{\sqrt{n~{}{\bm{\Theta}}_{jj}{\bm{\Theta}}_{kk}}}\bigg{|}\sum_{i=1}^{n}{\bm{\Theta}}_{j}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k}-\mathbf{e}_{k})\bigg{|}.$

Remark here $\breve{T}_{E}$ corresponds to the term $T_{\bm{Y}}$ in
Proposition C.2 and $\max_{e\in E}|Z_{e}|$ corresponds to the term
$T_{\bm{Z}}$ in Proposition C.1. And $c(\alpha,N_{0j})$ is the corresponding
Gaussian maxima quantile $q(\alpha;T_{\bm{Z}})$ (which does not need to be
computed). Since $\mathbb{P}(T_{\bm{Z}}>q(\alpha;T_{\bm{Z}}))=\alpha$, we
immediately have have $\mathbb{E}[W_{j}(\alpha)]=\mathbb{P}\left({\max_{e\in
N_{0j}}\sqrt{n}|Z_{e}|\geq c(\alpha,N_{0j})}\right)=\alpha$.

Now we replace $I_{j_{1}}(\alpha),I_{j_{2}}(\alpha)$ in $\mathrm{III}_{3}$ by
$W_{j_{1}}(\alpha),W_{j_{2}}(\alpha)$ and define $\mathrm{III}^{\prime}_{3}$
as

$\mathrm{III}^{\prime}_{3}:=\sum_{m=1}^{\lambda_{d}}\frac{\sum_{j_{1},j_{2}\in\mathcal{H}_{0},j_{1}\neq
j_{2}}\mathrm{Cov}\left({W_{j_{1}}(\alpha_{m})},{W_{j_{2}}(\alpha_{m})}\right)}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}}.$
(A.21)

To bound the difference between $\mathrm{III}_{3}$ and
$\mathrm{III}^{\prime}_{3}$, we first note $\operatorname{\rm
Cov}(I_{j_{1}}(\alpha),I_{j_{2}}(\alpha))=\mathbb{E}[I_{j_{1}}(\alpha)I_{j_{2}}(\alpha)]-\mathbb{E}[I_{j_{1}}(\alpha)]\mathbb{E}[I_{j_{2}}(\alpha)]$
then separately deal with the term
$\left|\mathbb{E}\left[{I_{j_{1}}(\alpha)I_{j_{2}}(\alpha)}\right]-\mathbb{E}\left[{W_{j_{1}}(\alpha)W_{j_{2}}(\alpha)}\right]\right|$
and the term
$\left|\mathbb{E}\left[{I_{j_{1}}(\alpha)}\right]\mathbb{E}\left[{I_{j_{2}}(\alpha)}\right]-\mathbb{E}\left[{W_{j_{1}}(\alpha)}\right]\mathbb{E}\left[{G_{j_{2}}(\alpha)}\right]\right|$.

By Lemma A.5, we have up to some constant factor,

$\frac{\left|\mathbb{E}\left[{I_{j_{1}}(\alpha)I_{j_{2}}(\alpha)}\right]-\mathbb{E}\left[{W_{j_{1}}(\alpha)W_{j_{2}}(\alpha)}\right]\right|}{\alpha^{2}}\leq\frac{\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})}{\alpha}.$

Applying the same strategy to the term
$\mathbb{E}\left[{I_{j_{1}}(\alpha)}\right]\mathbb{E}\left[{I_{j_{2}}(\alpha)}\right]$,
we obtain

$\frac{\left|\mathbb{E}\left[{I_{j_{1}}(\alpha)}\right]\mathbb{E}\left[{I_{j_{2}}(\alpha)}\right]-\mathbb{E}\left[{W_{j_{1}}(\alpha)}\right]\mathbb{E}\left[{G_{j_{2}}(\alpha)}\right]\right|}{\alpha^{2}}\leq\frac{\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})}{\alpha}.$

Combining the above two inequalities, and noting the definition of
$\mathrm{III}^{\prime}_{3}$ in (A.21), we derive the following bound on the
difference between $\mathrm{III}_{3}$ and $\mathrm{III}^{\prime}_{3}$,

$\left|\mathrm{III}_{3}-\mathrm{III}^{\prime}_{3}\right|\leq\frac{1}{\epsilon^{2}}\sum_{m=1}^{\lambda_{d}}\frac{\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})}{\alpha_{m}}\leq\frac{C^{\prime}t_{\lambda_{d}}}{\rho\epsilon^{2}h_{d}}\cdot\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L}).$

where the second inequality holds due to the fact $\alpha_{m}\geq\alpha_{L}$
$\forall~{}1\leq m\leq\lambda_{d}$ and $t_{\lambda_{d}}=(\lambda_{d}-1)h_{d}$,
the definition $\alpha_{L}=q|\mathcal{B}|/d$ in the proof of Theorem 5.2, and
the definition $\rho=|\mathcal{B}|/d$ in Section 5.

The above bound on $\left|\mathrm{III}_{3}-\mathrm{III}^{\prime}_{3}\right|$,
when combined with Lemma A.4, immediately establishes

$\displaystyle\mathrm{III}_{3}$ $\displaystyle\leq$
$\displaystyle\frac{C^{\prime}t_{\lambda_{d}}}{\rho\epsilon^{2}h_{d}}\cdot\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})+\frac{C^{\prime\prime}t_{\lambda_{d}}}{\rho\epsilon^{2}d_{0}h_{d}}\left(1+C_{\bm{\Theta}}\frac{|S|\log
d}{d_{0}p}\right)$ $\displaystyle\leq$
$\displaystyle\frac{C^{\prime\prime\prime}t_{\lambda_{d}}}{\rho\epsilon^{2}d_{0}h_{d}}\left(1+\eta(d,n,\zeta_{1},\zeta_{2},\alpha_{L})d_{0}+\frac{|S|\log
d}{d_{0}p}\right),$

for some constant $C^{\prime\prime\prime}$. ∎

######  Lemma A.4.

Recalling the term $\mathrm{III}^{\prime}_{3}$ from (A.21) in the proof of
Lemma A.3, we have

$\mathrm{III}^{\prime}_{3}=\sum_{m=1}^{\lambda_{d}}\frac{\sum_{j_{1},j_{2}\in\mathcal{H}_{0},j_{1}\neq
j_{2}}\operatorname{\rm
Cov}(W_{j_{1}}(\alpha_{m}),W_{j_{2}}(\alpha_{m}))}{\epsilon^{2}d_{0}^{2}\alpha_{m}^{2}}\leq\frac{C^{\prime\prime}t_{\lambda_{d}}}{\rho\epsilon^{2}d_{0}h_{d}}\left(1+C_{\bm{\Theta}}\frac{|S|\log
d}{d_{0}p}\right).$

###### Proof of Lemma A.4.

Similarly as in the proof of Lemma A.3, we define $(Z_{e})_{e\in
N_{0j_{1}}\cup N_{0j_{2}}}$ to be jointly Gaussian such that this
$(|N_{0j_{1}}|+|N_{0j_{2}}|)$-dimensional Gaussian random vector shares the
same mean vector and covariance matrix as the term
$(\frac{1}{\sqrt{n~{}{\bm{\Theta}}_{jj}{\bm{\Theta}}_{kk}}}\sum_{i=1}^{n}{\bm{\Theta}}_{j}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k}-\mathbf{e}_{k}))_{(j,k)\in
N_{0j_{1}}\cup N_{0j_{2}}}$. Note that the two sub-vectors $(Z_{e})_{e\in
N_{0j_{1}}}$ and $(Z_{e})_{e\in N_{0j_{1}}}$ are generally dependent. Then we
define $(Z^{\prime}_{e})_{e\in N_{0j_{1}}},(Z^{\prime}_{e})_{e\in N_{0j_{2}}}$
to be two Gaussian random vectors such that

$(Z^{\prime}_{e})_{e\in N_{0j_{1}}}\stackrel{{\scriptstyle
d}}{{=}}(Z_{e})_{e\in N_{0j_{1}}},~{}~{}(Z^{\prime}_{e})_{e\in
N_{0j_{2}}}\stackrel{{\scriptstyle d}}{{=}}(Z_{e})_{e\in
N_{0j_{2}}}~{}~{}\text{ and }(Z^{\prime}_{e})_{e\in
N_{0j_{1}}}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}(Z^{\prime}_{e})_{e\in N_{0j_{2}}}.$ (A.22)

Recalling the definition of $W_{j}(\alpha)$ in (A.20):
$W_{j}(\alpha)=\mathds{1}{(\max_{e\in N_{0j}}|Z_{e}|\geq c(\alpha,N_{0j}))}$,
we thus have the following,

$\displaystyle\mathrm{IV}_{j_{1}j_{2}}(\alpha)$ $\displaystyle:=$
$\displaystyle\frac{\left|\operatorname{\rm
Cov}(W_{j_{1}}(\alpha_{m}),W_{j_{2}}(\alpha_{m}))\right|}{\alpha^{2}}$
$\displaystyle=$
$\displaystyle\frac{\left|\mathbb{E}\left[{W_{j_{1}}(\alpha)W_{j_{2}}(\alpha)}\right]-\mathbb{E}\left[{W_{j_{1}}(\alpha)}\right]\mathbb{E}\left[{W_{j_{2}}(\alpha)}\right]\right|}{\alpha^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{\alpha^{2}}\Big{|}\mathbb{P}(\max_{e\in
N_{0j_{1}}}|Z_{e}|\geq c(\alpha,N_{0j_{1}}),\max_{e\in N_{0j_{2}}}|Z_{e}|\geq
c(\alpha,N_{0j_{2}}))-$ $\displaystyle\mathbb{P}(\max_{e\in
N_{0j_{1}}}|Z^{\prime}_{e}|\geq c(\alpha,N_{0j_{1}}),\max_{e\in
N_{0j_{2}}}|Z^{\prime}_{e}|\geq c(\alpha,N_{0j_{2}}))\Big{|}$ (A.24)
$\displaystyle=$
$\displaystyle\frac{1}{\alpha^{2}}\Big{|}\mathbb{P}(\max_{e\in
N_{0j_{1}}}|Z_{e}|\geq t,\max_{e\in N_{0j_{2}}}|Z_{e}|\geq
t)-\mathbb{P}(\max_{e\in N_{0j_{1}}}|Z^{\prime}_{e}|\geq t,\max_{e\in
N_{0j_{2}}}|Z^{\prime}_{e}|\geq t)\Big{|}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha^{2}}\Big{|}\mathbb{P}(\max_{e\in N_{0j_{1}}\cup
N_{0j_{2}}}|Z_{e}|\geq t)-\mathbb{P}(\max_{e\in N_{0j_{1}}\cup
N_{0j_{2}}}|Z^{\prime}_{e}|\geq t)\Big{|},$

where the third equality follows due to the construction of $(Z_{e})_{e\in
N_{0j_{1}}\cup N_{0j_{2}}},(Z^{\prime}_{e})_{e\in N_{0j_{1}}\cup N_{0j_{2}}}$.
Note that in the fourth equality, we assume
$c(\alpha,N_{0j_{1}})=c(\alpha,N_{0j_{2}}):=t$ without loss of generality,
since we can rescale one of the maximum statistic by rescaling the Gaussian
random vectors. Remark that the scaling will not break down the application of
Theorem 3.2, which will be explained in detail later in this proof. The last
inequality holds by (A.22) and the fact that $\mathbb{P}\left({A\cap
B}\right)=\mathbb{P}\left({A}\right)+\mathbb{P}\left({B}\right)-\mathbb{P}\left({A\cup
B}\right)$.

Notice that we can apply the Cramér-type Gaussian comparison bound with
$\ell_{0}$ norm to control (A.24). Specifically, we first figure out the
difference between the covariance matrices of $(Z_{e})_{e\in N_{0j_{1}}\cup
N_{0j_{2}}}$ and $(Z^{\prime}_{e})_{e\in N_{0j_{1}}\cup N_{0j_{2}}}$. Denote
the covariance matrices by $\bm{\Sigma}^{Z}$ and $\bm{\Sigma}^{Z^{\prime}}$
respectively. As these two Gaussian random vectors have two sub-vectors, we
write their covariance matrices in a block form

$\bm{\Sigma}^{Z}=\left(\begin{array}[]{cc}\bm{\Sigma}^{Z}_{11}&\bm{\Sigma}^{Z}_{12}\\\
\bm{\Sigma}^{Z}_{21}&\bm{\Sigma}^{Z}_{22}\end{array}\right),~{}~{}\bm{\Sigma}^{Z^{\prime}}=\left(\begin{array}[]{cc}\bm{\Sigma}^{Z^{\prime}}_{11}&\bm{O}\\\
\bm{O}&\bm{\Sigma}^{Z^{\prime}}_{22}\end{array}\right).$

where $\bm{\Sigma}^{Z^{\prime}}$ is block diagonal due to (A.22). Note that we
also have $\bm{\Sigma}^{Z}_{11}=\bm{\Sigma}^{Z^{\prime}}_{11}$ and
$\bm{\Sigma}^{Z}_{22}=\bm{\Sigma}^{Z^{\prime}}_{22}$. Then we have

$\bm{\Sigma}^{Z}-\bm{\Sigma}^{Z^{\prime}}=\left(\begin{array}[]{cc}\bm{O}&\bm{\Sigma}^{Z}_{12}\\\
\bm{\Sigma}^{Z}_{21}&\bm{O}\end{array}\right).$ (A.25)

Throughout the following proof, we assume $\bm{\Theta}_{jj}=1,j\in[d]$ without
loss of generality, since the standardized version is considered in
$\breve{T}_{E}$ (A.5). Recall that $(Z_{e})_{e\in N_{0j_{1}}\cup N_{0j_{2}}}$
shares the same covariance structure as $(Y_{e})_{e\in N_{0j_{1}}\cup
N_{0j_{2}}}$ where $Y_{e}$ (with $e=(j,k)$) is defined as

$Y_{e}:=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}{\bm{\Theta}}_{j}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k}-\mathbf{e}_{k}).$

Then we are ready to calculate the covariance matrix $\bm{\Sigma}^{Z}$.
Specifically, we compute the entries in each block. Regarding the block
$\bm{\Sigma}^{Z}_{11}$, for any $k,k^{\prime}\in N_{0j_{1}}$ where
$N_{0j_{1}}=\\{k:\bm{\Theta}_{j_{1}k}=0\\}$, we have the corresponding
$(k,k^{\prime})$ entry in $\bm{\Sigma}^{Z}_{11}$ equals

$\operatorname{\rm
Cov}({\bm{\Theta}}_{j_{1}}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k}-\mathbf{e}_{k}),{\bm{\Theta}}_{j_{1}}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k^{\prime}}-\mathbf{e}_{k^{\prime}}))=\bm{\Theta}_{j_{1}j_{1}}\bm{\Theta}_{kk^{\prime}}+\bm{\Theta}_{j_{1}k}\bm{\Theta}_{j_{1}k^{\prime}}=\bm{\Theta}_{kk^{\prime}},$
(A.26)

by applying Isserlis’ theorem (Isserlis, 1918) and noting
$\bm{\Theta}_{j_{1}k}=\bm{\Theta}_{j_{1}k^{\prime}}=0$. Similar results hold
for the block $\bm{\Sigma}^{Z}_{22}$. Regarding the block
$\bm{\Sigma}^{Z}_{12}$, consider $k_{1}\in N_{0j_{1}},k_{2}\in N_{0j_{2}}$,
then we have the corresponding $(k_{1},k_{2})$ entry in the block equals

$\operatorname{\rm
Cov}({\bm{\Theta}}_{j_{1}}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k_{1}}-\mathbf{e}_{k_{1}}),{\bm{\Theta}}_{j_{2}}^{\top}(\bm{X}_{i}\bm{X}_{i}^{\top}{\bm{\Theta}}_{k_{2}}-\mathbf{e}_{k_{2}}))=\bm{\Theta}_{j_{1}j_{2}}\bm{\Theta}_{k_{1}k_{2}}+\bm{\Theta}_{j_{1}k_{2}}\bm{\Theta}_{j_{2}k_{1}}.$
(A.27)

Now we have fully characterized the covariance matrix $\bm{\Sigma}^{Z}$ and
the covariance matrix difference in (A.25) for any
$j_{1},j_{2}\in\mathcal{H}_{0},j_{1}\neq j_{2}$. Specifically, we have
$||\bm{\Sigma}^{Z}-\bm{\Sigma}^{Z^{\prime}}||_{0}=||\bm{\Sigma}_{12}^{Z}||_{0}=\sum_{k_{1}\in
N_{0j_{1}},k_{2}\in
N_{0j_{2}}}\mathds{1}{(\bm{\Theta}_{j_{1}j_{2}}\bm{\Theta}_{k_{1}k_{2}}+\bm{\Theta}_{j_{1}k_{2}}\bm{\Theta}_{j_{2}k_{1}}\neq
0)}$. Based on whether $\bm{\Theta}_{j_{1}j_{2}}$ is zero or not, we consider
the following two cases then handle them separately:

  * •

Case 1: $\bm{\Theta}_{j_{1}j_{2}}=0$. If $k_{1}=k_{2}$, then we have the
covariance matrix entry (A.27) equal zero; If $k_{1}\neq k_{2}$, then (A.27)
is nonzero only if $\bm{\Theta}_{j_{1}k_{2}}\neq
0,\bm{\Theta}_{j_{2}k_{1}}\neq 0$ (i.e., $k_{2}\notin N_{0j_{1}}$,
$k_{1}\notin N_{0j_{2}}$). By the fact
$j_{1},j_{2}\in\mathcal{H}_{0},j_{1}\neq j_{2}$ and the definition of
$\mathcal{H}_{0}=\\{j:\|\bm{\Theta}_{j,-j}\|_{0}<k_{\tau}\\}$, we have
$\\#\\{(k_{1},k_{2}):k_{1}\neq k_{2},\bm{\Theta}_{j_{1}k_{2}}\neq
0,\bm{\Theta}_{j_{2}k_{1}}\neq 0\\}\leq k_{\tau}^{2}$. Hence
$||\bm{\Sigma}^{Z}-\bm{\Sigma}^{Z^{\prime}}||_{0}\leq k_{\tau}^{2}$.

  * •

Case 2: $\bm{\Theta}_{j_{1}j_{2}}\neq 0$. The covariance matrix entry (A.27)
is nonzero only if $\bm{\Theta}_{j_{1}k_{2}}\neq
0,\bm{\Theta}_{j_{2}k_{1}}\neq 0$ (i.e., $k_{2}\notin N_{0j_{1}}$,
$k_{1}\notin N_{0j_{2}}$) or $\bm{\Theta}_{k_{1}k_{2}}\neq 0$.

We start from the simpler case, i.e., Case 2 where
$\bm{\Theta}_{j_{1}j_{2}}\neq 0$. Simply, we obtain

$\mathrm{IV}_{j_{1}j_{2}}(\alpha)=\frac{|\operatorname{\rm
Cov}(W_{j_{1}}(\alpha),W_{j_{2}}(\alpha))|}{\alpha^{2}}\leq\frac{\operatorname{{\rm
Var}}(W_{j_{1}}(\alpha))}{\alpha^{2}}+\frac{\operatorname{{\rm
Var}}(W_{j_{2}}(\alpha))}{\alpha^{2}}\leq\frac{C}{\alpha},$

for some constant $C$ since
$\mathrm{Var}\left({W_{j}(\alpha)}\right)=\mathbb{E}\left[{W_{j}(\alpha)}\right](1-\mathbb{E}\left[{W_{j}(\alpha)}\right])=\alpha(1-\alpha)$
for $j=j_{1},j_{2}$. For a fixed $j_{1}$, we also know that
$|\\{j_{2}\in\mathcal{H}_{0}:j_{2}\neq j_{1},\bm{\Theta}_{j_{1}j_{2}}\neq
0\\}|<k_{\tau}$. Then we have

$\sum_{m=1}^{\lambda_{d}}\sum_{\bm{\Theta}_{j_{1}j_{2}}\neq
0}\frac{\mathrm{IV}_{j_{1}j_{2}}(\alpha_{m})}{\epsilon^{2}d^{2}_{0}}\leq\sum_{m=1}^{\lambda_{d}}\frac{d_{0}k_{\tau}}{\epsilon^{2}d_{0}^{2}}\cdot\frac{C}{\alpha_{m}}\leq\frac{1}{\epsilon^{2}d_{0}}\sum_{m=1}^{\lambda_{d}}\frac{C^{\prime}}{\alpha_{m}},$
(A.28)

where the last inequality holds due to the same derivations for
$\mathrm{III}_{1}$ in the proof of Lemma A.2.

Regarding Case $1$ where $\bm{\Theta}_{j_{1}j_{2}}=0$, we will give a more
careful treatment to $\mathrm{IV}_{j_{1}j_{2}}(\alpha)$ in (A.2). Due to the
discussion about Case 1, we have
$||\bm{\Sigma}^{Z}-\bm{\Sigma}^{Z^{\prime}}||_{0}\leq k_{\tau}^{2}$. This fact
will be utilized to derive a nice bound on $\mathrm{III^{\prime}_{3}}$.
Indeed, we can apply Theorem 3.2 to (A.24) (with $U$ and $V$ chosen to be
$Z_{e})_{e\in N_{0j_{1}}\cup N_{0j_{2}}}$ and $(Z^{\prime}_{e})_{e\in
N_{0j_{1}}\cup N_{0j_{2}}}$ respectively) and obtain

$\mathrm{IV}_{j_{1}j_{2}}(\alpha)\leq\frac{\log d}{\alpha
p}\left(\sum_{k_{1}\in N_{0j_{1}},k_{2}\in N_{0j_{2}},k_{1}\neq
k_{2}}\mathds{1}{(\bm{\Theta}_{j_{1}k_{2}}\bm{\Theta}_{j_{2}k_{1}}\neq
0)}\right).$ (A.29)

when $\bm{\Theta}_{j_{1}j_{2}}=0$ (i.e., under Case 1). Recall Theorem 3.2
assumes for Gaussian random vectors $U$ and $V$, there exists a disjoint
$\mathfrak{p}$-partition of nodes