# High-Dimensional Inference for Cluster-Based Graphical Models

Carson Eisenach Florentina Bunea Yang Ning Claudiu Dinicu Department of
Operations Research and Financial Engineering, Princeton University,
Princeton, NJ 08544, USA; e-mail<EMAIL_ADDRESS>of
Statistical Science, Cornell University, Ithaca, NY 14850, USA; e-mail:
<EMAIL_ADDRESS>of Statistical Science, Cornell University,
Ithaca, NY 14850, USA; e-mail<EMAIL_ADDRESS>of Statistical
Science, Cornell University, Ithaca, NY 14850, USA; e-mail<EMAIL_ADDRESS>
###### Abstract

Motivated by modern applications in which one constructs graphical models
based on a very large number of features, this paper introduces a new class of
cluster-based graphical models. Unlike standard graphical models, variable
clustering is applied as an initial step for reducing the dimension of the
feature space. We employ model assisted clustering, in which the clusters
contain features that are similar to the same unobserved latent variable. Two
different cluster-based Gaussian graphical models are considered: the latent
variable graph, corresponding to the graphical model associated with the
unobserved latent variables, and the cluster-average graph, corresponding to
the vector of features averaged over clusters. We derive estimates tailored to
these graphs, with the goal of pattern recovery under false discovery rate
(FDR) control. Our study reveals that likelihood based inference for the
latent graph is analytically intractable, and we develop alternative
estimation and inference strategies. We replace the likelihood of the data by
appropriate empirical risk functions that allow for valid inference in both
graphical models under study. Our main results are Berry-Esseen central limit
theorems for the proposed estimators, which are proved under weaker
assumptions than those employed in the existing literature on Gaussian
graphical model inference. As a corollary of the main results, we make
explicit the implications of the asymptotic approximations on graph recovery
under FDR control, and show when it can be controlled asymptotically. Our
analysis takes into account the uncertainty induced by the initial clustering
step. We find that in the model-assisted clustering framework, the errors
induced by clustering are asymptotically ignorable in the follow-up analysis,
under no further restrictions on the parameter space for which inference is
valid. The theoretical properties of the proposed procedures are verified on
simulated data and an fMRI data analysis.

Keyword: Berry-Esseen bound, Graphical model, High-dimensional inference,
Clustering, False discovery rate

##  1 Introduction

High dimensional graphical models have become increasingly popular, over the
last several decades, for understanding independence and conditional
independence relationships among components of high dimensional random
vectors. The challenges posed by the estimation and statistical analysis of a
graphical model with many more nodes than the the number of observations led
to renewed interest in these models, and to a very large body of literature,
including Meinshausen and Bühlmann (2006); Yuan and Lin (2007); Friedman et
al. (2008); Verzelen (2008); Lam and Fan (2009); Rothman et al. (2008); Peng
et al. (2009); Ravikumar et al. (2011a); Yuan (2010); Cai et al. (2011); Sun
and Zhang (2012); Liu et al. (2012); Xue and Zou (2012); Ning and Liu (2013);
Jankova and van de Geer (2014); Cai et al. (2016); Janková and van de Geer
(2017); Fan et al. (2017), to give only an incomplete list. However, in
practice, when the dimension is very large and the sample size is small, the
dependency among variables may become very weak, if it exists, and difficult
to detect without auxiliary information. Moreover, when the number of
variables is in the tens of thousands, it is difficult to form an opinion of
their overall dependency structure, even at the visual level, from a graph
estimated via a graphical model.

A solution to these issues is to employ an initial dimension reduction
procedure on the very high dimensional vector. For example, in neuroscience
applications, a typical functional magnetic resonance image (fMRI) consists of
blood-oxygen-level-dependent (BOLD) activities measured at 200,000+ voxels of
the brain, over a period of time. Instead of analyzing voxel-level data
directly, scientists routinely cluster voxels into several regions of interest
(ROI) with homogeneous functions based on the domain knowledge, and then carry
out the analysis at the ROI-level. In this example, the group structure of
variables may boost the dependency signals, in the context of graphical
models. Similar pre-processing steps are used in other application domains,
such as genomics, finance and economics.

Motivated by a rich set of applications, we consider variable clustering as
the initial dimension reduction step applied to the observed vector
$\bm{X}=:(X_{1},\ldots,X_{d})\in\mathbb{R}^{d}$. To the best of our knowledge,
very little is known about the effect of clustering on down-stream analysis
and, consequently, on the induced graphical models. Our contribution is the
provision of a framework that allows for such an analysis. We introduce
cluster-based graphical models and show how they can be estimated, under FDR
control.

The model is built on the assumption that the observed variables
$\bm{X}=(X_{1},\ldots,X_{d})\in\mathbb{R}^{d}$ can be partitioned into $K$
unknown clusters $G^{*}=\\{G^{*}_{1},\ldots,G^{*}_{K}\\}$ such that variables
in the same cluster share the same behavior. Following the intuition behind
the widely-used $K$-means type procedures, we define a population-level
cluster as a group of variables that are noise corrupted versions of a hub-
like variable, which is not directly observable, and treated as a latent
factor. Formally, we assume there exists a latent random vector
$\bm{Z}\in\mathbb{R}^{K}$, with mean zero and covariance matrix
$\operatorname{\rm Cov}(\bm{Z})=\mathbf{C}^{*}$, such that

$\bm{X}=\mathbf{A}\bm{Z}+\bm{E},$ (1.1)

for a zero mean error vector $\bm{E}$ with independent entries. The entries of
the $d\times K$ matrix $\mathbf{A}$ are given by $A_{jk}=\mathbb{I}\\{j\in
G^{*}_{k}\\}$. A cluster of variables consist in those components of
$\mathbf{X}$ with indices in the same $G^{*}_{k}$. We denote
$\operatorname{\rm Cov}(\bm{E})=\bm{\Gamma}^{*}$, a diagonal matrix with
entries $\Gamma^{*}_{jj}=\gamma_{j}^{*}$ for any $1\leq j\leq d$. We also
assume that the mean-zero noise $\bm{E}$ is independent of $\bm{Z}$. The
clusters are uniquely defined by model (1.1) provided that the smallest
cluster contains at least two variables and that $\mathbf{C}^{*}$ is strictly
positive definite, as shown in Bunea et al. (2018), and this result holds
irrespective of distributional assumptions. To keep the presentation focused,
in this work we assume that $\bm{Z}\sim N(0,\mathbf{C}^{*})$ and $\bm{E}\sim
N(0,\bm{\Gamma}^{*})$, which implies $\bm{X}\sim N(0,\bm{\Sigma}^{*})$ with
$\bm{\Sigma}^{*}=\mathbf{A}\mathbf{C}^{*}\mathbf{A}^{T}+\bm{\Gamma}^{*}$. In
this context we consider two related, but different, graphical models:

  * (i)

The latent variable graph, associated with the sparsity pattern of the
precision matrix

$\bm{\Theta}^{*}:=\mathbf{C}^{*-1}$ (1.2)

of the Gaussian vector $\bm{Z}\in\mathbb{R}^{K}$. The latent variable graph
encodes conditional independencies (CI’s) among the unobservable latent
variables $\bm{Z}$.

  * (ii)

The cluster-average graph, associated with the sparsity pattern of the
precision matrix

$\bm{\Omega}^{*}:=\bm{S}^{*-1},$ (1.3)

where $\bm{S}^{*}$ is the covariance matrix of
$\bar{\bm{X}}\in\mathbb{R}^{K}$, and
$\bar{\bm{X}}=:(\bar{X}_{1},\ldots,\bar{X}_{K})$ is the within cluster average
given by $\bar{X}_{k}=:\frac{1}{|G^{*}_{k}|}\sum_{i\in G^{*}_{k}}X_{i}$. The
cluster-average graph encodes CI’s among averages of observable random
variables. In particular, we have

$\bm{S}^{*}=\mathbf{C}^{*}+\bar{\bm{\Gamma}}^{*},$

where
$\bar{\bm{\Gamma}}^{*}=\textrm{diag}(\bar{\gamma}^{*}_{1},...,\bar{\gamma}^{*}_{K})$
with $\bar{\gamma}^{*}_{k}=\frac{1}{|G_{k}^{*}|^{2}}\sum_{j\in
G_{k}^{*}}\gamma_{j}^{*}$.

Although both graphs correspond to vectors of dimension $K$, possibly much
lower than $d$, they are in general different, as the sparsity patterns of
$\bm{\Theta}^{*}$ and $\bm{\Omega}^{*}$ will typically differ, and have
different interpretations. It would be therefore misleading to use one as a
proxy for the other when drawing conclusions in practice. For instance, in the
neuroscience example, if we interpret each latent variable as the function of
a ROI, the latent variable graph encodes the CI relationships between
functions, which may be one question of scientific interest. The relationship
between $\bm{\Theta}^{*}$ and $\bm{\Omega}^{*}$ shows that this question will
not be typically answered by studying the cluster-average graph, although that
may be tempting to do. On the other hand, the cluster-average graph may be of
independent interest, as it encodes CI’s relations among the average signals
within each ROI.

Since the two cluster-based graphical models introduced above can both be of
interest in a large array of applications, we provide inferential tools for
each of them in this work. We assume that we observe $n$ i.i.d. copies
$\bm{X}_{1},\ldots,\bm{X}_{n}$ of $\bm{X}$. The focus of our work is on post-
clustering and post-regularization inference for these two sparse precision
matrices, with the ultimate goal of sparsity pattern recovery under FDR
control.

Inference for the entries of a Gaussian precision matrix has received a large
amount of attention in the past few years, most notably post-regularization
inference, for instance Ren et al. (2013); Zhang and Zhang (2014); Jankova and
van de Geer (2014); Gu et al. (2015); Barber and Kolar (2015); Janková and van
de Geer (2017); Javanmard and Montanari (2013); van de Geer et al. (2013);
Ning and Liu (2017); Cai et al. (2017); Neykov et al. (2018). These works
generalize to the high-dimensional setup the classical ideas of one-step
estimation (Bickel, 1975), by first constructing a sparse estimator of the
precision matrix, via regularization, and then building de-sparsified updates
that are asymptotically normal. The effect of the initial regularization step
is controlled in the second step, and inference after regularization becomes
valid.

Whereas we will also consider a similar estimation strategy, we differ from
the existing literature in important ways. In our work, we add another layer
of data-dependent dimension reduction, via clustering, and provide a framework
within which the variability induced by clustering can be controlled. Even
after controlling for the clustering variability, we note that the existing
procedures for estimation and, especially, post-regularization inference in
Gaussian graphical models are not immediately applicable to our problem for
the following reasons: (1) They are developed for variables that can be
observed directly. From this perspective, they could, in principle, be applied
to the cluster-average graph, but are not directly extendable to the latent
graph; (2) To the best of our knowledge, all existing methods for precision
matrix inference require the largest eigenvalue of the corresponding
covariance matrix to be upper bounded by a constant. Such an assumption
implies, in turn, that the Euclidean norm of each row of the covariance matrix
is bounded, which reduces significantly the parameter space for which
inference is valid. The assumption holds, for instance, when the number of
variables is bounded, or when the entries of each row are appropriately small.

To overcome these limitations, we take a different approach in this work, that
allows us to lift unpleasant technical conditions associated with other
procedures, while maintaining the validity of inference for both the latent
and the average graph. We summarize our contributions below.

1. Post clustering inference does not impose additional restrictions on the parameter space. We discuss, in Section 2.2, clustering methods tailored to model (1.1), where the number of clusters $K$ is unknown and is allowed to grow with $n$. Using the results of Bunea et al. (2018), these methods yield a partition $\widehat{G}=G^{*}$, with high probability, provided that $\lambda_{\text{min}}(\mathbf{C}^{*})>c$, for a small positive quantity $c$ made precise in Section 2.2. A lower bound on the smallest eigenvalue of the covariance matrix is the minimal condition under which inference in any general graphical model can be performed. Therefore, consistent clustering via model (1.1) does not require a further reduction of the parameter space for which the more standard post-regularized inference can be developed. Moreover, as Section 4 shows, asymptotic inference based on the estimated clusters reduces to asymptotic inference relative to the true clusters, $G^{*}$, without any need for data splitting. This fact is in contrast with the phenomenon encountered in post-model selection inference, for instance in variable selection in linear regression (Lockhart et al., 2014; Lee et al., 2013; Taylor et al., 2014). In that case, reducing inference to the consistently selected set of variables can only be justified over a reduced part of the parameter space (Bunea, 2004).

2. Methods for estimation tailored to high dimensional inference in cluster-based graphical models. We develop a new estimation strategy tailored to our final goal, that of constructing approximately Gaussian estimators for the entries of the precision matrices $\bm{\Theta}^{*}$ and $\bm{\Omega}^{*}$ given in (1.2) and (1.3) above. Although we work under the assumption that the data is Gaussian, likelihood based estimators may be unsatisfactory, as their analysis may require stringent assumptions, as explained above, or may become analytically intractable, as argued in Section 3.3, for the latent graph. We do propose a method that mimics very closely the principles underlying the construction of an efficient score function for estimation in the presence of high dimensional nuisance parameters, see for instance Bickel et al. (1993), but we do not base it on the corresponding likelihood-derived quantities. We explain the underlying principles in Section 3.1.

3\.  Berry-Esseen-type bounds for Gaussian approximations and FDR control. Our
goal is to estimate the cluster-average and latent graphs, respectively, and
to provide guarantees on the recovery error. Existing literature focuses on
constructing approximate confidence intervals for one or a finite number of
entries of a CI graph, or known linear functionals of such entries (Ren et
al., 2013; Zhang and Zhang, 2014; Jankova and van de Geer, 2014; Gu et al.,
2015; Barber and Kolar, 2015; Janková and van de Geer, 2017; Javanmard and
Montanari, 2013; van de Geer et al., 2013; Ning and Liu, 2017; Cai et al.,
2017; Neykov et al., 2018). In all these cases, deriving the asymptotic
limiting distribution of appropriate test statistics suffices. Our focus here
is on the estimation of the sparsity pattern, which can be equivalently viewed
as a multiple-testing problem. It is well known that the exact sparsity
pattern can be recovered, with high probability, only if the entries of each
precision matrix are above the minimax optimal noise level $O(\sqrt{\log
d/n})$ (Ravikumar et al., 2011a; Meinshausen and Bühlmann, 2006). Since our
aim is inference on the sparsity pattern without further restrictions on the
parameter space, the next best type of error that we can control is the False
Discovery Rate (FDR) (Benjamini and Hochberg, 1995). For this, we first derive
Gaussian approximations for the distribution of our estimates. We establish
Berry-Esseen type bounds on the difference between the cumulative distribution
function of our estimators and that of a standard Gaussian random variable
that are valid for each $K,d$ and $n$, and are presented in Theorems 4.3 and
4.4, respectively. In Section 4.3 we use these results for pattern recovery
under FDR control, and explain the effect of the asymptotic approximations on
this quantity.

The paper is organized as follows. Section 2 below contains a brief summary of
existing results on model-assisted clustering, via model (1.1). Section 3
describes the estimation procedures for the latent variable graph and the
cluster-average graph, respectively. In Section 4 we establish Berry-Esseen
type central limit theorems for the estimators derived in Section 3, and
provide bounds on the FDR associated with each graphical model under study,
respectively. Section 5 gives numerical results using both simulated and real
data sets.

##  2 Background

###  2.1 Notation

The following notation is adopted throughout this paper. Let $d$ denote the
ambient dimension, $n$ the sample size, $K$ the number of clusters and $m$ the
minimum cluster size. The matrix $\mathbf{C}^{*}$ denotes the population
covariance of the latent vector $\bm{Z}$. Likewise, the matrices
$\bm{\Gamma}^{*}$, $\bm{\Sigma}^{*}$, $\bm{\Theta}^{*}$, $\bm{S}^{*}$ and
$\bm{\Omega}^{*}$ denote population-level quantities.

For $\mathbf{v}=(v_{1},...,v_{d})^{T}\in\mathbb{R}^{d}$, and $1\leq
q\leq\infty$, we define
$\|{\mathbf{v}}\|_{q}=(\sum_{i=1}^{d}|v_{i}|^{q})^{1/q}$,
$\|{\mathbf{v}}\|_{0}=|\textrm{supp}(\mathbf{v})|$, where
$\textrm{supp}(\mathbf{v})=\\{j:v_{j}\neq 0\\}$ and $|A|$ is the cardinality
of a set $A$. Denote $\|{\mathbf{v}}\|_{\infty}=\max_{1\leq i\leq d}|v_{i}|$
and $\mathbf{v}^{\otimes 2}=\mathbf{v}\mathbf{v}^{T}$. Assume that
$\mathbf{v}$ can be partitioned as
$\mathbf{v}=(\mathbf{v}_{1},\mathbf{v}_{2})$. Let $\nabla f(\mathbf{v})$
denote the gradient of the function $f(\mathbf{v})$, and
$\nabla_{1}f(\mathbf{v})=\partial f(\mathbf{v})/\partial\mathbf{v}_{1}$.
Similarly, let $\nabla^{2}f(\mathbf{v})$ denote the Laplacian of the function
$f(\mathbf{v})$ and
$\nabla^{2}_{12}f(\mathbf{v})=\partial^{2}f(\mathbf{v})/\partial\mathbf{v}_{1}\partial\mathbf{v}_{2}$.
For a $d\times d$ matrix $\mathbf{M}=[M_{jk}]$, let
$\|{\mathbf{M}}\|_{\max}=\max_{jk}|M_{jk}|$,
$\|{\mathbf{M}}\|_{1}=\sum_{jk}|M_{jk}|$, and
$\|{\mathbf{M}}\|_{\infty}=\max_{k}\sum_{j}|M_{jk}|$. If the matrix
$\mathbf{M}$ is symmetric, then $\lambda_{\min}(\mathbf{M})$ and
$\lambda_{\max}(\mathbf{M})$ are the minimal and maximal eigenvalues of
$\mathbf{M}$. Let $[d]=\\{1,2,....,d\\}$. For any $j\in[d]$, we denote the
$j$th row and $j$th column of $\mathbf{M}$ as $\mathbf{M}_{j\cdot}$ and
$\mathbf{M}_{\cdot j}$, respectively. Similarly, let $\mathbf{M}_{-j,-k}$ be
the sub-matrix of $\mathbf{M}$ with the $j^{th}$ row and $k^{th}$ column
removed. Define $\mathbf{M}^{\otimes 2}=\mathbf{M}\otimes\mathbf{M}$. The
notation ${\mathcal{S}}^{d\times d}$ refers to the set of all real, symmetric
$d\times d$ matrices. Likewise, ${\mathcal{S}}^{d\times
d}_{+}\subset{\mathcal{S}}^{d\times d}$ is the positive semi-definite cone. We
use $\otimes$ and $\circ$ to denote the Kronecker and Hadamard product of two
matrices, respectively. Let $\mathbf{e}_{j}$ denote the vector of all zeros
except for a one in the $j^{th}$ position. The vector $\mathbf{1}$ is the
vector of all ones. $a\vee b=\max(a,b)$.

###  2.2 Model Assisted Variable Clustering

In this section we review existing results on variable clustering that will be
used throughout this paper. If model (1.1) is used as a background for
defining clusters of variables, then Bunea et al. (2018) showed that these
clusters are uniquely defined, up to label switching, as soon as
$m=:\min_{1\leq k\leq K}|G_{k}^{*}|\geq 2$ and the components of the latent
vector $\bm{Z}$ are different a.s. :

$\Delta(\mathbf{C}^{*})=:\min_{j<k}\mathbb{E}(Z_{j}-Z_{k})^{2}>0.$

Since

$\Delta(\mathbf{C}^{*})=\min_{j<k}(\mathbf{e}_{j}-\mathbf{e}_{k})^{T}\mathbf{C}^{*}(\mathbf{e}_{j}-\mathbf{e}_{k})\geq
2\lambda_{\min}\left(\mathbf{C}^{*}\right),$

the clusters are uniquely defined as soon as
$\lambda_{\min}\left(\mathbf{C}^{*}\right)>0$, which is the minimal condition
under which one can study properties of the inverse of $\mathbf{C}^{*}$.
Moreover, Bunea et al. (2018) developed two algorithms, PECOK and COD, that
are shown to recover the clusters exactly, from $n$ i.i.d. copies
$\bm{X}_{1},\ldots,\bm{X}_{n}$ of $\bm{X}$, as soon as

$\lambda_{\min}\left(\mathbf{C}^{*}\right)\geq c,$

for a positive quantity $c$ that approaches 0 as $n$ grows. Specifically, for
the COD procedure,

$c=O\left(\|\bm{\Sigma}^{*}\|_{\max}\sqrt{\log(d\vee n)/n}\right).$

On the other hand, for the PECOK procedure

$c=O\left(\|\bm{\Gamma}^{*}\|_{\max}\sqrt{K\log(d\vee n)/mn}\right),$

which can be much smaller when one has a few, balanced, clusters.

These values of $c$ are further shown to be minimax or near-minimax optimal
for cluster recovery. We refer to Theorems 3 and 4 in Bunea et al. (2018) for
the precise expressions and details. Under these minimal conditions on
$\lambda_{\min}\left(\mathbf{C}^{*}\right)$, the exact recovery of the
clusters holds with probability larger than $1-1/(d\vee n)$. This will allow
us to show, in Section 4, that inference in cluster-based graphical models is
not hampered by the clustering step.

For completeness, we outline the PECOK algorithm below, which consists in a
convex relaxation of the $K$-means algorithm, further tailored to estimation
of clusters $G^{*}=\\{G^{*}_{1},\ldots,G^{*}_{K}\\}$ defined via the
interpretable model (1.1). The PECOK algorithm consists in the following three
steps:

  1. 1.

Compute an estimator $\widetilde{\bm{\Gamma}}$ of the matrix
$\bm{\Gamma}^{*}$.

  2. 2.

Solve the semi-definite program (SDP)

$\widehat{\mathbf{B}}=\mathop{\mathrm{argmax}}_{\mathbf{B}\in\mathcal{D}}\langle\widehat{\bm{\Sigma}}-\widetilde{\bm{\Gamma}},\mathbf{B}\rangle,$
(2.1)

where $\widehat{\bm{\Sigma}}$ is the sample covariance matrix and

$\mathcal{D}:=\left\\{\mathbf{B}\in R^{d\times d}:\begin{array}[]{l}\bullet\
\bm{B}\succcurlyeq 0\ \ \text{(symmetric and positive semidefinite)}\\\
\bullet\ \sum_{a}B_{ab}=1,\ \forall b\\\ \bullet\ B_{ab}\geq 0,\ \forall
a,b\\\ \bullet\ \mathop{\text{tr}}\kern
0.86108pt(\mathbf{B})=K\end{array}\right\\}.$ (2.2)

  3. 3.

Compute $\widehat{G}$ by applying a clustering algorithm on the rows (or
equivalently columns) of $\widehat{\mathbf{B}}$.

The construction of an accurate estimator $\widetilde{\bm{\Gamma}}$ of
$\bm{\Gamma}^{*}$, before the cluster structure is known, is a crucial step
for guaranteeing the statistical optimality of the PECOK estimator. Its
construction is given in Bunea et al. (2018), and included in Appendix F, for
the convenience of the reader.

We will employ an efficient algorithm for solving (2.1). Standard black-box
SDP solvers, for a fixed precision, exhibit $\mathcal{O}(d^{7})$ running time
on (2.1), which is prohibitively expensive. Eisenach and Liu (2017) recently
introduced the FORCE algorithm, which requires worst case
$\mathcal{O}(d^{6}K^{-2})$ time to solve the SDP, and in practice often
performs the clustering rapidly.

The key idea behind the FORCE algorithm is that an optimal solution to (2.1)
can be attained by first transforming (2.1) into an eigenvalue problem, and
then using a first-order method. Iterations of the first-order method are
interleaved with a dual step to round the current iterate to an integer
solution of the clustering problem, and then searches for an optimality
certificate. By using knowledge of both the primal and the dual SDPs, FORCE is
able to find the solution much faster than a standard SDP solver. We refer to
Eisenach and Liu (2017) for the detailed algorithm.

##  3 Estimation of Cluster-based Graphical Models

In this section, we propose a unified estimation approach, that utilizes
similar loss functions for estimation and inference in the cluster-average and
the latent variable graphs. We first describe our general principle, and then
apply it to the two graphical models, respectively.

###  3.1 One-step Estimators for High-Dimensional Inference

Assume that we observe $n$ i.i.d. realizations $\bm{X}_{1},...,\bm{X}_{n}$ of
$\bm{X}\in\mathbb{R}^{d}$. Let $Q(\bm{\beta},\bm{X})$ denote a known function
of $\bm{\beta}$ and $\bm{X}$, where $\bm{\beta}$ is a $q$-dimensional unknown
parameter that parametrizes the distribution of $\bm{X}$. We define the target
parameter $\bm{\beta}^{*}$ as

$\bm{\beta}^{*}=\mathop{\mathrm{argmin}}\mathbb{E}(Q(\bm{\beta},\bm{X})).$

Let us partition $\bm{\beta}$ as $\bm{\beta}=(\theta,\bm{\gamma})$, where
$\theta\in\mathbb{R}$ is a univariate parameter of interest, and
$\bm{\gamma}\in\mathbb{R}^{q-1}$ is a nuisance parameter. Our goal is to
construct a $n^{1/2}$-consistent and asymptotically normal estimator for
$\theta$ in high-dimensional models with $q=\textrm{dim}(\bm{\beta})\gg n$. In
this case, the dimension of the nuisance parameter $\gamma$ is large, which
makes the inference on $\theta$ challenging. We start from the empirical risk
function over $n$ observations defined as

$Q_{n}(\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}Q(\bm{\beta},\bm{X}_{i}).$ (3.1)

One standard instance of $Q_{n}$ is the negative log-likelihood function of
the data. In this work, we will conduct inference based on an alternative loss
function, as the analysis of the log-likelihood may require unpleasant
technical conditions that we would like to avoid, as discussed in Sections
4.1. However, we mimic as much as possible the likelihood principles, in order
to aid the understanding of the construction below. For these reasons we will
refer to $Q_{n}(\bm{\beta})$ as the negative pseudo-likelihood function. For
now, we will assume that $Q_{n}(\bm{\beta})$ is given, and a detailed
discussion of its respective choice for inference in the latent variable graph
and the cluster-average graph will be given in the following two subsections.
The pseudo-information matrix for one observation is defined as
$\mathbf{I}=\mathbb{E}(\nabla^{2}Q(\bm{\beta}^{*},\bm{X}_{i}))$, which can be
further partitioned as

$\mathbf{I}=\begin{bmatrix}\mathbf{I}_{11}&\mathbf{I}_{12}\\\
\mathbf{I}_{21}&\mathbf{I}_{22}\\\ \end{bmatrix},$ (3.2)

relative to the partition of $\bm{\beta}=(\theta,\bm{\gamma})$. When $Q_{n}$
is the negative log-likelihood function, and the dimension of the parameter is
independent of $n$, then $h(\bm{\theta};\bm{\gamma})$ given by (3.3) is called
the efficient score function for $\theta$, and classical theory shows that it
admits solutions that are consistent, asymptotically normal and attain the
information bound given by the reciprocal of (3.4) (Van der Vaart, 1998;
Bickel et al., 1993).

With these goals in mind, we similarly define the corresponding pseudo-score
function for estimating $\theta$ in the presence of the nuisance parameter
$\bm{\gamma}$ as

$h(\theta;\bm{\gamma})=\nabla_{1}Q_{n}(\bm{\beta})-\mathbf{I}_{12}\mathbf{I}_{22}^{-1}\nabla_{2}Q_{n}(\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}\Big{(}\nabla_{1}Q(\bm{\beta},\bm{X}_{i})-\mathbf{I}_{12}\mathbf{I}_{22}^{-1}\nabla_{2}Q(\bm{\beta},\bm{X}_{i})\Big{)},$
(3.3)

and define the pseudo information of $\theta$, in the presence of the nuisance
parameter $\bm{\gamma}$, as

$\mathbf{I}_{1\mid
2}=\mathbf{I}_{11}-\mathbf{I}_{12}\mathbf{I}_{22}^{-1}\mathbf{I}_{21}.$ (3.4)

When the dimension of $\bm{\gamma}$ is fixed, one can easily estimate
$\mathbf{I}_{12}$ and $\mathbf{I}_{22}$ in (3.3) by their sample versions
$\widehat{\mathbf{I}}_{12}$ and $\widehat{\mathbf{I}}_{22}$. However, such
simple procedure fails when the dimension of $\bm{\gamma}$ is greater than the
sample size, as $\widehat{\mathbf{I}}_{22}$ is rank deficient. To overcome
this difficulty, rather than estimating $\mathbf{I}_{12}$ and
$\mathbf{I}_{22}^{-1}$ separately, we directly estimate

$\mathbf{w}^{T}=\mathbf{I}_{12}\mathbf{I}_{22}^{-1}$ (3.5)

by

$\widehat{\mathbf{w}}=\mathop{\mathrm{argmin}}\|{\mathbf{w}}\|_{1},~{}~{}~{}~{}\textrm{s.t.}~{}~{}~{}~{}\|{\nabla^{2}_{12}Q_{n}(\widehat{\bm{\beta}})-\mathbf{w}^{T}\nabla^{2}_{22}Q_{n}(\widehat{\bm{\beta}})}\|_{\infty}\leq\lambda^{\prime},$
(3.6)

where $\lambda^{\prime}$ is a non-negative tuning parameter, and
$\widehat{\bm{\beta}}=(\widehat{\theta},\widehat{\bm{\gamma}})$ is an initial
estimator, which is usually defined case by case, for a given model. Then, we
can plug $\widehat{\mathbf{w}}$ and $\widehat{\bm{\gamma}}$ into the pseudo-
score function, which gives

$\widehat{h}(\theta,\widehat{\bm{\gamma}})=\nabla_{1}Q_{n}(\theta,\widehat{\bm{\gamma}})-\widehat{\mathbf{w}}^{T}\nabla_{2}Q_{n}(\theta,\widehat{\bm{\gamma}}).$
(3.7)

Following the Z-estimation principle (Van der Vaart, 1998; Bickel et al.,
1993), one could define the final estimator of $\theta$ as the solution of the
pseudo-score function $\widehat{h}(\theta,\widehat{\bm{\gamma}})$. However, in
many examples, the pseudo-score function
$\widehat{h}(\theta,\widehat{\bm{\gamma}})$ may have multiple solutions and it
becomes unclear which root serves as a consistent estimator; see Small et al.
(2000) for further discussion in the general estimating function context. To
bypass this issue, we consider the following simple one-step estimation
approach. Given the initial estimator $\widehat{\theta}$ from the partition of
$\widehat{\bm{\beta}}$, we perform a Newton-Raphson update based on the
pseudo-score function $\widehat{h}(\theta,\widehat{\bm{\gamma}})$, to obtain
$\widetilde{\theta}$, which is classically referred to as a one-step estimator
by Bickel (1975). Specifically, we construct

$\widetilde{\theta}=\widehat{\theta}-\widehat{\mathbf{I}}_{1|2}^{-1}\widehat{h}(\widehat{\theta},\widehat{\bm{\gamma}}),$
(3.8)

where $\widehat{\mathbf{I}}_{1|2}$ is an estimator of the partial information
matrix $\mathbf{I}_{1\mid 2}$. In Sections 3.2 and 3.3 below we show that,
under appropriate conditions, the one-step estimator $\widetilde{\theta}$
constructed relative to the empirical risk functions $Q_{n}$ defined in (3.10)
and (3.25), respectively, satisfies

$n^{1/2}(\widetilde{\theta}-\theta^{*})=-\mathbf{I}^{-1}_{1\mid
2}n^{1/2}h(\bm{\beta}^{*})+o_{p}(1).$ (3.9)

By applying the central limit theory to $h(\bm{\beta}^{*})$, we can establish
the asymptotic normality of $\widetilde{\theta}$ in Theorems 4.3 and 4.4. When
$Q_{n}(\bm{\beta})$ is the negative log-likelihood of the data, this approach
has been successfully used in Ning and Liu (2017) and, moreover, the estimator
$\widetilde{\theta}$ is asymptotically equivalent to the de-biased estimator
in Zhang and Zhang (2014); van de Geer et al. (2013). As will be explained in
the following subsections, the analysis based on the log-likelihood becomes
intractable for the latent graphical model and requires stringent technical
conditions for the cluster-average graphical model. To overcome this
difficulty, we employ the pseudo score functions relative to the empirical
risk functions $Q_{n}(\bm{\beta})$ defined in (3.10) and (3.25). The resulting
one-step estimator still attains the information bound established in the
literature, and more importantly requires weaker technical assumptions than
the existing methods. Moreover, in addition to (3.9), we also derive
explicitly the speed at which the normal approximation is attained.

###  3.2 Estimation of the Cluster-Average Graph

Recall that we assume $\bm{Z}\sim N(0,\mathbf{C}^{*})$ and $\bm{E}\sim
N(0,\bm{\Gamma}^{*})$, which implies $\bm{X}\sim N(0,\bm{\Sigma}^{*})$ with
$\bm{\Sigma}^{*}=\mathbf{A}\mathbf{C}^{*}\mathbf{A}^{T}+\bm{\Gamma}^{*}$. The
within-cluster average
$\bar{\bm{X}}=:(\bar{X}_{1},\ldots,\bar{X}_{K})\in\mathbb{R}^{K}$ is given by
$\bar{X}_{k}=:\frac{1}{|G^{*}_{k}|}\sum_{i\in G^{*}_{k}}X_{i}$, corresponding
to the population level clusters. Because $\bm{X}\sim N(0,\bm{\Sigma}^{*})$,
we can verify that $\bar{\bm{X}}\sim N(0,\bm{S}^{*})$, where

$\bm{S}^{*}=\mathbf{C}^{*}+\bar{\bm{\Gamma}}^{*},$

and
$\bar{\bm{\Gamma}}^{*}=\textrm{diag}(\bar{\gamma}^{*}_{1},...,\bar{\gamma}^{*}_{K})$
with $\bar{\gamma}^{*}_{k}=\frac{1}{|G_{k}^{*}|^{2}}\sum_{j\in
G_{k}^{*}}\gamma_{j}^{*}$. Recall that the precision matrix of $\bar{\bm{X}}$
is

$\bm{\Omega}^{*}={\bm{S}^{*}}^{-1}=(\mathbf{C}^{*}+\bar{\bm{\Gamma}}^{*})^{-1}.$

In this section we give the construction of the estimators of the cluster-
average graph corresponding to $\bar{\bm{X}}$. Specifically, we use the
generic strategy outlined in the previous section in order to construct
$n^{1/2}$-consistent and asymptotically normal estimators for each component
$\bm{\Omega}^{*}_{t,k}$ of the precision matrix $\bm{\Omega}^{*}$, for $1\leq
t<k\leq K$. For the estimation of each entry, the remaining $K(K+1)/2-1$
parameters in $\bm{\Omega}^{*}$ are treated as nuisance parameters.

Since we observe $n$ i.i.d. samples of $\bm{X}\in\mathbb{R}^{p}$, if the
clusters and their number were known, then we implicitly observe $n$ i.i.d.
samples of $\bar{\bm{X}}\in\mathbb{R}^{K}$. To explain our method, we first
assume that clustering is given, and then explain how to lift this assumption.

Following our general principle, we would naturally tend to choose the
negative log-likelihood function of the cluster-averages
$(\bar{\bm{X}}_{1},...,\bar{\bm{X}}_{n})$ as the empirical risk function
$Q_{n}(\bm{\beta})$ in (3.1). Along this line, Jankova and van de Geer (2014)
proposed the de-biased estimator for Gaussian graphical models. However, the
inference requires the irrepresentable condition (Ravikumar et al., 2011b) on
$\bm{S}^{*}$, which can be restrictive. The alternative methods proposed by
Ren et al. (2013); Janková and van de Geer (2017) imposed the condition that
the largest eigenvalue of $\bm{S}^{*}$ is bounded. These technical conditions
on $\bm{S}^{*}$ are difficult to justify and can be avoided by using our
approach. We propose to estimate each sparse row of $\bm{\Omega}^{*}$ as
explained below.

Let $\bar{\bm{S}}=n^{-1}\sum_{i=1}^{n}\bar{\bm{X}}_{i}\bar{\bm{X}}_{i}^{T}$
denote the sample covariance matrix of $\bar{\bm{X}}_{i}$. When $K$ is small,
the maximum likelihood estimator of $\bm{\Omega}^{*}$ is $\bar{\bm{S}}^{-1}$,
which can be viewed as the solution of the following equation
$\bar{\bm{S}}\bm{\Omega}-\mathbf{I}_{K}=0$. Thus, in the low dimensional
setting, this equation defines the maximum likelihood estimator. Since we are
only interested in $\bm{\Omega}^{*}_{t,k}$, we can extract the $k$th column
from the left hand side of the above equation, and use it as the pseudo-score
function $\bm{U}_{n}(\bm{\Omega}_{\cdot k})=\bar{\bm{S}}\bm{\Omega}_{\cdot
k}-\mathbf{e}_{k}$. To apply the inference strategy in Section 3.1, we need to
construct a valid empirical risk function $Q_{n}(\bm{\Omega}_{\cdot k})$ such
that $\nabla Q_{n}(\bm{\Omega}_{\cdot k})=\bm{U}_{n}(\bm{\Omega}_{\cdot k})$.

Simple algebra shows that a possible choice is

$Q_{n}(\bm{\Omega}_{\cdot k})=\frac{1}{2}\bm{\Omega}_{\cdot
k}^{T}\bar{\bm{S}}\bm{\Omega}_{\cdot k}-\mathbf{e}_{k}^{T}\bm{\Omega}_{\cdot
k}=\frac{1}{n}\sum_{i=1}^{n}(\frac{1}{2}\bm{\Omega}_{\cdot
k}^{T}\bar{\bm{X}}_{i}\bar{\bm{X}}_{i}^{T}\bm{\Omega}_{\cdot
k}-\mathbf{e}_{k}^{T}\bm{\Omega}_{\cdot k}),$ (3.10)

which we view in the sequel as the empirical risk corresponding to the
population level risk

$\mathbb{E}Q(\bm{\Omega}_{\cdot k},\bar{\bm{X}})=\frac{1}{2}\bm{\Omega}_{\cdot
k}^{T}\bm{S}^{*}\bm{\Omega}_{\cdot k}-\mathbf{e}_{k}^{T}\bm{\Omega}_{\cdot
k},$ (3.11)

based on the loss function

$Q(\bm{\Omega}_{\cdot k},\bar{\bm{X}})=:\frac{1}{2}\bm{\Omega}_{\cdot
k}^{T}\bar{\bm{X}}\bar{\bm{X}}^{T}\bm{\Omega}_{\cdot
k}-\mathbf{e}_{k}^{T}\bm{\Omega}_{\cdot k}.$ (3.12)

Since

$\nabla\mathbb{E}Q(\bm{\Omega}^{*}_{\cdot
k},\bar{\bm{X}})=\bm{S}^{*}\bm{\Omega}^{*}_{\cdot k}-\mathbf{e}_{k}=0$ (3.13)

and

$\nabla^{2}\mathbb{E}Q(\bm{\Omega}_{\cdot k}^{*},\bar{\bm{X}})=\bm{S}^{*},$
(3.14)

then the theoretical risk $\mathbb{E}Q(\bm{\Omega}^{*}_{\cdot
k},\bar{\bm{X}})$ has the rows $\bm{\Omega}^{*}_{\cdot k}$ of the target
precision matrix $\bm{\Omega}^{*}$ as the unique minimizers, as desired,
provided that $\bm{S}^{*}$ is positive definite, an assumption we make in
Section 4.1.

We note that the choice of the empirical risk $Q_{n}(\cdot)$ and that of the
corresponding pseudo-score $\bm{U}_{n}(\cdot)$ is not unique. We chose the
particular form (3.10) because it is quadratic in $\bm{\Omega}_{\cdot k}$,
which greatly simplifies the theoretical analysis and leads to weaker
technical assumptions. Moreover, the property (3.13) is the same as that of
the score function corresponding to the negative log-likelihood function,
supporting our terminology.

We use the general strategy presented in Section 3.1 to construct estimators
that employ the empirical risk $Q_{n}(\cdot)$ defined by (3.10) above. We
first recall that $Q_{n}(\cdot)$ depends on the unknown cluster structure
$G^{*}$ via $\bar{\bm{X}}_{i}$. We note that in general the estimated group
$\widehat{G}_{k}$ may differ from $G_{k}^{*}$ by a label permutation. For
notational simplicity, we ignore this label permutation issue and treat
$\widehat{G}_{k}$ as an estimate of $G_{k}^{*}$ (rather than $G_{j}^{*}$ for
some $j\neq k$). To define the estimator of $\bm{\Omega}^{*}_{t,k}$, we first
replace $\bar{\bm{X}}_{i}$ by $\widehat{\bm{X}}_{i}$ and denote
$\widehat{\bm{S}}=n^{-1}\sum_{i=1}^{n}\widehat{\bm{X}}_{i}\widehat{\bm{X}}_{i}^{T}$,
where
$\widehat{X}_{ik}=\frac{1}{|\widehat{G}_{k}|}\sum_{j\in\widehat{G}_{k}}X_{ij}$.

Let $(t,k)$ be arbitrary, fixed. Replacing $\bar{\bm{S}}$ by
$\widehat{\bm{S}}$ in $Q_{n}(\cdot)$, we follow Section 3.1 to define the
pseudo-score function

$h(\bm{\Omega}_{\cdot
k})=\mathbf{v}_{t}^{*T}(\widehat{\bm{S}}\bm{\Omega}_{\cdot
k}-\mathbf{e}_{k}),$ (3.15)

where ${\mathbf{v}}^{*}_{t}$ is a $K$-dimensional vector with
$(\mathbf{v}_{t}^{*})_{t}=1$ and
$(\mathbf{v}_{t}^{*})_{-t}=-\mathbf{w}^{*}_{t}$ with
$\mathbf{w}^{*}_{t}=(\bm{S}^{*}_{-t,-t})^{-1}\bm{S}^{*}_{-t,t}$ consistent
with the definition in (3.5) above. To make inference based on
$h(\bm{\Omega}_{\cdot k})$, we need to further estimate $\mathbf{w}^{*}_{t}$
and $\bm{\Omega}^{*}_{\cdot k}$. Following (3.6), an estimate of
$\mathbf{w}^{*}_{t}$ is given by

$\widehat{\mathbf{w}}_{t}=\mathop{\mathrm{argmin}}\|\mathbf{w}\|_{1},~{}~{}\textrm{s.t}~{}~{}\|\widehat{\bm{S}}_{t,-t}-\mathbf{w}^{T}\widehat{\bm{S}}_{-t,-t}\|_{\infty}\leq\lambda^{\prime},$
(3.16)

where $\lambda^{\prime}$ is a tuning parameter. Then we can define
$\widehat{\mathbf{v}}_{t}$ accordingly, and

$\widehat{h}(\bm{\Omega}_{\cdot
k})=\widehat{\mathbf{v}}_{t}^{T}(\widehat{\bm{S}}\bm{\Omega}_{\cdot
k}-\mathbf{e}_{k}).$ (3.17)

Recall that the construction of the one-step estimator (3.8) requires an
initial estimator of $\bm{\Omega}^{*}_{\cdot k}$. To be concrete, we consider
the following initial estimator of $\bm{\Omega}^{*}_{\cdot k}$,

$\widehat{\bm{\Omega}}_{\cdot
k}=\mathop{\mathrm{argmin}}\|\bm{\beta}\|_{1},~{}~{}\textrm{s.t}~{}~{}\|\widehat{\bm{S}}\bm{\beta}-\mathbf{e}_{k}\|_{\max}\leq\lambda,$
(3.18)

where $\lambda$ is a tuning parameter. This estimator has the same form as the
CLIME estimator for the $k$-th column of $\bm{\Omega}$ (Cai et al., 2011).
However, unlike the CLIME estimator which requires
$\lambda\asymp\|\bm{\Omega}^{*}_{\cdot k}\|_{1}\sqrt{\log K/n}$, in our
Theorem 4.3 we assume $\lambda=C\sqrt{\log(K\vee n)/n}$, where $C$ only
depends on the minimum eigenvalue of $\mathbf{C}^{*}$ and the largest diagonal
entries of $\mathbf{C}^{*}$ and $\bm{\Gamma}^{*}$ which are assumed bounded by
constants in Assumptions 4.1 and 4.2. With this choice of $\lambda$, we show
in Lemma A.8 in Appendix A that

$\|\widehat{\bm{\Omega}}_{\cdot k}-\bm{\Omega}_{\cdot k}^{*}\|_{1}\lesssim
s_{1}\sqrt{\frac{\log(K\vee n)}{n}},$ (3.19)

with high probability, where $s_{1}$ is the sparsity level of
$\bm{\Omega}_{\cdot k}^{*}$. As a comparison, Theorem 6 in Cai et al. (2011)
only implies $\|\widehat{\bm{\Omega}}_{\cdot k}-\bm{\Omega}_{\cdot
k}^{*}\|_{1}\lesssim s_{1}\|\bm{\Omega}_{\cdot
k}^{*}\|_{1}^{2}\sqrt{\frac{\log K}{n}}$. For many sparse matrices, the
$\ell_{1}$ norm of a column, $\|\bm{\Omega}_{\cdot k}^{*}\|_{1}$, can grow to
infinity with $K$ or $s_{1}$, and thus (3.19) gives a faster rate. This is
possible to obtain under the Gaussian assumption on $\bm{X}$, which can be
easily relaxed to sub-Gaussian. In this case, and when
$\lambda_{\min}(C^{*})>c$, Lemma A.7 is instrumental in showing that the extra
$\|\bm{\Omega}_{\cdot k}^{*}\|_{1}^{2}$ factor in the rate of the original
CLIME estimator can be avoided, whereas if only the marginal components of
$\bm{X}$ are assumed to be sub-Gaussian, as in Cai et al. (2011), it may be
unavoidable, without further conditions on $\bm{\Omega}^{*}$.

Based on the block matrix inverse formula, we can show that the partial pseudo
information matrix reduces to $\mathbf{I}_{1|2}=1/\Omega_{t,t}^{*}$. Finally,
the one-step estimator is defined as

$\widetilde{\Omega}_{t,k}=\widehat{\Omega}_{t,k}-\widehat{h}(\widehat{\bm{\Omega}}_{\cdot
k})\widehat{\Omega}_{t,t},$ (3.20)

in accordance with (3.8). In Section 4, we will show that under mild
regularity conditions
$n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega_{t,k}^{*})\leadsto N(0,s_{tk}^{2})$,
where $s_{tk}^{2}=\Omega_{t,k}^{*2}+\Omega^{*}_{t,t}\Omega^{*}_{k,k}$. Let
$\widehat{s}_{tk}^{2}=\widehat{\Omega}_{t,k}^{2}+\widehat{\Omega}_{t,t}\widehat{\Omega}_{k,k}$
be a consistent estimator of the asymptotic variance. Then, a
$(1-\alpha)\times 100\%$ confidence interval for $\Omega_{t,k}$ is

$[\widetilde{\Omega}_{t,k}-z_{1-\alpha/2}\widehat{s}_{tk}/n^{1/2},\widetilde{\Omega}_{t,k}+z_{1-\alpha/2}\widehat{s}_{tk}/n^{1/2}],$

where $z_{\alpha}$ is the $\alpha$-quantile of a standard normal distribution.
Equivalently, we can use the scaled test statistics $\widetilde{\Omega}_{t,k}$
to construct a test for $H_{0}:\Omega^{*}_{t,k}=0$ versus
$H_{1}:\Omega^{*}_{t,k}\neq 0$ with $\alpha$ significance level. Namely, the
null hypothesis is rejected if and only if the above $(1-\alpha)\times 100\%$
confidence interval does not contain $0$. We will employ such tests in Section
4.

###  3.3 Latent Variable Graph

Recall that the structure of the latent variable graph is encoded by the
sparsity pattern of $\bm{\Theta}^{*}={\mathbf{C}^{*}}^{-1}$, which is
generally different from the cluster-average group as ${\mathbf{C}^{*}}^{-1}$
and $\bm{\Omega}^{*}=(\mathbf{C}^{*}+\bar{\bm{\Gamma}}^{*})^{-1}$ may have
different sparsity patterns. In this section, we focus on the inference on the
component $\Theta_{t,k}^{*}$, for some $1\leq t<k\leq K$. Similar to the
cluster-average graph, we first discuss the likelihood approach. The negative
log-likelihood corresponding to model (1.1) indexed by the parameter
$(\bm{\Theta},\bm{\Gamma})$ is

$\ell(\bm{\Theta},\bm{\Gamma})=-\textrm{tr}(\widehat{\bm{\Sigma}}(\mathbf{A}\bm{\Theta}^{-1}\mathbf{A}^{T}+\bm{\Gamma})^{-1})+\log\det((\widehat{\bm{\Sigma}}(\mathbf{A}\bm{\Theta}^{-1}\mathbf{A}^{T}+\bm{\Gamma})^{-1}),$

where $\widehat{\bm{\Sigma}}=n^{-1}\sum_{i=1}^{n}\bm{X}_{i}\bm{X}^{T}_{i}$.
After tedious algebra, we show that the Fisher information matrix for
$(\bm{\Theta},\bm{\Gamma})$ is given by

$\displaystyle\mathbf{I}$
$\displaystyle=\begin{bmatrix}(\mathbf{M}^{*}\mathbf{A}^{T}\bm{\Gamma}^{*-1}\bm{\Sigma}^{*}\bm{\Gamma}^{*-1}\mathbf{A}\mathbf{M}^{*})^{\otimes
2}&(\mathbf{M}^{*}\mathbf{A}^{T}\bm{\Gamma}^{*-1}\mathbf{F}^{*T})^{\otimes
2}\mathbf{D}_{d}\\\
\mathbf{D}_{d}^{T}(\mathbf{F}^{*}\bm{\Sigma}^{*-1}\bm{\Gamma}^{*-1}\mathbf{A}\mathbf{M}^{*})^{\otimes
2}&\mathbf{D}_{d}^{T}(\mathbf{F}^{*}\bm{\Sigma}^{*-1}\mathbf{F}^{*T})^{\otimes
2}\mathbf{D}_{d}\end{bmatrix},$ (3.21)

where $\mathbf{D}_{d}=(\mathbf{I}_{d}\otimes
1_{d}^{T})\circ(1_{d}^{T}\otimes\mathbf{I}_{d})$,
$\mathbf{M}^{*}=(\bm{\Theta}^{*}+\mathbf{A}^{T}\bm{\Gamma}^{*-1}\mathbf{A})^{-1}$
and
$\mathbf{F}^{*}=\mathbf{I}_{d}-\mathbf{A}\mathbf{M}^{*}\mathbf{A}^{T}\bm{\Gamma}^{*-1}$.
As seen in Section 3.1, the inference based on the likelihood or equivalently
efficient score function (3.3) requires the estimation of
$\mathbf{I}_{12}\mathbf{I}_{22}^{-1}$ which, given the complicated structure
of the information matrix (3.21), becomes analytically intractable.

A solution to this problem is inference based on an empirical risk function
similar to (3.10), but tailored to the latent variable graph. With a slight
abuse of notation, and reasoning as in (3.13) and (3.14), we notice that, for
each $k$,

$\mathbb{E}Q(\bm{\Theta}_{\cdot k},\bm{X})=\frac{1}{2}\bm{\Theta}_{\cdot
k}^{T}\mathbf{C}^{*}\bm{\Theta}_{\cdot k}-\mathbf{e}_{k}^{T}\bm{\Theta}_{\cdot
k},$ (3.22)

has the target $\bm{\Theta}^{*}_{\cdot k}$ as a unique minimizer, where the
loss function $Q(\bm{\Theta}_{\cdot k},\bm{X})$ is defined as

$Q(\bm{\Theta}_{\cdot k},\bm{X})=\frac{1}{2}\bm{\Theta}_{\cdot
k}^{T}\bar{\mathbf{C}}\bm{\Theta}_{\cdot
k}-\mathbf{e}_{k}^{T}\bm{\Theta}_{\cdot k},$ (3.23)

and the matrix $\bar{\mathbf{C}}:=(\bar{C}_{jk})_{j,k}$ has entries

$\bar{C}_{jk}=\frac{1}{|G^{*}_{j}||G^{*}_{k}|}\sum_{a\in G^{*}_{j},b\in
G^{*}_{k}}(X_{a}X_{b}-\bar{\Gamma}_{ab}),$ (3.24)

and $\bar{\Gamma}_{ab}=0$ if $a\neq b$ and
$\bar{\Gamma}_{aa}=X_{a}X_{a}-\frac{1}{|G^{*}_{k}|-1}\sum_{a\in
G^{*}_{k},a\neq j}X_{a}X_{j}$. Since
$\mathbb{E}(\bar{\mathbf{C}})=\mathbf{C}^{*}$, the risk relative to the loss
function in (3.23) is indeed (3.22), and the empirical risk is

$Q_{n}(\bm{\Theta}_{\cdot
k})=\frac{1}{n}\sum_{i=1}^{n}(\frac{1}{2}\bm{\Theta}_{\cdot
k}^{T}\bar{\mathbf{C}}^{(i)}\bm{\Theta}_{\cdot
k}-\mathbf{e}_{k}^{T}\bm{\Theta}_{\cdot k}),$ (3.25)

where $\bar{\mathbf{C}}^{(i)}$ is obtained by replacing $\bm{X}$ in
$\bar{C}_{jk}$ by $\bm{X}_{i}$. Similar to the cluster-average graph,
$Q_{n}(\bm{\Theta}_{\cdot k})$ also depends on the unknown cluster structure.
We estimate $G^{*}_{k}$ by $\widehat{G}_{k}$, and define
$\widehat{\bm{\Gamma}}=(\widehat{\Gamma}_{ab})$, where
$\widehat{\Gamma}_{ab}=0$ if $a\neq b$ and
$\widehat{\Gamma}_{aa}=\frac{1}{n}\sum_{i=1}^{n}(X_{ia}X_{ia}-\frac{1}{|\widehat{G}_{k}|-1}\sum_{a\in\widehat{G}_{k},a\neq
j}X_{ia}X_{ij})$, and $\widehat{\mathbf{C}}=(\widehat{C}_{jk})$ where
$\widehat{C}_{jk}=\frac{1}{n}\sum_{i=1}^{n}(\frac{1}{|\widehat{G}_{j}||\widehat{G}_{k}|}\sum_{a\in\widehat{G}_{j},b\in\widehat{G}_{k}}(X_{ia}X_{ib}-\widehat{\Gamma}_{ab}))$.
We replace $\bar{\mathbf{C}}$ by $\widehat{\mathbf{C}}$ in (3.25) above and
follow exactly the strategy of Section 3.2, with $\widehat{\bm{S}}$ replaced
by $\widehat{\mathbf{C}}$, to construct the corresponding pseudo-score
function $\widehat{h}(\bm{\Theta}_{\cdot k})$, similarly to (3.17), and the
initial estimator $\widehat{\bm{\Theta}}_{\cdot k}$, similarly to (3.18). We
combine these quantities, following the general strategy (3.8), as above, to
obtain the final one-step estimator of $\Theta_{t,k}^{*}$, defined as

$\widetilde{\Theta}_{t,k}=\widehat{\Theta}_{t,k}-\widehat{h}(\widehat{\bm{\Theta}}_{\cdot
k})\widehat{\Theta}_{t,t},$ (3.26)

after observing that, in this case,
$\mathbf{I}_{1|2}=1/{\bm{\Theta}}^{*}_{t,t}$.

Although the form of this estimator is similar to (3.20), derived for the
cluster-average graph, the study of the asymptotic normality of
$\widetilde{\Theta}_{t,k}$ reveals that its asymptotic variance is much more
involved, and will be discussed in detail in Section 4.2.2.

##  4 Main Theoretical Results

###  4.1 Assumptions

In this section we state the two assumptions under which all our results are
proved.

######  Assumption 4.1.

The covariance matrix $\mathbf{C}^{*}$ of $\bm{Z}$ satisfies:
$c_{1}\leq\lambda_{\min}\left(\mathbf{C}^{*}\right)$ and
$\max_{t}C_{t,t}^{*}\leq c_{2}$, for some absolute constants $c_{1},c_{2}>0$.

######  Assumption 4.2.

The matrix $\bm{\Gamma}^{*}$ satisfies: $\max_{1\leq i\leq
d}\gamma^{*}_{i}\leq c_{3}$ for some absolute constant $c_{3}>0$, where
$\gamma_{i}^{*}$ are the entries of the diagonal matrix $\bm{\Gamma}^{*}$.

Assumptions 4.1 and 4.2 are minimal conditions for inference on precision
matrices. Furthermore, they imply the conditions needed for clustering
consistency derived in Bunea et al. (2018) and discussed in Section 2.2, for
$n$ large enough. The latter only require that
$\lambda_{\min}\left(\mathbf{C}^{*}\right)$ is bounded from below by a
sequence that converges to zero, as soon as $\|\bm{\Sigma}^{*}\|_{\max}$ and
$\|\bm{\Gamma}^{*}\|_{\max}$ are bounded. This is strengthened by our
assumptions. In particular, a constant lower bound on $\mathbf{C}^{*}$ is
standard in any inference on graphical models and needed to show the
asymptotic normality of the estimator introduced above (Ren et al., 2013;
Jankova and van de Geer, 2014; Janková and van de Geer, 2017).

###  4.2 Asymptotic Normality via Berry-Esseen-type Bounds

####  4.2.1 Results for the Cluster-Average Graph

In the section, we show that the estimators $\widetilde{\Omega}_{t,k}$ given
by (3.20) are asymptotically normal, for all $t<k$ . We define the sparsity of
the cluster-average graph as $s_{1}\in\mathbb{N}$ such that

$\max_{1\leq j\leq K}\sum_{k=1}^{K}\mathbb{I}(\Omega^{*}_{j,k}\neq 0)\leq
s_{1}.$

Recall that the estimators (3.16) and (3.18) depend on the tuning parameters
$\lambda$ and $\lambda^{\prime}$. In the following theorem, we choose
$\lambda\asymp\lambda^{\prime}\asymp\sqrt{\frac{\log(K\vee n)}{n}}$. For
notational simplicity, we use $C$ to denote a generic constant, the value of
which may change from line to line.

######  Theorem 4.3.

If Assumptions 4.1 and 4.2 hold, we have

$\max_{1\leq t<k\leq
K}\sup_{x\in\mathbb{R}}\Big{|}\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{\widehat{s}_{tk}}\leq
x\Big{)}-\Phi(x)\Big{|}\leq\frac{C}{(d\vee n)^{3}}+\frac{Cs_{1}\log(K\vee
n)}{n^{1/2}}+\frac{C}{(K\vee n)^{3}}$ (4.1)

where
$\widehat{s}_{tk}^{2}=\widehat{\Omega}_{t,k}^{2}+\widehat{\Omega}_{t,t}\widehat{\Omega}_{k,k}$
and $C$ is a positive constant.

Theorem 4.3, proved in Appendix A, gives the rate of the normal approximation
of the distribution of the scaled and centered entries
$\widetilde{\Omega}_{t,k}$. The right hand side in (4.1) is non-asymptotic and
is valid for each $K$, $n$ and $d$. Its first, small, term is the price to pay
for having first used the data for clustering, and it is dominated by the
other two terms. From this perspective, the clustering step is the least
taxing, as long as we can ensure its consistency, which in turn can be
guaranteed under the minimal assumptions 4.1 and 4.2 already needed for the
remaining steps.

The second, and dominant, term regards the normal approximation of the
distribution of

${n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}.$ (4.2)

Specifically, as an intermediate step, Proposition A.2 in Appendix A shows
that the difference between the c.d.f. of (4.2), scaled by
$s_{tk}=\sqrt{\Omega_{t,k}^{*2}+\Omega^{*}_{t,t}\Omega^{*}_{k,k}}$, and that
of a standard Gaussian random variable is bounded by $\frac{s_{1}\log(K\vee
n)}{n^{1/2}}$. Therefore, asymptotic normality holds as soon as this quantity
converges to zero, which agrees with the weakest sparsity conditions for
Gaussian graphical model inference in the literature (Ren et al., 2013;
Janková and van de Geer, 2017). In addition, the asymptotic variance
$s_{tk}^{2}$ agrees with the minimum variance bound in Gaussian graphical
models (Janková and van de Geer, 2017). Thus, inference based on the empirical
risk function (3.10) does not lead to any asymptotic efficiency loss. Unlike
the previous works, we do not require the bounded operator norm condition,
$\lambda_{\max}(\bm{S}^{*})\leq C$. This condition is avoided in our analysis
by using a more convenient empirical risk function (3.10), as opposed to the
log-likelihood in Jankova and van de Geer (2014), and a CLIME-type initial
estimator (3.18) satisfying (3.19), as opposed to the node-wise Lasso
estimator in Janková and van de Geer (2017).

The last term in the normal approximation is $O(\frac{1}{(K\vee n)^{3}})$
which is dominated by the second one, and is associated with the replacement
of the theoretical variance $s_{tk}^{2}$ by the estimate
$\widehat{s}_{tk}^{2}$. Finally, we note that the powers of the first and the
third term in the right hand side of (4.1) can be replaced by $2+\delta$, for
any $\delta>0$, and a change in this power also changes the associated
constant $C$ in the term $\frac{Cs_{1}\log(K\vee n)}{n^{1/2}}$. As shown in
Proposition 4.5, to obtain valid FDR control, we need $K^{2}/(K\vee
n)^{2+\delta}=o(1)$, which holds for any $\delta>0$. For simplicity, we choose
$\delta=1$ which gives the power 3.

####  4.2.2 Results for the Latent Variable Graph

In this section we show that the estimators $\widetilde{\Theta}_{t,k}$ given
by (3.26) are asymptotically normal, for all $t<k$. We define the sparsity of
the latent graph as $s_{0}\in\mathbb{N}$ such that

$\max_{1\leq j\leq K}\sum_{k=1}^{K}\mathbb{I}(\Theta_{j,k}^{*}\neq 0)\leq
s_{0}.$

Inference for the estimator $\widetilde{\Theta}_{tk}$ follows the general
approach outlined in Section 3.1. We prove in Proposition B.1 in Appendix B
that

$n^{1/2}(\widetilde{\Theta}_{t,k}-\Theta_{t,k}^{*})=\frac{1}{n^{1/2}}\sum_{i=1}^{n}\Theta_{t,t}^{*}\mathbf{v}_{t}^{*T}(\bar{\mathbf{C}}^{(i)}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k})+o_{p}(1),$ (4.3)

where ${\mathbf{v}}^{*}_{t}$ is a $K$-dimensional vector with
$(\mathbf{v}_{t}^{*})_{t}=1$ and
$(\mathbf{v}_{t}^{*})_{-t}=-\mathbf{w}^{*}_{t}$ with
$\mathbf{w}^{*}_{t}=(\mathbf{C}^{*}_{-t,-t})^{-1}\mathbf{C}^{*}_{-t,t}$. and
$\bar{\mathbf{C}}^{(i)}$ is defined in (3.24). The terms of the sum in display
(4.3) are mean zero random variables, and their variance is

$\sigma^{2}_{tk}=\mathbb{E}(\Theta_{t,t}^{*}\mathbf{v}_{t}^{*T}(\bar{\mathbf{C}}^{(i)}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k}))^{2},$

which does not have an explicit closed form, unlike the asymptotic variance of
the estimates of the entries of $\bm{\Omega}^{*}$ above. However, we show in
Proposition B.7 in Appendix B that $\sigma^{2}_{tk}$ admits an approximation
that is easy to estimate:

$\Big{|}\sigma^{2}_{tk}-[(\Theta_{t,k}^{*})^{2}+\Theta_{k,k}^{*}\Theta_{t,t}^{*}]\Big{|}\lesssim\frac{s_{0}}{m},$

where $m=\min_{1\leq k\leq K}|G^{*}_{k}|$. Guided by this approximation, we
estimate $\sigma^{2}_{tk}$ by

$\widehat{\sigma}^{2}_{tk}=\widehat{\Theta}_{t,k}^{2}+\widehat{\Theta}_{k,k}\widehat{\Theta}_{t,t}.$

When all clusters have the equal size, we obtain $K=d/m$. Thus the
$O(\frac{s_{0}}{m})$ terms can be ignored asymptotically in the sense that
$\frac{s_{0}}{m}=\frac{s_{0}K}{d}\leq\frac{K^{2}}{d}=o(1)$, when the clusters
are approximately balanced, and their number satisfies $K^{2}=o(d)$. This is a
reasonable assumption in most applied clustering problems. We note that the
estimator $\widehat{\sigma}^{2}_{tk}$ may be inconsistent when the size of
some clusters is too small. However, we recall that our ultimate goal is to
use these estimators for recovering the sparsity pattern of $\bm{\Omega}^{*}$
under FDR control. To evaluate the sensitivity of our overall procedure to the
size of the smallest cluster, we conduct simulation studies in Section 5. The
results shows that the proposed method works well as soon as $m>4$.

The following theorem gives the Berry–Esseen normal approximation bound for
the estimators of the entries of the precision matrix corresponding to the
latent variable graph.

######  Theorem 4.4.

If Assumptions 4.1 and 4.2 hold, then

$\max_{1\leq t<k\leq
K}\sup_{x\in\mathbb{R}}\Big{|}\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Theta}_{t,k}-\Theta^{*}_{t,k})}{\widehat{\sigma}_{tk}}\leq
x\Big{)}-\Phi(x)\Big{|}\leq\frac{C}{(d\vee n)^{3}}+\frac{C}{(K\vee
n)^{3}}+\frac{Cs_{0}\log(K\vee n)}{n^{1/2}}+\frac{Cs_{0}}{m},$ (4.4)

where $C$ is a positive constant.

Compared to the average graph, the Berry-Esseen bound in (4.4) contains an
additional $\mathcal{O}(\frac{s_{0}}{m})$ term, stemming from the
approximation of the analytically intractable asymptotic variance by an
estimable quantity. The proof is deferred to Appendix B.

###  4.3 Post-clustering FDR Control

Given the edge-wise inferential results for the cluster-average and latent
variable graphs established above, we explain in this section how to combine
them to control graph-wise inferential uncertainty. We view the problem of
recovering the sparsity pattern of a precision matrix as one of multiple
hypotheses testing, by setting:

$\displaystyle\mathbf{H}_{0;tk}:\Omega^{*}_{t,k}=0\quad{\rm
vs.}\quad\mathbf{H}_{1;tk}:\Omega^{*}_{t,k}\neq 0\quad\text{for all }1\leq
t<k\leq K,$ (4.5)

for the cluster-average graph, and

$\displaystyle\mathbf{H}_{0;tk}^{{}^{\prime}}:\Theta_{t,k}^{*}=0\quad{\rm
vs.}\quad\mathbf{H}_{1;tk}^{{}^{\prime}}:\Theta_{t,k}^{*}\neq 0\quad\text{for
all }1\leq t<k\leq K.$ (4.6)

for the latent variable graph. In the following, we describe our procedure for
the estimation of the cluster-average graph. The same procedure applies to the
latent variable graph.

Define the set of true null hypotheses,
$\mathcal{H}_{0}\coloneqq\\{(t,k):\,1\leq t<k\leq K,\textrm{ such that
}\Omega^{*}_{t,k}=0\\}$, as the set of indices $(t,k)$ for which there is no
edge between the nodes $t$ and $k$. To control the error incurred by multiple
testing, we focus on the false discovery rate (FDR), which is the average
number of Type I errors relative to the total number of discoveries (Benjamini
and Hochberg, 1995). Recall that $\widetilde{\Omega}_{t,k}$ is a consistent
and asymptotically normal estimator of $\Omega_{t,k}$. We consider the natural
test statistic
$\widetilde{W}_{t,k}=n^{1/2}\widetilde{\Omega}_{t,k}/\widehat{s}_{tk}$ for
$\mathbf{H}_{0;tk}$, where
$\widehat{s}^{2}_{tk}=\widehat{\Omega}_{t,k}^{2}+\widehat{\Omega}_{t,t}\widehat{\Omega}_{k,k}$.
Given a cutoff $\tau>0$, the total number of discoveries is

$R_{\tau}:=\sum_{1\leq t<k\leq K}\mathbb{I}[|\widetilde{W}_{t,k}|>\tau].$

Similarly, the number of false positives or false discoveries is given by

$V_{\tau}:=\sum_{(t,k)\in\mathcal{H}_{0}}\mathbb{I}[|\widetilde{W}_{t,k}|>\tau].$

The FDR is formally defined as the expected ratio of $V_{\tau}$ over
$R_{\tau}$,

$\mathrm{FDR}(\tau):=\mathbb{E}\left[\frac{V_{\tau}}{R_{\tau}}\mathbb{I}[R_{\tau}>0]\right],$

where the indictor function is included to remove the trivial case
$R_{\tau}=0$.

Our goal is to find a data-dependent cutoff $\tau$ such that
$\mathrm{FDR}(\tau)\leq\alpha+o(1)$ for any given $0<\alpha<1$. This is the
best one can hope for when, as in our case, the distribution of the test
statistics $\widetilde{W}_{t,k}$ is only available asymptotically. The Berry-
Esseen type bounds derived in Theorem 4.3 and 4.4 above allow us to quantify
precisely the price that needs to be paid for the asymptotic approximation,
and become instrumental for understanding asymptotic FDR control.

In addition, the test statistics $\widetilde{W}_{t,k}$ for different
hypotheses are dependent. To allow for the dependence, instead of the standard
B-H procedure (Benjamini and Hochberg, 1995), we consider the more flexible
B-Y procedure by Benjamini and Yekutieli (2001). The resulting FDR procedure
is as follows: reject all hypotheses such that
$|\widetilde{W}_{t,k}|\geq\widehat{\tau}$, where

$\widehat{\tau}:=\min\left\\{\tau>0:\tau\geq\Phi^{-1}\left(1-\frac{\alpha
R_{\tau}}{2N_{BY}|\mathcal{H}|}\right)\right\\}\text{ and
}N_{BY}=\sum_{i=1}^{|\mathcal{H}|}\frac{1}{i},$ (4.7)

where $|\mathcal{H}|=K(K-1)/2$ is the total number of hypotheses. Our next
result shows when the FDR based on our test statistics is guaranteed to be no
greater than $\alpha$, asymptotically.

######  Proposition 4.5.

  1. 1.

Assume that the conditions in Theorem 4.3 hold. For any $0<\alpha<1$, we have

$\mathrm{FDR}(\widehat{\tau})\leq\alpha+2|\mathcal{H}_{0}|b_{n},$ (4.8)

where $b_{n}=\frac{C}{(d\vee n)^{3}}+\frac{C}{(K\vee
n)^{3}}+\frac{Cs_{1}\log(K\vee n)}{n^{1/2}}$, and $|\mathcal{H}_{0}|$ is the
number of true null hypotheses in (4.5).

  2. 2.

Assume that the conditions in Theorem 4.4 hold. If we define the test
statistic as
$\widetilde{W}_{t,k}=n^{1/2}\widetilde{\Theta}_{t,k}/\widehat{\sigma}_{tk}$,
and $\widehat{\tau}$ as in (4.7), we have

$\mathrm{FDR}(\widehat{\tau})\leq\alpha+2|\mathcal{H}^{\prime}_{0}|c_{n},$
(4.9)

where $c_{n}=\frac{C}{(d\vee n)^{3}}+\frac{C}{(K\vee
n)^{3}}+\frac{Cs_{0}\log(K\vee n)}{n^{1/2}}+\frac{Cs_{0}}{m}$, and
$|\mathcal{H}^{\prime}_{0}|$ is the number of true null hypotheses in (4.6).

Thus, this theorem implies that our method can control the FDR asymptotically,
in the sense that $\mathrm{FDR}(\widehat{\tau})\leq\alpha+o_{p}(1)$, provided
that $s_{1}|\mathcal{H}_{0}|\log(K\vee n)=o(n^{1/2})$, for the average graph,
and $s_{0}|\mathcal{H}^{{}^{\prime}}_{0}|\log(K\vee n)=o(n^{1/2})$ for the
latent graph. Gaussian graphical model estimation under FDR control was
recently studied by Liu (2013). Their approach is based on the following
Cramer-type moderate deviation result using our terminology,

$\max_{(t,k)\in\mathcal{H}_{0}}\sup_{0\leq t\leq 2\sqrt{\log
K}}\Big{|}\frac{\mathbb{P}(\widehat{T}_{t,k}\geq
t)}{2-2\Phi(t)}-1\Big{|}=o(1),$ (4.10)

where $\widehat{T}_{t,k}$ is test statistic they proposed for estimation of
the Gaussian graphical model structure. The above result (4.10) controls the
relative error of the Gaussian approximation within the moderate deviation
regime $[0,2\sqrt{\log K}]$, whereas our result is based on the control of the
absolute error via the Berry-Esseen-type Gaussian approximation. One of the
main advantages of their result is that the number of clusters is allowed to
be $K=o(n^{r})$, where $r$ is a constant that can be greater than 1. However,
to prove (4.10), they required that the number of strong signals tends to
infinity, that is $|\\{(t,k):\Omega_{t,k}^{*}\geq C\sqrt{\log
K/n}\\}|\rightarrow\infty$, which reduces significantly the parameter space
for which inference is valid. In contrast, the aim of this work is the study
of pattern recovery under no conditions on the signal strength of the entries
of the target precision matrices, as in practice it is difficult to assess
whether these conditions are met.

The overall message conveyed by Proposition 4.5 is that, in the absence of any
signal strength assumptions, cluster-based graphical models can still be
recovered, under FDR control, provided that the number of clusters $K$ is not
very high relative to $n$, and provided that the clusters are not very small.
This further stresses the importance of an initial dimension reduction step in
high-dimensional graphical model estimation. For instance, results similar to
those of Theorem 4.3 can be derived along the same lines for the estimation of
the sparsity pattern of $\bm{\Sigma}^{-1}$, for a generic, unstructured,
covariance matrix of $\bm{X}$, where one replaces $K$ by $d$ throughout, and
$s_{1}$ is replaced by $s$, the number of non-zero entries in the $d\times d$
matrix $\bm{\Sigma}^{-1}$. Then, the analogue of (4.9) of Proposition 4.5
shows that FDR control in generic graphical models, based on asymptotic
approximations of $p$-values, cannot generally be guaranteed if $d>n$. Our
work shows that extra structural assumptions, for instance those motivated by
clustering, do alleviate this problem. The simulation study presented in the
next section provides further support to our findings.

##  5 Numerical Results

This section contains simulations and a real data analysis that illustrate the
finite sample performance of the inferential procedures developed in the
previous sections for the latent variable graph and cluster-average graph,
respectively.

###  5.1 Synthetic Datasets

In this subsection, we demonstrate the effectiveness of the FDR control
procedures on synthetic datasets. We consider two settings $(n,d)=(800,400)$
and $(500,1000)$, and in each setting we vary the value of $K$ and $m$. The
error variable $\bm{E}$ is sampled from the multivariate normal distribution
with covariance $\bm{\Gamma}^{*}$ whose entries $\gamma_{i}^{*}$ are generated
from $U[0.25,0.5]$. Recall that the latent variable $\bm{Z}$ follows from
$\bm{Z}\sim N(0,\mathbf{C}^{*})$. We consider three different models to
generate the graph structure of $\bm{Z}$. Once the graph structure is
determined, the corresponding adjacency matrix $\mathbf{W}$ is found, and the
precision matrix $\bm{\Theta}^{*}=\mathbf{C}^{*-1}$ is taken as
$\bm{\Theta}^{*}=c\mathbf{W}+(|\lambda_{\min}\left(\mathbf{W}\right)|+0.2)\mathbf{I}$,
where $c=0.3$ when $d=400$ and $c=0.5$ when $d=1000$. Finally, we assign the
cluster labels for all variables so that all clusters have approximately equal
size, which gives us the matrix $\mathbf{A}$. Given $\mathbf{A},\bm{Z}$ and
$\bm{E}$, we can generate $\bm{X}$ according to the model (1.1).

We consider the following three generating models for the graph structure of
$\bm{Z}$:

  * •

Scale-Free Graph –  The Scale-Free model is a generative model for network
data, whose degree distribution follows a power law. To be concrete, we
generate the graph one node at a time, starting with a 2 node chain. For nodes
$3,\dots,K$, node $t$ is added and one edge is added between $t$ and one of
the $t-1$ previous nodes. Denoting by $k_{i}$ the current degree of node $i$
in the graph, the probability that node $t$ and node $i$ are connected is
$p_{i}=k_{i}/(\sum_{i}k_{i})$. The number of edges in the resulting graph is
always $K$.

  * •

Hub Graph –  The $K$ nodes of the graph are partitioned evenly into groups of
size $N$. Within each group, one node is selected to be the group hub, and an
edge is added between it and the remainder of its group. $N$ is either 5 or 6
depending upon the choice of $K$. The number of edges in the graph is
$K(N-1)/N$, so for $K=100$ with $N=5$, the number of edges in the resulting
graph is 80.

  * •

Band3 Graph –  This model generates a graph with a Toeplitz adjacency matrix.
There is an edge between node $i$ and node $j$ if and only if $|i-j|\leq B$,
where we set $B=3$ in this scenario. In general, the number of edges in a band
graph with $K$ nodes is given by $BK-\frac{3}{2}B^{2}+\frac{5}{2}B$. So, for
$K=100$ and $B=3$, the number of edges in the graph is 294.

Recall that $\bar{\bm{X}}\sim N(0,\bm{S}^{*})$, where $\bm{S}^{*}$ is defined
in (3.2). To determine the structure of the average graph, we numerically
compute $\bm{S}^{*-1}$ and threshold the matrix at $10^{-8}$.

We examine the empirical FDR of our procedures on some synthetic datasets. The
following protocol is followed in all the experiments:

  1. 1.

Generate the graph structure of $\bm{Z}$ as specified above.

  2. 2.

Simulate $n$ observations from our model (1.1).

  3. 3.

Estimate the cluster partition $\widehat{G}$. For computational convenience,
we apply the FORCE algorithm (Eisenach and Liu, 2017) when $d=400$ and the COD
algorithm (Bunea et al., 2018) when $d=1000$.

  4. 4.

Construct the test statistic $\widetilde{W}_{t,k}$ defined in Section 4.3. The
regularization parameters $\lambda$ and $\lambda^{\prime}$ are chosen by
5-fold cross validation.

  5. 5.

Find the FDR cutoff (4.7) at level $\alpha$, where we consider three cases
$\alpha=0.05,0.1,0.2$.

The simulation is repeated 100 times. To compare with our Benjamini-Yekutieli
based FDR procedure, we also report the empirical FDR based on the more
classical Benjamini-Hochberg procedure. That is, we apply the same procedures
1-4, but in step 5 we replace the FDR cutoff in (4.7) with the Benjamini-
Hochberg (B-H) cutoff, i.e.,

$\widehat{\tau}_{BH}:=\min\left\\{\tau>0:\tau\geq\Phi^{-1}\left(1-\frac{\alpha
R_{\tau}}{2|\mathcal{H}|}\right)\right\\}.$

Table 1 compares the empirical FDR based on our method with the B-H procedure
under different $m,K$ settings when $d=400$. When $m$ is relatively large
(e.g., $m=20$), both methods can control FDR on average, although our method
is relatively more conservative. As expected, the FDR control problem becomes
more challenging for large $K$ and small $m$. In this case the graph contains
more nodes and each cluster contains fewer variables. We observe that when
$m=5$ our method can still control FDR reasonably well but the B-H method
produces empirical FDR far beyond the nominal level, especially for the hub
graph. The inferior performance of the B-H procedure is due to the fact that
the dependence among the test statistics is not accounted in the B-H method.
Finally, we examine the empirical power of the FDR procedure under each
scenario, which is defined as

$\textrm{Average}\left[\sum_{(t,k)\in\mathcal{H}_{1}}\frac{\mathbb{I}[\widetilde{W}_{t,k}\geq\widehat{\tau}]}{|\mathcal{H}_{1}|}\right],$

where $\mathcal{H}_{1}$ is the set of alternative hypothesis. Table 2 gives
the empirical power of our FDR procedure, and the B-H procedure when $d=400$.
It shows that our procedure and the B-H procedure have very high power in all
scenarios. The same phenomenon is observed when $d$ is large, i.e., $d=1000$;
see Tables 3 and 4. In summary, the proposed procedure can identify most
signals in the graph with well controlled FDR.

|  |  |  | $\mathbf{\bm{\alpha}=0.05}$ | $\mathbf{\bm{\alpha}=0.1}$ | $\mathbf{\bm{\alpha}=0.2}$
---|---|---|---|---|---|---
|  | $K$ | $m$ | Scalefree | Band3 | Hub | Scalefree | Band3 | Hub | Scalefree | Band3 | Hub
B-Y Based Procedure |  Latent | 80 | 5 | 1.16% | 1.42% | 5.99% | 2.01% | 2.64% | 7.60% | 4.01% | 4.73% | 10.33%
66 | 6 | 0.93% | 1.03% | 1.08% | 1.99% | 1.98% | 1.96% | 3.73% | 3.68% | 4.00%
50 | 8 | 1.16% | 0.99% | 1.09% | 2.20% | 1.77% | 1.98% | 3.90% | 3.49% | 3.73%
20 | 20 | 0.99% | 0.88% | 1.30% | 1.81% | 1.66% | 2.26% | 3.75% | 3.38% | 4.71%
Grp. Av. | 80 | 5 | 1.15% | 1.40% | 6.29% | 2.14% | 2.67% | 7.89% | 4.00% | 4.70% | 10.43%
66 | 6 | 0.91% | 1.04% | 1.04% | 1.85% | 2.00% | 2.14% | 3.70% | 3.71% | 3.83%
50 | 8 | 1.16% | 1.00% | 1.11% | 2.22% | 1.80% | 1.94% | 3.98% | 3.49% | 3.61%
20 | 20 | 0.94% | 0.88% | 1.23% | 1.81% | 1.68% | 2.26% | 3.65% | 3.38% | 4.71%
B-H Based Procedure |  Latent | 80 | 5 | 8.23% | 9.24% | 15.60% | 15.01% | 16.39% | 24.00% | 28.16% | 28.97% | 38.58%
66 | 6 | 7.18% | 7.38% | 7.31% | 14.12% | 14.01% | 13.78% | 26.56% | 25.84% | 27.28%
50 | 8 | 6.94% | 6.75% | 6.69% | 13.23% | 12.73% | 12.74% | 25.20% | 23.62% | 26.05%
20 | 20 | 5.43% | 4.54% | 6.38% | 11.09% | 8.47% | 10.71% | 21.46% | 17.15% | 20.75%
Grp. Av. | 80 | 5 | 8.42% | 9.26% | 15.77% | 15.25% | 16.43% | 24.15% | 28.03% | 29.02% | 38.57%
66 | 6 | 7.21% | 7.37% | 7.45% | 13.99% | 13.92% | 13.81% | 26.19% | 25.77% | 27.17%
50 | 8 | 6.82% | 6.78% | 6.67% | 13.23% | 12.67% | 12.89% | 25.10% | 23.61% | 25.91%
20 | 20 | 5.38% | 4.51% | 6.49% | 11.01% | 8.52% | 10.66% | 21.58% | 17.09% | 20.71%
Table 1: Averaged empirical FDR for synthetic data experiments with $d=400$ and $n=800$.  |  |  |  | $\mathbf{\bm{\alpha}=0.05}$ | $\mathbf{\bm{\alpha}=0.1}$ | $\mathbf{\bm{\alpha}=0.2}$
---|---|---|---|---|---|---
|  | $K$ | $m$ | Scalefree | Band3 | Hub | Scalefree | Band3 | Hub | Scalefree | Band3 | Hub
B-Y Based Procedure |  Latent | 80 | 5 | 88.22% | 99.99% | 99.00% | 91.10% | 100.00% | 99.04% | 93.38% | 100.00% | 99.11%
66 | 6 | 93.62% | 100% | 100% | 95.42% | 100% | 100% | 96.86% | 100% | 100%
50 | 8 | 97.18% | 100% | 100% | 97.90% | 100% | 100% | 98.47% | 100% | 100%
20 | 20 | 100% | 100% | 100% | 100% | 100% | 100% | 100% | 100% | 100%
Grp. Av. | 80 | 5 | 88.27% | 99.99% | 98.93% | 91.06% | 100.00% | 98.98% | 93.37% | 100.00% | 99.07%
66 | 6 | 93.62% | 100% | 100% | 95.40% | 100% | 100% | 96.91% | 100% | 100%
50 | 8 | 97.18% | 100% | 100% | 97.86% | 100% | 100% | 98.47% | 100% | 100%
20 | 20 | 100% | 100% | 100% | 100% | 100% | 100% | 100% | 100% | 100%
B-H Based Procedure |  Latent | 80 | 5 | 95.75% | 100.00% | 99.18% | 97.15% | 100.00% | 99.23% | 98.48% | 100.00% | 99.26%
66 | 6 | 97.85% | 100% | 100% | 98.78% | 100% | 100% | 99.34% | 100% | 100%
50 | 8 | 99.27% | 100% | 100% | 99.57% | 100% | 100% | 99.80% | 100% | 100%
20 | 20 | 100% | 100% | 100% | 100% | 100% | 100% | 100% | 100% | 100%
Grp. Av. | 80 | 5 | 95.82% | 100.00% | 99.14% | 97.14% | 100.00% | 99.18% | 98.53% | 100.00% | 99.23%
66 | 6 | 97.88% | 100% | 100% | 98.78% | 100% | 100% | 99.34% | 100% | 100%
50 | 8 | 99.24% | 100% | 100% | 99.59% | 100% | 100% | 99.78% | 100% | 100%
20 | 20 | 100% | 100% | 100% | 100% | 100% | 100% | 100% | 100% | 100%
Table 2: Averaged FDR power for synthetic data experiments with $d=400$ and $n=800$.  |  |  |  | $\mathbf{\bm{\alpha}=0.05}$ | $\mathbf{\bm{\alpha}=0.1}$ | $\mathbf{\bm{\alpha}=0.2}$
---|---|---|---|---|---|---
|  | $K$ | $m$ | Scalefree | Band3 | Hub | Scalefree | Band3 | Hub | Scalefree | Band3 | Hub
B-Y Based Procedure |  Latent | 80 | 12 | 1.40% | 1.52% | 1.78% | 2.78% | 2.81% | 3.29% | 5.71% | 5.09% | 6.30%
66 | 15 | 1.23% | 1.29% | 1.31% | 2.62% | 2.39% | 2.69% | 5.08% | 4.33% | 4.84%
50 | 20 | 1.05% | 0.94% | 1.26% | 2.27% | 1.85% | 2.79% | 4.28% | 3.65% | 4.58%
20 | 50 | 1.04% | 1.03% | 1.05% | 1.81% | 1.96% | 2.20% | 3.46% | 3.26% | 4.19%
Grp. Av. | 80 | 12 | 1.55% | 1.57% | 1.78% | 2.78% | 2.80% | 3.26% | 5.57% | 5.13% | 6.13%
66 | 15 | 1.30% | 1.30% | 1.24% | 2.64% | 2.35% | 2.64% | 5.07% | 4.34% | 4.89%
50 | 20 | 1.09% | 0.93% | 1.26% | 2.32% | 1.89% | 2.75% | 4.29% | 3.64% | 4.51%
20 | 50 | 1.09% | 1.01% | 1.05% | 1.91% | 1.96% | 2.14% | 3.50% | 3.28% | 4.25%
B-H Based Procedure |  Latent | 80 | 12 | 10.99% | 10.04% | 12.00% | 20.04% | 18.18% | 21.47% | 35.50% | 31.46% | 37.24%
66 | 15 | 9.28% | 8.20% | 9.14% | 16.88% | 14.94% | 16.82% | 31.53% | 27.75% | 31.81%
50 | 20 | 7.73% | 6.91% | 7.96% | 14.74% | 13.13% | 15.33% | 27.03% | 23.92% | 28.72%
20 | 50 | 5.00% | 4.85% | 6.38% | 11.01% | 9.38% | 11.60% | 21.49% | 17.00% | 23.15%
Grp. Av. | 80 | 12 | 10.99% | 10.00% | 12.00% | 20.29% | 18.20% | 21.38% | 35.56% | 31.48% | 37.10%
66 | 15 | 9.28% | 8.24% | 8.99% | 17.03% | 14.87% | 16.87% | 31.49% | 27.82% | 31.63%
50 | 20 | 7.60% | 6.93% | 7.94% | 14.61% | 13.13% | 15.40% | 27.12% | 23.94% | 28.71%
20 | 50 | 4.95% | 4.91% | 6.38% | 10.97% | 9.35% | 11.60% | 21.36% | 16.94% | 23.41%
Table 3: Averaged empirical FDR for synthetic data experiments with $d=1000$ and $n=500$.  |  |  |  | $\mathbf{\bm{\alpha}=0.05}$ | $\mathbf{\bm{\alpha}=0.1}$ | $\mathbf{\bm{\alpha}=0.2}$
---|---|---|---|---|---|---
|  | $K$ | $m$ | Scalefree | Band3 | Hub | Scalefree | Band3 | Hub | Scalefree | Band3 | Hub
B-Y Based Procedure |  Latent | 80 | 12 | 71.42% | 99.95% | 100.00% | 77.14% | 99.97% | 100.00% | 82.29% | 99.98% | 100.00%
66 | 15 | 82.98% | 99.98% | 100.00% | 86.32% | 99.99% | 100.00% | 90.02% | 100.00% | 100.00%
50 | 20 | 90.45% | 99.99% | 100.00% | 93.04% | 99.99% | 100.00% | 94.86% | 100.00% | 100.00%
20 | 50 | 99.95% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00%
Grp. Av. | 80 | 12 | 71.54% | 99.95% | 100.00% | 77.13% | 99.97% | 100.00% | 82.25% | 99.98% | 100.00%
66 | 15 | 82.74% | 99.98% | 100.00% | 86.26% | 99.99% | 100.00% | 90.08% | 100.00% | 100.00%
50 | 20 | 90.45% | 99.99% | 100.00% | 93.00% | 99.99% | 100.00% | 94.76% | 100.00% | 100.00%
20 | 50 | 99.95% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00%
B-H Based Procedure |  Latent | 80 | 12 | 87.15% | 99.99% | 100.00% | 91.33% | 100.00% | 100.00% | 94.66% | 100.00% | 100.00%
66 | 15 | 93.57% | 100.00% | 100.00% | 95.63% | 100.00% | 100.00% | 97.38% | 100.00% | 100.00%
50 | 20 | 96.27% | 100.00% | 100.00% | 97.78% | 100.00% | 100.00% | 98.94% | 100.00% | 100.00%
20 | 50 | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00%
Grp. Av. | 80 | 12 | 87.22% | 99.99% | 100.00% | 91.33% | 100.00% | 100.00% | 94.58% | 100.00% | 100.00%
66 | 15 | 93.52% | 100.00% | 100.00% | 95.65% | 100.00% | 100.00% | 97.35% | 100.00% | 100.00%
50 | 20 | 96.33% | 100.00% | 100.00% | 97.78% | 100.00% | 100.00% | 98.94% | 100.00% | 100.00%
20 | 50 | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00% | 100.00%
Table 4: Averaged FDR power for synthetic data experiments with $d=1000$ and
$n=500$.

###  5.2 fMRI Dataset

Power et al. (2011) finds that the human brain can be divided into regions of
interest (ROIs) which can be further organized into functional networks. We
apply our inferential procedures to these regions of interest in publicly
available resting-state fMRI data from the Neuro-bureau pre-processed
repository (Bellec et al., 2015). Specifically, we use the data from patient
1018959, session 1 in the KKI dataset. This fMRI data was pre-processed using
the Athena pipeline and mapped to T1 MNI152 coordinate space. We choose this
dataset to make our experiments easily reproducible, as the data are available
pre-processed using standard alignment and denoising procedures.

Using the T1 MNI152 coordinates, we extract the 264 ROIs identified in Power
et al. (2011), which gives us $d=264$ mean activities across $n=148$ time
periods. We apply the FORCE algorithm to cluster the 264 ROIs, and obtain an
estimate of $K=53$ and the corresponding clusters. Using the FDR control
procedures with $\alpha=0.01$, we obtain the networks shown in Figures 1 and
2. For clarification purpose, in these two figures, we only display the nodes
and connections corresponding to the 10 largest groups. The groups are colored
according to the functional network the majority of their nodes belong to as
given in Power et al. (2011). In the latent graph, the nodes $Z_{5}$,
$Z_{29}$, $Z_{30}$, and $Z_{32}$ are highly connected. This finding is
consistent with Power et al. (2011) where the authors note that the graph of
the observable variables is highly connected within a functional group. By
contrast, the average graph shows completely different patterns. For instance,
$\bar{X}_{5}$, $\bar{X}_{29}$, $\bar{X}_{30}$, and $\bar{X}_{32}$ which belong
to the same functional group are not connected in the cluster-averages graph.
If one uses the cluster-averages graph to interpret the dependence structure
of the functional network, the scientific results could be misleading.

Figure 1: Recovered latent graph structure between 10 largest clusters in fMRI
data with FDR level $\alpha=0.01$ colored according to their functions. Figure
2: Recovered cluster-averages graph structure between 10 largest clusters in
fMRI data with FDR level $\alpha=0.01$.

## Acknowledgement

Florentina Bunea was partially supported by NSF-DMS 712709.

## References

  * Barber and Kolar (2015) Barber, R. F. and Kolar, M. (2015).  Rocket: Robust confidence intervals via kendall’s tau for transelliptical graphical models.  arXiv:1502.07641.
  * Bellec et al. (2015) Bellec, P., Chu, C., Chouinard-Decorte, F., Benhajali, Y., Margulies, D. S. and Craddock, R. C. (2015).  The Neuro Bureau ADHD-200 Preprocessed repository.  NeuroImage .
  * 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) 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 29 1165–1188.
  * Bickel (1975) Bickel, P. J. (1975).  One-step huber estimates in the linear model.  Journal of the American Statistical Association 70 428–434.
  * Bickel et al. (1993) Bickel, P. J., Klaassen, C. A., Bickel, P. J., Ritov, Y., Klaassen, J., Wellner, J. A. and Ritov, Y. (1993).  Efficient and adaptive estimation for semiparametric models.  Johns Hopkins University Press Baltimore.
  * Bunea (2004) Bunea, F. (2004).  Consistent covariate selection and post model selection inference in semiparametric regression.  Annals of Statistics 898–927.
  * Bunea et al. (2018) Bunea, F., Giraud, C., Luo, X., Royer, M. and Verzelen, N. (2018).  Model Assisted Variable Clustering: Minimax-optimal Recovery and Algorithms.  arXiv:1508.01939.
  * 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. (2017) Cai, T. T., Guo, Z. et al. (2017).  Confidence intervals for high-dimensional linear regression: Minimax rates and adaptivity.  The Annals of Statistics 45 615–646.
  * Cai et al. (2016) Cai, T. T., Ren, Z. and Zhou, H. H. (2016).  Estimating structured high-dimensional covariance and precision matrices: Optimal rates and adaptive estimation.  Electronic Journal of Statistics 10 1–59.
  * Eisenach and Liu (2017) Eisenach, C. and Liu, H. (2017).  Efficient, Certifiably Optimal High-Dimensional Clustering.  arXiv:1806.00530.
  * Fan et al. (2017) Fan, J., Liu, H., Ning, Y. and Zou, H. (2017).  High dimensional semiparametric latent graphical model for mixed data.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79 405–421.
  * Friedman et al. (2008) Friedman, J., Hastie, T. and Tibshirani, R. (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics 9 432–441.
  * Giraud (2014) Giraud, C. (2014).  Introduction to high-dimensional statistics, vol. 138.  CRC Press.
  * Gu et al. (2015) Gu, Q., Cao, Y., Ning, Y. and Liu, H. (2015).  Local and global inference for high dimensional gaussian copula graphical models.  arXiv:1502.02347.
  * Jankova and van de Geer (2014) Jankova, J. and van de Geer, S. (2014).  Confidence intervals for high-dimensional inverse covariance estimation.  arXiv:1403.6752.
  * 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.
  * Javanmard and Montanari (2013) Javanmard, A. and Montanari, A. (2013).  Confidence intervals and hypothesis testing for high-dimensional regression.  arXiv:1306.3171.
  * Lam and Fan (2009) Lam, C. and Fan, J. (2009).  Sparsistency and rates of convergence in large covariance matrix estimation.  Annals of Statistics 37 42–54.
  * Lee et al. (2013) Lee, J. D., Sun, D. L., Sun, Y. and Taylor, J. E. (2013).  Exact inference after model selection via the lasso.  arXiv:1311.6238.
  * Liu et al. (2012) Liu, H., Han, F., Yuan, M., Lafferty, J. and Wasserman, L. (2012).  High-dimensional semiparametric gaussian copula graphical models.  The Annals of Statistics 40 2293–2326.
  * Liu (2013) Liu, W. (2013).  Gaussian graphical model estimation with false discovery rate control.  The Annals of Statistics 2948–2978.
  * Lockhart et al. (2014) Lockhart, R., Taylor, J., Tibshirani, R. J., Tibshirani, R. et al. (2014).  A significance test for the lasso.  The Annals of Statistics 42 413–468.
  * Meinshausen and Bühlmann (2006) Meinshausen, N. and Bühlmann, P. (2006).  High dimensional graphs and variable selection with the lasso.  Annals of Statistics 34 1436–1462.
  * Neykov et al. (2018) Neykov, M., Ning, Y., Liu, J. S. and Liu, H. (2018).  A unified theory of confidence regions and testing for high dimensional estimating equations.  Statistical Science .
  * Ning and Liu (2013) Ning, Y. and Liu, H. (2013).  High-dimensional semiparametric bigraphical models.  Biometrika 100 655–670.
  * Ning and Liu (2017) Ning, Y. and Liu, H. (2017).  A general theory of hypothesis tests and confidence regions for sparse high dimensional models.  The Annals of Statistics 45 158–195.
  * 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. (2011) Power, J. D., Cohen, A. L., Nelson, S. M., Wig, G. S., Barnes, K. A., Church, J. A., Vogel, A. C., Laumann, T. O., Miezin, F. M., Schlaggar, B. L. and Petersen, S. E. (2011).  Functional network organization of the human brain.  Neuron 17 665–678.
  * Ravikumar et al. (2011a) Ravikumar, P., Wainwright, M., Raskutti, G. and Yu, B. (2011a).  High-dimensional covariance estimation by minimizing $\ell_{1}$-penalized log-determinant divergence.  Electronic Journal of Statistics 5 935–980.
  * Ravikumar et al. (2011b) Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B. (2011b).  High-dimensional covariance estimation by minimizing $l\\_1$-penalized log-determinant divergence.  Electronic Journal of Statistics 5 935–980.
  * Ren et al. (2013) Ren, Z., Sun, T., Zhang, C.-H. and Zhou, H. H. (2013).  Asymptotic normality and optimalities in estimation of large gaussian graphical model.  arXiv:1309.6024.
  * 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.
  * Small et al. (2000) Small, C. G., Wang, J. and Yang, Z. (2000).  Eliminating multiple root problems in estimation.  Statistical Science 15 313–341.
  * Sun and Zhang (2012) Sun, T. and Zhang, C.-H. (2012).  Sparse matrix inversion with scaled lasso.  Tech. rep.
  * Taylor et al. (2014) Taylor, J., Lockhart, R., Tibshirani, R. J. and Tibshirani, R. (2014).  Post-selection adaptive inference for least angle regression and the lasso.  arXiv:1401.3889.
  * van de Geer et al. (2013) van de Geer, S., Bühlmann, P. and Ritov, Y. (2013).  On asymptotically optimal confidence regions and tests for high-dimensional models.  arXiv:1303.0518.
  * Van der Vaart (1998) Van der Vaart, A. V. (1998).  Asymptotic statistics.  Cambridge University Press, Cambridge, UK.
  * Verzelen (2008) Verzelen, N. (2008).  Gaussian graphical models and Model selection.  Ph.D. thesis, Université Paris Sud-Paris XI.
  * Xue and Zou (2012) Xue, L. and Zou, H. (2012).  Regularized rank-based estimation of high-dimensional nonparanormal graphical models.  The Annals of Statistics 40 2541–2571.
  * Yuan (2010) Yuan, M. (2010).  High dimensional inverse covariance matrix estimation via linear programming.  Journal of Machine Learning Research 11 2261–2286.
  * Yuan and Lin (2007) Yuan, M. and Lin, Y. (2007).  Model selection and estimation in the gaussian graphical model.  Biometrika 94 19–35.
  * Zhang and Zhang (2014) Zhang, C.-H. and Zhang, S. S. (2014).  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76 217–242.

##  Appendix A Proofs Regarding Estimation in the Cluster-Average Graph

In this section we provide the proofs of the results needed for establishing
the asymptotic normality of the estimators of the entries of
$\bm{\Omega}^{*}$. The results below make use of the fact that consistent
clustering is possible, under our assumptions, as stated below.

######  Lemma A.1.

Let $\mathcal{E}=:\\{\widehat{G}=G^{*}\\}$, for $\widehat{G}$ estimated by
either the COD or the PECOK algorithm of Bunea et al. (2018). Then, under
Assumptions 4.1 and 4.2, we have

$\mathbb{P}(\mathcal{E})\geq 1-\frac{C}{(d\vee n)^{3}}.$

The conclusion of this Lemma is proved in Theorem 3, for the COD algorithm,
and Theorem 4, for the PECOK algorithm, of Bunea et al. (2018).

Lemma A.1 allows us to replace $\widehat{G}$ by $G^{*}$ in all the results
below, while incurring a small error, of order $O\left(\frac{C}{(d\vee
n)^{3}}\right)$, which will be shown to be dominated by other error bounds.

###  A.1 Main Proofs for the Cluster-Average Graph Estimators

###### Proof of Theorem 4.3.

The proof relies crucially on Proposition A.2, stated and proved after the end
of this proof. Assuming that this result has been obtained, the proof of
Theorem 4.3 follows the standard steps explained below.

Denote $\mathcal{E}^{\prime}=\\{\max_{1\leq t<k\leq
K}|\widehat{s}_{tk}/s_{tk}-1|\leq r\\}$, where $r=C\sqrt{\frac{s_{1}\log(K\vee
n)}{n}}$, and $\bar{\mathcal{E}}^{\prime}$ is the complement of the event
$\mathcal{E}^{\prime}$. We first consider the bound

$\displaystyle\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{\widehat{s}_{tk}}\leq
x\Big{)}-\Phi(x)$
$\displaystyle\leq\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{\widehat{s}_{tk}}\leq
x,\mathcal{E}^{\prime},\mathcal{E}\Big{)}-\Phi(x)+\mathbb{P}(\bar{\mathcal{E}}^{\prime})+\mathbb{P}(\bar{\mathcal{E}})$
$\displaystyle=\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{s_{tk}}\leq
x\frac{\widehat{s}_{tk}}{s_{tk}},\mathcal{E}^{\prime},\mathcal{E}\Big{)}-\Phi(x)+\mathbb{P}(\bar{\mathcal{E}}^{\prime})+\mathbb{P}(\bar{\mathcal{E}}).$

Proposition A.3 below implies $\mathbb{P}(\bar{\mathcal{E}}^{\prime})\leq
C(K\vee n)^{-3}$ and Lemma A.1 above implies
$\mathbb{P}(\bar{\mathcal{E}})\leq C(d\vee n)^{-3}$ for some constant $C$. In
addition, for $x\geq 0$,

$\displaystyle\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{s_{tk}}\leq
x\frac{\widehat{s}_{tk}}{s_{tk}},\mathcal{E}^{\prime},\mathcal{E}\Big{)}-\Phi(x)\leq\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{s_{tk}}\leq
x(1+r),\mathcal{E}\Big{)}-\Phi(x)$
$\displaystyle=\Big{\\{}\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{s_{tk}}\leq
x(1+r),\mathcal{E}\Big{)}-\Phi(x(1+r))\Big{\\}}+\Big{\\{}\Phi(x(1+r))-\Phi(x)\Big{\\}}.$
(A.1)

For the first term, Proposition A.2 implies

$\max_{1\leq t<k\leq
K}\sup_{x\in\mathbb{R}}\Big{|}\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{s_{tk}}\leq
x(1+r),\mathcal{E}\Big{)}-\Phi(x(1+r))\Big{|}\lesssim\frac{s_{1}\log(K\vee
n)}{n^{1/2}}.$

By the mean value theorem, the second term
$\Phi(x(1+r))-\Phi(x)=\phi(x(1+tr))xr$, for some $t\in[0,1]$. It is easily
seen that $\sup_{x\in\mathbb{R}}\sup_{t\in[0,1]}|\phi(x(1+tr))x|\leq C$ for
some constant $C$. Plugging it into (A.1), we obtain

$\max_{1\leq t<k\leq
K}\sup_{x\in\mathbb{R}}\Big{\\{}\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{s_{tk}}\leq
x\frac{\widehat{s}_{tk}}{s_{tk}},\mathcal{E}^{\prime},\mathcal{E}\Big{)}-\Phi(x)\Big{\\}}\lesssim\frac{s_{1}\log(K\vee
n)}{n^{1/2}}+r\lesssim\frac{s_{1}\log(K\vee n)}{n^{1/2}}.$ (A.2)

When $x<0$, similar to (A.1), the bound is

$\displaystyle\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{s_{tk}}\leq
x\frac{\widehat{s}_{tk}}{s_{tk}},\mathcal{E}^{\prime},\mathcal{E}\Big{)}-\Phi(x)$
$\displaystyle\leq\Big{\\{}\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{s_{tk}}\leq
x(1-r),\mathcal{E}\Big{)}-\Phi(x(1-r))\Big{\\}}+\Big{\\{}\Phi(x(1-r))-\Phi(x)\Big{\\}}.$

Thus (A.2) holds for $x<0$ as well. Combining these results, we obtain

$\max_{1\leq t<k\leq
K}\sup_{x\in\mathbb{R}}\Big{\\{}\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{\widehat{s}_{tk}}\leq
x\Big{)}-\Phi(x)\Big{\\}}\lesssim\frac{s_{1}\log(K\vee
n)}{n^{1/2}}+\frac{1}{(K\vee n)^{3}}+\frac{1}{(d\vee n)^{3}}.$

Following the similar argument, we can also derive

$\max_{1\leq t<k\leq
K}\sup_{x\in\mathbb{R}}\Big{\\{}\Phi(x)-\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{\widehat{s}_{tk}}\leq
x\Big{)}\Big{\\}}\lesssim\frac{s_{1}\log(K\vee n)}{n^{1/2}}+\frac{1}{(K\vee
n)^{3}}+\frac{1}{(d\vee n)^{3}}.$

This completes the proof. ∎

######  Proposition A.2 (Asymptotic Normality of $\widetilde{\Omega}_{t,k}$).

Under the same conditions as in Theorem 4.3, we have

$\max_{1\leq t<k\leq
K}\sup_{x\in\mathbb{R}}\Big{|}\mathbb{P}\Big{(}\frac{n^{1/2}(\widetilde{\Omega}_{t,k}-\Omega^{*}_{t,k})}{s_{tk}}\leq
x,\mathcal{E}\Big{)}-\Phi(x)\Big{|}\leq\frac{C}{(K\vee
n)^{3}}+\frac{Cs_{1}\log(K\vee n)}{n^{1/2}}.$

###### Proof.

The proof is done in two steps. In Step 1, we show that intersected with the
event $\mathcal{E}$,

$n^{1/2}|(\widetilde{\Omega}_{t,k}-\Omega_{t,k}^{*})/s_{t,k}+\Omega_{t,t}^{*}h(\bm{\Omega}_{\cdot
k}^{*})/s_{t,k}|\leq\frac{s_{1}\log(K\vee n)}{n^{1/2}},$ (A.3)

with probability at least $1-C/(K\vee n)^{3}$ and then use Lemma A.9 to obtain
the result. To prove (A.3), we decompose it as

$\displaystyle
n^{1/2}|(\widetilde{\Omega}_{t,k}-\Omega_{t,k}^{*})+\Omega_{t,t}^{*}h(\bm{\Omega}_{\cdot
k}^{*})|$
$\displaystyle=n^{1/2}|(\widehat{\Omega}_{t,k}-\Omega_{t,k}^{*})-\widehat{\Omega}_{t,t}\widehat{h}(\widehat{\bm{\Omega}}_{\cdot
k})+\Omega_{t,t}^{*}h(\bm{\Omega}_{\cdot k}^{*})|$
$\displaystyle\leq\underbrace{n^{1/2}|(\widehat{\Omega}_{t,k}-\Omega_{t,k}^{*})-\Omega_{t,t}^{*}(h(\widehat{\bm{\Omega}}_{\cdot
k})-h(\bm{\Omega}_{\cdot k}^{*}))|}_{\text{I.1}}$
$\displaystyle\quad+\underbrace{n^{1/2}|\Omega_{t,t}^{*}(\widehat{h}(\widehat{\bm{\Omega}}_{\cdot
k})-h(\widehat{\bm{\Omega}}_{\cdot
k}))|}_{\text{I.2}}+\underbrace{n^{1/2}|(\widehat{\Omega}_{t,t}-\Omega_{t,t}^{*})\widehat{h}(\widehat{\bm{\Omega}}_{\cdot
k})|}_{\text{I.3}}.$

In the following, we study these three terms separately. Recall that
$h(\bm{\Omega}_{\cdot
k})=\mathbf{v}_{t}^{*T}(\widehat{\bm{S}}\bm{\Omega}_{\cdot
k}-\mathbf{e}_{k}),$ and $\widehat{h}(\bm{\Omega}_{\cdot
k})=\widehat{\mathbf{v}}_{t}^{T}(\widehat{\bm{S}}\bm{\Omega}_{\cdot
k}-\mathbf{e}_{k})$, where ${\mathbf{v}}^{*}_{t}$ is a $K$-dimensional vector
with $(\mathbf{v}_{t}^{*})_{t}=1$ and
$(\mathbf{v}_{t}^{*})_{-t}=-\mathbf{w}^{*}_{t}$ with
$\mathbf{w}^{*}_{t}=(\bm{S}^{*}_{-t,-t})^{-1}\bm{S}^{*}_{-t,t}$. Term I.1
reduces to

$\displaystyle|I.1|$
$\displaystyle=n^{1/2}|(\widehat{\Omega}_{t,k}-\Omega_{t,k}^{*})-\Omega_{t,t}^{*}\mathbf{v}_{t}^{*T}\widehat{\bm{S}}(\widehat{\bm{\Omega}}_{\cdot
k}-\bm{\Omega}_{\cdot k}^{*})|$ $\displaystyle\leq
n^{1/2}|(\widehat{\Omega}_{t,k}-\Omega_{t,k}^{*})(1-\Omega_{t,t}^{*}(\widehat{\bm{S}}_{t,t}-\mathbf{w}^{*T}_{t}\widehat{\bm{S}}_{-t,t}))|$
(A.4)
$\displaystyle+n^{1/2}\Omega_{t,t}^{*}|(\widehat{\bm{S}}_{t,-t}-\mathbf{w}^{*T}_{t}\widehat{\bm{S}}_{-t,-t})(\widehat{\bm{\Omega}}_{-t,k}-\bm{\Omega}^{*}_{-t,k})|.$
(A.5)

Note that
$1/\Omega_{t,t}^{*}=\bm{S}^{*}_{t,t}-\mathbf{w}^{*T}_{t}\bm{S}^{*}_{-t,t}$.
The term in (A.4) can be bounded by

$\displaystyle
n^{1/2}|(\widehat{\Omega}_{t,k}-\Omega_{t,k}^{*})\Omega_{t,t}^{*}(\widehat{\bm{S}}_{t,-t}-\bm{S}_{t,-t}^{*})|+n^{1/2}|(\widehat{\Omega}_{t,k}-\Omega_{t,k}^{*})\Omega_{t,t}^{*}\mathbf{w}_{t}^{*T}(\widehat{\bm{S}}_{-t,t}-\bm{S}_{-t,t}^{*})|$
$\displaystyle\leq n^{1/2}|\Omega_{t,t}^{*}|\|\widehat{\bm{\Omega}}_{\cdot
k}-\bm{\Omega}_{\cdot
k}^{*}\|_{1}\max(\|\widehat{\bm{S}}-\bm{S}^{*}\|_{\max},\|\mathbf{w}_{t}^{*T}(\widehat{\bm{S}}_{-t\cdot}-\bm{S}_{-t\cdot}^{*})\|_{\infty})$
$\displaystyle\leq\frac{Cs_{1}\log(K\vee n)}{n^{1/2}},$

with probability at least $1-(K\vee n)^{-3}$, by
$\lambda_{\max}(\bm{\Omega}^{*})\leq C$ and the concentration and error bound
results in Lemmas A.6, A.7, A.8. The term in (A.5) can be bounded by

$n^{1/2}\Omega^{*}_{t,t}\|\widehat{\bm{S}}_{t,-t}-\mathbf{w}^{*T}_{t}\widehat{\bm{S}}_{-t,-t}\|_{\infty}\|\widehat{\bm{\Omega}}_{-t,k}-\bm{\Omega}^{*}_{-t,k}\|_{1}\leq\frac{Cs_{1}\log(K\vee
n)}{n^{1/2}},$

with probability at least $1-(K\vee n)^{-3}$ again by Lemmas A.7, A.8. Thus,
$|\text{I.1}|\leq\frac{s_{1}\log(K\vee n)}{n^{1/2}}$ with probability at least
$1-(K\vee n)^{-3}$.

For term I.2, we have

$\displaystyle|I.2|$
$\displaystyle=n^{1/2}\Omega_{t,t}^{*}|(\widehat{\mathbf{v}}_{t}-\mathbf{v}^{*})^{T}(\widehat{\bm{S}}\widehat{\bm{\Omega}}_{\cdot
k}-\mathbf{e}_{k})|$ $\displaystyle\leq
n^{1/2}\Omega_{t,t}^{*}\|\widehat{\mathbf{v}}_{t}-\mathbf{v}^{*}\|_{1}\|\widehat{\bm{S}}\widehat{\bm{\Omega}}_{\cdot
k}-\mathbf{e}_{k}\|_{\infty}\leq\frac{Cs_{1}\log(K\vee n)}{n^{1/2}},$

with probability at least $1-(K\vee n)^{-3}$ by Lemma A.8 and the constraint
of the CLIME-type estimator.

To control term I.3, first we observe that

$\displaystyle|\widehat{h}(\widehat{\bm{\Omega}}_{\cdot k})|$
$\displaystyle=|\widehat{\mathbf{v}}_{t}^{T}(\widehat{\bm{S}}\widehat{\bm{\Omega}}_{\cdot
k}-\mathbf{e}_{k})|$
$\displaystyle\leq|\mathbf{v}_{t}^{*T}(\widehat{\bm{S}}\bm{\Omega}^{*}_{\cdot
k}-\mathbf{e}_{k})|+|\mathbf{v}_{t}^{*T}\widehat{\bm{S}}(\widehat{\bm{\Omega}}_{\cdot
k}-\bm{\Omega}^{*}_{\cdot
k})|+|(\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*})^{T}(\widehat{\bm{S}}\widehat{\bm{\Omega}}_{\cdot
k}-\mathbf{e}_{k})|$
$\displaystyle\leq\|\mathbf{v}^{*}_{t}\|_{1}\|\widehat{\bm{S}}\bm{\Omega}^{*}_{\cdot
k}-\mathbf{e}_{k}\|_{\infty}+\|\widehat{\bm{S}}_{t,-t}-\mathbf{w}_{t}^{*T}\widehat{\bm{S}}_{-t,-t}\|_{\infty}\|\widehat{\bm{\Omega}}_{\cdot
k}-\bm{\Omega}^{*}_{\cdot
k}\|_{1}+\|\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*}\|_{1}\|\widehat{\bm{S}}\widehat{\bm{\Omega}}_{\cdot
k}-\mathbf{e}_{k}\|_{\infty}.$

As shown in the proof of Lemma A.9, $\|\mathbf{v}^{*}_{t}\|_{1}\leq
s_{1}^{1/2}\|\mathbf{v}^{*}_{t}\|_{2}\leq Cs_{1}^{1/2}$. The rest of the
bounds on the above terms follows easily from Lemmas A.7, A.8. Thus, we have
$|\widehat{h}(\widehat{\bm{\Omega}}_{\cdot k})|\leq C(s_{1}\log(K\vee
n)/n)^{1/2}$ with high probability. Since

$|\widehat{\Omega}_{t,t}-\Omega_{t,t}^{*}|\leq C(s_{1}\log(K\vee n)/n)^{1/2},$

by Lemma A.8, we obtain that $|\text{I.3}|\leq\frac{s_{1}\log(K\vee
n)}{n^{1/2}}$ with probability at least $1-(K\vee n)^{-3}$. It is easily seen
that $\Omega_{t,t}^{*}\geq\frac{1}{S^{*}_{t,t}}\geq C>0$, see Remark A.4. This
implies that $s^{2}_{t,k}=\Omega_{t,t}^{*}\Omega_{k,k}^{*}+\Omega_{t,k}^{*2}$
is lower bounded by a positive constant. The proof of (A.3) is complete.

In step 2, we need to verify that

$\max_{1\leq t<k\leq
K}\sup_{x\in\mathbb{R}}\Big{|}\mathbb{P}\Big{(}\frac{n^{1/2}\Omega_{t,t}^{*}h(\bm{\Omega}_{\cdot
k}^{*})}{s_{tk}}\leq x\Big{)}-\Phi(x)\Big{|}\leq\frac{C}{(K\vee
n)^{3}}+\frac{1}{n^{1/2}},$

which has been done in Lemma A.9. Thus, combining with result (A.3), we can
use the same simple union bound in the proof of Theorem 4.3 to obtain the
desired result. ∎

######  Proposition A.3 (Convergence Rate of the Variance Estimator).

Under the same conditions as in Theorem 4.3, we have

$\max_{1\leq t<k\leq K}|\widehat{s}_{t,k}^{2}-s_{t,k}^{2}|\leq
C\sqrt{\frac{s_{1}\log(K\vee n)}{n}},~{}~{}\max_{1\leq t<k\leq
K}\Big{|}\frac{\widehat{s}_{t,k}}{s_{t,k}}-1\Big{|}\leq
C\sqrt{\frac{s_{1}\log(K\vee n)}{n}},$

with probability at least $1-(K\vee n)^{-3}$.

###### Proof.

By Lemma A.8, under the event $\widehat{G}=G^{*}$, we have

$\max_{1\leq t,k\leq
K}|\widehat{\Omega}_{t,k}-\Omega^{*}_{t,k}|\leq\max_{1\leq k\leq
K}\|\widehat{\bm{\Omega}}_{\cdot k}-\bm{\Omega}^{*}_{\cdot k}\|_{2}\leq
C_{1}\sqrt{\frac{s_{1}\log(K\vee n)}{n}},$

with probability at least $1-\frac{C_{4}}{(K\vee n)^{3}}$. Under this event,

$\displaystyle\max_{1\leq t<k\leq K}|\widehat{s}_{t,k}^{2}-s_{t,k}^{2}|$
$\displaystyle=\max_{1\leq t<k\leq
K}|\widehat{\Omega}_{t,k}^{2}+\widehat{\Omega}_{t,t}\widehat{\Omega}_{k,k}-(\Omega_{t,k}^{*2}+\Omega_{t,t}^{*}\Omega_{k,k}^{*})|$
$\displaystyle\leq\max_{1\leq t<k\leq
K}|(\widehat{\Omega}_{t,k}-\Omega_{t,k}^{*})(\widehat{\Omega}_{t,k}+\Omega_{t,k}^{*})|+\widehat{\Omega}_{t,t}|\widehat{\Omega}_{k,k}-\Omega_{k,k}^{*}|+\Omega_{k,k}^{*}|\widehat{\Omega}_{t,t}-\Omega_{t,t}^{*}|$
$\displaystyle\leq C_{1}\sqrt{\frac{s_{1}\log(K\vee
n)}{n}}(4\|\bm{\Omega}^{*}\|_{\max}+\delta)\leq C\sqrt{\frac{s_{1}\log(K\vee
n)}{n}},$

for some constant $\delta>0$ since
$\|\bm{\Omega}^{*}\|_{\max}\leq\lambda_{\max}(\bm{\Omega}^{*})\leq C$. It is
easily seen that $\Omega_{t,t}^{*}\geq\frac{1}{S^{*}_{t,t}}\geq C>0$, see
Remark A.4. This implies that
$s^{2}_{t,k}=\Omega_{t,t}^{*}\Omega_{k,k}^{*}+\Omega_{t,k}^{*2}$ is lower
bounded by a positive constant. Thus,

$\max_{1\leq t<k\leq
K}\Big{|}\frac{\widehat{s}_{t,k}}{s_{t,k}}-1\Big{|}\leq\max_{1\leq t<k\leq
K}\Big{|}\frac{\widehat{s}^{2}_{t,k}-s^{2}_{t,k}}{s_{t,k}(s_{t,k}+\widehat{s}_{t,k})}\Big{|}\leq
C\sqrt{\frac{s_{1}\log(K\vee n)}{n}}.$

∎

###  A.2 Key Lemmas for Estimators of the Cluster-Average Graph

######  Remark A.4.

While Assumptions 4.1 and 4.2 are made for $\mathbf{C}^{*}$, they do imply
that $c_{1}\leq\lambda_{\min}\left(\bm{S}^{*}\right)$ and
$\max_{t}S_{t,t}^{*}\leq c_{2}+c_{3}$ holds for $\bm{S}^{*}$. Furthermore
Lemma D.3 implies the same restricted eigenvalue condition on
$\lambda_{\min}\left(\bm{S}^{*}\right)$ as on $\mathbf{C}^{*}$.

######  Remark A.5.

In the following proof, we always assume the event
$\mathcal{E}=\\{\widehat{G}=G^{*}\\}$ holds. Using a similar argument to that
used in the proof of Theorem 4.3, the following bounds will hold with
probability at least $1-\frac{C}{(K\vee n)^{3}}$.

######  Lemma A.6 (Consistency of $\widehat{\bm{S}}$).

If Assumptions 4.1 and 4.2 hold, then with probability greater than
$1-\frac{C}{(K\vee n)^{3}}$,

$\|\widehat{\bm{S}}-\bm{S}^{*}\|_{\max}\leq C\sqrt{\frac{\log(K\vee n)}{n}}$

for some constant $C$ dependent only on $c_{1}$, $c_{2}$, and $c_{3}$ from
Assumptions 4.1 and 4.2.

###### Proof of Lemma A.6.

The proof follows from the proof of Theorem 1 in Cai et al. (2011). ∎

######  Lemma A.7 (Concentration of Gradient and Hessian).

If Assumptions 4.1 and 4.2 hold, then with probability greater than
$1-\frac{C}{(K\vee n)^{3}}$, we have that

  * (a)

$\max_{1\leq k\leq K}\|\widehat{\bm{S}}\bm{\Omega}_{\cdot
k}^{*}-\mathbf{e}_{k}\|_{\infty}\leq C_{1}\sqrt{\frac{\log(K\vee n)}{n}}$,

  * (b)

$\max_{1\leq k\leq
K}\|\widehat{\bm{S}}_{t,-t}-\mathbf{w}_{t}^{*T}\widehat{\bm{S}}_{-t,-t}\|_{\infty}\leq
C_{2}\sqrt{\frac{\log(K\vee n)}{n}}$, and

  * (c)

$\max_{1\leq k\leq
K}\|\mathbf{w}_{t}^{*T}(\widehat{\bm{S}}_{-t,-t}-\bm{S}_{-t,-t}^{*})\|_{\infty}\leq
C_{2}\sqrt{\frac{\log(K\vee n)}{n}}$,

for some constants $C_{1}$ and $C_{2}$ dependent only on $c_{1}$, $c_{2}$, and
$c_{3}$ from Assumptions 4.1 and 4.2.

###### Proof of Lemma A.7.

We start from the decomposition

$\displaystyle\widehat{\bm{S}}-\bm{S}^{*}$
$\displaystyle=(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Sigma}}-\bm{\Sigma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}.$

Denoting by $\mathbf{B}^{*}=\mathbf{A}^{*T}\mathbf{A}^{*}$, we can write

$\displaystyle(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Sigma}}-\bm{\Sigma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbf{B}^{*-1}\mathbf{A}^{*T}\Big{\\{}(\mathbf{A}^{*}\bm{Z}_{i}+\bm{E}_{i})(\mathbf{A}^{*}\bm{Z}_{i}+\bm{E}_{i})^{T}-\mathbf{A}^{*}\mathbf{C}^{*}\mathbf{A}^{*T}-\bm{\Gamma}^{*}\Big{\\}}\mathbf{A}^{*}\mathbf{B}^{*-1}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\bm{Z}_{i}\bm{Z}_{i}^{T}-\mathbf{C}^{*}+\frac{1}{n}\sum_{i=1}^{n}\bm{Z}_{i}\bm{E}_{i}^{T}\mathbf{A}^{*}\mathbf{B}^{*-1}+\frac{1}{n}\sum_{i=1}^{n}\mathbf{B}^{*-1}\mathbf{A}^{*T}\bm{E}_{i}\bm{Z}_{i}^{T}$
$\displaystyle~{}~{}~{}+\frac{1}{n}\sum_{i=1}^{n}\mathbf{B}^{*-1}\mathbf{A}^{*T}(\bm{E}_{i}\bm{E}_{i}^{T}-\bm{\Gamma}^{*})\mathbf{A}^{*}\mathbf{B}^{*-1}.$
(A.6)

Note that $\widehat{\bm{S}}\bm{\Omega}_{\cdot
k}^{*}-\mathbf{e}_{k}=(\widehat{\bm{S}}-\bm{S}^{*})\bm{\Omega}_{\cdot k}^{*}$,
and $||\bm{\Omega}_{\cdot
k}^{*}||_{2}\leq\lambda_{\max}\left(\bm{\Omega}^{*}\right)\leq\lambda_{\min}\left(\bm{S}^{*}\right)^{-1}$.
As seen in Remark A.4, we can bound the smallest eigenvalue of $\bm{S}^{*}$
from below by $\lambda_{\min}\left(\mathbf{C}^{*}\right)$. We therefore obtain
that $||\bm{\Omega}_{\cdot k}^{*}||_{2}\leq c_{1}+c_{3}$. By using the
triangle inequality, we can apply Lemma C.1, Lemma C.2 and Lemma C.3 to bound
(A.6). Combining these results, we obtain

$\max_{1\leq k\leq K}\|\widehat{\bm{S}}\bm{\Omega}_{\cdot
k}^{*}-\mathbf{e}_{k}\|_{\infty}\leq C_{1}\sqrt{\frac{\log(K\vee n)}{n}}$

with probability at least $1-\frac{C}{(K\vee n)^{3}}$, concluding the proof of
claim (a).

For the remaining two claims, we can rewrite

$\mathbf{w}_{t}^{*}=\left(S_{t,t}^{*}-\bm{S}^{*T}_{-t,t}(\bm{S}_{-t,-t}^{*})^{-1}\bm{S}_{-t,t}^{*}\right)\bm{\Omega}_{-t,t}^{*}=\frac{1}{\Omega_{t,t}^{*}}\bm{\Omega}_{-t,t}^{*}$

by the block matrix inverse formula. Using Lemma D.4, it follows that
$||\mathbf{w}^{*}_{t}||_{2}\leq\lambda_{\min}\left(\bm{\Omega}^{*}\right)\max_{t}S_{t,t}^{*}$.
Then we can see that

$\displaystyle\max_{1\leq k\leq
K}||\widehat{\bm{S}}_{t,-t}-\mathbf{w}_{t}^{*T}\widehat{\bm{S}}_{-t,-t}||_{\infty}$
$\displaystyle=\max_{1\leq k\leq
K}||(\widehat{\bm{S}}_{t,-t}-\bm{S}_{t,-t}^{*})-\mathbf{w}_{t}^{*T}(\widehat{\bm{S}}_{-t,-t}-\bm{S}_{-t,-t}^{*})||_{\infty}$
$\displaystyle\leq\underbrace{\max_{1\leq k\leq
K}||(\widehat{\bm{S}}_{t,-t}-\bm{S}_{t,-t}^{*})||}_{\text{(i)}}+\underbrace{\max_{i\neq
t}|\mathbf{w}_{t}^{*T}(\widehat{\bm{S}}_{-t,i}-\bm{S}_{-t,i}^{*})|}_{\text{(ii)}}.$

By using Lemma A.6, (i) is bounded with high probability. Similar to the the
proof of Theorem 1 in Cai et al. (2011), we can show that $(ii)\leq
C_{2}\sqrt{\frac{\log(K\vee n)}{n}}$, with probability at least
$1-\frac{C}{(K\vee n)^{3}}$. Part (c) is the same as the term (ii), concluding
the proof. ∎

######  Lemma A.8 (Consistency of Initial Estimator).

If Assumptions 4.1 and 4.2 hold, then,

  * (a)

$\max_{1\leq k\leq K}\|\widehat{\bm{\Omega}}_{\cdot k}-\bm{\Omega}_{\cdot
k}^{*}\|_{1}\leq C_{1}s_{1}\sqrt{\frac{\log(K\vee n)}{n}}$, $\max_{1\leq k\leq
K}\|\widehat{\bm{\Omega}}_{\cdot k}-\bm{\Omega}_{\cdot k}^{*}\|_{2}\leq
C_{1}\sqrt{\frac{s_{1}\log(K\vee n)}{n}}$,

  * (b)

$\max_{1\leq t\leq K}\|\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*}\|_{1}\leq
C_{2}s_{1}\sqrt{\frac{\log(K\vee n)}{n}}$, and

  * (c)

$\max_{1\leq k\leq t\leq
K}|(\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*})^{T}\widehat{\bm{S}}(\widehat{\bm{\Omega}}_{\cdot
k}-\bm{\Omega}_{\cdot k}^{*})|\leq C_{3}\frac{s_{1}\log(K\vee n)}{n}$,

with probability at least $1-\frac{C_{4}}{(K\vee n)^{3}}$. $C_{1}$, $C_{2}$,
$C_{3}$, $C_{4}$ and $C_{5}$ are constants, dependent only upon $c_{0}$,
$c_{1}$, $c_{2}$, and $c_{3}$ from Assumptions 4.1 and 4.2.

###### Proof of Lemma A.8.

The proof follows from the same argument as in the proof of Lemma B.5 with
Lemma B.4 replaced by Lemma A.7.

∎

######  Lemma A.9 (CLT for the Pseudo-Score Function).

Recall that $s_{tk}^{2}=\Omega_{t,k}^{*2}+\Omega^{*}_{t,t}\Omega^{*}_{k,k}$.
Let $F_{n}$ denote the CDF of
$n^{1/2}\mathbf{v}_{t}^{*T}(\widehat{\bm{S}}\bm{\Omega}_{\cdot
k}^{*}-\mathbf{e}_{k})/(s_{tk}/\Omega^{*}_{t,t})$. If Assumptions 4.1 and 4.2
hold, then

$\max_{1\leq t<k\leq K}\sup_{x\in\mathbb{R}}|F_{n}(x)-\Phi(x)|\leq
C(n^{-1/2}+(d\vee n)^{-3}),$

where $C$ is a constant dependent only upon $c_{0}$, $c_{1}$, and $c_{2}$.

###### Proof of Lemma A.9.

We have the similar bound that
$F_{n}(x)-\Phi(x)\leq\widetilde{F}_{n}(x)-\Phi(x)+\mathbb{P}(\bar{\mathcal{E}})$,
where $\mathcal{E}$ is the event $\widehat{G}=G^{*}$, and
$\widetilde{F}_{n}(x)$ is the CDF of
$n^{-1/2}\sum_{i=1}^{n}\mathbf{v}_{t}^{*T}(\bar{\bm{X}}_{i}\bar{\bm{X}}_{i}^{T}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k})/(s_{tk}/\Theta^{*}_{tt})$. Lemma A.1 implies that
$\mathbb{P}(\bar{\mathcal{E}})\leq(d\vee n)^{-3}$. We can similarly verify the
Lyapunov condition to control $\widetilde{F}_{n}(x)-\Phi(x)$. Finally,
$\mathbb{E}[\mathbf{v}_{t}^{*T}(\bar{\bm{X}}_{i}\bar{\bm{X}}_{i}^{T}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k})^{2}]=(s_{tk}/\Theta^{*}_{tt})^{2}$ by applying the
Isserlis’ theorem $\operatorname{{\rm
Var}}(\bm{v}_{1}^{T}\bar{\bm{X}}\bar{\bm{X}}^{T}\bm{v}_{2})=(\bm{v}_{1}^{T}\bm{S}^{*}\bm{v}_{1})(\bm{v}_{2}^{T}\bm{S}^{*}\bm{v}_{2})+(\bm{v}_{1}^{T}\bm{S}^{*}\bm{v}_{2})^{2}$,
for any vector $\bm{v}_{1}$ and $\bm{v}_{2}$. ∎

###### Proof of Proposition 4.5.

The proof follows verbatim that of Theorem 8.5 in Giraud (2014), with the
exception of the fact that exact $p$-values are replaced by approximate
$p$-values, including the rate of approximation. We included the full proof
for the convenience of the reader. By the definition of the FDR,

$\mathrm{FDR}(\widehat{\tau})=\mathbb{E}\Big{[}\frac{\sum_{(t,k)\in\mathcal{H}_{0}}\mathbb{I}[|\widetilde{W}_{t,k}|>\widehat{\tau}]\mathbb{I}[R_{\widehat{\tau}}>0]}{R_{\widehat{\tau}}}\Big{]}=\sum_{(t,k)\in\mathcal{H}_{0}}\mathbb{E}\Big{[}\frac{\mathbb{I}[|\widetilde{W}_{t,k}|>\widehat{\tau}]\mathbb{I}[R_{\widehat{\tau}}>0]}{R_{\widehat{\tau}}}\Big{]}.$
(A.7)

To handle the $R_{\widehat{\tau}}$ in the denominator, we use the identity
$1=\sum_{i=R_{\widehat{\tau}}}^{\infty}\frac{R_{\widehat{\tau}}}{i(i+1)}$.
This implies

$1/R_{\widehat{\tau}}=\sum_{i=R_{\widehat{\tau}}}^{\infty}\frac{1}{i(i+1)}=\sum_{i=1}^{\infty}\frac{\mathbb{I}[i\geq
R_{\widehat{\tau}}]}{i(i+1)}$

Plugging this into (A.7) and bringing the expectation inside the summation
gives that

$\displaystyle\mathrm{FDR}(\widehat{\tau})$
$\displaystyle=\sum_{(t,k)\in\mathcal{H}_{0}}\sum_{i=1}^{\infty}\frac{1}{i(i+1)}\mathbb{E}\Big{[}\mathbb{I}[|\widetilde{W}_{t,k}|>\widehat{\tau}]\mathbb{I}[R_{\widehat{\tau}}>0]\mathbb{I}[i\geq
R_{\widehat{\tau}}]\Big{]}$
$\displaystyle\leq\sum_{(t,k)\in\mathcal{H}_{0}}\sum_{i=1}^{\infty}\frac{1}{i(i+1)}\mathbb{E}\Big{[}\mathbb{I}[|\widetilde{W}_{t,k}|>\Phi^{-1}\left(1-\frac{\alpha
R_{\widehat{\tau}}}{2N_{BY}|\mathcal{H}|}\right)]\mathbb{I}[R_{\widehat{\tau}}>0]\mathbb{I}[i\geq
R_{\widehat{\tau}}]\Big{]}$
$\displaystyle\leq\sum_{(t,k)\in\mathcal{H}_{0}}\sum_{i=1}^{\infty}\frac{1}{i(i+1)}\mathbb{E}\Big{[}\mathbb{I}[|\widetilde{W}_{t,k}|>\Phi^{-1}\left(1-\frac{\alpha(i\wedge|\mathcal{H}|)}{2N_{BY}|\mathcal{H}|}\right)]\Big{]},$
(A.8)

where the second line follows from the definition of the FDR cutoff and the
last inequality holds since $R_{\widehat{\tau}}\leq(i\wedge|\mathcal{H}|)$.
The Berry-Esseen bound in Theorem 4.3 implies that

$\mathbb{P}(|\widetilde{W}_{t,k}|>\Phi^{-1}\left(1-\frac{\alpha(i\wedge|\mathcal{H}|)}{2N_{BY}|\mathcal{H}|}\right))\leq\frac{\alpha(i\wedge|\mathcal{H}|)}{N_{BY}|\mathcal{H}|}+2b_{n}.$

Thus, it follows that

$\displaystyle\mathrm{FDR}(\widehat{\tau})$
$\displaystyle\leq\sum_{(t,k)\in\mathcal{H}_{0}}\sum_{i=1}^{\infty}\frac{1}{i(i+1)}\Big{(}\frac{\alpha(i\wedge|\mathcal{H}|)}{N_{BY}|\mathcal{H}|}+2b_{n}\Big{)}$
$\displaystyle=\alpha\frac{|\mathcal{H}_{0}|}{|\mathcal{H}|}\Big{(}\sum_{i=1}^{|\mathcal{H}|}\frac{i}{i(i+1)N_{BY}}+\sum_{i=|\mathcal{H}|+1}^{\infty}\frac{|\mathcal{H}|}{i(i+1)N_{BY}}\Big{)}+2|\mathcal{H}_{0}|b_{n}$
$\displaystyle=\alpha\frac{|\mathcal{H}_{0}|}{|\mathcal{H}|}\Big{(}\sum_{i=1}^{|\mathcal{H}|}\frac{1}{i+1}\frac{1}{N_{BY}}+\frac{|\mathcal{H}|}{|\mathcal{H}|+1}\frac{1}{N_{BY}}\Big{)}+2|\mathcal{H}_{0}|b_{n}$
$\displaystyle\leq\alpha+2|\mathcal{H}_{0}|b_{n},$

where the last step follows from $\frac{|\mathcal{H}_{0}|}{|\mathcal{H}|}\leq
1$ and the definition of $N_{BY}$. This completes the proof. ∎

##  Appendix B Proofs Regarding Estimation in the Latent Variable Graph

This section contains the proofs that establish the asymptotic normality of
the estimator of $\bm{\Theta}^{*}$.

###  B.1 Main Proofs for Estimators of the Latent Variable Graph

###### Proof of Theorem 4.4.

The proof follows exactly the line of the proof of Theorem 4.3, but invokes
Proposition B.2, instead of Proposition A.3, and Proposition B.1, instead of
Proposition A.2, as one needs to establish different intermediate results,
specifically tailored to estimation of the latent graph. ∎

######  Proposition B.1 (Asymptotic Normality of $\widetilde{\Theta}_{t,k}$).

Under the same conditions as in Theorem 4.4, we get

$\max_{1\leq t<k\leq
K}\sup_{x\in\mathbb{R}}\Big{|}\mathbb{P}\left(\frac{n^{1/2}(\widetilde{\Theta}_{t,k}-\Theta_{t,k}^{*})}{\sigma_{t,k}}<x,\mathcal{E}\right)-\Phi(x)\Big{|}\leq\frac{C}{(K\vee
n)^{3}}+\frac{Cs_{0}\log(K\vee n)}{n^{1/2}}.$

###### Proof.

The proof follows all the steps of Proposition A.2, with Lemma A.6 replaced by
Lemma B.3, A.7 by B.4, A.8 by B.5 and A.9 by B.6.

∎

######  Proposition B.2 (Convergence Rate of the Variance Estimator).

Under the same conditions as in Theorem 4.4, we get

$\max_{1\leq t<k\leq K}|\widehat{\sigma}_{t,k}^{2}-\sigma_{t,k}^{2}|\leq
C\sqrt{\frac{s_{0}\log(K\vee n)}{n}}+\frac{Cs_{0}}{m},~{}~{}\max_{1\leq
t<k\leq K}\Big{|}\frac{\widehat{\sigma}_{t,k}}{\sigma_{t,k}}-1\Big{|}\leq
C\sqrt{\frac{s_{0}\log(K\vee n)}{n}}+\frac{Cs_{0}}{m},$

with probability at least $1-(K\vee n)^{-3}$.

###### Proof.

Similar to the proof of Proposition A.3, we can prove that

$\max_{1\leq t<k\leq
K}|\widehat{\sigma}_{t,k}^{2}-(\Theta_{t,k}^{*2}+\Theta^{*}_{t,t}\Theta^{*}_{k,k})|\leq
C\sqrt{\frac{s_{0}\log(K\vee n)}{n}},$

with probability at least $1-(K\vee n)^{-3}$. Then, by Lemma B.7, we obtain

$\max_{1\leq t<k\leq K}|\widehat{\sigma}_{t,k}^{2}-\sigma_{t,k}^{2}|\leq
C\sqrt{\frac{s_{0}\log(K\vee n)}{n}}+\frac{Cs_{0}}{m}.$

The second statement can be similarly derived. ∎

###  B.2 Key Lemmas for Estimators of the Latent Graph

######  Lemma B.3 (Consistency of $\widehat{\mathbf{C}}$).

If Assumptions 4.1 and 4.2 hold, then with probability greater than
$1-\frac{C}{(K\vee n)^{3}}$,

$\|\widehat{\mathbf{C}}-\mathbf{C}^{*}\|_{\max}\leq C\sqrt{\frac{\log(K\vee
n)}{n}},$

for some constant $C$ dependent only on $c_{1}$, $c_{2}$, and $c_{3}$ from
Assumptions 4.1 and 4.2.

###### Proof of Lemma B.3.

Let $\mathcal{E}$ denote the event that $\widehat{G}=G^{*}$. Under the event
$\mathcal{E}$, we have

$\widehat{\mathbf{C}}-\mathbf{C}^{*}=(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Sigma}}-\bm{\Sigma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}-(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Gamma}}-\bm{\Gamma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}.$
(B.1)

We can write, for $i\in G^{*}_{k}$,

$\displaystyle(\widehat{\bm{\Gamma}}-\bm{\Gamma}^{*})_{i,i}$
$\displaystyle=\widehat{\gamma}_{i}-\gamma_{i}^{*}$
$\displaystyle=\widehat{\Sigma}_{i,i}-\frac{1}{|G^{*}_{k}|-1}\sum_{j\in
G^{*}_{k},j\neq i}\widehat{\Sigma}_{i,j}-\gamma_{i}^{*}$
$\displaystyle=\widehat{\Sigma}_{i,i}-\frac{1}{|G^{*}_{k}|-1}\sum_{j\in
G^{*}_{k},j\neq
i}\widehat{\Sigma}_{i,j}-\Sigma_{i,i}^{*}+\frac{1}{|G^{*}_{k}|-1}\sum_{j\in
G^{*}_{k},j\neq i}\Sigma_{i,j}^{*}$
$\displaystyle=\widehat{\Sigma}_{i,i}-\Sigma_{i,i}^{*}-\frac{1}{|G^{*}_{k}|-1}\sum_{j\in
G^{*}_{k},j\neq i}\left[\widehat{\Sigma}_{i,j}-\Sigma_{i,j}^{*}\right],$

which implies that $||\widehat{\bm{\Gamma}}-\bm{\Gamma}^{*}||_{\max}\leq
2||\widehat{\bm{\Sigma}}-\bm{\Sigma}^{*}||_{\max}$. Therefore from Lemma D.1
we see that

$||(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Gamma}}-\bm{\Gamma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}||_{\max}\leq\frac{2}{m}||(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Sigma}}-\bm{\Sigma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}||_{\max},$

demonstrating that it will suffice to bound the first term in (B.1), which we
now do.

Let $\mathbf{B}^{*}=\mathbf{A}^{*T}\mathbf{A}^{*}$. For the first term in
(B.1), we have

$\displaystyle(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Sigma}}-\bm{\Sigma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbf{B}^{*-1}\mathbf{A}^{*T}\Big{\\{}(\mathbf{A}^{*}\bm{Z}_{i}+\bm{E}_{i})(\mathbf{A}^{*}\bm{Z}_{i}+\bm{E}_{i})^{T}-\mathbf{A}^{*}\mathbf{C}^{*}\mathbf{A}^{*T}-\bm{\Gamma}^{*}\Big{\\}}\mathbf{A}^{*}\mathbf{B}^{*-1}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\bm{Z}_{i}\bm{Z}_{i}^{T}-\mathbf{C}^{*}+\frac{1}{n}\sum_{i=1}^{n}\bm{Z}_{i}\bm{E}_{i}^{T}\mathbf{A}^{*}\mathbf{B}^{*-1}+\frac{1}{n}\sum_{i=1}^{n}\mathbf{B}^{*-1}\mathbf{A}^{*T}\bm{E}_{i}\bm{Z}_{i}^{T}$
$\displaystyle~{}~{}~{}+\frac{1}{n}\sum_{i=1}^{n}\mathbf{B}^{*-1}\mathbf{A}^{*T}(\bm{E}_{i}\bm{E}_{i}^{T}-\bm{\Gamma}^{*})\mathbf{A}^{*}\mathbf{B}^{*-1}.$
(B.2)

Using the triangle inequality, we can apply Lemma C.1, Lemma C.2 and Lemma C.3
to bound B.2. Combining these results with that $\mathbb{P}(\mathcal{E})\geq
1-c_{0}/(d\vee n)^{3}$, we obtain

$\|\widehat{\mathbf{C}}-\mathbf{C}^{*}\|_{\max}\leq C\sqrt{\frac{\log(K\vee
n)}{n}}$

with probability at least $1-\frac{C}{(K\vee n)}$ for some constant $C$
dependent only on $c_{1}$, $c_{2}$, and $c_{3}$ from Assumptions 4.1 and 4.2.
∎

######  Lemma B.4 (Concentration of Gradient and Hessian of the Loss
Function).

If Assumptions 4.1 and 4.2 hold, then with probability greater than
$1-\frac{C}{(K\vee n)^{3}}$, we have that

  * (a)

$\max_{1\leq k\leq K}\|\widehat{\mathbf{C}}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k}\|_{\infty}\leq C_{1}\sqrt{\frac{\log(K\vee n)}{n}}$,

  * (b)

$\max_{1\leq k\leq
K}\|\widehat{\mathbf{C}}_{t,-t}-\mathbf{w}_{t}^{*T}\widehat{\mathbf{C}}_{-t,-t}\|_{\infty}\leq
C_{2}\sqrt{\frac{\log(K\vee n)}{n}}$, and

  * (c)

$\max_{1\leq k\leq
K}\|\mathbf{w}_{t}^{*T}(\widehat{\mathbf{C}}_{-t,-t}-\mathbf{C}^{*}_{-t,-t})\|_{\infty}\leq
C_{2}\sqrt{\frac{\log(K\vee n)}{n}}$,

for absolute constants $C_{1}$ and $C_{2}$ dependent only on $c_{1}$, $c_{2}$,
and $c_{3}$ from Assumptions 4.1 and 4.2.

###### Proof of Lemma B.4.

Let $\mathcal{E}$ denote the event that $\widehat{G}=G^{*}$. Note that
$\widehat{\mathbf{C}}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k}=(\widehat{\mathbf{C}}-\mathbf{C}^{*})\bm{\Theta}^{*}_{\cdot
k}$. Under $\mathcal{E}$, following the decomposition (B.2), we can similarly
show that

$\displaystyle(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Sigma}}-\bm{\Sigma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\bm{\Theta}^{*}_{\cdot
k}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\bm{Z}_{i}\bm{Z}_{i}^{T}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{C}^{*}\bm{\Theta}^{*}_{\cdot
k}+\frac{1}{n}\sum_{i=1}^{n}\bm{Z}_{i}\bm{E}_{i}^{T}\mathbf{A}^{*}\mathbf{B}^{*-1}\bm{\Theta}^{*}_{\cdot
k}+\frac{1}{n}\sum_{i=1}^{n}\mathbf{B}^{*-1}\mathbf{A}^{*T}\bm{E}_{i}\bm{Z}_{i}^{T}\bm{\Theta}^{*}_{\cdot
k}$
$\displaystyle~{}~{}~{}+\frac{1}{n}\sum_{i=1}^{n}\mathbf{B}^{*-1}\mathbf{A}^{*T}(\bm{E}_{i}\bm{E}_{i}^{T}-\bm{\Gamma}^{*})\mathbf{A}^{*}\mathbf{B}^{*-1}\bm{\Theta}^{*}_{\cdot
k}.$ (B.3)

As in the proof of Lemma B.3, we have that

$||(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Gamma}}-\bm{\Gamma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\bm{\Theta}^{*}_{\cdot
k}||_{\infty}\leq\frac{2}{m}||(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\mathbf{A}^{*T}(\widehat{\bm{\Sigma}}-\bm{\Sigma}^{*})\mathbf{A}^{*}(\mathbf{A}^{*T}\mathbf{A}^{*})^{-1}\bm{\Theta}^{*}_{\cdot
k}||_{\infty},$

demonstrating that again it will suffice to bound the first term in B.3.

Note that $||\bm{\Theta}_{\cdot
k}^{*}||_{2}\leq\lambda_{\max}\left(\bm{\Theta}^{*}\right)\leq c_{1}^{-1}$.
Therefore, by using the triangle inequality, we can apply Lemma C.1, Lemma C.2
and Lemma C.3 to bound the first term in B.3. Combining these results and
$\mathbb{P}(\mathcal{E})\geq 1-C/(d\vee n)^{3}$, we obtain

$\max_{1\leq k\leq K}\|\widehat{\mathbf{C}}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k}\|_{\infty}\leq C_{1}\sqrt{\frac{\log(K\vee n)}{n}},$

with probability at least $1-\frac{C}{(K\vee n)^{3}}$ for some constant
$C_{1}$ dependent only on $c_{1}$, $c_{2}$, and $c_{3}$ from Assumptions 4.1
and 4.2.

For the remaining two claims, we can rewrite

$\mathbf{w}_{t}^{*}=\left(C_{t,t}^{*}-\mathbf{C}^{*T}_{-t,t}(\mathbf{C}_{-t,-t}^{*})^{-1}\mathbf{C}_{-t,t}^{*}\right)\bm{\Theta}_{-t,t}^{*}=\frac{1}{\Theta_{t,t}^{*}}\bm{\Theta}_{-t,t}^{*}$

by the block matrix inverse formula. Using Lemma D.4, it follows that
$||\mathbf{w}^{*}_{t}||_{2}\leq\lambda_{\max}\left(\bm{\Theta}^{*}\right)\max_{t}C_{t,t}^{*}$.
Then we see that

$\displaystyle\max_{1\leq k\leq
K}||\widehat{\mathbf{C}}_{t,-t}-\mathbf{w}_{t}^{*T}\widehat{\mathbf{C}}_{-t,-t}||_{\infty}$
$\displaystyle=\max_{1\leq k\leq
K}||(\widehat{\mathbf{C}}_{t,-t}-\mathbf{C}^{*}_{t,-t})-\mathbf{w}_{t}^{*T}(\widehat{\mathbf{C}}_{-t,-t}-\mathbf{C}^{*}_{-t,-t})||_{\infty}$
$\displaystyle\leq\underbrace{\max_{1\leq k\leq
K}||(\widehat{\mathbf{C}}_{t,-t}-\mathbf{C}^{*}_{t,-t})||}_{\text{(i)}}+\underbrace{\max_{i\neq
t}|\mathbf{w}_{t}^{*T}(\widehat{\mathbf{C}}_{-t,i}-\mathbf{C}^{*}_{-t,i})|}_{\text{(ii)}}$

Clearly, using Lemma B.3, (i) is bounded with high probability. Likewise,
Lemma B.3 demonstrates that
$\widehat{\mathbf{C}}_{-t,i}-\mathbf{C}^{*}_{-t,i}$ is a sub-exponential
random vector with parameters dependent only on $c_{1}$, $c_{2}$, and $c_{3}$
from Assumptions 4.1 and 4.2. Thus
$\mathbf{w}_{t}^{*T}(\widehat{\mathbf{C}}_{-t,i}-\mathbf{C}^{*}_{-t,i})$ is
sub-exponential and because
$||\mathbf{w}^{*}_{t}||_{2}\leq\lambda_{\max}\left(\bm{\Theta}^{*}\right)\max_{t}C_{t,t}^{*}$,
we obtain that

$\max_{1\leq k\leq
K}\|\widehat{\mathbf{C}}_{t,-t}-\mathbf{w}_{t}^{*T}\widehat{\mathbf{C}}_{-t,-t}\|_{\infty}\leq
C_{2}\sqrt{\frac{\log(K\vee n)}{n}}$

with probability at least $1-\frac{C}{(K\vee n)^{3}}$ for some constant
$C_{2}$. $C_{2}$ is dependent only on $c_{1}$, $c_{2}$, and $c_{3}$ from
Assumptions 4.1 and 4.2. The final result is bounded by the previous one,
concluding the proof. ∎

######  Lemma B.5 (Consistency of Initial Estimators).

If Assumptions 4.1 and 4.2 hold, then

  * (a)

$\max_{1\leq k\leq K}\|\widehat{\Theta}_{\cdot k}-\Theta_{\cdot
k}^{*}\|_{1}\leq C_{1}s_{0}\sqrt{\frac{\log(K\vee n)}{n}}$, $\max_{1\leq k\leq
K}\|\widehat{\Theta}_{\cdot k}-\Theta_{\cdot k}^{*}\|_{2}\leq
C_{1}\sqrt{\frac{s_{0}\log(K\vee n)}{n}}$,

  * (b)

$\max_{1\leq t\leq K}\|\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*}\|_{1}\leq
C_{2}s_{0}\sqrt{\frac{\log(K\vee n)}{n}}$, and

  * (c)

$\max_{1\leq k\leq t\leq
K}|(\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*})^{T}\widehat{\mathbf{C}}(\widehat{\bm{\Theta}}_{\cdot
k}-\bm{\Theta}_{\cdot k}^{*})|\leq C_{3}\frac{s_{0}\log(K\vee n)}{n}$,

with probability at least $1-\frac{C_{4}}{(K\vee n)^{3}}$. $C_{1}$, $C_{2}$,
$C_{3}$, $C_{4}$ are constants, dependent only upon $c_{0}$, $c_{1}$, $c_{2}$,
and $c_{3}$ from Assumptions 4.1 and 4.2.

###### Proof of Lemma B.5.

Below, the constants $C_{a}$, $C_{a}^{\prime}$, $C_{b}$, $C_{b}^{\prime}$,
$C_{b}^{\prime\prime}$, $C_{c}$ and $C_{c}^{\prime}$ will depend only upon
$c_{0}$, $c_{1}$, $c_{2}$, and $c_{3}$ from Assumptions 4.1 and 4.2.

We first prove part (b). The proof of part (a) is similar. Let
$\widehat{\bm{\Delta}}=\widehat{\mathbf{w}}_{t}-\mathbf{w}_{t}^{*}$, noting
that we can consider $\mathbf{w}_{t}$ instead of $\mathbf{v}_{t}$ as the
$t^{th}$ entries in both the estimated and true value are 1. By $S$ denote the
support of $\mathbf{w}_{t}^{*}$. $\mathbf{w}_{t}^{*}$ is $s_{0}$-sparse
because $\mathbf{w}_{t}^{*}$ is a multiple of $\bm{\Theta}_{-t,k}^{*}$, which
we know to be $s_{0}$-sparse.

By Lemma B.4, there exists $C_{b}$ such that for $\lambda\geq
C_{b}\sqrt{\frac{\log(K\vee n)}{n}}$, $\mathbf{w}_{t}^{*}$ is feasible for
(3.16) with probability at least $1-\frac{C_{b}^{\prime}}{(K\vee n)^{3}}$.
Assuming $\mathbf{w}_{t}^{*}$ is feasible, then it follows by definition that
$||(\mathbf{w}_{t}^{*})_{S}||_{1}\geq||(\widehat{\mathbf{w}}_{t})_{S}||_{1}+||(\widehat{\mathbf{w}}_{t})_{\bar{S}}||_{1}$.
This in turn implies by the triangle inequality that
$||\widehat{\bm{\Delta}}_{S}||_{1}\geq||\widehat{\bm{\Delta}}_{\bar{S}}||_{1}$.
Letting $\lambda=C_{b}\sqrt{\frac{\log(K\vee n)}{n}}$, it follows from the
triangle inequality that

$||\widehat{\mathbf{C}}_{-t,-t}\widehat{\bm{\Delta}}||_{\infty}\leq||\widehat{\mathbf{w}}_{t}^{T}\widehat{\mathbf{C}}_{-t,-t}-\mathbf{C}_{t,-t}\||_{\infty}+||\mathbf{w}_{t}^{*T}\widehat{\mathbf{C}}_{-t,-t}-\mathbf{C}_{t,-t}\||_{\infty}\leq
2C_{b}\sqrt{\frac{\log(K\vee n)}{n}}.$

In addition, note that $||\widehat{\bm{\Delta}}||_{1}\leq
2||\widehat{\bm{\Delta}}_{S}||_{1}\leq
2\sqrt{s_{0}}||\widehat{\bm{\Delta}}_{S}||_{2}\leq
2\sqrt{s_{0}}||\widehat{\bm{\Delta}}||_{2}$. Therefore combining with the
above, this gives

$\widehat{\bm{\Delta}}^{T}\widehat{\mathbf{C}}_{-t,-t}\bm{\Delta}\leq||\widehat{\bm{\Delta}}||_{1}||\widehat{\mathbf{C}}_{-t,-t}\widehat{\bm{\Delta}}||_{\infty}\leq
2C\sqrt{\frac{\log(K\vee n)}{n}}||\widehat{\bm{\Delta}}||_{1}\leq
4C_{b}\sqrt{\frac{s_{0}\log(K\vee n)}{n}}||\widehat{\bm{\Delta}}||_{2}.$

From Lemma D.3,
$\widehat{\bm{\Delta}}^{T}\widehat{\mathbf{C}}_{-t,-t}\widehat{\bm{\Delta}}\geq\frac{4c_{1}}{3}||\widehat{\bm{\Delta}}||_{2}^{2}$
with probability at least $1-\frac{C_{b}^{\prime\prime}}{(K\vee n)^{3}}$.
Therefore

$||\widehat{\bm{\Delta}}||_{2}\leq\frac{16Cc_{1}}{3}\sqrt{\frac{s_{0}\log(K\vee
n)}{n}}\quad\quad\text{ and
}\quad\quad||\widehat{\bm{\Delta}}||_{1}\leq\frac{32C_{b}s_{0}c_{1}}{3}\sqrt{\frac{\log(K\vee
n)}{n}},$

with probability at least
$1-\frac{\max\\{C_{b}^{\prime},C_{b}^{\prime\prime}\\}}{(K\vee n)^{3}}$.

To obtain part (c), first we apply Holder’s inequality and the triangle
inequality which give

$\displaystyle\max_{1\leq k\leq t\leq
K}|(\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*})^{T}\widehat{\mathbf{C}}(\widehat{\bm{\Theta}}_{\cdot
k}-\bm{\Theta}_{\cdot k}^{*})|$ $\displaystyle\leq\max_{1\leq k\leq t\leq
K}||\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*}||_{1}||\widehat{\mathbf{C}}(\widehat{\bm{\Theta}}_{\cdot
k}-\bm{\Theta}_{\cdot k}^{*})||_{\infty}$ $\displaystyle\leq\max_{1\leq k\leq
t\leq
K}||\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*}||_{1}\left(||\widehat{\mathbf{C}}\widehat{\bm{\Theta}}_{\cdot
k}-\bm{e}_{k}||_{\infty}+||\widehat{\mathbf{C}}\bm{\Theta}_{\cdot
k}^{*}-\bm{e}_{k}||_{\infty}\right).$ (B.4)

With choice of $\lambda$ as above, the KKT conditions give that
$||\widehat{\mathbf{C}}\widehat{\bm{\Theta}}_{\cdot
k}-\bm{e}_{k}||_{\infty}\leq C_{b}\sqrt{\frac{\log(K\vee n)}{n}}$. From Lemma
B.4, we have that $||\widehat{\mathbf{C}}\bm{\Theta}_{\cdot
k}^{*}-\bm{e}_{k}||_{\infty}\leq C\sqrt{\frac{\log(K\vee n)}{n}}$ with
probability at least $1-\frac{C}{(K\vee n)^{3}}$. Using part (b), we get that
$\max_{1\leq t\leq K}\|\widehat{\mathbf{v}}_{t}-\mathbf{v}_{t}^{*}\|_{1}\leq
Cs_{0}\sqrt{\frac{\log(K\vee n)}{n}}$ with probability at least
$1-\frac{C}{(K\vee n)^{3}}$. The desired result now follows from (B.4). ∎

######  Lemma B.6 (CLT for the Pseudo-Score Function).

Recall that
$\sigma^{2}_{tk}=\mathbb{E}(\Theta_{tt}^{*}\mathbf{v}_{t}^{*T}(\bar{\mathbf{C}}^{(i)}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k}))^{2}$ with $\bar{\mathbf{C}}^{(i)}$ defined in (3.24). Let
$F_{n}$ denote the CDF of
$n^{-1/2}\mathbf{v}_{t}^{*T}(\widehat{\mathbf{C}}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k})/(\sigma_{tk}/\Theta^{*}_{tt})$. If Assumptions 4.1 and 4.2
hold, then we have

$\max_{1\leq t<k\leq K}\sup_{x\in\mathbb{R}}|F_{n}(x)-\Phi(x)|\leq
C(n^{-1/2}+(d\vee n)^{-3}).$

where $C$ is a constant dependent only upon $c_{0}$, $c_{1}$, and $c_{2}$.

###### Proof of Lemma B.6.

Denote by $\mathcal{E}$ the event that $\widehat{G}=G^{*}$. We have

$F_{n}(x)-\Phi(x)\leq\widetilde{F}_{n}(x)-\Phi(x)+\mathbb{P}(\bar{\mathcal{E}}),$

where $\widetilde{F}_{n}(x)$ is the CDF of
$n^{-1/2}\sum_{i=1}^{n}\mathbf{v}_{t}^{*T}(\bar{\mathbf{C}}^{(i)}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k})/(\sigma_{tk}/\Theta^{*}_{tt})$. To control
$\widetilde{F}_{n}(x)-\Phi(x)$, we now verify the Lyapunov condition. As in
the proof of Lemma B.4, we can write

$\mathbf{v}_{t}^{*T}(\bar{\mathbf{C}}^{(i)}\bm{\Theta}^{*}_{\cdot
k}-\mathbf{e}_{k})=\mathbf{v}_{t}^{*T}(\bar{\mathbf{C}}^{(i)}-\mathbf{C}^{*})\bm{\Theta}^{*}_{\cdot
k}.$

From Lemmas C.1 \- C.3 we see that the entries in
$\bm{Q}_{i}=(\bar{\mathbf{C}}^{(i)}-\mathbf{C}^{*})\bm{\Theta}^{*}_{\cdot k}$
are sub-exponential with parameters $\alpha=C_{1}$ and $\nu=C_{2}$ which
depend only upon $\lambda_{\max}\left(\bm{\Theta}^{*}\right)$,
$\max_{k}\gamma_{k}^{*}$, and $\max_{t}C_{t,t}^{*}$.

Recall the definition of $\bm{v}_{t}^{*}$: $(\bm{v}_{t}^{*})_{t}=1$ and
$(\bm{v}_{t}^{*})_{-t}=-\mathbf{w}^{*}_{t}=-(\mathbf{C}_{-t,-t}^{*})^{-1}\mathbf{C}_{-t,t}^{*}$.
By the block matrix inverse formula, we can rewrite

$\mathbf{w}_{t}^{*}=\left(C_{t,t}^{*}-\mathbf{C}^{*T}_{-t,t}(\mathbf{C}_{-t,-t}^{*})^{-1}\mathbf{C}_{-t,t}^{*}\right)\bm{\Theta}_{-t,t}^{*}=\frac{1}{\Theta_{t,t}^{*}}\bm{\Theta}_{-t,t}^{*}.$

Using Lemma D.4, it follows that
$||\mathbf{w}^{*}_{t}||_{2}\leq\lambda_{\max}\left(\bm{\Theta}^{*}\right)\max_{t}C_{t,t}^{*}$
and
$||\mathbf{v}^{*}_{t}||_{2}\leq\lambda_{\max}\left(\bm{\Theta}^{*}\right)\max_{t}C_{t,t}^{*}+1$.
From Lemma E.3 and the above, $\mathbf{v}_{t}^{*T}\bm{Q}_{i}$ is sub-