# Gaussian Graphical Model Selection
for Huge Data via Minipatch Learning

Tianyi Yao and Minjie Wang and Genevera I. Allen , Department of Statistics,
Rice University, Houston, TXSchool of Statistics, University of Minnesota,
Minneapolis, MNDepartments of Electrical and Computer Engineering, Statistics,
and Computer Science, Rice University, Houston, TXJan and Dan Duncan
Neurological Research Institute, Baylor College of Medicine, Houston, TX

###### Abstract

Gaussian graphical models are essential unsupervised learning techniques to
estimate conditional dependence relationships between sets of nodes. While
graphical model selection is a well-studied problem with many popular
techniques, there are typically three key practical challenges: i) many
existing methods become computationally intractable in huge-data settings with
tens of thousands of nodes; ii) the need for separate data-driven tuning
hyperparameter selection procedures considerably adds to the computational
burden; iii) the statistical accuracy of selected edges often deteriorates as
the dimension and/or the complexity of the underlying graph structures
increase. We tackle these problems by proposing the Minipatch Graph (MPGraph)
estimator. Our approach builds upon insights from the latent variable
graphical model problem and utilizes ensembles of thresholded graph estimators
fit to tiny, random subsets of both the observations and the nodes, termed
minipatches. As estimates are fit on small problems, our approach is
computationally fast with integrated stability-based hyperparameter tuning.
Additionally, we prove that under certain conditions our MPGraph algorithm
achieves finite-sample graph selection consistency. We compare our approach to
state-of-the-art computational approaches to Gaussian graphical model
selection including the BigQUIC algorithm, and empirically demonstrate that
our approach is not only more accurate but also extensively faster for huge
graph selection problems.

##  1 Introduction

Gaussian graphical models are popular unsupervised learning approaches for
estimating and uncovering conditional dependence structure between variables
or nodes (Lauritzen, 1996). Graphical models have found wide application in
systems biology, imaging, finance, and natural language processing, among many
others. In several of these domains, we do not know the inherent graph
structure and seek to learn this from data, a problem called graphical model
selection. With larger and larger data sets in a variety of domains, we are
interested in performing graphical model selection on data with a huge number
of nodes. Consider examples from biomedicine. In neuroscience, Gaussian
graphical models are often used to learn the functional connectivity between
brain regions or neurons (Friston, 2011; Yatsenko et al., 2015; Chang et al.,
2019). In fMRI data, the number of brain regions and hence graph nodes can be
on the order of tens- to hundreds-of-thousands (Friston, 2011) whereas in
calcium imaging, the number of neurons or graph nodes in some newer
technologies can measure in the tens-of-thousands (Stringer and Pachitariu,
2019; Vinci et al., 2019).

Many existing Gaussian graphical model selection techniques are
computationally too cumbersome to be applied in such large data settings.
Further, these techniques require hyperparamter tuning to determine the
sparsity or edge set of the graph; this process often creates an
insurmountable computational burden for large data. Finally, the statistical
accuracy of Gaussian graphical model selection methods is known to degrade
with a large number of nodes, creating further challenges. In this paper, our
goal is to develop a radically different type of computational method for
learning the structure of Gaussian graphical models that not only yields
statistically more accurate estimation, but is also more scalable for huge
data.

Formally, let $X=(X_{1},X_{2},\ldots,X_{M})$ be a $M$-dimensional random
vector following a multivariate Gaussian distribution
$\mathcal{N}(\mathbf{0},\bm{\mathbf{\Sigma}})$ with covariance matrix
$\bm{\mathbf{\Sigma}}$ and inverse covariance matrix
$\bm{\mathbf{\Theta}}=\bm{\mathbf{\Sigma}}^{-1}$. Graphical models are often
denoted by $G=(V,E)$, in which the node set $V$ represents the collection of
random variable $X$ and edge set $E$ characterizes the conditional dependence
relationship between these random variables. An absence of edge between
$X_{i}$ and $X_{j}$ means that $X_{i}$ and $X_{j}$ are conditionally
independent; for Gaussian graphical models, this corresponds to zero entries
in the precision matrix $\Theta_{ij}=\Theta_{ji}=0$. One important objective
of Gaussian graphical model selection or learning the graph structure is to
infer the edge set $E$, or the sparsity pattern of $\bm{\mathbf{\Theta}}$,
from the observed data.

###  1.1 Related Works

One popular approach to Gaussian graphical model selection problem is the
sparse inverse covariance estimator introduced by (Yuan and Lin, 2007;
Banerjee et al., 2008). Since its inception, a plethora of optimization
algorithms (Duchi et al., 2008; d’Aspremont et al., 2008; Rolfs et al., 2012)
have been proposed to solve this problem, including the widely-used graphical
Lasso procedure (Friedman et al., 2008). However, as noted in Hsieh et al.
(2013), many of these solvers do not scale to problems with thousands or tens
of thousands of nodes, which is commonly observed in modern applications. The
BigQUIC algorithm (Hsieh et al., 2013) is arguably the state-of-the-art
computational approach to the sparse inverse covariance estimation. While
numerical studies (Hsieh et al., 2013) have shown that it can solve the
estimation problem with a million nodes within a day, these experiments were
conducted in somewhat ideal conditions where the single optimal tuning
hyperparameter value is assumed to be known a priori, which is rare in
practice. The statistical and computational performance of the BigQUIC
algorithm in more realistic settings with data-driven hyperparameter tuning
have been inadequately investigated.

A related line of work seeks to lower the computational cost of estimating
sparse Gaussian graphical models by soft-thresholding the sample covariance
matrix prior to directly solving the sparse inverse covariance estimation
problem (Fattahi and Sojoudi, 2019; Zhang et al., 2020). However, these
approaches generally require imposing specific sparsity constraints such as
assuming the underlying graph has a chordal structure, which are often too
restrictive to hold in real-world applications of Gaussian graphical model
selection.

In addition, Meinshausen et al. (2006) proposed the neighborhood selection
procedure, which separately fits the Lasso regression estimator to each nodes
using all remaining nodes as predictors and then infers the edge set by
aggregating the neighbors selected from all nodewise regressions. Liu and Wang
(2017) further extended the neighborhood selection approach by replacing the
Lasso estimator with the more tuning-insensitive SQRT-Lasso estimator (Belloni
et al., 2011), which yields the TIGER method for estimating high-dimensional
Gaussian graphical models. Another line of work proposed to employ the Dantzig
selector for solving each nodewise sparse regression subproblem, resulting in
the CLIME method (Cai et al., 2011). Despite generally faster computation due
to the easily parallelizable nature of the nodewise regressions, these
Gaussian graphical model selection strategies can still become computationally
challenging as the number of nodes $M$ increases because they need to solve
$M$-dimensional sparse regression problems for a total of $M$ times even for a
single value of tuning hyperparameter.

Importantly, most of the aforementioned computational approaches to the
Gaussian graphical model selection problem have at least one tuning
hyperparameter that controls how many edges and which edges are selected, thus
significantly impacting the statistical accuracy of the estimated sets of
edges (Müller et al., 2016). Even though many of the existing methods were
shown to enjoy certain guarantees for the accuracy of the edges that they
select, such theoretical guarantees are often contingent on some oracle tuning
hyperparameter choices that cannot be used in practice, as astutely pointed
out by Liu and Wang (2017). In realistic applications of Gaussian graphical
model selection, these existing approaches need to be employed in conjunction
with data-driven tuning hyperparameter selection procedures such as cross-
validation (Efron, 1982; Wasserman and Roeder, 2009), extended BIC (Chen and
Chen, 2008; Foygel and Drton, 2010), StARS (Liu et al., 2010), and B-StARS
(Müller et al., 2016). However, the statistical accuracy of selected edges
often degrades as the dimension and/or the complexity of the underlying graph
structures increase when the tuning hyperparameters are chosen via these data-
dependent ways without oracle tuning information. Additionally, these data-
driven tuning hyperparameter selection strategies usually require solving the
graph selection problem at least once for each candidate tuning hyperparameter
value for sometimes hundreds of them, thus considerably adding to the overall
computational burden.

#### Contributions

We summarize our main contributions as follows: We develop a fast
computational approach to learn the structure of Gaussian graphical models
named Minipatch Graph (MPGraph) that ensembles thresholded graph estimators
trained on tiny, random subsets of both observations and nodes, termed
minipatches (see Sec. 2). Leveraging such ensembles, our approach not only
selects edges with improved statistical accuracy, but is also computationally
fast with integrated stability-based hyperparameter tuning. Additionally, we
theoretically analyze MPGraph showing that it attains finite-sample graph
selection consistency under certain conditions. We further empirically
demonstrate the practical advantages of our approach through extensive
comparative studies, showing that MPGraph dominates state-of-the-art
computational approaches to Gaussian graphical model selection in terms of
both edge selection accuracy and computational time for huge graph selection
problems. The major innovation of our work lies in the idea to split up huge
graph selection problems into smaller and more computationally manageable
units (minipatches) and then ensemble the results back together to obtain a
consistent graph estimate.

##  2 Minipatch Graphical Model Selection

###  2.1 Review: Minipatch Learning

The idea of using random subsets of observations and/or features for model
training has appeared in many areas of machine learning including random
forests (Breiman, 2001), stochastic optimization (Hardt et al., 2016), dropout
training (Srivastava et al., 2014). Recently, a line of work in feature
selection (Yao and Allen, 2021) and ensemble learning (Yao et al., 2021;
Toghani and Allen, 2021) coined the term “minipatches” to denote tiny, random
subsets of both observations and features of the data. Following the notations
in Yao et al. (2021), given observed data matrix
$\mathbf{X}\in\mathbb{R}^{N\times M}$ that comprises $N$ observations each
having $M$ features, a minipatch is obtained by simultaneously subsampling $n$
rows (observations) and $m$ columns (features) without replacement using some
form of randomization, usually with $n\ll N$ and $m\ll M$.

###  2.2 Minipatch Graph (MPGraph)

Inspired by the minipatch feature selection approaches (Yao and Allen, 2021)
that achieve major computational savings through the use of minipatches, we
seek to leverage the concept of minipatch learning for tackling both the
statistical and computational challenges with Gaussian graphical model
selection in huge-data settings.

At first glance, directly applying the idea of minipatch learning to Gaussian
graphical model selection seems straightforward as one might ask whether we
can simply train separate graph estimators on many minipatches and then
aggregate all these subgraph estimates to form a complete graph. However, such
simplistic approach is actually problematic due to the complex dependency
structure between the nodes/variables in the data. Concretely, suppose that we
take a random minipatch that consists of a subset of nodes from the original
data set and try to estimate the conditional dependence graph between these
nodes by fitting a graph estimator to this minipatch. Then, the remaining
nodes outside of this minipatch effectively become latent unobserved nodes
with respect to the minipatch, thus inducing many false positive edges of
small magnitude in this subgraph estimate (Vinci et al., 2019). Estimating
graphs for each minipatch thus leads to the well-known latent variable
graphical model problem (Chandrasekaran et al., 2012).

The main idea of our approach is then to quickly solve the latent variable
graphical model problem to yield a consistent graph structural estimate for
each minipatch and then ensemble the results. To this end, we leverage a
simple and computationally fast approach recently proposed by Wang and Allen
(2021) to solve each latent variable problem, which does not require low-
rankness for the latent components. They apply a hard thresholding operator to
the graphical Lasso (GLasso) estimate and show that this is sparsistent for
Gaussian graphical model selection in the presence of latent variables under
fairly weak assumptions. Putting everything together, we propose a new graph
selection method, MPGraph, which utilizes minipatch learning and thresholded
graph estimation techniques to select statistically accurate edges in a
computationally efficient manner.

Our proposed MPGraph framework is summarized in Algorithm 1. Following the
notational conventions, $\mathbf{X}_{I_{k},F_{k}}$ denotes the submatrix of
$\mathbf{X}$ containing its rows indexed by $I_{k}$ and its columns indexed by
$F_{k}$. Similarly, $\bm{\mathbf{\Theta}}_{i,j}$ denotes the element the at
the $i$th row and $j$th column. Also, $F_{k,a}$ represents the $a$th element
in the vector $F_{k}$. For simplicity, $[N]$ denotes the set
$\\{1,2,\ldots,N\\}$. MPGraph fits base thresholded graph estimators to many
tiny, random minipatches. On the $k$th minipatch, the node sampling indicator
$D^{(k)}_{ij}$ is set to 1 if and only if both nodes $i$ and $j$ are sampled
into this minipatch; similarly, the edge selection indicator $S^{(k)}_{ij}$ is
set to 1 if there is an edge between node $i$ and $j$, as determined by the
base thresholded graph estimator. After $K$ iterations, MPGraph computes the
edge selection frequencies by taking an ensemble of edge selection events over
all minipatches. We define the selection frequency of the edge between node
$i$ and $j$, $\hat{\Pi}^{(K)}_{ij}$, to be the number of times both nodes are
sampled together and have an estimated edge between them as determined by the
base graph estimator divided by the number of times both are sampled together
into minipatches. MPGraph finally produces a set of stable edges
$\hat{E}^{\text{stable}}$ whose selection frequencies are above a user-
specific threshold $\pi_{\text{thr}}\in(0,1)$. We discuss the choice of
$\pi_{\text{thr}}$ in Sec. 2.4. Because the computations on each minipatch can
be done independently, the MPGraph framework is embarrassingly parallelizable,
which can yield enormous computational savings.

Algorithm 1 MPGraph

Input: $\mathbf{X}\in\mathbb{R}^{N\times M}$, $n$, $m$,
$\pi_{\text{thr}}\in(0,1)$

for $k=1,2,\ldots,K$ do // In parallel

1) Sample a minipatch: subsample $n$ observations $I_{k}\subset[N]$ and $m$
features $F_{k}\subset[M]$ uniformly at random without replacement to get a
minipatch $\mathbf{X}_{I_{k},F_{k}}\in\mathbb{R}^{n\times m}$

2) Fit base thresholded graph estimator to $\mathbf{X}_{I_{k},F_{k}}$ to
obtain estimated subgraph
$\widetilde{\bm{\mathbf{\Theta}}}^{(k)}\in\mathbb{R}^{m\times m}$

3) Map indices: for each pair of nodes in the minipatch $\\{(i,j)\in
F_{k}\times F_{k}:i<j\\}$, set
$(i^{\prime},j^{\prime})=\\{(a,b)\in[m]\times[m]:F_{k,a}=i,F_{k,b}=j\\}$

4) Update edge selection indicator $S_{ij}^{(k)}$ and node sampling indicator
$D_{ij}^{(k)}$, $\forall\>1\leq i<j\leq M$:  $\displaystyle S_{ij}^{(k)}$
$\displaystyle=\mathds{1}(i\in F_{k},j\in
F_{k},\widetilde{\bm{\mathbf{\Theta}}}^{(k)}_{i^{\prime},j^{\prime}}\neq
0)\>\>\text{and}\>\>D_{ij}^{(k)}=\mathds{1}(i\in F_{k},j\in F_{k})$

end for

Compute edge selection frequencies $\hat{\Pi}_{ij}^{(K)},\forall\>1\leq
i<j\leq M$:
$\displaystyle\hat{\Pi}^{(K)}_{ij}=\frac{\sum_{k=1}^{K}S_{ij}^{(k)}}{\text{max}(1,\sum_{k=1}^{K}D_{ij}^{(k)})}$

Output: $\hat{E}^{\text{stable}}=\big{\\{}1\leq i<j\leq
M:\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}}\big{\\}}$

As mentioned, we adopt the hard thresholded GLasso (TGLasso) approach from
Wang and Allen (2021) as the latent variable graphical model estimator in Step
2 of Algorithm 1. We describe this in detail in Algorithm 2. We save the
investigations of other types of internal graph estimation strategies for
future work. Given a minipatch, Algorithm 2 fits the GLasso at a small amount
of regularization $\lambda_{0}$ to get an initial subgraph estimate
$\hat{\bm{\mathbf{\Theta}}}^{\lambda_{0}}_{k}$, as suggested by Wang and Allen
(2021). After that, it hard-thresholds the initial estimate at a sequence of
threshold values $\tau_{k}$ and eventually outputs the best thresholded
subgraph estimate $\widetilde{\bm{\mathbf{\Theta}}}^{(k)}$ according to the
eBIC criterion. $\gamma$ is a user-specific parameter for controlling the
amount of regularization from the eBIC criterion. We found setting $\gamma$ to
a default value of $0.5$ works well in practice.

Algorithm 2 Example Base Thresholded Graph Estimator for Step 2) in Algorithm
1

Input: $\mathbf{X}_{I_{k},F_{k}}$, $\gamma$

1) Fit GLasso to minipatch $\mathbf{X}_{I_{k},F_{k}}$ at regularization
parameter $\lambda_{0}=\sqrt{\log{m}/n}$ to get
$\hat{\bm{\mathbf{\Theta}}}^{\lambda_{0}}_{k}$

2) Extract the maximum off-diagonal element: $\theta_{k}=\max_{1\leq i<j\leq
m}|[\hat{\bm{\mathbf{\Theta}}}^{\lambda_{0}}_{k}]_{ij}|$

for $\tau_{k}=0.1\theta_{k},0.2\theta_{k},\ldots,\theta_{k}$ do

1) Hard-threshold $\hat{\bm{\mathbf{\Theta}}}^{\lambda_{0}}_{k}$ element-wise
at $\tau_{k}$:
$\widetilde{\bm{\mathbf{\Theta}}}^{(k)}_{\tau_{k}}=T_{\tau_{k}}(\hat{\bm{\mathbf{\Theta}}}^{\lambda_{0}}_{k})$

2) Record the number of edges in
$\widetilde{\bm{\mathbf{\Theta}}}^{(k)}_{\tau_{k}}$ as $e_{\tau_{k}}$

end for

3) Choose the best threshold level using eBIC criterion:
$\displaystyle\tau_{k}^{*}=\operatorname*{arg\,min}_{\tau_{k}\in\\{0.1\theta_{k},0.2\theta_{k},\ldots,\theta_{k}\\}}$
$\displaystyle-2\ell(\widetilde{\bm{\mathbf{\Theta}}}^{(k)}_{\tau_{k}})+e_{\tau_{k}}\log
n+4e_{\tau_{k}}\gamma\log m$ where $\ell$ is the log-likelihood of the
multivariate Gaussian distribution

Output:
$\widetilde{\bm{\mathbf{\Theta}}}^{(k)}=\widetilde{\bm{\mathbf{\Theta}}}^{(k)}_{\tau_{k}^{*}}\in\mathbb{R}^{m\times
m}$

#### Relation to the StARS Procedure

The Stability Approach to Regularization Selection (StARS) (Liu et al., 2010)
is a popular data-driven approach to choosing tuning hyperparameter values for
a variety of Gaussian graphical model estimators. Similar to our MPGraph
method, StARS repeatedly fits graph estimators to $K$ random subsamples of the
observations but includes all the features and then computes a so-called edge
instability score by aggregating edge selection events from all subsamples. In
contrast, our MPGraph approach employs tiny subsamples of both observations
and features, thus necessitating latent variable graphical model estimates but
also dramatically reducing the overall computational costs. Secondly, StARS
seeks to find a sparse, stable graph and hence solves the graph selection
problem repeatedly at a sequence of tuning hyperparameter values on each
subsample. On the other hand, our MPGraph approach performs graph sparsity
tuning on each minipatch separately, thus dramatically reducing computational
costs. In the end, our method produces a set of stable edges by simply
aggregating over the minipatches with no additional tuning or computation.

###  2.3 Theoretical Analysis

We have presented a new, fast computational method for learning the structure
of Gaussian graphical models. But since our method is novel and does not solve
known Gaussian graphical model selection problems (Meinshausen et al., 2006;
Yuan and Lin, 2007; Cai et al., 2011), one naturally asks: Does MPGraph
correctly recover the edge structure of the Gaussian graphical model? In other
words, is our method sparsistent or model selection consistent? In this
section, we establish such theoretical results.

Recall that we employ the thresholded GLasso on each minipatch to obtain a
graph estimate. Wang and Allen (2021) show that this method is graph selection
consistent in the presence of latent variables under certain conditions on the
true graph and corresponding covariance matrix as well as the Schur complement
of the latent variables. Then, if we can show that all of these conditions are
satisfied for all possible random minipatches, our MPGraph approach will also
be graph selection consistent. This is the approach we take.

First, denote the true edge set by
$E^{*}=E(\Theta^{*})=\\{(i,j):\Theta^{*}_{ij}\neq 0\>\text{and}\>i\neq j\\}$
and $s=\operatorname*{arg\,max}_{k}s_{k}$ where
$s_{k}=|E\left(\Theta_{F_{k}}^{*}\right)|$ refers to the total number of non-
zero edges for the $k$th minipatch $F_{k}$. Also let $\varphi_{\max}()$ and
$\varphi_{\min}()$ denote the maximum and minimum eigenvalue respectively.
Denote the elementwise $\ell_{\infty}$-norm of a matrix by
$\|\mathbf{A}\|_{\infty}=\max_{i\neq j}|A_{ij}|$. Consider the following
assumptions:

######  Assumption 1 (restate = assumptionMPGone, name = )

$X_{i}\text{ be i.i.d. }\mathcal{N}\left(\mathbf{0},\Sigma^{*}\right)$.

######  Assumption 2 (restate = assumptionMPGtwo, name = )

$\varphi_{\min}\left(\Sigma^{*}\right)\geq\underline{\kappa}>0,$ or
equivalently $\varphi_{\max}\left(\operatorname{\Theta}^{*}\right)\leq
1/\underline{\kappa}$.

######  Assumption 3 (restate = assumptionMPGthree, name = )

$\varphi_{\max}\left(\Sigma^{*}\right)\leq\overline{\kappa}$.

######  Assumption 4 (restate = assumptionMPGfour, name = )

Minimum edge strength:

$\theta_{\min}:=\min_{(i,j)\in
E\left(\operatorname{\Theta}^{*}\right)}\left|\operatorname{\Theta}_{ij}^{*}\right|>c_{1}\sqrt{\frac{s\log
m}{n}}.$

######  Assumption 5 (restate = assumptionMPGfiveraw, name = )

Maximum effect of unsampled nodes:

$\displaystyle\max\limits_{F_{k}}\|\left(\operatorname{\Theta}_{F_{k}}^{*}-\operatorname{\Theta}_{F_{k}F_{k}^{c}}^{*}\left(\operatorname{\Theta}_{F_{k}^{c}}^{*}\right)^{-1}\operatorname{\Theta}_{F_{k}^{c}F_{k}}^{*}\right)^{-1}-\left(\operatorname{\Theta}_{F_{k}}^{*}\right)^{-1}\|_{\infty}$
$\displaystyle\leq\|\hat{\Sigma}_{F_{k}}-\Sigma_{F_{k}}^{*}\|_{\infty}.$

Now, let $R_{ij}$ denote the number of times any pair of nodes $i,j$ are
sampled together into minipatches. Define constant $\tau$ to satisfy
$m^{\tau}=b_{1}\exp(-b_{2}n\lambda^{2})$ where $b_{1}$ and $b_{2}$ depends on
$\overline{\kappa}$ and $\underline{\kappa}$; additionally, define
$w_{1}=1-m^{\tau}$ and $w_{2}=m^{\tau}$. Given these, we can state our main
result:

###### Theorem 1

Let Assumptions 1-5 be satisfied and let $n$ grow proportionally with $N$.
Then, the minipatch graph selection estimator, MPGraph, with
$\lambda\asymp\sqrt{\frac{\log m}{n}}$, is graph selection consistent with
high probability:

$\displaystyle\mathbb{P}\left(\hat{E}^{\text{stable}}=E^{*}\right)$
$\displaystyle\geq
1-\sum_{l=1}^{2}\exp\left\\{-\min_{i,j}R_{ij}\cdot\frac{(w_{l}-\pi_{\text{thr}})^{2}}{2m^{\tau}(1-m^{\tau})}+2\ln
m\right\\}$ $\displaystyle\to 1,\hskip 17.07164pt\text{ as }K\to\infty,\text{
or as }N\to\infty.$

Theorem 1 establishes the graph selection consistency of our MPGraph method.
We provide a detailed proof in the Supplementary Materials, but pause to
briefly outline our approach. We show that Assumptions 1-5 imply all the
assumptions required from Wang and Allen (2021) for graph selection
consistency in the presence of latent variables are satisfied for all possible
minipatches. Specifically, we show that the stated eigenvalue conditions on
$\Sigma^{*}$ imply the necessary restrictions on the covariance matrix.
Assumptions 4 and 5 are the key assumptions ensuring the overall graph
selection consistency of our MPGraph approach. Intuitively, Assumption 4
requires the minimum edge strength of the true graph to be sufficiently large
in magnitude such that the thresholding procedure will not remove true edges
on a minipatch. Further, Assumption 5 suggests that the effect of unsampled
(latent) variables should not exceed the estimation error of the covariance
matrix. This assumption ensures that the false edges induced by latent
variables are small in magnitude, thus allowing the thresholding procedure to
successfully remove them. In practice, we have some control to ensure
Assumption 5 holds by making sure the number of nodes $m=|F_{k}|$ in all
minipatches is not too small (as investigated further in Section 3.1).
Additionally, it is worth noting that we do not require low-rankness for the
unsampled (latent) component because we employ the thresholded GLasso (Wang
and Allen, 2021) approach on each minipatch.

Further, Theorem 1 suggests that the sample complexity required for our
MPGraph estimator is $n=\mathcal{O}(s\log m)$ for each minipatch. We expect
that the overall sample complexity is of the same order as neighborhood
selection and verify this empirically in Figure1A. Overall, this theoretical
result provides the reassurance that one can correctly estimate the graph
structure with high probability using ensembles of minipatches.

###  2.4 Practical Considerations

Our MPGraph is a general meta-algorithm with three tuning hyperparameters:
minipatch size ($n$ and $m$) and threshold $\pi_{\text{thr}}$. In practice,
our method is quite robust to the choice of these hyperparameters. We have
additional empirical experiments in the Supplementary Materials investigating
how edge selection accuracy and computational time of MPGraph vary for various
minipatch sizes. As a general rule of thumb, we found taking $m$ to be
$5\%\sim 10\%$ of $M$ and then picking $n$ relative to $m$ so that it
surpasses the sample complexity of the base graph estimator used in the meta-
algorithm strikes a good balance between statistical accuracy and
computational time. In addition, we found fixing the threshold
$\pi_{\text{thr}}$ to $0.5$ and the total number of iterations $K=1000$ works
well for many problems in practice. Note that the internal latent variable
graph selection estimator also has tuning parameters, in our case the
threshold level $\tau$, but we strongly suggest tuning these separately for
each minipatch (as we outline in Algorithm 2) and not all together across
minipatches.

##  3 Empirical Studies

In this section, we follow standard computational benchmarks (Liu et al.,
2010; Hsieh et al., 2013; Müller et al., 2016) to compare our proposed MPGraph
method with a wide range of existing Gaussian graphical model selection
strategies in terms of both edge selection accuracy and computational time. In
particular, we consider the following widely-used groundtruth graph structures
for generating synthetic benchmark data sets:

  * •

Scenario 1: Chain graph: The groundtruth precision matrix
$\bm{\mathbf{\Theta}}\in\mathbb{S}^{M\times M}_{+}$ is set to be
$\Theta_{i,i-1}=\Theta_{i-1,i}=0.6$ and $\Theta_{i,i}=1.25$ for all
$i=1,\ldots,M$.

  * •

Scenario 2: Erdős-Renyi graph: First generate an adjacency matrix
$\mathbf{A}\in\\{0,1\\}^{M\times M}$ whose edge patterns follow a Erdős-Renyi
graph with edge probability $p$ such that the total number of undirected edges
is $(M-1)$. Then the groundtruth precision matrix $\bm{\mathbf{\Theta}}$ is
set to be $\Theta_{i,j}=\Theta_{j,i}=\rho_{ij}$ if $A_{i,j}=1$ and
$\Theta_{i,j}=\Theta_{j,i}=0$ otherwise, where
$\rho_{ij}\sim\text{uniform}([-0.6,-0.3]\cup[0.3,0.6])$. Lastly, set
$\Theta_{i,i}=1$ for all $i=1,\ldots,M$.

  * •

Scenario 3: Small-World graph: First generate an adjacency matrix
$\mathbf{A}\in\\{0,1\\}^{M\times M}$ whose edge patterns follow a Watts-
Strogatz small-world graph with $2$ nearest neighbors and edge rewiring
probability of $0.5$. Then the groundtruth precision matrix
$\bm{\mathbf{\Theta}}$ is generated in the same way as in the Erdős-Renyi
graph scenario.

For a given groundtruth precision matrix $\bm{\mathbf{\Theta}}$, we generate
the data matrix $\mathbf{X}\in\mathbb{R}^{N\times M}$ whose rows are
independently drawn from a $M$-variate Gaussian distribution
$\mathcal{N}(\mathbf{0},\bm{\mathbf{\Theta}}^{-1})$. The groundtruth edge set
is $E=\\{1\leq i<j\leq M:\Theta_{i,j}\neq 0\\}$. We evaluate the edge
selection accuracy of each method in terms of the F1 Score, which takes values
between $0$ and $1$ with $1$ signifying perfect match between the estimated
edge set and the groundtruth edge set $E$. We consider a range of low-
dimensional and high-dimensional experimental scenarios with various $N$ and
$M$ sizes in Sec. 3.1-3.3. All comparisons were conducted on a VM with $12$
vCPUs (Intel Cascade Lake) with $240$ GB of memory.

Figure 1: Validation of Theoretical Results. (A) Probability of exact edge-set
recovery versus total sample size $N$ for chain graph simulations with varying
number of nodes $M$. Each point represents the average over $100$ trials. (B)
Edge selection accuracy (F1 Score) comparisons for small-world graph
simulations with fixed $N=500$ and varying $M$. Note that oracle parameter
tuning is used for all methods. (C) Edge selection accuracy of MPGraph using
various minipatch sizes (i.e. $m/M$) for the same simulations in (B).

###  3.1 Validation of Theoretical Results

Before conducting thorough comparative studies in Sec. 3.2-3.3, we want to
first empirically illustrate and verify the theoretical properties of MPGraph.
Following similar experiments in Ravikumar et al. (2011), we show the
probability of exact edge-set recovery versus sample size $N$ for chain graph
simulations with varying number of nodes $M$ in Figure 1A. We see that the
recovery probability of MPGraph follows that of neighhorhood selection (NS)
closely, and reaches 1 at the same $N$ as NS. This supports our theoretical
analysis as it suggests the overall sample complexity of MPGraph is of the
same order as NS for consistent graph selection. In addition, we compare the
edge selection accuracy of MPGraph to standard GGM selection procedures using
oracle tuning approaches (i.e. assume $|E|$ is known) in Figure 1B. We see
that even though accuracy of all methods degrades as the graph dimension $M$
increases, MPGraph generally achieves better F1 Score even in high-dimensional
settings where exact edge-set recovery might not be attainable.

In Figure 1C, we verify that the number of nodes in a minipatch $m$ is the
determining factor for how well MPGraph recovers the true graph because it
dictates whether Assumptions 4 and 5 would be largely satisfied in practice.
For the challenging small-world graph situations in Figure 1C, Assumption 5 is
likely not satisfied when $m$ is too small, hence performance of MPGraph
degrades precipitously. But with a sufficiently large $m$ (e.g. at least $5\%$
of $M$), we can achieve better performance. More results are in the
Supplementary Materials. Thus, we recommend choosing $m$ to be $5\%\sim 10\%$
of $M$ and then picking $n$ relative to $m$ so that it surpasses the sample
complexity of the base graph estimator used on the minipatches, which would
allow ample room for Assumption 4 and 5 to be satisfied for reasonable graphs
in practice.

###  3.2 Comparative Empirical Studies

In this section, our goal is to empirically compare our proposed MPGraph
method to a wide range of existing computational approaches and state-of-the-
art solvers implemented in widely-used software packages for Gaussian
graphical model selection. The approaches we consider include the sparse
inverse covariance estimation (GL), neighborhood selection (NS), CLIME, and
TIGER. In particular, for the GL approach, we consider the QUIC algorithm
(Hsieh et al., 2014) (skggm (Laska and Narayan, 2017)), the BigQUIC algorithm
(BigQuic) (Hsieh et al., 2013), and various coordinate-descent-based
algorithms (glasso (Friedman et al., 2008), huge (Zhao et al., 2012), pulsar
(Müller et al., 2016), sklearn (Pedregosa et al., 2011)). In addition, we use
both huge and pulsar for NS, pulsar for TIGER, and flare (Li et al., 2012) for
CLIME. Our MPGraph algorithm is implemented in Python.

In order to mimic realistic use cases of Gaussian graphical model selection in
practice, we do not assume the oracle optimal tuning hyperparameter is known.
Instead, widely-used tuning hyperparameter selection schemes from the
aforementioned software packages are used to choose the optimal hyperparameter
values for each method in a data-driven manner. Specifically, a geometric
sequence of candidate hyperparameter values $\Lambda$ of length $100$ in the
range of $(\delta\lambda_{\text{max}},\lambda_{\text{max}})$ are used for all
software packages, where the ratio $\delta\in(0,1)$ is set to the value
recommended by the respective software package. Additionally, we use the
default minipatch size (i.e. $m/M=5\%$ and $m/n=0.8$) and default settings
(i.e. $\pi_{\text{thr}}=0.5$ and $K=1000$) for MPGraph for comparisons in this
section. Last but not least, parallelism is enabled for all methods whose
software implementations include this functionality.

Figure 2: Edge Selection Accuracy (F1 Score) and Computational Time from
Simulation Scenarios 1-3 for a Variety of Dimensions. Parallelism is enabled
for all methods whose software packages include this functionality, as
indicated by the (P) after the method name. Our MPGraph method outperforms all
competitors in terms of edge selection accuracy across all four data sets.
MPGraph is one of the computationally fastest methods.

In Figure 2, we display the performance of all methods on the three graph
types for a variety of dimensions. While some competing methods such as GL
BigQuic StARS and NS pulsar B-StARS attain similar edge selection accuracy to
MPGraph in scenarios with simpler graph structure (Figure 2A) or low-
dimensional settings (Figure 2C), their accuracy deteriorates considerably
when the complexity of the underlying graph structure and/or the dimension of
the problem increase (Figure 2B, 2D). On the other hand, the edge selection
accuracy of MPGraph does not degrade as much when the underlying graphical
model selection problem becomes more challenging. In fact, MPGraph achieves
the best edge selection accuracy across all four data sets, outperforming most
of the competitors by a large margin. In the meantime, MPGraph is one of the
computationally fastest methods, with computational time on the order of
seconds to a minute as opposed to hours required by many competitors,
especially in the more challenging scenarios (Figure 2B, 2D).

Table 1: Experimental Results for Large-Scale Gaussian Graphical Model
Selection Problems. The number of selected edges, true positive rate (TPR),
precision, F1 Score, and computational time in seconds are reported for each
method. Note that parallelism is enabled for all methods whose software
implementations include this functionality, as indicated by the (P) after the
method name. The best F1 Score and runtime for each scenario are bold-faced.
Our MPGraph method achieves the best edge selection accuracy (F1 Score) across
all scenarios, and it is the computationally the fastest, especially for huge
graph selection problems with ten thousand nodes.

| Method | # Selected Edges | TPR | Precision | F1 Score | Time (S)
---|---|---|---|---|---|---
Small-World $N=500$ $M=1000$ | MPGraph (P) | 1373 | 0.872 | 0.635 | 0.735 | 31.542
GL huge eBIC | 4872 | 0.996 | 0.204 | 0.339 | 135.399
NS pulsar BStARS (P) | 2767 | 0.999 | 0.361 | 0.530 | 41.388
GL BigQuic eBIC (P) | 1697 | 0.822 | 0.484 | 0.610 | 76.097
GL glasso eBIC | 3910 | 0.956 | 0.245 | 0.389 | 3217.082
GL skggm MA (P) | 1009 | 0.665 | 0.659 | 0.662 | 3794.389
Erdos-Renyi $N=500$ $M=1000$ | MPGraph (P) | 1058 | 0.831 | 0.776 | 0.803 | 39.863
GL huge eBIC | 7079 | 0.994 | 0.139 | 0.243 | 120.357
NS pulsar BStARS (P) | 3022 | 0.999 | 0.327 | 0.492 | 41.488
GL BigQuic eBIC (P) | 2186 | 0.744 | 0.336 | 0.463 | 73.211
GL glasso eBIC | 4592 | 0.921 | 0.198 | 0.326 | 1716.447
GL skggm MA (P) | 763 | 0.484 | 0.626 | 0.546 | 4010.957
Chain Graph $N=500$ $M=1000$ | MPGraph (P) | 1039 | 0.993 | 0.955 | 0.974 | 26.001
GL huge eBIC | 3719 | 1.000 | 0.269 | 0.423 | 154.775
NS pulsar BStARS (P) | 1785 | 1.000 | 0.560 | 0.718 | 42.757
GL BigQuic eBIC (P) | 2074 | 1.000 | 0.482 | 0.650 | 106.036
GL glasso eBIC | 2929 | 1.000 | 0.341 | 0.509 | 1878.841
GL skggm MA (P) | 1394 | 0.998 | 0.715 | 0.833 | 1437.701
Chain Graph $N=2500$ $M=5000$ | MPGraph (P) | 4659 | 0.920 | 0.987 | 0.952 | 901.814
GL huge eBIC | 14686 | 1.000 | 0.340 | 0.508 | 28731.445
NS pulsar BStARS (P) | 17716 | 1.000 | 0.282 | 0.440 | 3062.926
GL BigQuic eBIC (P) | 11490 | 1.000 | 0.435 | 0.606 | 5815.430
GL glasso eBIC | - | - | - | - |  $>8$ hours
GL skggm MA (P) | 7597 | 1.000 | 0.658 | 0.794 | 17011.503
Chain Graph $N=5000$ $M=10000$ | MPGraph (P) | 7840 | 0.780 | 0.995 | 0.875 | 5183.570
GL huge eBIC | - | - | - | - |  $>8$ hours
NS pulsar BStARS (P) | 61372 | 1.000 | 0.163 | 0.280 | 26737.656
GL BigQuic eBIC (P) | 24321 | 1.000 | 0.411 | 0.583 | 27157.235
GL glasso eBIC | - | - | - | - |  $>8$ hours
GL skggm MA (P) | - | - | - | - |  $>8$ hours

###  3.3 Large-Scale Comparative Studies

In this section, we consider the problem of estimating the structure of large-
scale Gaussian graphical models with the number of nodes ranging from
thousands to tens-of-thousands, similar to the experiments in Hsieh et al.
(2014). Problems of such scale are often encountered in neuroscience and
genomics, among many others. To save computational resources, we only include
the most accurate and/or the fastest competitors from Section 3.2 for the
large-scale comparisons. In addition, following similar practices in Hsieh et
al. (2014), we stop any method whose runtime exceeds 8 hours. The same data-
driven hyperparameter tuning schemes as in Sec. 3.2 are used for these
methods.

The experimental results are summarized in Table 1. We see that our MPGraph
method consistently obtains the highest edge selection accuracy across all
scenarios, often leading the next best performing method by a large margin. At
the same time, our MPGraph is computationally the fastest across the board,
and it is extensively faster than the state-of-the-art competitors for huge
graph selection problems, as shown in the chain graph with $M=5000$ and
$M=10000$ scenarios.

##  4 Conclusions

In this work, we have developed a fast computational approach to learn the
structure of Gaussian graphical models named MPGraph that leverages both
minipatch learning and thresholded graph estimators. We empirically
demonstrate that our approach is not only more accurate but also extensively
faster for huge graph selection problems. Additionally, we theoretically
analyze MPGraph showing that it attains finite-sample graph selection
consistency under certain conditions, lending theoretical support to the
strong empirical performance of our method. Future methodological work might
involve developing adaptive sampling schemes, following ideas from Yao and
Allen (2021). In future theoretical work, one could perhaps extend our graph
selection consistency results to provide bounds on the number of falsely
selected edges or perhaps conduct inference and control the false discovery
rate of the selected edges. We would also like to explore memory-efficient
implementation of our approach, possibly considering approximate computing
schemes such as hashing to store intermediate results more efficiently.
Further, we plan to study whether our approach can be extended to other large-
scale graphical models such as the Ising model and the exponential family
graphical model.

## Acknowledgements

The authors acknowledge support from NSF DMS-1554821, NSF NeuroNex-1707400,
and NIH 1R01GM140468.

Gaussian Graphical Model Selection for Huge Data via Minipatch Learning:
Supplementary Materials

Tianyi Yao and Minjie Wang and Genevera I. Allen

##  A Theoretical Analysis

In this section, we prove the theoretical properties of our MPGraph approach
established in Section 2.3. In particular, we show the graphical model
selection consistency of our method.

We first show that our estimate is graph selection consistent on each
minipatch; hence our estimate is consistent aggregating over all minipatches.
For each minipatch, we treat the sampled nodes as observed nodes, and
unsampled nodes as unobserved latent variables, leading to the problem of
graph selection in the presence of latent variables.

Denote $F_{k}$ as the minipatch set for features. The inverse covariance
matrix $\operatorname{\Theta}=\Sigma^{-1}$ can be written as
$\operatorname{\Theta}=\left(\begin{array}[]{cc}\operatorname{\Theta}_{F_{k}}&\operatorname{\Theta}_{F_{k}F_{k}^{c}}\\\
\operatorname{\Theta}_{F_{k}^{c}F_{k}}&\operatorname{\Theta}_{F_{k}^{c}}\end{array}\right).$
The marginal concentration matrix $(\Sigma^{*}_{F_{k}})^{-1}$ associated with
the sampled features $X_{F_{k}}$, is given by the Schur complement:

$\left(\Sigma_{F_{k}}^{*}\right)^{-1}=\operatorname{\Theta}_{F_{k}}^{*}-\operatorname{\Theta}_{F_{k}F_{k}^{c}}^{*}\left(\operatorname{\Theta}_{F_{k}^{c}}^{*}\right)^{-1}\operatorname{\Theta}_{F_{k}^{c}F_{k}}^{*},$

We restate the assumptions we have made in Section 2.3. These assumptions
guarantee that the assumptions in Wang and Allen (2021) are satisfied for each
minipatch. First, denote the true edge set by
$E^{*}=E(\Theta^{*})=\\{(i,j):\Theta^{*}_{ij}\neq 0\>\text{and}\>i\neq j\\}$
and $s=\operatorname*{arg\,max}_{k}s_{k}$ where
$s_{k}=|E\left(\Theta_{F_{k}}^{*}\right)|$ refers to the total number of non-
zero edges for the $k$th minipatch $F_{k}$. Also let $\varphi_{\max}()$ and
$\varphi_{\min}()$ denote the maximum and minimum eigenvalue respectively.
Denote the elementwise $\ell_{\infty}$-norm of a matrix by
$\|\mathbf{A}\|_{\infty}=\max_{i\neq j}|A_{ij}|$. Consider the following
assumptions:

######  Assumption 1 (restate = assumptionMPGone, name = )

$X_{i}\text{ be i.i.d. }\mathcal{N}\left(\mathbf{0},\Sigma^{*}\right)$.

######  Assumption 2 (restate = assumptionMPGtwo, name = )

$\varphi_{\min}\left(\Sigma^{*}\right)\geq\underline{\kappa}>0,$ or
equivalently $\varphi_{\max}\left(\operatorname{\Theta}^{*}\right)\leq
1/\underline{\kappa}$.

######  Assumption 3 (restate = assumptionMPGthree, name = )

$\varphi_{\max}\left(\Sigma^{*}\right)\leq\overline{\kappa}$.

######  Assumption 4 (restate = assumptionMPGfour, name = )

Minimum edge strength:

$\theta_{\min}:=\min_{(i,j)\in
E\left(\operatorname{\Theta}^{*}\right)}\left|\operatorname{\Theta}_{ij}^{*}\right|>c_{1}\sqrt{\frac{s\log
m}{n}}.$

######  Assumption 5 (restate = assumptionMPGfiveraw, name = )

Maximum effect of unsampled nodes:

$\max\limits_{F_{k}}\|\left(\operatorname{\Theta}_{F_{k}}^{*}-\operatorname{\Theta}_{F_{k}F_{k}^{c}}^{*}\left(\operatorname{\Theta}_{F_{k}^{c}}^{*}\right)^{-1}\operatorname{\Theta}_{F_{k}^{c}F_{k}}^{*}\right)^{-1}-\left(\operatorname{\Theta}_{F_{k}}^{*}\right)^{-1}\|_{\infty}\leq\|\hat{\Sigma}_{F_{k}}-\Sigma_{F_{k}}^{*}\|_{\infty}$

.

First, we have the following lemma on the boundedness of the eigenvalues of
the covariance matrix of each minipatch.

###### Lemma 1

If Assumption 2 and 3 hold true, then
$\varphi_{\min}\left(\Sigma_{F_{k}}^{*}\right)\geq\underline{\kappa}>0,$ and
$\varphi_{\max}\left(\Sigma_{F_{k}}^{*}\right)\leq\overline{\kappa}$ for all
minipatches $F_{k}$.

Proof of Lemma 1: Denote minimal eigenvalue $\varphi_{\text{min}}$ and maximal
eigenvalue $\varphi_{\max}$. Then we have

$\displaystyle\varphi_{\text{min }}(\Sigma^{*})$
$\displaystyle=\min_{v^{T}v=1}v^{T}\Sigma^{*}v$
$\displaystyle\varphi_{\operatorname{max}}(\Sigma^{*})$
$\displaystyle=\max_{v^{T}v=1}v^{T}\Sigma^{*}v.$

Note that the matrix $X_{F_{k}}^{T}X_{F_{k}}$ is a submatrix of $X^{T}X$.
Namely, $X_{F_{k}}^{T}X_{F_{k}}$ consists of the rows and columns of $X^{T}X$
indexed by the elements of $F_{k}$. Define the projection:

$\widetilde{.}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{m},y=\left(y_{1},\ldots,y_{p}\right)^{T}\mapsto\tilde{y}=\left(y_{F_{1}},\ldots,y_{F_{k}}\right)^{T}.$

Moreover, define the set of column vectors:

$V_{F_{k}}=\left\\{x=\left(x_{1},\ldots,x_{p}\right)^{T}\in\mathbb{R}^{p}\mid
x^{T}x=1\text{ and }\forall j\notin F_{k}:x_{j}=0\right\\}.$

Note that

$\widetilde{V}_{F_{k}}=\left\\{y\in\mathbb{R}^{m}\mid y^{T}y=1\right\\}.$

Furthermore, by construction, we have:

$\forall v\in
V_{F_{k}}:v^{T}X^{T}Xv=(\tilde{v})^{T}X_{F_{k}}^{T}X_{F_{k}}\tilde{v}.$

Combining all these, we get

$\displaystyle\varphi_{\min}(\Sigma^{*}_{F_{k}})=\varphi_{\min}\left(X_{F_{k}}^{T}X_{F_{k}}\right)$
$\displaystyle=\min_{\left\\{y\in\mathbb{R}^{m}\mid
y^{T}y=1\right\\}}y^{T}X_{F_{k}}^{T}X_{F_{k}}y$
$\displaystyle=\min_{\tilde{v}\in\widetilde{V}_{F_{k}}}(\tilde{v})^{T}X_{F_{k}}^{T}X_{F_{k}}\tilde{v}$
$\displaystyle=\min_{v\in V_{F_{k}}}v^{T}X^{T}Xv$
$\displaystyle\geq\min_{\left\\{v\in\mathbb{R}^{p}\mid
v^{T}v=1\right\\}}v^{T}X^{T}Xv$
$\displaystyle=\varphi_{\min}\left(X^{T}X\right).$

and similarly,

$\displaystyle\varphi_{\max}(\Sigma^{*}_{F_{k}})=\varphi_{\max}\left(X_{F_{k}}^{T}X_{F_{k}}\right)$
$\displaystyle=\max_{\left\\{y\in\mathbb{R}^{m}\mid
y^{T}y=1\right\\}}y^{T}X_{F_{k}}^{T}X_{F_{k}}y$ $\displaystyle=\max_{v\in
V_{F_{k}}}(\tilde{v})^{T}X_{F_{k}}^{T}X_{F_{k}}\tilde{v}$
$\displaystyle=\max_{v\in V_{F_{k}}}v^{T}X^{T}Xv$
$\displaystyle\leq\max_{\left\\{v\in\mathbb{R}^{p}\mid
v^{r}v=1\right\\}}v^{T}X^{T}Xv$
$\displaystyle=\varphi_{\max}\left(X^{T}X\right),$

as desired.

$\hfill\blacksquare$

We now prove graph selection consistency by first proving that we have graph
selection consistency for each minipatch.

###### Lemma 2

Let Assumptions 1-5 be satisfied. For each minipatch $k$, denote the feature
set $F_{k}$. Then for each minipatch $F_{k}$, the minipatch graph selection
estimator with $\lambda\asymp\sqrt{\frac{\log m}{n}}$ is graph selection
consistent:

$\displaystyle\mathbb{P}\left(\operatorname{sign}(\widetilde{\operatorname{\Theta}}^{(k)}_{ij})=\operatorname{sign}\left(\operatorname{\Theta}_{ij}^{*}\right),\forall\widetilde{\operatorname{\Theta}}^{(k)}_{ij}\in\widetilde{\operatorname{\Theta}}^{(k)},i,j\in
F_{k}\right)\geq 1-b_{1}\exp(-b_{2}n\lambda^{2})\geq 1-m^{\tau}\to 1,\hskip
17.07164pt\forall i\neq j.$

Proof of Lemma 2: Lemma 1 suggests that, under Assumption 2 and 3,
$\varphi_{\min}\left(\Sigma_{F_{k}}^{*}\right)\geq\underline{\kappa}>0,$ and
$\varphi_{\max}\left(\Sigma_{F_{k}}^{*}\right)\leq\overline{\kappa}$ for all
minipatches $F_{k}$. Note that Assumption 5 is a direct extension of
Assumption 5 of Wang and Allen (2021). We replace Assumption 5 of Wang and
Allen (2021) with Assumption 5 and the thresholded graphical lasso estimator
is still graph selection consistent by noting that the quantity $\mathrm{I}$
in the proof of Lemma 2 of Wang and Allen (2021) now becomes:

$\displaystyle\mathrm{I}$ $\displaystyle\leq\max_{i\neq
j}\left|\hat{\sigma}_{ij}-\sigma_{sij}\right|\cdot\left\|\Delta^{-}\right\|_{1}\leq\max_{i\neq
j}\left|(\hat{\sigma}_{ij}-\sigma_{0ij})+(\sigma_{0ij}-\sigma_{sij})\right|\cdot\left\|\Delta^{-}\right\|_{1}$
$\displaystyle\leq 2\max_{i\neq
j}\left|\hat{\sigma}_{ij}-\sigma_{0ij}\right|\cdot\left\|\Delta^{-}\right\|_{1}\leq
2C_{1}\sqrt{\frac{\log m}{n}}\left\|\Delta^{-}\right\|_{1},$

where the second last inequality holds true due to Assumption 5. Here, in our
setup, $\hat{\sigma}_{ij}$ refers to the entries of $\hat{\Sigma}_{F_{k}}$,
$\sigma_{0ij}$ refers to the entries of $\Sigma^{*}_{F_{k}}$, $\sigma_{sij}$
refers to the entries of $(\Theta^{*}_{F_{k}})^{-1}$ and
$\Delta=\widehat{\Theta}_{F_{k}}-\Theta^{*}_{F_{k}}$. The rest of the proof
follows Wang and Allen (2021) paper. Therefore, Assumption 5 still guarantees
that the thresholded graphical lasso estimator is graph selection consistent
for the $k$th minipatch $F_{k}$.

Hence, all the assumptions required for Theorem 3 in the work of Wang and
Allen (2021) are satisfied. Therefore, we have graph selection consistency for
all minipatches $F_{k}$, i.e.,

$\displaystyle\mathbb{P}\left(\operatorname{sign}(\widetilde{\operatorname{\Theta}}^{(k)}_{ij})=\operatorname{sign}\left(\operatorname{\Theta}_{ij}^{*}\right),\forall\widetilde{\operatorname{\Theta}}^{(k)}_{ij}\in\widetilde{\operatorname{\Theta}}^{(k)},i,j\in
F_{k}\right)\geq 1-b_{1}\exp(-b_{2}n\lambda^{2})\geq 1-m^{\tau}\to 1,\hskip
17.07164pt\forall i\neq j.$

for all minipatches $F_{k}$.

$\hfill\blacksquare$

Now, let $R_{ij}$ denote the number of times any pair of nodes $i,j$ are
sampled together into minipatches. Define constant $\tau$ (as in Lemma 2) to
satisfy $m^{\tau}=b_{1}\exp(-b_{2}n\lambda^{2})$ where $b_{1}$ and $b_{2}$
depends on $\overline{\kappa}$ and $\underline{\kappa}$; additionally, define
$w_{1}=1-m^{\tau}$ and $w_{2}=m^{\tau}$. Given these, we have the following
theorem on the graphical model selection consistency of our MPGraph approach.

###### Theorem 1

Let Assumptions 1-5 be satisfied and let $n$ grow proportionally with $N$.
Then, the minipatch graph selection estimator, MPGraph, with
$\lambda\asymp\sqrt{\frac{\log m}{n}}$, is graph selection consistent with
high probability:

$\displaystyle\mathbb{P}\left(\hat{E}^{\text{stable}}=E^{*}\right)$
$\displaystyle\geq 1-\exp\left\\{-\min
R_{ij}\cdot\frac{((1-m^{\tau})-\pi_{\text{thr}})^{2}}{2m^{\tau}(1-m^{\tau})}+{\color[rgb]{0,0,0}2\ln
m}\right\\}$ $\displaystyle-\exp\left\\{-\min
R_{ij}\cdot\frac{(\pi_{\text{thr}}-m^{\tau})^{2}}{2m^{\tau}(1-m^{\tau})}+{\color[rgb]{0,0,0}2\ln
m}\right\\}$
$\displaystyle=1-\sum_{l=1}^{2}\exp\left\\{-\min_{i,j}R_{ij}\cdot\frac{(w_{l}-\pi_{\text{thr}})^{2}}{2m^{\tau}(1-m^{\tau})}+2\ln
m\right\\}$ $\displaystyle\to 1,\hskip 17.07164pt\text{ as }K\to\infty,\text{
or as }N\to\infty.$

Proof of Theorem 1: By definition, $\hat{E}^{\text{stable}}=\big{\\{}1\leq
i<j\leq M:\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}}\big{\\}}$. Meanwhile,

$\displaystyle\mathbb{P}(\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}})$
$\displaystyle=\mathbb{P}\left(\frac{\sum_{k=1}^{K}S_{ij}^{(k)}}{\text{max}(1,\sum_{k=1}^{K}D_{ij}^{(k)})}\geq\pi_{\text{thr}}\right)$
$\displaystyle=P\left(\sum_{k=1}^{K}S_{ij}^{(k)}\geq\pi_{\text{thr}}\cdot{\text{max}(1,\sum_{k=1}^{K}D_{ij}^{(k)})}\right).$

We first consider the case of true edges, $(i,j)\in E^{*}$; the case for null
edges will be discussed later. Denote $R_{ij}=\sum_{k=1}^{K}D_{ij}^{(k)}$ as
the number of times both nodes $i$ and $j$ are sampled together into
minipatches. Without loss of generality, we assume that each pair of nodes is
sampled at least once during the MPGraph procedure.

From Lemma 2, we know that, for a true edge $(i,j)\in E^{*}$, $S_{ij}^{(k)}$
is a Bernoulli random variable with success probability at least $1-m^{\tau}$,
where $\tau$ depends on $n$, $\overline{\kappa}$ and $\underline{\kappa}$;
hence, $\sum_{k=1}^{K}S_{ij}^{(k)}\overset{i.i.d.}{\sim}$
Binomial$(R_{ij},1-m^{\tau})$. Here, without loss of generality, we consider
the case when the lower bound of Lemma 2 is attained, i.e., success
probability equals $1-m^{\tau}$; results when the success probability is
greater than $1-m^{\tau}$ still hold.

Now we have:

$\displaystyle\mathbb{P}(\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}})=P\left(\sum_{k=1}^{K}S_{ij}^{(k)}\geq\pi_{\text{thr}}R_{ij}\right)$
$\displaystyle=\sum_{k=\pi_{\text{thr}}R_{ij}}^{R_{ij}}\binom{R_{ij}}{k}(1-m^{\tau})^{k}(m^{\tau})^{R_{ij}-k}.$

We can approximate the summation above using Gaussian tail bound provided that
$R_{ij}\to\infty$, and therefore, $K\to\infty$. By the Central Limit Theorem,
$\sum_{k=1}^{K}S_{ij}^{(k)}\sim
N(R_{ij}(1-m^{\tau}),R_{ij}m^{\tau}(1-m^{\tau}))$.

Using Gaussian tail bound, we have:

$\displaystyle\mathbb{P}(\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}})$
$\displaystyle=\mathbb{P}\left(\sum_{k=1}^{K}S_{ij}^{(k)}\geq\pi_{\text{thr}}R_{ij}\right)=1-\mathbb{P}\left(\sum_{k=1}^{K}S_{ij}^{(k)}<\pi_{\text{thr}}R_{ij}\right)$
$\displaystyle=1-\mathbb{P}\left(\sum_{k=1}^{K}S_{ij}^{(k)}-R_{ij}(1-m^{\tau})<\pi_{\text{thr}}R_{ij}-R_{ij}(1-m^{\tau})\right)$
$\displaystyle\geq
1-\exp\left\\{-\frac{(R_{ij}(1-m^{\tau})-\pi_{\text{thr}}R_{ij})^{2}}{2R_{ij}m^{\tau}(1-m^{\tau})}\right\\}$
$\displaystyle=1-\exp\left\\{-R_{ij}\frac{((1-m^{\tau})-\pi_{\text{thr}})^{2}}{2m^{\tau}(1-m^{\tau})}\right\\}.$

Or equivalently,
$\mathbb{P}(\hat{\Pi}^{(K)}_{ij}\leq\pi_{\text{thr}})\leq\exp\left\\{-R_{ij}\frac{((1-m^{\tau})-\pi_{\text{thr}})^{2}}{2m^{\tau}(1-m^{\tau})}\right\\}$.
Here, we require that $1-m^{\tau}\geq\pi_{\text{thr}}$. By the union bound,
for all true edges, we have:

$\displaystyle\mathbb{P}\left(\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}},\forall(i,j)\in
E^{*}\right)$
$\displaystyle=1-\mathbb{P}\left(\hat{\Pi}^{(K)}_{ij}\leq\pi_{\text{thr}},\exists(i,j)\in
E^{*}\right)$ $\displaystyle\geq 1-m^{2}\exp\left\\{-\min
R_{ij}\cdot\frac{((1-m^{\tau})-\pi_{\text{thr}})^{2}}{2m^{\tau}(1-m^{\tau})}\right\\}$
$\displaystyle=1-\exp\left\\{-\min
R_{ij}\cdot\frac{((1-m^{\tau})-\pi_{\text{thr}})^{2}}{2m^{\tau}(1-m^{\tau})}+{\color[rgb]{0,0,0}2\ln
m}\right\\}.$

Note $\min R_{ij}$ refers to the minimum number of times any pair of nodes are
sampled together into minipatches. We know that when $m\to\infty$,
$m^{\tau}(1-m^{\tau})\to 0$ for $\tau<0$. Further, when $\tau<-1$, we have
$m^{-\tau}>\ln m$ and $-m^{-\tau}+\ln m\to-\infty$. Therefore, $-\min
R_{ij}\cdot\frac{((1-m^{\tau})-\pi_{\text{thr}})^{2}}{2m^{\tau}(1-m^{\tau})}+2\ln
m\to-\infty$ when $R_{ij}\to\infty$ and we have
$\mathbb{P}(\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}},\forall(i,j)\in E^{*})\to
1$.

Similarly, for a null edge $(i,j)\in{E^{*}}^{c}$, $S_{ij}^{(k)}$ is a
Bernoulli random variable with success probability $m^{\tau}$; hence,
$\sum_{k=1}^{K}S_{ij}^{(k)}\overset{i.i.d.}{\sim}$
Binomial$(R_{ij},m^{\tau})$. When $K\to\infty$, by the Central Limit Theorem,
$\sum_{k=1}^{K}S_{ij}^{(k)}\sim N(R_{ij}m^{\tau},R_{ij}m^{\tau}(1-m^{\tau}))$.

Using Gaussian tail bound, we have:

$\displaystyle\mathbb{P}(\hat{\Pi}^{(K)}_{ij}<\pi_{\text{thr}})$
$\displaystyle=\mathbb{P}\left(\sum_{k=1}^{K}S_{ij}^{(k)}<\pi_{\text{thr}}R_{ij}\right)=1-\mathbb{P}\left(\sum_{k=1}^{K}S_{ij}^{(k)}\geq\pi_{\text{thr}}R_{ij}\right)$
$\displaystyle=1-\mathbb{P}\left(\sum_{k=1}^{K}S_{ij}^{(k)}-R_{ij}m^{\tau}\geq\pi_{\text{thr}}R_{ij}-R_{ij}m^{\tau}\right)$
$\displaystyle\geq
1-\exp\left\\{-\frac{(\pi_{\text{thr}}R_{ij}-R_{ij}m^{\tau})^{2}}{2R_{ij}m^{\tau}(1-m^{\tau})}\right\\}$
$\displaystyle=1-\exp\left\\{-R_{ij}\frac{(\pi_{\text{thr}}-m^{\tau})^{2}}{2m^{\tau}(1-m^{\tau})}\right\\}.$

Here, we require that $\pi_{\text{thr}}\geq m^{\tau}$. Similarly, we have:

$\displaystyle\mathbb{P}\left(\hat{\Pi}^{(K)}_{ij}<\pi_{\text{thr}},\forall(i,j)\in{E^{*}}^{c}\right)$
$\displaystyle=1-\mathbb{P}\left(\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}},\exists(i,j)\in{E^{*}}^{c}\right)$
$\displaystyle\geq 1-m^{2}\exp\left\\{-\min
R_{ij}\cdot\frac{(\pi_{\text{thr}}-m^{\tau})^{2}}{2m^{\tau}(1-m^{\tau})}\right\\}$
$\displaystyle=1-\exp\left\\{-\min
R_{ij}\cdot\frac{(\pi_{\text{thr}}-m^{\tau})^{2}}{2m^{\tau}(1-m^{\tau})}+{\color[rgb]{0,0,0}2\ln
m}\right\\}.$

Similarly, we have
$\mathbb{P}(\hat{\Pi}^{(K)}_{ij}<\pi_{\text{thr}},\forall(i,j)\in{E^{*}}^{c})\to
1$.

Combining these two, we have:

$\displaystyle\mathbb{P}\left(\hat{E}^{\text{stable}}=E^{*}\right)$
$\displaystyle=1-P\left(\hat{\Pi}^{(K)}_{ij}<\pi_{\text{thr}},\exists(i,j)\in
E^{*}\text{ or
}\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}},\exists(i,j)\in{E^{*}}^{c}\right)$
$\displaystyle\geq 1-\exp\left\\{-\min
R_{ij}\cdot\frac{((1-m^{\tau})-\pi_{\text{thr}})^{2}}{2m^{\tau}(1-m^{\tau})}+{\color[rgb]{0,0,0}2\ln
m}\right\\}-\exp\left\\{-\min
R_{ij}\cdot\frac{(\pi_{\text{thr}}-m^{\tau})^{2}}{2m^{\tau}(1-m^{\tau})}+{\color[rgb]{0,0,0}2\ln
m}\right\\}$ $\displaystyle\to 1,\hskip 17.07164pt\text{ as }K\to\infty.$

$\hfill\blacksquare$

On the other hand, we can show some other theoretical properties of our
MPGraph such as the familywise error rate. Following the theoretical results
of Yao and Allen (2021), we have the familywise error rate result below.

###### Proposition 1

Assume that Assumptions 1-5 hold. Then as $K\rightarrow\infty$, MPGraph
controls the familywise error rate (FWER):

$\mathbb{P}\left[\left|\hat{E}^{\text{stable }}\cap{E^{*}}^{c}\right|\geq
1\right]\leq m^{\tau}|{E^{*}}^{c}|/\pi_{\text{thr}},$

where $\tau$ refers to the constant in Lemma 2 with
$m^{\tau}=b_{1}\exp(-b_{2}n\lambda^{2})$ where $b_{1}$ and $b_{2}$ depends on
$\overline{\kappa}$ and $\underline{\kappa}$.

##  B Additional Empirical Studies to Investigate Effects of Minipatch Size

Figure 3: Performance of MPGraph for Small-World Graph Simulations with Fixed
$N=500$ and Varying Dimensionality $M$. (A) Edge selection accuracy of MPGraph
using various minipatch sizes (i.e. $m/M$) with oracle tuning approaches (i.e.
assume total number of true edges $|E|$ is known). (B) Edge selection accuracy
of MPGraph using various minipatch sizes (i.e. $m/M$) for the same simulations
in (A), but with data-driven tuning approaches.

In Figure 3A, we empirically demonstrate that the number of nodes in a
minipatch $m$ is the determining factor for how well MPGraph recovers the true
graph. For the challenging small-world graph situations in Figure 3A, the
performance of MPGraph degrades precipitously when $m$ is too small. But with
a sufficiently large $m$ (e.g. at least $5\%$ of $M$), MPGraph achieves strong
performance. Because the total number of true edges $|E|$ is usually unknown
in realistic use cases of graphical model selection in practice, we also show
the edge selection accuracy of MPGraph with data-driven tuning (i.e.
$\hat{E}^{\text{stable}}=\big{\\{}1\leq i<j\leq
M:\hat{\Pi}^{(K)}_{ij}\geq\pi_{\text{thr}}\big{\\}}$ with
$\pi_{\text{thr}}=0.5$) in Figure 3B. As we can see, the performance of
MPGraph with data-driven tuning is close to that of MPGraph with oracle
tuning.

We have more empirical studies investigating how minipatch sizes impact on
edge selection accuracy and computational time of MPGraph. Following the same
simulation procedure outlined in Sec. 3 of the main paper, we generate data
$\mathbf{X}\in\mathbb{R}^{3000\times 500}$ from the groundtruth precision
matrix that has the Erdős-Renyi graph structure. To investigate the effects of
minipatch size on performance, we run MPGraph on a grid of $n$ and $m$ values.
In particular, we use a sequence of $m$ values such that
$m/M\in\\{0.01,0.05,\ldots,0.4\\}$ and then choose $n$ relative to $m$ so that
$m/n\in\\{0.1,\ldots,0.8\\}$. All other parameters are set to the default
values (i.e. $\pi_{\text{thr}}=0.5$ and $K=1000$). Edge selection accuracy in
terms of the F1 Scores and computational time are reported for these minipatch
sizes in Figure 4. We see that MPGraph has stable edge selection accuracy for
a sensible range of $n$ and $m$ values. As a general rule of thumb, we
recommend taking $m$ to be $5\%\sim 10\%$ of $M$ and then picking $n$ relative
to $m$ so that it surpasses the sample complexity of the base thresholded
graph estimator used on the minipatches, which would strike a good balance
between statistical accuracy and computational time in practice.

Figure 4: Effects of Minipatch Size. We demonstrate how edge selection
accuracy and computational time of MPGraph change with different minipatch
sizes in terms of $m/M$ and $m/n$, where $M$ is the total number of nodes. (A)
Edge selection accuracy in terms of F1 Score (with data-driven tuning); (B)
Computational time on $\log_{10}(\text{second})$ scale. We see that our method
has stable edge selection accuracy for a sensible range of $n$ and $m$ values.

## References

  * Banerjee et al. (2008) O. Banerjee, L. E. Ghaoui, and A. d’Aspremont.  Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data.  _Journal of Machine learning research_ , 9(Mar):485–516, 2008.
  * Belloni et al. (2011) A. Belloni, V. Chernozhukov, and L. Wang.  Square-root lasso: pivotal recovery of sparse signals via conic programming.  _Biometrika_ , 98(4):791–806, 2011.
  * Breiman (2001) L. Breiman.  Random forests.  _Machine Learning_ , 45(1):5–32, 2001.
  * Cai et al. (2011) T. Cai, W. Liu, and X. Luo.  A constrained $\ell_{1}$ minimization approach to sparse precision matrix estimation.  _Journal of the American Statistical Association_ , 106(494):594–607, 2011.
  * Chandrasekaran et al. (2012) V. Chandrasekaran, P. A. Parrilo, and A. S. Willsky.  Latent variable graphical model selection via convex optimization.  _Ann. Statist._ , 40(4):1935–1967, 08 2012.  doi: 10.1214/11-AOS949.  URL https://doi.org/10.1214/11-AOS949.
  * Chang et al. (2019) A. Chang, T. Yao, and G. I. Allen.  Graphical models and dynamic latent factors for modeling functional brain connectivity.  In _2019 IEEE Data Science Workshop (DSW)_ , pages 57–63. IEEE, 2019\.
  * Chen and Chen (2008) J. Chen and Z. Chen.  Extended bayesian information criteria for model selection with large model spaces.  _Biometrika_ , 95(3):759–771, 2008.
  * d’Aspremont et al. (2008) A. d’Aspremont, O. Banerjee, and L. El Ghaoui.  First-order methods for sparse covariance selection.  _SIAM Journal on Matrix Analysis and Applications_ , 30(1):56–66, 2008.
  * Duchi et al. (2008) J. Duchi, S. Gould, and D. Koller.  Projected subgradient methods for learning sparse gaussians.  In _Proceedings of the Twenty-Fourth Conference on Uncertainty in Artificial Intelligence_ , UAI’08, page 153–160, 2008.
  * Efron (1982) B. Efron.  _The jackknife, the bootstrap and other resampling plans_.  SIAM, 1982.
  * Fattahi and Sojoudi (2019) S. Fattahi and S. Sojoudi.  Graphical lasso and thresholding: Equivalence and closed-form solutions.  _The Journal of Machine Learning Research_ , 20(1):364–407, 2019.
  * Foygel and Drton (2010) R. Foygel and M. Drton.  Extended bayesian information criteria for gaussian graphical models.  In _Advances in neural information processing systems_ , pages 604–612, 2010.
  * Friedman et al. (2008) J. Friedman, T. Hastie, and R. Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  _Biostatistics_ , 9(3):432–441, 2008.
  * Friston (2011) K. J. Friston.  Functional and effective connectivity: a review.  _Brain connectivity_ , 1(1):13–36, 2011.
  * Hardt et al. (2016) M. Hardt, B. Recht, and Y. Singer.  Train faster, generalize better: Stability of stochastic gradient descent.  In _Proceedings of the 33rd International Conference on Machine Learning - Volume 48_ , ICML’16, page 1225–1234, 2016.
  * Hsieh et al. (2013) C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, P. K. Ravikumar, and R. Poldrack.  Big & quic: Sparse inverse covariance estimation for a million variables.  In _Advances in Neural Information Processing Systems_ , volume 26, 2013.
  * Hsieh et al. (2014) C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, and P. Ravikumar.  Quic: quadratic approximation for sparse inverse covariance estimation.  _Journal of Machine Learning Research_ , 15(1):2911–2947, 2014.
  * Laska and Narayan (2017) J. Laska and M. Narayan.  skggm 0.2.7: A scikit-learn compatible package for Gaussian and related Graphical Models, 2017.
  * Lauritzen (1996) S. L. Lauritzen.  _Graphical models_ , volume 17.  Clarendon Press, 1996.
  * Li et al. (2012) X. Li, T. Zhao, X. Yuan, and H. Liu.  An r package flare for high dimensional linear regression and precision matrix estimation.  _R Package Vigette_ , 2012.
  * Liu and Wang (2017) H. Liu and L. Wang.  TIGER: A tuning-insensitive approach for optimally estimating Gaussian graphical models.  _Electronic Journal of Statistics_ , 11(1):241 – 294, 2017.
  * Liu et al. (2010) H. Liu, K. Roeder, and L. Wasserman.  Stability approach to regularization selection (stars) for high dimensional graphical models.  In _Advances in neural information processing systems_ , pages 1432–1440, 2010.
  * Meinshausen et al. (2006) N. Meinshausen, P. Bühlmann, et al.  High-dimensional graphs and variable selection with the lasso.  _The annals of statistics_ , 34(3):1436–1462, 2006.
  * Müller et al. (2016) C. L. Müller, R. Bonneau, and Z. Kurtz.  Generalized stability approach for regularized graphical models.  _arXiv preprint arXiv:1605.07072_ , 2016.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Ravikumar et al. (2011) P. Ravikumar, M. J. Wainwright, G. Raskutti, B. Yu, et al.  High-dimensional covariance estimation by minimizing $\ell_{1}$-penalized log-determinant divergence.  _Electronic Journal of Statistics_ , 5:935–980, 2011.
  * Rolfs et al. (2012) B. Rolfs, B. Rajaratnam, D. Guillot, I. Wong, and A. Maleki.  Iterative thresholding algorithm for sparse inverse covariance estimation.  In _Advances in Neural Information Processing Systems_ , volume 25, 2012.
  * Srivastava et al. (2014) N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _Journal of Machine Learning Research_ , 15(56):1929–1958, 2014.
  * Stringer and Pachitariu (2019) C. Stringer and M. Pachitariu.  Computational processing of neural recordings from calcium imaging data.  _Current opinion in neurobiology_ , 55:22–31, 2019.
  * Toghani and Allen (2021) M. Toghani and G. I. Allen.  MP-Boost: Minipatch boosting via adaptive feature and observation sampling.  In _2021 IEEE International Conference on Big Data and Smart Computing_. IEEE, 2021.
  * Vinci et al. (2019) G. Vinci, G. Dasarathy, and G. I. Allen.  Graph quilting: graphical model selection from partially observed covariances.  _arXiv preprint arXiv:1912.05573_ , 2019.
  * Wang and Allen (2021) M. Wang and G. I. Allen.  Thresholded graphical lasso adjusts for latent variables: Application to functional neural connectivity.  _arXiv preprint arXiv:2104.06389_ , 2021.
  * Wasserman and Roeder (2009) L. Wasserman and K. Roeder.  High dimensional variable selection.  _Annals of statistics_ , 37(5A):2178, 2009.
  * Yao and Allen (2021) T. Yao and G. I. Allen.  Feature selection for huge data via minipatch learning.  _arXiv preprint arXiv:2010.08529_ , 2021.
  * Yao et al. (2021) T. Yao, D. LeJeune, H. Javadi, R. Baraniuk, and G. I. Allen.  Minipatch learning as implicit ridge-like regularization.  In _2021 IEEE International Conference on Big Data and Smart Computing_. IEEE, 2021.
  * Yatsenko et al. (2015) D. Yatsenko, K. Josić, A. S. Ecker, E. Froudarakis, R. J. Cotton, and A. S. Tolias.  Improved estimation and interpretation of correlations in neural circuits.  _PLoS Comput Biol_ , 11(3):e1004083, 2015.
  * Yuan and Lin (2007) M. Yuan and Y. Lin.  Model selection and estimation in the gaussian graphical model.  _Biometrika_ , 94(1):19–35, 2007.
  * Zhang et al. (2020) J. Zhang, M. Wang, Q. Li, S. Wang, X. Chang, and B. Wang.  Quadratic sparse gaussian graphical model estimation method for massive variables.  In _Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence, IJCAI-20_ , pages 2964–2972, 2020.
  * Zhao et al. (2012) T. Zhao, H. Liu, K. Roeder, J. Lafferty, and L. Wasserman.  The huge package for high-dimensional undirected graph estimation in r.  _The Journal of Machine Learning Research_ , 13(1):1059–1062, 2012.