# Nonparanormal Graph Quilting with Applications to Calcium Imaging

Andersen Chang∗1 Lili Zheng∗2 Gautam Dasarthy3 Genevera I. Allen2,4,5,6

(1Department of Neuroscience, Baylor College of Medicine
2Department of Electrical and Computer Engineering, Rice University
3School of Electrical, Computer and Energy Engineering, Arizona State
University
4Department of Computer Science, Rice University
5Department of Statistics, Rice University
6Jan and Dan Duncan Neurological Research Institute, Texas Children’s Hospital
)

###### Abstract

Probabilistic graphical models have become an important unsupervised learning
tool for detecting network structures for a variety of problems, including the
estimation of functional neuronal connectivity from two-photon calcium imaging
data. However, in the context of calcium imaging, technological limitations
only allow for partially overlapping layers of neurons in a brain region of
interest to be jointly recorded. In this case, graph estimation for the full
data requires inference for edge selection when many pairs of neurons have no
simultaneous observations. This leads to the Graph Quilting problem, which
seeks to estimate a graph in the presence of block-missingness in the
empirical covariance matrix. Solutions for the Graph Quilting problem have
previously been studied for Gaussian graphical models; however, neural
activity data from calcium imaging are often non-Gaussian, thereby requiring a
more flexible modeling approach. Thus, in our work, we study two approaches
for nonparanormal Graph Quilting based on the Gaussian copula graphical model,
namely a maximum likelihood procedure and a low-rank based framework. We
provide theoretical guarantees on edge recovery for the former approach under
similar conditions to those previously developed for the Gaussian setting, and
we investigate the empirical performance of both methods using simulations as
well as real data calcium imaging data. Our approaches yield more
scientifically meaningful functional connectivity estimates compared to
existing Gaussian graph quilting methods for this calcium imaging data set.

11footnotetext: Co-first authors on this manuscript.

Index terms— Graphical models; Graph quilting; Nonparanormal graphical models;
Rank-based correlation; Functional connectivity; Covariance completion

##  1 Introduction

Probabilistic graphical models are a popular unsupervised learning technique
for inference and sparse edge selection in network estimation, and are an
important tool for understanding dependency structures in high-dimensional
data. Graphical modeling approaches have been employed for data analysis in a
wide variety of areas, including neuroscience (Yatsenko et al.,, 2015;
Carrillo-Reid et al.,, 2021; Subramaniyan et al.,, 2018), genomics (Allen and
Liu,, 2013; Hartemink et al.,, 2000), and sensor networks (Chen et al.,, 2016;
Dasarathy et al.,, 2016). One particular research problem where graphical
models are applied is in the study of functional neuronal connectivity,
defined as statistical relationships between the activities of neurons in the
brain, from two-photon calcium imaging data (Horwitz,, 2003; Fingelkurts et
al.,, 2005). Functional neuronal connectivity is of particular interest in the
realm of neuroscience as a mechanism for describing how neuronal circuits in
the brain are organized and for finding patterns in neuronal activity that
underlies how information is passed between different regions of the brain
(Feldt et al.,, 2011); it may also serve well as a proxy for deriving the
structural synaptic connectivity between individual neurons in the brain, as
well as provide insight into the relationship between the two (Deco et al.,,
2014).

Due to technological limitations, in many calcium imaging experiments, the
full set of neurons in a brain region of interest are never simultaneously
observed. Instead, scans of functional activity are recorded in sequential,
partially overlapping layers, each containing only a subset of the population
of neurons (Grienberger and Konnerth,, 2012; Berens et al.,, 2017;
Pnevmatikakis et al.,, 2016). This data collection scheme leads to block-
missingness in the ensuing computation of the empirical covariance matrix for
the functional recording data on all neurons, as the joint activity of many
pairs are never observed; we demonstrate this visually in Figures 1(a) and
1(b). Therefore, in order to estimate a graphical model for functional
neuronal connectivity for the full set of observed neurons, the edge structure
in the missing portion must be inferred using the information from the
existing contemporaneous joint observations. The estimation of a graph in the
presence of block missing entries in the covariance matrix is known as the
Graph Quilting problem (Vinci et al.,, 2019).

Previously, several approaches have been developed in order to address the
Graph Quilting problem. (Vinci et al.,, 2019) originally proposed the Maximum
Determinant (MADGQ) algorithm, which first finds an $\ell_{1}$-regularized
maximum likelihood estimate of the graph on the observed portion of the
covariance matrix with the constraint that no edges are affiliated with
unobserved elements, then applies thresholding and Schur complements on the
result in order to identify graph edges and a minimal superset of edges.
Later, (Chang et al.,, 2022) introduced the Low Rank Graph Quilting (LRGQ)
approach, which utilizes a two-step procedure of covariance imputation under
the assumption of the existence of a low-rank representation of the underlying
covariance matrix, followed by the application of the graphical Lasso. Both of
the aforementioned Graph Quilting procedures have shown positive results for
graph imputation and edge recovery in the presence of block-wise missingness
in the entries of the full covariance matrices. However, their scope is
currently restricted to the case where the full underlying data follow a
multivariate Gaussian distribution. In particular, as the sample covariance is
the sufficient statistic for Gaussian graphical models, these methods and
their theory are both developed with a focus on accurately estimating or
imputing the covariance in the Gaussian setting. The empirical studies in
these papers for validating their methods also focused only on data fit to
standard Gaussian graphical models. Another related work, (Massa and Chiogna,,
2013), also proposes and empirically studies an approach of combining multiple
graphical models of subsets of variables to a full graph, which are in similar
spirit to the MADGQ algorithm in (Vinci et al.,, 2019), with a focus on the
Gaussian setting.

(a)

(b)

(c)

(d)

Figure 1: (a): An example of a typical schema which requires Graph Quilting.
Here, we have four subsets of features in a full feature set, observed across
different sessions; each rectangle represents the features and observations in
a particular block. (b): The corresponding incomplete empirical covariance
matrix for the same four patches of nodes from (a). The parts of the
covariance matrix not covered by any block are never jointly observed; graph
edges in this part of the covariance must be inferred from existing entries.
(c): Example of raw fluorescence trace for one neuron recorded in a calcium
imaging experiment, showing neuronal activity over time. (d): Resulting
empirical density of distribution of trace values.

Functional activity data from calcium imaging, on the other hand, tends to be
highly non-Gaussian. In particular, as shown in Figures 1(c) and 1(d), the
distribution of the fluorescence traces, which represent the firing activity
of each individual neuron from calcium imaging data, is heavily right skewed
with extreme positive outliers. To address this problem in the realm of
graphical models, approaches which assume that the data follow a parametric
distribution other than the Gaussian have been developed. For example, many
have explored elliptical distributions (Vogel and Fried,, 2011; Finegold and
Drton,, 2011) and general exponential family distributions (Yang et al.,,
2015, 2018) when the Gaussian assumption may not be appropriate. For the
specific case of estimating functional neuronal connectivity, the most common
alternative to a Gaussian-based method is the Poisson graphical model (Yang et
al.,, 2013), which assumes that the number of spikes of each neuron across
time bins follows a Poisson distribution (Xue et al.,, 2014; Vinci et al.,,
2018). The Poisson distribution is also the basis for a variety of other
approaches for estimating functional connectivity outside of the graphical
modeling paradigm, such as the linear-nonlinear model (Pillow et al.,, 2008;
Stevenson et al.,, 2008) and time series modeling with inter-spike intervals
(Masud and Borisyuk,, 2011). More recently, (Chang and Allen,, 2021) proposed
a different class of graphical models specifically for functional neuronal
connectivity based on the Subbotin distribution in order to capture
conditional dependencies between extreme values in the activity traces, which
directly represent neuronal activity. Nonlinear correlation methods, which
employ information theory metrics such as mutual information and joint
entropy, have also been used to derive functional neuronal connectivity from
calcium imaging data (Garofalo et al.,, 2009; Stetter et al.,, 2012).

While the aforementioned non-Gaussian functional neuronal connectivity models
have shown encouraging results for the calcium imaging application, they may
not be naturally suitable for the Graph Quilting problem. In order to leverage
existing ideas (Vinci et al.,, 2019; Chang et al.,, 2022) for solving this
problem, we still hope to operate on some pairwise similarity matrix like the
empirical covariance calculated from the observations rather than the raw data
itself. Therefore, we consider one particular alternative to the Gaussian
graphical model that can be applied to the Graph Quilting problem in the non-
Gaussian setting and which has been used for the analysis of functional
neuronal connectivity from calcium imaging data: the nonparanormal graphical
model (Liu et al.,, 2009; Dobra and Lenkoski,, 2011). This model assumes a
Gaussian copula for the joint distribution of all features while allowing
arbitrary marginal distributions for each feature. Although some prior methods
(Liu et al.,, 2009) estimate a univariate transform for each feature before
graph estimation (Zhao et al.,, 2012), it has also been shown that
nonparanormal graphical models can be estimated by applying the graphical
Lasso (Yuan and Lin,, 2007) on a transformed rank-based correlation matrix
(Liu et al.,, 2012; Harris and Drton,, 2013; He et al.,, 2017; Xue and Zou,,
2012); thus, this particular procedure could be utilized in the Graph Quilting
setting.

In this paper, we study two potential Graph Quilting techniques for
nonparanormal graphical models for graph estimation with non-Gaussian data
under block-missingness in the pairwise observation set. Our methods utilize
the MADGQ and LRGQ approaches of (Vinci et al.,, 2019) and (Chang et al.,,
2022), as applied to rank-based correlation matrices; we call the
nonparanormal adaptations MADGQ-NPN and LRGQ-NPN, respectively. While these
previous works have studied the nonparanormal graphical model using rank-based
correlation matrices, such as those mentioned above, ours is the first to
consider the performance of these methods in the presence of block-missingness
in the empirical covariance matrix. In particular, because the missingness
pattern is highly structured rather than random and the number of observed
entries is relatively sparse, the performance of nonparanormal Graph Quilting
will not follow directly from prior results on nonparanormal graphical models.
Therefore, we explore below the nonparanormal Graph Quilting approaches both
from a theoretical and empirical perspective in order to ascertain whether
they are appropriate for the Graph Quilting setting.

The rest of this paper is structured as follows. In Section 2, we describe the
MADGQ-NPN and LRGQ-NPN procedures for nonparanormal Graph Quilting. We then
show in Section 3 the conditions under which the MADGQ-NPN achieves exact
graph recovery for observed node pairs and minimal superset recovery for
unobserved node pairs. In Section 4, we present simulation studies which
compare the performances of the MADGQ-NPN procedure and one of the LRGQ-NPN
algorithms on data from non-Gaussian parametric distributions. Lastly, in
Section 5, we investigate the efficacy of the nonparanormal Graph Quilting
methods for estimating functional neuronal connectivity networks on real-world
calcium imaging data sets in comparison to each other as well as to the Graph
Quilting methods with Gaussian assumptions.

##  2 Nonparanormal Graph Quilting

###  2.1 Problem Set-up

We consider the nonparanormal graphical model (Liu et al.,, 2009), where each
sample vector $X_{i}\in\mathbb{R}^{p}$ follows a nonparanormal distribution
$\mathrm{NPN}_{p}(f,\Sigma)$, formally defined in Definition 1.

######  Definition 1.

Let $f=(f_{1},\dots,f_{p})$ be an ordered list of monotone univariate
functions and let $\Sigma\in\mathbb{R}^{p\times p}$ be a positive-definite
correlation matrix. Then a random vector $X=(X_{1},\dots,X_{p})^{\top}$
follows the nonparanormal distribution $\mathrm{NPN}_{p}(f,\Sigma)$ if
$(f_{1}(X_{1}),\dots,f_{p}(X_{p}))^{\top}\sim\mathcal{N}(0,\Sigma)$.

Define the precision matrix of the latent Gaussian vector
$(f_{1}(X_{1}),\dots,f_{p}(X_{p}))$ as $\Theta=\Sigma^{-1}$; from this, we
consider the following the graph structure that encodes the conditional
dependence relationship among $X$: $\mathcal{G}=(V,E)$,
$V=[p]$,$E=\\{(j,k):j,k\in[p],\Theta_{j,k}\neq 0\\}$. The primary interest is
to recover the unknown graph structure, or equivalently, the non-zero patterns
in the off-diagonal elements of $\Theta$. When one observes i.i.d. samples
$X_{1},\dots,X_{n}$ from this nonparanormal model, rank-based methods (Liu et
al.,, 2012; Harris and Drton,, 2013) have been proposed to learn the graph
structure with selection consistency guarantees.

In the Graph Quilting setting, however, we do not have access to all features
for each sample. Instead, we observe $K$ partially overlapping blocks
$V_{1},\dots,V_{K}$, each of size $p_{k}=|V_{k}|<p$. The corresponding
observed data matrix is denoted by $X^{(k)}\in\mathbb{R}^{n_{k}\times p_{k}}$
where $n_{k}$ is the sample size for block $k$, and our goal is to learn the
graph structure of all $p$ features from $\\{X^{(k)}\\}_{k=1}^{K}$. We define
the jointly observed feature pairs as $O=\\{(i,j):\exists 1\leq k\leq K,i,j\in
V_{k}\\}\subset[p]\times[p]$, and let $O^{c}=[p]\times[p]\backslash O$ denote
the feature pairs that has no joint measurement. Thus, we need to infer the
graph structure of the missing portion $O^{c}$ of the corresponding covariance
matrix from the known pairwise observations in $O$.

###  2.2 Nonparanormal Graph Quilting Methods

In this section, we extend two prior graph quilting approaches, the MADGQ
(MAximum Determinant graph-quilting) approach (Vinci et al.,, 2019) and the
LRGQ (Low-rank Graph Quilting) approach (Chang et al.,, 2022), from the
Gaussian graphical model setting to the nonparanormal setting. We first note
that both approaches take an estimate of the covariance/correlation
$\Sigma_{i,j}$ of the Gaussian variable pairs $(i,j)$ in
$O\subset[p]\times[p]$ as the input. In the nonparanormal setting, however, we
need the covariance of the latent Gaussian variables
$f_{1}(X_{1}),\dots,f_{p}(X_{p})$, which can be estimated using rank-based
correlations. For our work, we consider two rank-based correlations:
Spearman’s rho and Kendall’s tau. For each block $1\leq k\leq K$, let
$r^{(k)}_{i,j}$ be the rank of $X^{(k)}_{i,j}$ among
$X^{(k)}_{1,j},\dots,X^{(k)}_{n_{k},j}$. Also, let
$\bar{r}^{(k)}_{j}=\frac{1}{n_{k}}\sum_{i=1}^{n_{k}}r^{(k)}_{i,j}=\frac{n_{k}+1}{2}$.
From these, we compute Spearman’s rho
$\widehat{\rho}^{(k)}\in\mathbb{R}^{p_{k}\times p_{k}}$ and Kendall’s tau
$\widehat{\tau}^{(k)}\in\mathbb{R}^{p_{k}\times p_{k}}$ correlations as
follows:

$\begin{split}\widehat{\rho}^{(k)}_{j,l}&=\frac{\sum_{i=1}^{n_{k}}(r^{(k)}_{i,j}-\bar{r}^{(k)}_{j})(r^{(k)}_{i,l}-\bar{r}^{(k)}_{l})}{\sqrt{\sum_{i=1}^{n_{k}}(r^{(k)}_{i,j}-\bar{r}^{(k)}_{j})^{2}\sum_{i=1}^{n_{k}}(r^{(k)}_{i,l}-\bar{r}^{(k)}_{l})^{2}}},\\\
\widehat{\tau}^{(k)}_{j,l}&=\frac{2}{n_{k}(n_{k}-1)}\sum_{1\leq
i<i^{\prime}\leq
n_{k}}\mathrm{sign}((r^{(k)}_{i,j}-r^{(k)}_{i^{\prime},j})(r^{(k)}_{i,l}-r^{(k)}_{i^{\prime},l})).\end{split}$
(1)

To obtain a $p\times p$ correlation matrix, we combine the $K$ rank
correlations $\widehat{\rho}^{(k)}$, $\widehat{\tau}^{(k)}$ together.
Specifically, for any index $j\in[p]$, if $j\in V_{k}$, let $j_{k}$ be its
corresponding index in $V_{k}$. That is, $(V_{k})_{j_{k}}=j$. Then we formally
define $\widehat{\rho},\widehat{\tau}\in\mathbb{R}^{p\times p}$ as follow: for
$(j,l)\in O,$

$\widehat{\rho}_{j,l}=\frac{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}\widehat{\rho}^{(k)}_{j_{k},l_{k}}}{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}},\quad\widehat{\tau}_{j,l}=\frac{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}\widehat{\tau}^{(k)}_{j_{k},l_{k}}}{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}};$ (2)

otherwise, $\widehat{\rho}_{j,l}=\widehat{\tau}_{j,l}=0$. As in (Liu et al.,,
2012), since both Spearman’s rho and Kendall’s tau rank correlations are
biased for the population correlation, we apply elementwise sine function
transformations to $\widehat{\rho}$ and $\widehat{\tau}$: for any
$(j,l)\in[p]\times[p]$,

$\widehat{\Sigma}^{(\rho)}_{j,l}=2\sin\left(\frac{\pi}{6}\widehat{\rho}_{j,l}\right),\quad\widehat{\Sigma}^{(\tau)}_{j,l}=\sin\left(\frac{\pi}{2}\widehat{\tau}_{j,l}\right).$
(3)

Both initial observed correlation matrices take the value of zero in $O^{c}$,
and $\widehat{\Sigma}^{(\rho)}_{O},\widehat{\Sigma}^{(\tau)}_{O}$ serve as
estimates for the population correlation $\Sigma_{O}$. In the following, we
describe the MADGQ-NPN and LRGQ-NPN approaches based on
$\widehat{\Sigma}^{(\rho)}_{O}$ or $\widehat{\Sigma}^{(\tau)}_{O}$.

#### MADGQ-NPN approach

This approach is based on a partially observed likelihood, following the MADGQ
method in (Vinci et al.,, 2019). In a nutshell, the algorithm consists of two
steps: (i) first estimating the edge set within $O$ by minimizing an
$\ell_{1}$-regularized log-likelihood loss under the constraint of having no
edges in $O^{c}$, followed by a thresholding step that eliminates the bias
effect of this artificial constraint; (ii) then using the Schur complement to
estimate a superset of the edges in $O^{c}$, based on potential graph
distortions in $O$ caused by edges in $O^{c}$. The full procedure is
summarized in Algorithm 1. Similar to the original MADGQ method, Algorithm 1
involves three tuning parameters: the $\ell_{1}$ regularization parameter
$\Lambda\in\mathbb{R}^{p\times p}$ in (4), and two thresholding parameters
$\tau_{1}$ and $\tau_{2}$ for obtaining edge sets in $O$ and $O^{c}$. For
$\Lambda$, one can either simply let $\Lambda_{j,l}=\lambda$ for all $1\leq
j,\,l\leq p$, or choose $\Lambda_{j,l}=C_{0}\sqrt{\frac{\log p}{n_{j,l}}}$
where $n_{j,l}$ is the joint sample size of node $j,\,l$ together. The
specific value of $\lambda$ or the scaling factor $C_{0}$ can then be chosen
based on the extended Bayesian information criterion (Foygel and Drton,, 2010;
Gao et al.,, 2012), which is commonly used in the graphical model literature
and computationally efficient. For the first thresholding parameter $\tau_{1}$
that defines the edge set $\widehat{E}_{O}\subset O$, we can use the stability
selection approach (Liu et al.,, 2010) by examining the stability of
$\widehat{E}_{O}$ when running step 1 and 2 on randomly subsampled data. For
the second thresholding parameter $\tau_{2}$, motivated by our theoretical
results presented in Section 3, we can let $\tau_{2}=c\sqrt{\frac{\log
p}{n_{k}}}$ for a small constant $c>0$. Throughout our empirical studies in
this paper, we use $c=0.05$ which turns out to work reasonably well in
different settings.

Under certain assumptions on the graph structure and edge strength, (Vinci et
al.,, 2019) shows that, in the Gaussian setting, the original MADGQ approach
with sample covariance as the input is guaranteed to recover the graph
structure in $O$ and find a minimum edge superset in $O^{c}$ (see Definition
2). As we will show in Section 3, such results can also be extended to our
Algorithm 1 in the nonparanormal setting.

Input: Block indices $\\{V_{k}\\}_{k=1}^{K}$, Data sets
$\\{X^{(k)}\in\mathbb{R}^{n_{k}\times p_{k}}\\}_{k=1}^{K}$, thresholding
parameters $\tau_{1},\,\tau_{2}>0$

  1. 1.

Compute the rank-based correlation matrix $\widehat{\Sigma}_{O}$ as
$\widehat{\Sigma}^{(\tau)}_{O}$ or $\widehat{\Sigma}^{(\rho)}_{O}$ in (3).

  2. 2.

Solve the MAD${}_{\text{GQlasso}}$-NPN problem

$\widehat{\widetilde{\Theta}}=\operatorname*{arg\,max}_{\Theta\succ
0,\Theta_{O^{c}}=0}\log\det\Theta-\langle\widehat{\Sigma}_{O},\Theta_{O}\rangle+\|\Lambda\circ\Theta\|_{1,\text{off}}.$
(4)

  3. 3.

Find the edge set in $O$: $\widehat{E}_{O}=\\{(i,j)\in O:i\neq
j,|\widehat{\widetilde{\Theta}}_{i,j}|>\tau_{1}\\}$

  4. 4.

For $k=1,\dots,K$, compute the Schur complement

$\widehat{\widetilde{\Theta}}^{(k)}=\widehat{\widetilde{\Theta}}_{V_{k},V_{k}}-\widehat{\widetilde{\Theta}}_{V_{k},V_{k}^{c}}\widehat{\widetilde{\Theta}}_{V_{k}^{c},V_{k}^{c}}^{-1}\widehat{\widetilde{\Theta}}_{V_{k}^{c},V_{k}}.$
(5)

  5. 5.

Find the node set

$\widehat{W}_{\tau_{1},\tau_{2}}=\\{i\in V:\forall k~{}s.t.~{}i\in
V_{k},~{}\exists j\neq
i,~{}\tau_{2}<|\widehat{\widetilde{\Theta}}^{(k)}_{i,j}|<\tau_{1}\\}$ (6)

  6. 6.

Obtain a superset of edges in $O^{c}$:
$\widehat{E}_{O^{c}}=O^{c}\cap(\widehat{W}_{\tau_{1},\tau_{2}}\times\widehat{W}_{\tau_{1},\tau_{2}})$.

Output: $\widehat{E}=\widehat{E}_{O}\cup\widehat{E}_{O^{c}}$

Algorithm 1 MADGQ-NPN Algorithm

#### LRGQ-NPN approach

Another approach we consider is extending the LRGQ method proposed by (Chang
et al.,, 2022), which makes use of the potential low-rankness in covariance
matrices. This approach is motivated by the commonly seen approximate low-
rankness in many neuroscience data sets and exhibits promising performance in
these applications. The LRGQ approach is a two-step procedure that first
completes the covariance estimate from the partially observed sample
covariance in $O$ using appropriate low-rank matrix completion methods, then
uses the completed covariance matrix as input to the graphical Lasso algorithm
to produce a graph estimate. For the imputation step, (Chang et al.,, 2022)
introduces three different imputation methods, including a block-wise SVD
method (BSVDgq), a nuclear norm minimization method (NNMgq), and a non-convex
gradient descent method for low-rank factorization (LRFgq). In the
nonparanormal setting, we propose to apply the aforementioned two-step
approach on a rank-based correlation matrix instead of the sample covariance;
we summarize the full procedure in Algorithm 2. Below, throughout our
empirical investigation for the LRGQ-NPN approach, we focus on the BSVDgq
approach in (Chang et al.,, 2022) for the imputation step. More details on its
implementation are included in the Appendix. Algorithm 2 also involves two
tuning parameters: the rank $r$ and the regularization parameter $\Lambda$.
Similar to (Chang et al.,, 2022) and other prior works on low-rank matrix
completion, one can choose the appropriate rank $r$ using the Bayesian
information criterion (BIC) (Burnham and Anderson,, 2004). While for
$\Lambda$, same as the MADGQ procedure, we can let each entry $\Lambda_{j,k}$
be the same or scale proportionally to $\sqrt{\frac{\log p}{n_{j,k}}}$. The
specific values of $\Lambda$ can then be chosen based on the stability
criterion (Liu et al.,, 2010).

As has been pointed out in many prior works on graphical models (Ravikumar et
al.,, 2011; Liu et al.,, 2012), the consistency of graphical Lasso hinges on
sufficiently accurate covariance estimation in terms of each entry, which
implies that the two-step procedure is guaranteed to give a consistent graph
selection provided that the imputation step leads to small
$\|\widehat{\Sigma}^{\mathrm{LR}}-\Sigma\|_{\infty}$. The BSVDgq method in
particular has been shown in (Chang et al.,, 2022) to achieve sufficient
imputation accuracy that eventually leads to graph selection consistency in
the Gaussian setting. However, such theoretical results are based on delicate
spectral analysis of sample covariance matrices, which is extremely
challenging to be extended to the non-linear rank-based correlations. We leave
the theoretical investigation for the LRGQ-NPN approach as future work and
instead focus on empirical validation here instead.

Input: Block indices $\\{V_{k}\\}_{k=1}^{K}$, Data sets
$\\{X^{(k)}\in\mathbb{R}^{n_{k}\times p_{k}}\\}_{k=1}^{K}$, rank of full
covariance matrix $r$.

  1. 1.

Compute the rank-based correlation matrix $\widehat{\Sigma}_{O}$ as
$\widehat{\Sigma}^{(\tau)}_{O}$ or $\widehat{\Sigma}^{(\rho)}_{O}$ in (3).

  2. 2.

Obtain imputed covariance matrix $\widehat{\Sigma}^{\mathrm{LR}}$ using low-
rank covariance completion methods.

  3. 3.

Apply the graphical Lasso to the imputed full covariance matrix
$\widehat{\Sigma}^{\mathrm{LR}}$ in order to obtain the estimated graph

$\widehat{\Theta}^{\mathrm{LR}}=\operatorname*{arg\,min}_{\Theta\succ
0}\,\log\det\Theta-\langle\widehat{\Sigma}^{\mathrm{LR}},\Theta\rangle+\|\Lambda\circ\Theta\|_{1,\text{off}}$

Output: $\widehat{E}=\\{(i,j):i\neq
j\in[p],\widehat{\Theta}^{\mathrm{LR}}_{i,j}\neq 0\\}$

Algorithm 2 LRGQ-NPN Algorithm

##  3 Edge Recovery of MADGQ-NPN

Although the idea of using rank-based correlations to substitute the sample
covariance / Pearson correlation matrix is straightforward, it is not clear if
this idea would succeed in the nonparanormal graph quilting setting. In this
section, we examine the theoretical properties of the MADGQ-NPN approach
(Algorithm 1), giving an affirmative answer to this question by showing
similar theoretical guarantees as those in (Vinci et al.,, 2019) established
for Gaussian data. Specifically, we will show that under similar assumptions,
Algorithm 1 exactly recovers the edge set in $O$ and a minimal superset of
edges in $O^{c}$. Recall the MAD${}_{\text{GQlasso}}$-NPN solution and its
Schur complements defined in Algorithm 1, here let’s first define their
population versions as follows:

$\widetilde{\Theta}=\operatorname*{arg\,max}_{\Theta\succ
0,\Theta_{O^{c}}=0}\log\det\Theta-\sum_{(i,j)\in O}\Theta_{i,j}\Sigma_{i,j},$
(7)

$\widetilde{\Sigma}=\widetilde{\Theta}^{-1}$, and for $k=1,\dots,K$,

$\widetilde{\Theta}^{(k)}=\widetilde{\Theta}_{V_{k},V_{k}}-\widetilde{\Theta}_{V_{k},V_{k}^{c}}\widetilde{\Theta}_{V_{k}^{c},V_{k}^{c}}^{-1}\widetilde{\Theta}_{V_{k}^{c},V_{k}}.$
(8)

(Vinci et al.,, 2019) has established nice theoretical properties for graph
selection if one has access to $\widetilde{\Theta}$ and
$\widetilde{\Theta}^{(k)}$. These results form the basis of our theory, and
much of our analysis is devoted to showing the proximity of our finite sample
solutions $\widehat{\widetilde{\Theta}}$ in (4) to its population counterpart
$\widetilde{\Theta}$, and $\widehat{\tilde{\Theta}}^{(k)}$ in (5) to
$\widetilde{\Theta}^{(k)}$. In the following, we will follow the terminology
and assumptions developed for the population theory in (Vinci et al.,, 2019).

As has been shown in (Vinci et al.,, 2019), hard thresholding
$\widetilde{\Theta}_{O}$ can lead to the graph recovery in $O$ under certain
assumptions on graph signals. While for the edges in $O^{c}$, the exact edge
set $E_{O^{c}}$ is not identifiable in the graph quilting setting. However, it
is possible to recover a superset of $E_{O^{c}}$ based on the distortion in
the Schur complements of $\widetilde{\Theta}_{V_{k},V_{k}}$ for each block
$k$, created by the out-of-block edges. Let the distortion created by the out-
of-block edges in the block $k$ be
$\delta^{(k)}_{i,j}=\Theta_{i,j}-\widetilde{\Theta}^{(k)}_{i,j}$. Given the
observed covariance $\Sigma_{O}$ and these distortion information, we define
the following minimal superset of $E_{O^{c}}$ as the smallest possible set
that include all possible edges which cannot be ruled out without further
information.

######  Definition 2 (Minimal Superset of $E_{O^{c}}$).

Let

$\mathcal{D}_{\rm off}(\Sigma,O)=\big{\\{}(i,j,k):\delta^{(k)}_{ij}\neq
0,i\neq j\big{\\}}$ (9)

be the set of known distortions over the off-diagonal elements of the Schur
complements $\tilde{\Theta}^{(1)},...,\tilde{\Theta}^{(K)}$, and let

$\mathcal{A}_{\rm off}(\Sigma,O):=\left\\{\Sigma^{\prime}\succ
0:~{}\Sigma^{\prime}_{O}=\Sigma_{O},~{}\mathcal{D}_{\rm
off}(\Sigma^{\prime},O)=\mathcal{D}_{\rm off}(\Sigma,O)\right\\}$ (10)

be the set of all positive definite covariance matrices that agree with the
observed $\Sigma_{O}$ and distortions $\mathcal{D}_{\rm off}(\Sigma,O)$. A set
$\mathcal{S}_{\rm off}\subseteq O^{c}$ is the minimal superset of $E_{O^{c}}$
with respect to $\Sigma_{O}$ and $\mathcal{D}_{\rm off}(\Sigma,O)$ if it
satisfies the following properties:

  1. 1.

$\forall\Sigma^{\prime}\in\mathcal{A}(\Sigma,O,Q)$ we have
$E_{O^{c}}^{\prime}\subseteq\mathcal{S}_{\rm off}$;

  2. 2.

$\forall\mathcal{S}^{\prime}\subsetneq\mathcal{S}_{\rm off}$,
$\exists\Sigma^{\prime}\in\mathcal{A}_{\rm off}(\Sigma,O)$ such that
$E^{\prime}_{O^{c}}\cap(\mathcal{S}_{\rm
off}\setminus\mathcal{S}^{\prime})\neq\emptyset$.

We also define the following quantities that will be useful in our theory. Let
$\nu:=\min_{(i,j)\in E_{O}}|\Theta_{i,j}|$ be the minimum signal in $E_{O}$,
$\delta:=\max_{(i,j)\in O,i\neq j}|\Theta_{i,j}-\widetilde{\Theta}_{i,j}|$ be
the maximum distortion induced by constraining $\widetilde{\Theta}_{O^{c}}=0$.
Also let
$\psi:=\min_{(i,j,k):0<|\tilde{\Theta}_{ij}^{(k)}|<\delta}\min\\{|\tilde{\Theta}_{i,j}^{(k)}|,\delta-|\tilde{\Theta}_{i,j}^{(k)}|\\}$,
$d=\max_{i}\|\Theta_{i,:}\|_{0}$,
$\widetilde{d}=\max_{i}\|\widetilde{\Theta}_{i,:}\|_{0}$,
$\widetilde{s}=\|\widetilde{\Theta}\|_{0,\mathrm{off}}$ be the number of non-
zero off-diagonal elements of $\widetilde{\Theta}$. Some other technical
quantities include
$\widetilde{\kappa}=\frac{\lambda_{\max}(\widetilde{\Theta})}{\lambda_{\min}(\widetilde{\Theta})}$,
$\kappa_{\widetilde{\Sigma}}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widetilde{\Sigma}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}=\max_{j}\sum_{k=1}^{p}|\widetilde{\Sigma}_{j,k}|$,
$\kappa_{\widetilde{\Gamma}}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\widetilde{\Gamma}_{S,S})^{-1}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}$,
where $S$ is the support set of $\widetilde{\Theta}$. Let $H_{i}=\\{j:(i,j)\in
O^{c}\\}$, $N_{H_{i}}(i)=N(i)\cap H_{i}$ be the neighborhood of $i$ in
$H_{i}$. We say that two nodes $i$ and $j$ are $V$-connected if there is a
path in $V$ connecting $i$ and $j$. We also require the following assumptions:

######  Assumption 1 (Weak distortion compared to signal).

We assume that the maximum off-diagonal distortion of the MADGQ solution is
smaller than half the signal strength in the original precision matrix:
$\delta<\frac{\nu}{2}$.

As has been proven in Theorem 3.1 in (Vinci et al.,, 2019), Assumption 1 can
be satisfied as long as $\|\Theta_{O^{c}}\|_{\infty}$ is bounded above by a
function of the edge weights within $\Theta_{O}$; this assumption is likely to
hold as long as the pairwise observation set includes all edges with strong
signal strength.

######  Assumption 2.

For every node $i\in V$ with $N_{H_{i}}(i)\neq\emptyset$, we have that for
every $k$ such that $i\in V_{k}$, there exists at least one node $j\in
V_{k}\setminus\\{i\\}$ that is $(H_{i}\cup\\{j\\})$-connected to some node in
$N_{H_{i}}(i)$.

Assumption 2 requires that if any node $i$ has an edge in $O^{c}$, then for
any block it belongs to, there exists a path that starts from $i$ and an edge
in $O^{c}$, and eventually back to this block. This assumption ensures that
any edge $(i,j)$ in $O^{c}$ would cause certain distortions in the blocks $i$
and $j$ belong to. Therefore, given the distortion set
$\mathcal{D}_{\mathrm{off}}(\Sigma,O)$, it is possible to identify the node
set with edges in $O^{c}$. However, we can’t directly compute the distortion
$\delta_{i,j}^{(k)}$, a discrepancy between the unknown $\Theta_{i,j}$ and
$\widetilde{\Theta}_{i,j}$. Instead, we can only assume that such discrepancy
leads to a small non-zero $\widetilde{\Theta}_{i,j}$. This motivates our
following assumption:

######  Assumption 3.

If $\delta_{i,\backslash i}^{(k)}\neq 0$, then there exists $j\neq i$ such
that $0<|\tilde{\Theta}_{ij}^{(k)}|<\delta$.

######  Assumption 4 (Incoherence condition).

Let $\Gamma=\widetilde{\Sigma}\otimes\widetilde{\Sigma}$,
$S=\\{(j,l):\widetilde{\Theta}_{j,l}\neq 0\\}$. We assume $\max_{e\in O\cap
S^{c}}\|\Gamma_{e,S}\Gamma_{S,S}^{-1}\|_{1}\leq 1-\alpha$ for some
$0<\alpha\leq 1$.

Similar incoherence conditions are commonly assumed in the literature of
graphical models (Ravikumar et al.,, 2011).

######  Assumption 5 (Sufficient block measurements).

The $K$ blocks cover all nodes: $\cup_{k=1}^{K}V_{k}=[p]$, and at least one
off-diagonal element: $|O|>p$.

This is a mild assumption on the block measurements, which is typically
satisfied by our motivating neuroscience applications.

######  Assumption 6 (Regularization parameter).

$\Lambda_{j,l}=\frac{C_{0}}{\alpha}\sqrt{\frac{\log p}{\min_{k}n_{k}}}$ for
all $(j,l)\in O$ and some universal constant $C_{0}>0$.

Assumptions 1-6 also appear in (Vinci et al.,, 2019) in the Gaussian graphical
model setting. Note that we do not require any additional assumptions when
extending the theoretical guarantees to the nonparanormal setting, while
accommodating for strictly weaker assumptions on the joint distribution.

######  Theorem 1.

Suppose that Assumptions 1-6 hold, and there exist at least one edge in the
graph encoded by $\Theta$ and $\widetilde{\Theta}$: $d,\widetilde{d}>2$. Then
we have the following guarantees for Algorithm 1 with probability at least
$1-\sum_{k}p_{k}^{-10}$:

  * •

_Exact recovery in $O$_. If

$n_{k}>\left[\frac{C_{0}}{4}\kappa_{\widetilde{\Gamma}}\left(1+\frac{8}{\alpha}\right)\left(\left(\frac{\nu}{2}-\delta\right)^{-1}+3\left(1+\frac{8}{\alpha}\right)(\kappa_{\widetilde{\Sigma}}+\kappa^{3}_{\widetilde{\Sigma}}\kappa_{\widetilde{\Gamma}})\widetilde{d}\right)\right]^{2}\log
p_{k},$ (11)

$\delta+\varepsilon_{1}\leq\tau_{1}<\nu-\delta-\varepsilon_{1}$, where
$\varepsilon_{1}=\frac{C_{0}}{4}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}$, then $\widehat{E}_{O}=E_{O}$.

  * •

_Minimal superset recovery in $O^{c}$_. If

$n_{k}>C_{0}\kappa_{\widetilde{\Gamma}}^{2}\left(1+\frac{8}{\alpha}\right)^{2}\left[\frac{9\widetilde{\kappa}^{4}}{4\psi^{2}}+\frac{1}{4\lambda_{\min}(\widetilde{\Theta})^{2}}\right]\min\\{p+\widetilde{s},\widetilde{d}^{2}\\}\log
p_{k},$ (12)

$\varepsilon_{2}\leq\tau_{2}<\psi-\varepsilon_{2}$,
$\delta-\varepsilon_{2}<\tau_{1}\leq\nu-\varepsilon_{2}$, where
$\varepsilon_{2}=\frac{3C_{0}}{4}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\widetilde{\kappa}^{2}\min\\{\sqrt{p+\widetilde{s}},\widetilde{d}\\}\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}$, then $\widehat{E}_{O^{c}}=\mathcal{S}_{\mathrm{off}}$.

The proof of Theorem 1 can be found in the Appendix. Theorem 1 suggests that
under the same population assumptions as the Gaussian graph quilting setting
in (Vinci et al.,, 2019), as long as the sample size for each block is
sufficiently large: $n_{k}=\Omega(\widetilde{d}^{2}\log p_{k})$, and the
thresholding parameters are appropriately chosen:
$\tau_{1}\in[\delta+C\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}},\nu-\delta-C\max_{k}\sqrt{\frac{\log p_{k}}{n_{k}}})$,
$\tau_{2}\asymp\widetilde{d}\max_{k}\sqrt{\frac{\log p_{k}}{n_{k}}}$, we can
achieve exact recovery for edges in $O$ and construct a minimal superset for
edges in $O^{c}$. This result is comparable to the main theory in (Vinci et
al.,, 2019) for Gaussian graphical model, although we are considering a
strictly broader distribution family. To prove Theorem 1, we first extend the
existing error bounds for rank-based correlation matrices (Liu et al.,, 2012)
in the full data setting to our modified correlation estimates defined in 2.2,
which is computed from $K$ semi-overlapping blocks of measurements. We then
utilize this error bound to show the proximity of
$\widehat{\widetilde{\Theta}}$ and its Schur complements
$\widehat{\widetilde{\Theta}}^{(k)}$ to their population counterparts, which
eventually leads to Theorem 1 when combined with the population theory
developed in (Vinci et al.,, 2019).

##  4 Simulation Studies

(a)

(b)

(c)

(d)

Figure 2: Performance of MADGQ-NPN, BSVDgq-NPN, and zero imputation on
simulated data from a Gamma distribution. (a) TPR, changing block size. (b)
FDP, changing block size. (c) TPR, number of blocks. (d) FDP, changing number
of blocks.

We now investigate the empirical performance of the nonparanormal Graph
Quilting procedures outlined in Section 2 on simulated data. For each of the
simulation trials, we create synthetic block observation patterns by randomly
ordering the features and assigning them to $K$ partially overlapping blocks
of size $o$. Data is generated for each block from multivariate Gaussian data
from a sparse inverse covariance matrix with a small-world structure and then
perform copula transform on each column to a non-Gaussian distribution. Rank-
based correlations are calculated entry-wise for each pair of features in the
observation set, with correlations averaged for pairs observed in multiple
blocks. Our goal is to recover the nonzero entries in the population sparse
inverse covariance matrix of the unobserved Gaussian variable. Below, we
compare the MADGQ-NPN approach and one of the LRGQ-NPN procedures, namely the
block SVD (BSVDgq-NPN) approach, along with a basic zero imputation procedure
in which the missing entries of the correlation matrix are imputed as 0 before
applying the graphical Lasso algorithm. The methods are evaluated using the
the true positive rate (TPR) and false discovery proportion (FDP) of their
respective resulting graph estimates as compared to the true underlying graph;
in particular, we run 50 replications on each set of simulation parameters,
with new random block assignment in each replication, and show the average and
standard deviation of the TPR and FDP of edge selection for each method.
Hyperparameter selection is performed via optimal F-1 score tuning with
respect to the true graph.

We first test the nonparanormal Graph Quilting methods with Spearman
correlation matrices calculated on data containing 100 features generated from
a Gamma distribution with shape parameter 5 and scale parameter 1. For each
block, 2000 observations are generated per feature. Figures 2(a) and 2(b) show
the TPR and FDP of each method compared to the true underlying graph when the
size of the block $o$ is 52, 56, 60, 64, and 68 while the number of blocks is
held constant at 2, and Figures 2(c) and 2(d) show the TPR and FDP for $K=$ 2,
3, 4, 5 and 6 blocks while keeping the total number of observed node pairs
across all blocks ($K\times o$) constant at 120. From these results, we see
that both methods achieve a consistently high true positive rate for edge
recovery, and generally outperforms the zero imputation method. Notably, even
though the low-rank assumption inherent to the BSVDgq-NPN method is not met
here, the method is still able to recover the true edges of the graph at a
decently high rate. Comparing the two nonparanormal Graph Quilting methods, we
see that the MADGQ-NPN approach consistently has a higher TPR and FDP compared
to BSVDgq-NPN; this result follows what we expect, as the former approach is
designed to construct a superset of the true edges for the unobserved entries
in the covariance matrix and thus will return an estimate with both more true
positives and false positives compared to BSVDgq-NPN. Across different
simulation parameters, we generally observe that TPR is higher and FDP is
lower for the estimates from both methods when there are fewer blocks and when
each block is larger, which we would expect to see as these conditions
effectively provide more samples for estimation. Additionally, the difference
in performance between the two nonparanormal Graph Quilting methods appear to
be fairly consistent across the different block sizes and number of blocks.

We also evaluate the nonparanormal Graph Quilting methods on data simulated
from a Cauchy distribution with mean 0 and scale parameter 3 with 2000
observations and 100 features, from which a Kendall correlation matrix is
calculated and used as the input to the graphical Lasso. For this particular
experiment, we consider the case where the underlying covariance matrix is
approximately low-rank; this is generated via the spiked covariance model
(Johnstone,, 2001), and we also enforce a small-world graph structure in its
inverse. Figures 3(a) and 3(b) compare the MADGQ-NPN and BSVDgq-NPN methods in
terms of TPR and FDR for edge selection for $K=2$ blocks and varying block
size $o$ is 52, 56, 60, 64, and 68, and Figures 3(c) and 3(d) compare the same
methods for $K=$ 2, 3, 4, 5 and 6 blocks with a constant total number of
observations ($K\times o$) of 120. Both nonparanormal Graph Quilting methods
perform well in terms of true positive rate for edge recovery here as well,
and both considerably outperform the zero imputation approach. As opposed to
the previous simulation study in which the MADGQ-NPN method outperformed the
BSVDgq-NPN method in terms of selecting true edges in the graph, we see in
this case that the latter outperforms the former for both the TPR and FDP
metrics. The relative performance of the two nonparanormal Graph Quilting
methods matches what we would expect from this particular simulation setting,
as the structure of the full true underlying covariance matrix more closely
match the model set-up of the BSVDgq-NPN method, and also aligns with
comparative results seen from (Chang et al.,, 2022) comparing Graph Quilting
methods under a Gaussian assumption. Even in this case, though, the MADGQ-NPN
method still recovers the true edges of the underlying graph reasonably well.
These results show that the choice of which of the nonparanormal Graph
Quilting approach to apply for a problem will depend whether a low-rank
assumption makes sense in the particular scientific context. Additionally, as
above, we observe that performance generally improves with larger blocks and
with a fewer total number of blocks.

(a)

(b)

(c)

(d)

Figure 3: Performance of MADGQ-NPN, BSVDgq-NPN, and zero imputation on
simulated data from a Cauchy distribution with low-rank covariance. (a) TPR,
changing block size. (b) FDP, changing block size. (c) TPR, changing number of
blocks. (d) FDP, changing number of blocks.

##  5 Calcium Imaging Example

We now study the nonparanormal Graph Quilting procedures on a real-world
calcium imaging data set in order to assess the applicability of these methods
for estimating functional neuronal connectivity. The data come from the Allen
Institute (de Vries et al.,, 2016) and contains functional activity recordings
for 227 neurons in a mouse V1 cortex during spontaneous activity across
approximately 9000 time points. For this analysis, we compare the performances
of the MADGQ-NPN and BSVDgq-NPN nonparanormal Graph Quilting methods to each
other, as well as to a zero imputation approach, in a similar fashion to the
methodology in Section 4. We also compare the two nonparanormal Graph Quilting
methods to their analogous Gaussian-based procedures.

(a)

(b)

(c)

(d)

Figure 4: Performance of MADGQ-NPN, BSVDgq-NPN, and zero imputation on Allen
Institute data. (a) TPR, changing block size. (b) FDP, changing block size.
(c) TPR, changing number of blocks. (d) FDP, changing number of blocks.

To do the former, we measure the performances of the MADGQ-NPN, BSVDgq-NPN,
and zero imputation procedure by how well the graph estimates from each
nonparanormal Graph Quilting method on rank-based correlation matrices with
synthetic block-missingness recovers the edge structure of the graph estimated
using the nonparanormal graphical model on the same rank-based correlation
matrix with all pairwise entries observed. Specifically, we calculate the
Spearman correlation matrix for all observed neurons using binned spike counts
and apply the graphical Lasso to the fully observed covariance, the result of
which we consider as the true underlying graph structure. We create artificial
block-missingness in the empirical covariance matrix by randomly assigning
features to $K$ partially overlapping blocks of size $o$ and masking all
pairwise entries in the covariance which are not contained in any block. The
masked Spearman correlation matrix is then used as input for the nonparanormal
Graph Quilting methods. We show the average and standard deviation of the TPR
and FDR for recovering the graph estimate on the full Spearman covariance
matrix for each method across 50 replications on each set of parameters, with
new random block assignments each replication. Hyperparameter selection is
performed using optimal F-1 score with respect to the graph estimated from the
graphical Lasso fit on the fully observed data with the rank-based correlation
matrix.

Figures 4(a) and 4(b) show the TPR and FDP of each method compared to the true
underlying graph when the size of the block $o$ is 130, 140, 150, 160, and 170
while the number of blocks is held constant at 2, and Figures 4(c) and 4(d)
show the TPR and FDP for $K=$ 2, 3, 4, 5 and 6 blocks while keeping the total
number of observations across all blocks ($K\times o$) constant at 300. In
general, both nonparanormal Graph Quilting methods are able to recover most of
edges as the graph estimated with the nonparanormal graphical model when all
features are observed simultaneously, which shows that both methods can
reliably recover the edges in functional neuronal connectivity networks
derived from calcium imaging data that would be present if all neurons were
observed simultaneously. However, we also observe a relatively high FDP for
both methods, which seems to show that both methods tend to slightly
overselect the total number of edges in the underlying full graph. The BSVDgq-
NPN method outperforms MADGQ-NPN for edge recovery in terms of both TPR and
FDP; this can likely be attributed to the approximate low-rank structure often
found in the empirical covariance matrices from calcium imaging data (Stringer
et al.,, 2019). Also, as seen in Section 4, all methods perform better with
fewer blocks and with larger block sizes.

(a)

(b)

(c)

(d)

Figure 5: Estimated functional connectivity networks from Graph Quilting
methods on Allen Institute data. Nodes are colored by visual angular tuning
category. (a) BSVDgq. (b) BSVDgq-NPN. (c) MADGQ. (d) MADGQ-NPN.

We then compare the functional neuronal connectivity graphs from the
nonparanormal Graph Quilting methods to those found using Graph Quilting with
a Gaussian assumption. Specifically, we use the MADGQ and BSVDgq Graph
Quilting methods to obtain graph estimates on the same data set with the same
artificial block-missingness pattern as above, and compare the selected edges
to functional connectivity graphs estimated by MADGQ-NPN and BSVDgq-NPN using
neural properties from provided metadata. We first assess each of the
estimated functional connectivity networks by the proportion of edges that
connect pairs of neurons with the same visual angular tuning category. In the
neuroscience literature, it has been hypothesized that neurons are tuned such
that they fire in the presence of specific stimuli (Sakia and Miyashita,,
1994) and that neurons with the similar tunings are more likely to be
functionally connected (Stevenson et al.,, 2012). Thus, we expect a sizable
proportion of the edges in the estimated functional connectivity networks to
link pairs of neurons within the same tuning category. For this particular
calcium imaging data set, neural angular tuning is categorized into 8
different bins, each comprising 45 degree intervals. Hyperparameter tuning for
the BSVDgq and BSVDgq-NPN methods is performed using the extended Bayesian
information criterion with respect to the original data; for the MADGQ and
MADGQ-NPN methods, in order to crate a fair comparison between methods, we
tune the thresholding parameters such that the number of edges is
approximately similar to the estimates from BSVDgq and BSVDgq-NPN.

Figures 5(a) and 5(b) show examples of estimated functional connectivity
networks from the BSVDgq and BSVDgq-NPN methods, and Figures 5(c) and 5(d)
show estimated functional connectivity networks from the MADGQ and MADGQ-NPN
methods, respectively. Structurally, we see that the graph estimates from the
nonparanormal Graph Quilting methods are much more likely to be comprised of
edges between neurons with the same angular tuning category. Specifically,
across different replications with synthetic block-missingness, the BSVDgq-NPN
graph estimates contain an average of 38.9% of edges that link neurons with
the same tuning category, compared to 21.4% of edges in the BSVDgq estimates.
Similarly, we see that 46.8% of edges in the MADGQ-NPN graph estimates are
between pairs of neurons in the same tuning bin, as opposed to just 28.7% in
the MADGQ graph estimates.

We also compare the recorded firing activity of one particular example neuron
and its selected edge neighbors in the the functional connectivity graphs
estimated from the BSVDgq and BSVDgq-NPN methods in terms of how closely the
neural firing patterns match one another. Specifically, for the problem of
functional neuronal connectivity, our goal is to find neurons with consistent
synchronous firing activity across time, which is represented by
contemporaneous large positive spikes in the fluorescence traces (Smetters et
al.,, 1996; Turaga et al.,, 2013). In Figures 6(a) and 6(b), we visualize the
fluorescence trace of the selected neuron, overlayed with the fluorescence
traces neurons which are edge neighbors unique to the BSVDgq graph; we also do
the same in Figures 6(c) and 6(d) for neurons that are edge neighbors unique
to the BSVDgq-NPN graph. The top 10 periods of spiking activity of the example
selected neuron are represented via blue dotted lines in the plots. From the
results, we see that the firing activity of the selected neuron seems to match
relatively closely with the edge neighbors selected only in the graph estimate
from the BSVDgq-NPN algorithm. On the other hand, the edge neighbors selected
only by the BSVDgq graph do not appear to have the same firing pattern.
Quantitatively, the top 10 firing times for the selected neuron and its edge
neighbors match 72.4% of the time for the BSVDgq-NPN functional connectivity
graph estimates, as opposed to just 24.5% for the BSVDgq functional
connectivity graphs. Overall, from this real-world calcium imaging data study,
we see that nonparanormal Graph Quilting provides more logical functional
connectivity estimates in the nueroscience context compared to the ordinary
Graph Quilting procedures.

(a)

(b)

(c)

(d)

Figure 6: Fluorescence traces of one particular neuron (in grey) and edge
neighbors (in red) from functional connectivity graphs estimated via BSVDgq
((a, b)) and BSVDgq-NPN ((c, d)).

##  6 Discussion

In this work, we have presented two potential approaches to nonparanormal
Graph Quilting, MADGQ-NPN and DSVDgq-NPN, which broaden the scope of Graph
Quilting procedures to be applicable in the nonparanormal graphical model
setting. We demonstrate theoretical properties of the MADGQ-NPN method,
showing criterion for exact edge recovery in the observed portion and minimum
superset recovery in the missing portion of the graph. Through our empirical
studies, we demonstrate that both nonparanormal Graph Quilting methods can be
effective for edge selection for non-Gaussian data depending on the structure
of the underlying covariance matrix. Through our real-world calcium imaging
data example, we show that the nonparanormal Graph Quilting methods can be
used to recover the same functional neuronal connectivity network edges from
calcium imaging data in the presence of non-simultaneous observations for the
full population of neurons as would be found if all neurons are observed
concurrently, and that these methods can be applied to estimate more
appropriate functional neuronal connectivity networks from calcium imaging
data compared to Gaussian methods.

There are many potential directions for future research that can be taken from
our work. While we have characterized the theoretical performance of the
MADGQ-NPN approach, we do not currently have any guarantees for the LRGQ-NPN
methods. Methodologically, extensions to the nonparanormal Graph Quilting
procedure to account for potential other data effects such as latent
variables, autocorrelations, or covariates could improve graph estimation and
edge selection accuracy. Also, the nonparanormal graphical models could be
applied to research problems in other fields where joint observations may be
missing, such as RNA-seq in genomics and signal processing in power systems.
In conclusion, our work has helped to extend graph inference for nonparanormal
graphical models in the presence of block-missingness in the observed
covariance matrix, with theoretical guarantees for performance and promising
empirical results for calcium imaging data.

## Acknowledgements

The authors gratefully acknowledge support by NSF NeuroNex-1707400, NIH
1R01GM140468, and NSF DMS-2210837.

## References

  * Allen and Liu, (2013) Allen, G. and Liu, Z. (2013).  A local poisson graphical model for inferring networks from sequencing data.  IEEE transactions on nanobioscience, 12(3):189–198.
  * Berens et al., (2017) Berens, P., Theis, L., Stone, J., Sofroniew, N., Tolias, A., Bethge, M., and Freeman, J. (2017).  Standardizing and benchmarking data analysis for calcium imaging.  Computational and Systems Neuroscience Meeting (COSYNE 2017), pages 66–67.
  * Burnham and Anderson, (2004) Burnham, K. P. and Anderson, D. R. (2004).  Multimodel inference: understanding aic and bic in model selection.  Sociological methods & research, 33(2):261–304.
  * Carrillo-Reid et al., (2021) Carrillo-Reid, L., Han, S., O’Neil, D., Taralova, E., Jebara, T., and Yuste, R. (2021).  Identification of pattern completion neurons in neuronal ensembles using probabilistic graphical models.  Journal of Neuroscience, 41(41):8577–8588.
  * Chang and Allen, (2021) Chang, A. and Allen, G. (2021).  Subbotin graphical models for extreme value dependencies with applications to functional neuronal connectivity.  arXiv preprint arXiv:2106.11554.
  * Chang et al., (2022) Chang, A., Zheng, L., and Allen, G. I. (2022).  Low-rank covariance completion for graph quilting with applications to functional connectivity.  arXiv preprint arXiv:2209.08273.
  * Chen et al., (2016) Chen, L., Wainwright, M., Cetin, M., and Willsky, A. (2016).  Data association based on optimization in graphical models with application to sensor networks.  Mathematical and computer modelling, 43(9-10):1114–1135.
  * Dasarathy et al., (2016) Dasarathy, G., Singh, A., Balcan, M., and Park, J. (2016).  Active learning algorithms for graphical model selection.  Artificial Intelligence and Statistics, pages 1356–1364.
  * de Vries et al., (2016) de Vries, S., Lecoq, J., Buice, M., Groblewski, P., Ocker, G., Oliver, M., Feng, D., Cain, N., Ledochowitsch, P., Millman, D., and Roll, K. (2016).  Dataset: Allen institute mindscope program.  Allen Institute. Dataset., https://observatory.brain-map.org/visualcoding/search/cell_list?experiment_container_id=511510989&sort_field=p_sg&sort_dir=asc.
  * Deco et al., (2014) Deco, G., McIntosh, A., Shen, K., Hutchison, R., Menon, R., Everling, S., Hagmann, P., and Jirsa, V. (2014).  Identification of optimal structural connectivity using functional connectivity and neural modeling.  Journal of Neuroscience, 34(23):7910–7916.
  * Dobra and Lenkoski, (2011) Dobra, A. and Lenkoski, A. (2011).  Copula Gaussian graphical models and their application to modeling functional disability data.  The Annals of Applied Statistics, 5(2A):969 – 993.
  * Feldt et al., (2011) Feldt, S., Bonifazi, P., and Cossart, R. (2011).  Dissecting functional connectivity of neuronal microcircuits: experimental and theoretical insights.  Trends in neurosciences, 34(5):225–236.
  * Finegold and Drton, (2011) Finegold, M. and Drton, M. (2011).  Robust graphical modeling of gene networks using classical and alternative t-distributions.  The Annals of Applied Statistics, pages 1057–1080.
  * Fingelkurts et al., (2005) Fingelkurts, A., Fingelkurts, A., and Kähkönen, S. (2005).  Functional connectivity in the brain—is it an elusive concept?  Neuroscience & Biobehavioral Reviews, 28(8):827–836.
  * Foygel and Drton, (2010) Foygel, R. and Drton, M. (2010).  Extended bayesian information criteria for gaussian graphical models.  Advances in neural information processing systems, 23.
  * Gao et al., (2012) Gao, X., Pu, D. Q., Wu, Y., and Xu, H. (2012).  Tuning parameter selection for penalized likelihood estimation of gaussian graphical model.  Statistica Sinica, pages 1123–1146.
  * Garofalo et al., (2009) Garofalo, M., Nieus, T., Massobrio, P., and Martinoia, S. (2009).  Evaluation of the performance of information theory-based methods and cross-correlation to estimate the functional connectivity in cortical networks.  PloS one, 4(8):e6482.
  * Grienberger and Konnerth, (2012) Grienberger, C. and Konnerth, A. (2012).  Imaging calcium in neurons.  Neuron, 73(5):862–885.
  * Harris and Drton, (2013) Harris, N. and Drton, M. (2013).  Pc algorithm for nonparanormal graphical models.  Journal of Machine Learning Research, 14(11):3365–3383.
  * Hartemink et al., (2000) Hartemink, A., Gifford, D., Jaakkola, T., and Young, R. (2000).  Using graphical models and genomic expression data to statistically validate models of genetic regulatory networks.  Biocomputing, pages 422–433.
  * He et al., (2017) He, Y., Zhang, X., Wang, P., and Zhang, L. (2017).  High dimensional gaussian copula graphical model with fdr control.  Computational Statistics & Data Analysis, 113:457–474.
  * Horwitz, (2003) Horwitz, B. (2003).  The elusive concept of brain connectivity.  Neuroimage, 19(2):466–470.
  * Johnstone, (2001) Johnstone, I. (2001).  On the distribution of the largest eigenvalue in principal components analysis.  The Annals of statistics, 29(2):295–327.
  * 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(4):2293–2326.
  * Liu et al., (2009) Liu, H., Lafferty, J., and Wasserman, L. (2009).  The nonparanormal: Semiparametric estimation of high dimensional undirected graphs.  Journal of Machine Learning Research, 10(10):2295–2328.
  * Liu et al., (2010) Liu, H., Roeder, K., and Wasserman, L. (2010).  Stability approach to regularization selection (stars) for high dimensional graphical models.  Advances in neural information processing systems, 23.
  * Massa and Chiogna, (2013) Massa, M. S. and Chiogna, M. (2013).  Effectiveness of combinations of gaussian graphical models for model building.  Journal of Statistical Computation and Simulation, 83(9):1602–1612.
  * Masud and Borisyuk, (2011) Masud, M. and Borisyuk, R. (2011).  Statistical technique for analysing functional connectivity of multiple spike trains.  Journal of Neuroscience Methods, 196(1):201–219.
  * Pillow et al., (2008) Pillow, J., Shlens, J., Paninski, L., Sher, A., Litke, A., Chichilnisky, E., and Simoncelli, E. (2008).  Spatio-temporal correlations and visual signalling in a complete neuronal population.].  Nature, 454(7207):995–999.
  * Pnevmatikakis et al., (2016) Pnevmatikakis, E., Soudry, D., Gao, Y., Machado, T., Merel, J., Pfau, D., Reardon, T., Mu, Y., Lacefield, C., Yang, W., and Ahrens, M. (2016).  Simultaneous denoising, deconvolution, and demixing of calcium imaging data.  Neuron, 89(2):285–299.
  * Ravikumar et al., (2011) Ravikumar, P., Wainwright, M. J., Raskutti, G., and Yu, B. (2011).  High-dimensional covariance estimation by minimizing $\ell_{1}$-penalized log-determinant divergence.  Electronic Journal of Statistics, 5:935–980.
  * Sakia and Miyashita, (1994) Sakia, K. and Miyashita, Y. (1994).  Neuronal tuning to learned complex forms in vision.  Neuroreport, 5:829–832.
  * Smetters et al., (1996) Smetters, D., Majewska, A., and Yuste, R. (1996).  Detecting action potentials in neuronal populations with calcium imaging.  Methods, 18(2):215–221.
  * Stetter et al., (2012) Stetter, O., Battaglia, D., Soriano, J., and Geisel, T. (2012).  Model-free reconstruction of excitatory neuronal connectivity from calcium imaging signals.  PLoS Computational Biology, 8:e1002653.
  * Stevenson et al., (2012) Stevenson, I., London, B., Oby, E., Sachs, N., Reimer, J., Englitz, B., David, S., Shamma, S., Blanche, T., Mizuseki, K., and Zandvakili, A. (2012).  Neuronal tuning to learned complex forms in vision.  Neuroreport, 8(11):e1002775.
  * Stevenson et al., (2008) Stevenson, I., Rebesco, J., Miller, L., and Körding, K. (2008).  Inferring functional connections between neurons.  Current opinion in neurobiology, 18(6):582–588.
  * Stringer et al., (2019) Stringer, C., Pachitariu, M., Steinmetz, N., Reddy, C., Carandini, M., and Harris, K. (2019).  Spontaneous behaviors drive multidimensional, brainwide activity.  Science, 364(6437):eaav789.
  * Subramaniyan et al., (2018) Subramaniyan, M., Ecker, A., Patel, S., Cotton, R., Bethge, M., Pitkow, X., Berens, P., and Tolias, A. (2018).  Faster processing of moving compared with flashed bars in awake macaque v1 provides a neural correlate of the flash lag illusion.  Journal of Neurophysiology, 120(5):2430–2452.
  * Turaga et al., (2013) Turaga, S., Buesing, L., Packer, A., Dalgleish, H., Pettit, N., Hausser, M., and Macke, J. (2013).  Inferring neural population dynamics from multiple partial recordings of the same neural circuit.  Advances in Neural Information Processing Systems,, 26.
  * Vinci et al., (2019) Vinci, G., Dasarathy, G., and Allen, G. I. (2019).  Graph quilting: graphical model selection from partially observed covariances.  arXiv preprint arXiv:1912.05573.
  * Vinci et al., (2018) Vinci, G., Ventura, V., Smith, M., and Kass, R. (2018).  Adjusted regularization in latent graphical models: Application to multiple-neuron spike count data.  The annals of applied statistics, 12(2):1068–1095.
  * Vogel and Fried, (2011) Vogel, D. and Fried, R. (2011).  Elliptical graphical modelling.  Biometrika, 98(4):935–951.
  * 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(5):2541–2571.
  * Xue et al., (2014) Xue, W., Kang, J., Bowman, F., Wager, T., and Guo, J. (2014).  Identifying functional co‐activation patterns in neuroimaging studies via poisson graphical models.  Biometrics, 70(4):812–822.
  * Yang et al., (2015) Yang, E., Ravikumar, P., Allen, G. I., and Liu, Z. (2015).  Graphical models via univariate exponential family distributions.  The Journal of Machine Learning Research, 16(1):3813–3847.
  * Yang et al., (2013) Yang, E., Ravikumar, P. K., Allen, G. I., and Liu, Z. (2013).  On poisson graphical models.  Advances in neural information processing systems, 26.
  * Yang et al., (2018) Yang, Z., Ning, Y., and Liu, H. (2018).  On semiparametric exponential family graphical models.  The Journal of Machine Learning Research, 19(1):2314–2372.
  * Yatsenko et al., (2015) Yatsenko, D., Josić, K., Ecker, A., Froudarakis, E., Cotton, R., and Tolias, A. (2015).  Improved estimation and interpretation of correlations in neural circuits.  PLoS computational biology, 11(3):e1004083.
  * Yuan and Lin, (2007) Yuan, M. and Lin, Y. (2007).  Model selection and estimation in the gaussian graphical model.  Biometrika, 94(1):19–35.
  * Zhao et al., (2012) Zhao, T., Liu, H., Roeder, K., Lafferty, J., and Wasserman, L. (2012).  The huge package for high-dimensional undirected graph estimation in r.  Journal of Machine Learning Research, 14(11):3365–3383.

##  Appendix A Algorithmic Details

In our LRGQ-NPN approach, one key step is to impute the covariance matrix
$\widehat{\Sigma}^{\mathrm{LR}}$ from partially computed
$\widehat{\Sigma}_{O}=\widehat{\Sigma}_{O}^{(\tau)}$ or
$\widehat{\Sigma}_{O}^{(\rho)}$. We consider the following BSVDgq algorithm
(Algorithm 3) proposed in (Chang et al.,, 2022) for this step.

Input: $\\{V_{k},k\in 1,\ldots K\\}$,
$\widehat{\Sigma}_{O}\in\mathbb{R}^{p\times p}$, $r>0$.

Initialize: $\widetilde{H}=\boldsymbol{0}_{p\times r}$, $A=V_{1}$,
$\widehat{q}=\text{median}(\\{\widehat{\Sigma}_{ii},1\leq i\leq p\\})$.

Find low-rank solution for first patch:

  1. 1.

Calculate SVD of
$(\widehat{\Sigma}-\widehat{q}\boldsymbol{I})_{V_{1},\,V_{1}}=\mathbf{W}\boldsymbol{\Lambda}\mathbf{U}^{T}.$

  2. 2.

Set $\mathbf{h}$ as indices of the largest $r$ diagonal elements of
$\boldsymbol{\Lambda}$.

  3. 3.

Set
$\mathbf{C}_{V_{1},\,:}=\mathbf{U}_{:,\mathbf{h}}\boldsymbol{\Lambda}_{\mathbf{h},\mathbf{h}}^{1/2}.$

For $s\in 2,\ldots K$

  1. 1.

Find low-rank solution for $s$-th patch:

    1. (a)

Calculate SVD of
$(\widehat{\Sigma}-\widehat{q}\boldsymbol{I})_{V_{s},\,V_{s}}=\mathbf{W}\boldsymbol{\Lambda}\mathbf{U}^{T}.$

    2. (b)

Set $\mathbf{h}$ as indices of the largest $r$ diagonal elements of
$\boldsymbol{\Lambda}$.

    3. (c)

Calculate
$\mathbf{D}=\mathbf{U}_{:,\mathbf{h}}\boldsymbol{\Lambda}_{\mathbf{h},\mathbf{h}}^{1/2}.$

  2. 2.

Merge with previous patches:

    1. (a)

Find overlaps $E=\\{A:a\in V_{s}\\},\,J=\\{j:V_{s}[j]\in A\\}$.

    2. (b)

Calculate SVD of
$(\mathbf{D}_{J}^{T}\mathbf{C}_{E,\,:})=\mathbf{W}\boldsymbol{\Lambda}\mathbf{U}^{T}.$

    3. (c)

Set
$\mathbf{M}=\mathbf{C}_{E,\,:},\mathbf{C}_{V_{s},\,:}=\mathbf{D}\mathbf{W}\mathbf{U}^{T}$

    4. (d)

Set $\mathbf{C}_{E,\,:}=\mathbf{M}$

* 3.

Update $A=\bigcup_{k=1}^{s}V_{s}$

return _$\widetilde{\Sigma}=\mathbf{C}\mathbf{C}^{\top}.$_

Algorithm 3 Spiked Block Singular Value Decomposition (Spiked BSVDgq) (Chang
et al.,, 2022)

##  Appendix B Proofs

For completeness, here we restate our MADGQ-NPN algorithm, its theory and
assumptions presented in the main paper. The rank-based correlation is defined
as follows:

$\widehat{\Sigma}^{(\rho)}_{j,l}=2\sin\left(\frac{\pi}{6}\widehat{\rho}_{j,l}\right),\quad\widehat{\Sigma}^{(\tau)}_{j,l}=\sin\left(\frac{\pi}{2}\widehat{\tau}_{j,l}\right),$
(13)

where

$\widehat{\rho}_{j,l}=\frac{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}\widehat{\rho}^{(k)}_{j_{k},l_{k}}}{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}},\quad\widehat{\tau}_{j,l}=\frac{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}\widehat{\tau}^{(k)}_{j_{k},l_{k}}}{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}},$ (14)

and

$\begin{split}\widehat{\rho}^{(k)}_{j,l}&=\frac{\sum_{i=1}^{n_{k}}(r^{(k)}_{i,j}-\bar{r}^{(k)}_{j})(r^{(k)}_{i,l}-\bar{r}^{(k)}_{l})}{\sqrt{\sum_{i=1}^{n_{k}}(r^{(k)}_{i,j}-\bar{r}^{(k)}_{j})^{2}\sum_{i=1}^{n_{k}}(r^{(k)}_{i,l}-\bar{r}^{(k)}_{l})^{2}}},\\\
\widehat{\tau}^{(k)}_{j,l}&=\frac{2}{n_{k}(n_{k}-1)}\sum_{1\leq
i<i^{\prime}\leq
n_{k}}\mathrm{sign}((r^{(k)}_{i,j}-r^{(k)}_{i^{\prime},j})(r^{(k)}_{i,l}-r^{(k)}_{i^{\prime},l})).\end{split}$
(15)

######  Assumption 1 (Weak distortion compared to signal).

We assume that the maximum off-diagonal distortion of the MADGQ solution is
smaller than half the signal strength in the original precision matrix:
$\delta<\frac{\nu}{2}$.

######  Assumption 2.

For every node $i\in V$ with $N_{H_{i}}(i)\neq\emptyset$, we have that for
every $k$ such that $i\in V_{k}$, there exists at least one node $j\in
V_{k}\setminus\\{i\\}$ that is $(H_{i}\cup\\{j\\})$-connected to some node in
$N_{H_{i}}(i)$.

######  Assumption 3.

If $\delta_{i,\backslash i}^{(k)}\neq 0$, then there exists $j\neq i$ such
that $0<|\tilde{\Theta}_{ij}^{(k)}|<\delta$.

######  Assumption 4 (Incoherence condition).

Let $\Gamma=\widetilde{\Sigma}\otimes\widetilde{\Sigma}$,
$S=\\{(j,l):\widetilde{\Theta}_{j,l}\neq 0\\}$. We assume $\max_{e\in O\cap
S^{c}}\|\Gamma_{e,S}\Gamma_{S,S}^{-1}\|_{1}\leq 1-\alpha$ for some
$0<\alpha\leq 1$.

######  Assumption 5 (Sufficient block measurements).

The $K$ blocks cover all nodes: $\cup_{k=1}^{K}V_{k}=[p]$, and at least one
off-diagonal element: $|O|>p$.

######  Assumption 6 (Regularization parameter).

$\Lambda_{j,l}=\frac{C_{0}}{\alpha}\sqrt{\frac{\log p}{\min_{k}n_{k}}}$ for
all $(j,l)\in O$ and some universal constant $C_{0}>0$.

######  Theorem 1.

Suppose that Assumptions 1-6 hold, and there exist at least one edge in the
graph encoded by $\Theta$ and $\widetilde{\Theta}$: $d,\widetilde{d}>2$. Then
we have the following guarantees for Algorithm 1, with probability at least
$1-\sum_{k}p_{k}^{-10}$:

  * •

_Exact recovery in $O$_. If

$n_{k}>\left[\frac{C_{0}}{4}\kappa_{\widetilde{\Gamma}}\left(1+\frac{8}{\alpha}\right)\left(\left(\frac{\nu}{2}-\delta\right)^{-1}+3\left(1+\frac{8}{\alpha}\right)(\kappa_{\widetilde{\Sigma}}+\kappa^{3}_{\widetilde{\Sigma}}\kappa_{\widetilde{\Gamma}})\widetilde{d}\right)\right]^{2}\log
p_{k},$ (16)

$\delta+\varepsilon_{1}\leq\tau_{1}<\nu-\delta-\varepsilon_{1}$, where
$\varepsilon_{1}=\frac{C_{0}}{4}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}$, then $\widehat{E}_{O}=E_{O}$.

  * •

_Minimal superset recovery in $O^{c}$_. If

$n_{k}>C_{0}\kappa_{\widetilde{\Gamma}}^{2}\left(1+\frac{8}{\alpha}\right)^{2}\left[\frac{9\widetilde{\kappa}^{4}}{4\psi^{2}}+\frac{1}{4\lambda_{\min}(\widetilde{\Theta})^{2}}\right]\min\\{p+\widetilde{s},\widetilde{d}^{2}\\}\log
p_{k},$ (17)

$\varepsilon_{2}\leq\tau_{2}<\psi-\varepsilon_{2}$,
$\delta-\varepsilon_{2}<\tau_{1}\leq\nu-\varepsilon_{2}$, where
$\varepsilon_{2}=\frac{3C_{0}}{4}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\widetilde{\kappa}^{2}\min\\{\sqrt{p+\widetilde{s}},\widetilde{d}\\}\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}$, then $\widehat{E}_{O^{c}}=\mathcal{S}_{\mathrm{off}}$.

###### Proof of Theorem 1.

First of all, we note that a population version of Theorem 1 has been proved
in (Vinci et al.,, 2019) (see Theorem 3.1 and Theorem 3.4), which shows the
exact edge recovery in $O$ and minimal superset recovery in $O^{c}$ based on
thresholding $\widetilde{\Theta}$ and $\widetilde{\Theta}^{(k)}$. Our proof of
Theorem 1 mainly involves applying the population results in (Vinci et al.,,
2019) and some new estimation error bounds for
$\|\widehat{\widetilde{\Theta}}-\widetilde{\Theta}\|_{\infty}$ and
$\|\widehat{\widetilde{\Theta}}^{(k)}-\widetilde{\Theta}^{(k)}\|_{\infty}$.

In particular, when Assumptions 1-3 hold, (Vinci et al.,, 2019) shows that

$E_{O}=\\{(i,j)\in O:|\widetilde{\Theta}_{i,j}|>t\\},\quad\forall\delta\leq
t\leq\nu-\delta,$ (18)

and for any $\delta\leq t\leq\nu$,

$\mathcal{S}_{\mathrm{off}}=O^{c}\cap(W_{t}\times W_{t}),W_{t}=\big{\\{}i\in
V:~{}\forall k~{}s.t.~{}i\in V_{k},~{}\exists j\neq
i,0<|\tilde{\Theta}^{(k)}_{ij}|<t\big{\\}}.$ (19)

On the other hand, the following two theorems show the proximity between
$\widehat{\widetilde{\Theta}}$ and $\widetilde{\Theta}$,
$\widehat{\widetilde{\Theta}}^{(k)}$ and $\widetilde{\Theta}^{(k)}$:

######  Theorem 2.

Suppose that Assumptions 4 and 5 hold,
$\Lambda_{j,l}=\frac{C_{0}}{\alpha}\sqrt{\frac{\log p}{\min_{k}n_{k}}}$ for
all $(j,l)\in O$ and some universal constant $C_{0}>0$, and there exist at
least one edge in the graph encoded by $\Theta$ and $\widetilde{\Theta}$:
$d,\widetilde{d}>2$. If for all $1\leq k\leq K$,

$n_{k}\geq\left[\frac{3C_{0}}{4}\left(1+\frac{8}{\alpha}\right)^{2}(\kappa_{\widetilde{\Sigma}}\kappa_{\widetilde{\Gamma}}+\kappa_{\widetilde{\Sigma}}^{3}\kappa_{\widetilde{\Gamma}}^{2})\right]^{2}\widetilde{d}^{2}\log
p_{k},$

then with probability at least $1-\sum_{k=1}^{K}p_{k}^{-10}$,
$\widehat{\widetilde{\Theta}}$ defined in (4) satisfies

$\|\widehat{\widetilde{\Theta}}-\widetilde{\Theta}\|_{\infty}\leq\frac{C_{0}}{4}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}.$

######  Theorem 3.

Suppose all conditions in Theorem 2 hold. In addition, if for $1\leq k\leq K$,

$n_{k}\geq\frac{C_{0}\kappa_{\widetilde{\Gamma}}^{2}(1+\frac{8}{\alpha})^{2}}{4\lambda_{\min}(\widetilde{\Theta})^{2}}\min\\{p+\widetilde{s},\widetilde{d}^{2}\\}\log
p_{k},$

then for any $1\leq k\leq K$,

$\|\widehat{\widetilde{\Theta}}^{(k)}-\widetilde{\Theta}^{(k)}\|_{\infty}\leq\frac{3C_{0}}{4}\widetilde{\kappa}^{2}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\min\\{\sqrt{p+\widetilde{s}},\widetilde{d}\\}\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}$

holds with probability at least $1-\sum_{k=1}^{K}p_{k}^{-10}$.

When the sample size condition (16) holds, the conditions in Theorem 2 are all
satisfied, and hence
$\|\widehat{\widetilde{\Theta}}-\widetilde{\Theta}\|_{\infty}\leq\varepsilon_{1}$
defined in Theorem 1. Thus when
$\delta+\varepsilon_{1}\leq\tau_{1}\leq\nu-\delta-\varepsilon_{1}$,

$\widehat{E}_{O}=\\{(i,j)\in O:i\neq
j,|\widehat{\widetilde{\Theta}}_{i,j}|>\tau_{1}\\}\subset\\{(i,j)\in O:i\neq
j,|\widetilde{\Theta}_{i,j}|>\delta\\}=E_{O},$

and

$\widehat{E}_{O}\supset\\{(i,j)\in O:i\neq
j,|\widetilde{\Theta}_{i,j}|>\tau_{1}+\varepsilon_{1}\\}\supset\\{(i,j)\in
O:i\neq j,|\widetilde{\Theta}_{i,j}|>\nu-\delta\\}=E_{O}.$

The sample size condition (16) also ensures that
$\nu-\delta-\varepsilon_{1}>\delta+\varepsilon_{1}$ so that the required range
for $\tau_{1}$ is not empty.

In addition, when the sample size condition (17) holds, all conditions in
Theorem 3 are satisfied, and the range for $\tau_{2}$ is not empty. Thus
$\|\widehat{\widetilde{\Theta}}^{(k)}-\widetilde{\Theta}^{(k)}\|_{\infty}\leq\varepsilon_{2}$.
When $\varepsilon_{2}\leq\tau_{2}<\psi-\varepsilon_{2}$,
$\delta-\varepsilon_{2}<\tau_{3}\leq\nu-\varepsilon_{2}$, the set
$\widehat{W}_{\tau_{2},\tau_{3}}$ defined in (6) satisfies

$\begin{split}\widehat{W}_{\tau_{2},\tau_{3}}&=\\{i\in V:\forall
k~{}s.t.~{}i\in V_{k},~{}\exists j\neq
i,~{}\tau_{2}<|\widehat{\widetilde{\Theta}}^{(k)}_{i,j}|<\tau_{3}\\}\\\
&\subset\\{i\in V:\forall k~{}s.t.~{}i\in V_{k},~{}\exists j\neq
i,~{}0<|\widetilde{\Theta}^{(k)}_{i,j}|<\nu\\}\\\ &=W_{\nu},\end{split}$

and

$\begin{split}\widehat{W}_{\tau_{2},\tau_{3}}&\supset\\{i\in V:\forall
k~{}s.t.~{}i\in V_{k},~{}\exists j\neq
i,~{}\psi\leq|\widetilde{\Theta}^{(k)}_{i,j}|\leq\delta-2\varepsilon\\}\\\
&=\\{i\in V:\forall k~{}s.t.~{}i\in V_{k},~{}\exists j\neq
i,~{}0<|\widetilde{\Theta}^{(k)}_{i,j}|<\delta\\}\\\ &=W_{\delta},\end{split}$

where the second line is due to Assumption 3. Therefore, combining these
results with (19), one has
$\widehat{E}_{O^{c}}=O^{c}\cap(\widehat{W}_{\tau_{2},\tau_{3}}\times\widehat{W}_{\tau_{2},\tau_{3}})=\mathcal{S}_{\mathrm{off}}$.
∎

###### Proof of Theorem 2.

Our proof follows similar ideas to the proof of Theorem 4.1 in (Vinci et al.,,
2019), while the main difference lies in the correlation matrix estimation. In
particular, note that the population MADGQ solution $\widetilde{\Theta}$ is
the same as the one in (Vinci et al.,, 2019), computed from the population
correlation matrix. While for the finite sample estimate
$\widehat{\widetilde{\Theta}}$ we defined in (4), it depends on the rank based
correlation $\widehat{\Sigma}_{O}=\widehat{\Sigma}^{(\tau)}_{O}$ or
$\widehat{\Sigma}^{(\rho)}_{O}$, instead of the sample covariance matrix.

However, one thing to note in the proof of Theorem 4.1 in (Vinci et al.,,
2019) is that we only care about the infinity-norm error bound
$\|\widehat{\Sigma}_{O}-\Sigma_{O}\|_{\infty}$ instead of its specific form.
The following lemma characterizes this error bound for both rank-based
correlation estimates.

######  Lemma 1.

There exists a universal constant $C>0$, such that as long as $n_{k}\geq C$
holds for all $1\leq k\leq K$, with probability at least
$1-\sum_{k=1}^{K}p_{k}^{-10}$, $\widehat{\Sigma}^{(\rho)}_{O}$ and
$\widehat{\Sigma}^{(\tau)}_{O}$ defined in (13) satisfy

$\|\widehat{\Sigma}^{(\rho)}_{O}-\Sigma_{O}\|_{\infty}\leq
C\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}},\quad\|\widehat{\Sigma}^{(\tau)}_{O}-\Sigma_{O}\|_{\infty}\leq
C\max_{k}\sqrt{\frac{\log p_{k}}{n_{k}}}.$

Lemma 1 is proved by applying Theorems 4.1 and 4.2 in (Liu et al.,, 2012) on
each block $k$ and a union bound over $K$ blocks.

Inspired by Lemma 1, we specifically change the upper bound
$\sigma(\bar{n},p^{b})$ in (Vinci et al.,, 2019) to $C\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}$, which is an upper bound for
$\|\widehat{\Sigma}_{O}-\Sigma_{O}\|_{\infty}$ with probability at least
$1-\sum_{k=1}^{K}p_{k}^{-10}$. Similar to (Vinci et al.,, 2019), we let
$\lambda_{j,l}=\frac{C_{0}}{\alpha}\max_{k}\sqrt{\frac{\log p_{k}}{n_{k}}}$
with $C_{0}=8C$. Recognizing the fact that our sample size condition in
Theorem 2 suggests

$C\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}\leq\left[6\left(1+\frac{8}{\alpha}\right)^{2}\widetilde{d}\max\\{\kappa_{\widetilde{\Sigma}}\kappa_{\widetilde{\Gamma}},\kappa_{\widetilde{\Sigma}}^{3}\kappa_{\widetilde{\Gamma}}^{2}\\}\right]^{-1},$

we can follow the same arguments as in (Vinci et al.,, 2019), and eventually
show that

$\\{(j,l):\widehat{\widetilde{\Theta}}_{j,l}\neq
0\\}\subset\\{(j,l):\widetilde{\Theta}_{j,l}\neq 0\\},$ (20)

and

$\|\widehat{\widetilde{\Theta}}-\widetilde{\Theta}\|_{\infty}\leq\frac{C_{0}}{4}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}},$

with probability at least $1-\sum_{k=1}^{K}p_{k}^{-10}$. ∎

###### Proof of Theorem 3.

To prove Theorem 3, we make use of Lemma C.2 and Lemma C.3 in (Vinci et al.,,
2019), which bounds the entrywise error bounds of Schur complements as a
function of the error of full matrices. For completeness, here we restate
these two lemmas.

######  Lemma 2 (Lemma C.2 in (Vinci et al.,, 2019)).

Let $X,Y\in\mathbb{R}^{p\times p}$ be positive definite matrices. Then, for
any nonempty set $A\subset[p]$,

$\left\|X/X_{A^{c}A^{c}}-Y/Y_{A^{c}A^{c}}\right\|_{\infty}\leq\frac{\lambda_{\rm
max}(X)}{\lambda_{\rm min}(X)}\frac{\lambda_{\rm max}(Y)}{\lambda_{\rm
min}(Y)}\left\|X-Y\right\|_{2}$ (21)

where $X/X_{A^{c}A^{c}}=X_{AA}-X_{AA^{c}}(X_{A^{c}A^{c}})^{-1}X_{A^{c}A}$ is
the Schur Complement of the block $A^{c}\times A^{c}$ of the matrix $X$,
$\lambda_{\rm min}(X)$ and $\lambda_{\rm max}(X)$ are the smallest and the
largest eigenvalues of $X$.

######  Lemma 3 (Lemma C.3 in (Vinci et al.,, 2019)).

For any $p\times p$ symmetric matrix $X$ with max row-degree smaller than or
equal to $d$, we have

$\|X\|_{2}\leq\min(\sqrt{\|X\|_{0}},d)\|X\|_{\infty}$ (22)

where $\|X\|_{2}$ is the spectral norm, $\|X\|_{0}:=|\\{(i,j):X_{ij}\neq
0\\}|$, and $\|X\|_{\infty}$ is the max norm.

As has been shown in the proof of Theorem 2, under appropriate conditions,
with probability at least $1-\sum_{k=1}^{K}p_{k}^{-10}$, the support set of
$\widehat{\widetilde{\Theta}}$ is a subset of the support set of
$\widetilde{\Theta}$; and
$\|\widehat{\widetilde{\Theta}}-\widetilde{\Theta}\|_{\infty}\leq\frac{C_{0}}{4}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}$. Therefore, by applying Lemma 2, we have

$\begin{split}\|\widehat{\widetilde{\Theta}}^{(k)}-\widetilde{\Theta}^{(k)}\|_{\infty}\leq\frac{\lambda_{\max}(\widetilde{\Theta})}{\lambda_{\min}(\widetilde{\Theta})}\frac{\lambda_{\max}(\widehat{\widetilde{\Theta}})}{\lambda_{\min}(\widehat{\widetilde{\Theta}})}\|\widehat{\widetilde{\Theta}}-\widetilde{\Theta}\|_{2},\end{split}$
(23)

and an application of Lemma 3 leads to

$\begin{split}\|\widehat{\widetilde{\Theta}}-\widetilde{\Theta}\|_{2}&\leq\min\\{\sqrt{p+\widetilde{s}},\widetilde{d}\\}\|\widehat{\widetilde{\Theta}}-\widetilde{\Theta}\|_{\infty}\\\
&\leq\frac{C_{0}}{4}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\min\\{\sqrt{p+\widetilde{s}},\widetilde{d}\\}\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}\\\
&\leq\frac{1}{2}\lambda_{\min}(\widetilde{\Theta}),\end{split}$

where the last line is due to the sample size condition in Theorem 3. Hence
$\lambda_{\max}(\widehat{\widetilde{\Theta}})\leq\lambda_{\max}(\widetilde{\Theta})+\|\widehat{\widetilde{\Theta}}^{(k)}-\widetilde{\Theta}^{(k)}\|_{2}\leq\frac{3}{2}\lambda_{\max}(\widetilde{\Theta})$,
and
$\lambda_{\min}(\widehat{\widetilde{\Theta}})\geq\lambda_{\min}(\widetilde{\Theta})-\|\widehat{\widetilde{\Theta}}^{(k)}-\widetilde{\Theta}^{(k)}\|_{2}\geq\frac{1}{2}\lambda_{\min}(\widetilde{\Theta})$.
Combining these bounds for $\lambda_{\max}(\widehat{\widetilde{\Theta}})$,
$\lambda_{\min}(\widehat{\widetilde{\Theta}})$,
$\|\widehat{\widetilde{\Theta}}^{(k)}-\widetilde{\Theta}^{(k)}\|_{2}$, with
(23), we have

$\|\widehat{\widetilde{\Theta}}^{(k)}-\widetilde{\Theta}^{(k)}\|_{2}\leq\frac{3C_{0}}{4}\kappa_{\widetilde{\Theta}}^{2}\kappa_{\widetilde{\Gamma}}(1+\frac{8}{\alpha})\min\\{\sqrt{p+\widetilde{s}},\widetilde{d}\\}\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}}.$

∎

###### Proof of Lemma 1.

Firstly, we note that our rank-based correlation matrices
$\widehat{\Sigma}^{(\rho)}$ and $\widehat{\Sigma}^{(\tau)}$ are computed in
three steps: we compute the rank statistics, Kendall’s tau and Spearman’s rho,
for each block as in (15); we then combine the results of $K$ blocks by taking
an average on the overlapping part; finally, two sine transformations are
applied entrywise on the combined rank statistics as in (13). In the
following, we will prove the error bounds following the same three steps, and
leverage similar ideas from the proof of Theorems 4.1 and 4.2 in (Liu et al.,,
2012).

In the first step, we note that for $1\leq k\leq K$,
$\widehat{\rho}^{(k)}_{j,l}$ and $\widehat{\tau}^{(k)}_{j,l}$ are computed as
the standard Spearman’s rho and Kendall’s tau from the data $X^{(k)}$ of the
$k$th block, and hence we can apply the Hoeffding’s inequality for
U-statistics as in the proof of Theorems 4.1 and 4.2 in (Liu et al.,, 2012):

$|\widehat{\rho}^{(k)}_{j,l}-\mathbb{E}\widehat{\rho}^{(k)}_{j,l}|\leq
C\sqrt{\frac{\log
p_{k}}{n_{k}}},\quad|\widehat{\tau}^{(k)}_{j,l}-\mathbb{E}\widehat{\tau}^{(k)}_{j,l}|\leq
C\sqrt{\frac{\log p_{k}}{n_{k}}},$ (24)

for $j,l\in[p_{k}]$ and a universal constant $C>0$, with probability at least
$1-p_{k}^{-10}$. Here, the exponent $-10$ can also be changed by any negative
constant $-c$, as long as $C>0$ in the upper bound in (24) is chosen
appropriately. By the definition of $\widehat{\rho}$ and $\widehat{\tau}$ in
(14) and a union bound over $1\leq k\leq K$, we know that with probability at
least $1-\sum_{k=1}^{K}p_{k}^{-10}$,

$|\widehat{\rho}_{j,l}-\mathbb{E}\widehat{\rho}_{j,l}|\leq
C\max_{k}\sqrt{\frac{\log
p_{k}}{n_{k}}},\quad|\widehat{\tau}_{j,l}-\mathbb{E}\widehat{\tau}_{j,l}|\leq
C\max_{k}\sqrt{\frac{\log p_{k}}{n_{k}}}.$

To transform the bounds above to the bounds for the estimated correlation
matrices, we note that
$\Sigma_{j,l}=\sin\left(\frac{\pi}{2}\mathbb{E}\widehat{\tau}_{j,l}\right)$,
$\widehat{\Sigma}^{(\tau)}_{j,l}=\sin\left(\frac{\pi}{2}\widehat{\tau}_{j,l}\right)$,
and the Lipschitz property of the sine function implies

$|\widehat{\Sigma}^{(\tau)}_{j,l}-\Sigma^{(\tau)}_{j,l}|\leq\frac{\pi}{2}|\widehat{\tau}_{j,l}-\mathbb{E}\widehat{\tau}_{j,l}|.$
(25)

While for $\widehat{\Sigma}^{(\rho)}$, as has been shown in the proof of
Theorem 4.1 in (Liu et al.,, 2012),

$\mathbb{E}\widehat{\rho}^{(k)}_{j_{k},l_{k}}=\frac{6}{\pi}\arcsin\left(\frac{\Sigma_{j,l}}{2}\right)+\frac{6}{\pi(n_{k}+1)}\left(\arcsin(\Sigma_{j,l})-3\arcsin\left(\frac{\Sigma_{j,l}}{2}\right)\right).$

Denote by $a^{(k)}_{j,l}$ the excess error term:

$a^{(k)}_{j,l}=\mathbb{E}\widehat{\rho}^{(k)}_{j_{k},l_{k}}-\frac{6}{\pi}\arcsin\left(\frac{\Sigma_{j,l}}{2}\right)=\frac{6}{\pi(n_{k}+1)}\left(\arcsin(\Sigma_{j,l})-3\arcsin\left(\frac{\Sigma_{j,l}}{2}\right)\right),$

then we know that $|a^{(k)}_{j,l}|\leq\frac{12}{n_{k}+1}$. Thus, by the
definition of $\widehat{\rho}$ in (14),

$\mathbb{E}\widehat{\rho}_{j,l}=\frac{6}{\pi}\arcsin\left(\frac{\Sigma_{j,l}}{2}\right)+\frac{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}a^{(k)}_{j,l}}{\sum_{k=1}^{K}{{\mathbbm{1}}}_{\\{j,l\in V_{k}\\}}},$

and

$\begin{split}|\widehat{\Sigma}^{(\rho)}_{j,l}-\Sigma^{(\rho)}_{j,l}|&=\left|2\sin\left(\frac{\pi}{6}\widehat{\rho}_{j,l}\right)-2\sin\left(\frac{\pi}{6}\mathbb{E}\widehat{\rho}_{j,l}-\frac{\sum_{k}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}a^{(k)}_{j,l}}{\sum_{k}{{\mathbbm{1}}}_{\\{j,l\in
V_{k}\\}}}\right)\right|\\\
&\leq\frac{\pi}{3}\left(|\widehat{\rho}_{j,l}-\mathbb{E}\widehat{\rho}_{j,l}|+\frac{12}{\min_{k}n_{k}+1}\right).\end{split}$
(26)

Therefore, combining (26), (25), (24), the proof of Lemma 1 is completed. ∎