{'fraction_non_alphanumeric': 0.03360008207653637, 'fraction_numerical': 0.013311788242536165, 'mean_word_length': 4.670302501454334, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gaussian Graphical Model Selection for Huge Data via Minipatch Learning

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.


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 , . . . , X M ) be a M -dimensional random vector following a multivariate Gaussian distribution N (0, Σ) with covariance matrix Σ and inverse covariance matrix Θ = Σ −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 Θ ij = Θ 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 Θ, from the observed data.

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;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.

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 ) and ensemble learning 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 , given observed data matrix X ∈ R N ×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 N and m M .

Minipatch Graph (MPGraph)
Inspired by the minipatch feature selection approaches ) 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, X I k ,F k denotes the submatrix of X containing its rows indexed by I k and its columns indexed by F k . Similarly, Θ i,j denotes the element the at the ith row and jth column. Also, F k,a represents the ath element in the vector F k . For simplicity, [N ] denotes the set {1, 2, . . . , N }. MPGraph fits base thresholded graph estimators to many tiny, random minipatches. On the kth minipatch, the node sampling indicator D 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,Π 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Ê stable whose selection frequencies are above a user-specific threshold π thr ∈ (0, 1). We discuss the choice of π 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: X ∈ R N ×M , n, m, π thr ∈ (0, 1) for k = 1, 2, . . . , K do // In parallel 1) Sample a minipatch: subsample n observations I k ⊂ [N ] and m features F k ⊂ [M ] uniformly at random without replacement to get a minipatch X I k ,F k ∈ R n×m 2) Fit base thresholded graph estimator to X I k ,F k to obtain estimated subgraph Θ 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 λ 0 to get an initial subgraph estimateΘ λ 0 k , as suggested by Wang and Allen (2021). After that, it hard-thresholds the initial estimate at a sequence of threshold values τ k and eventually outputs the best thresholded subgraph estimate Θ (k) according to the eBIC criterion. γ is a user-specific parameter for controlling the amount of regularization from the eBIC criterion. We found setting γ 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: τ k as e τ k end for 3) Choose the best threshold level using eBIC criterion: where is the log-likelihood of the multivariate Gaussian distribution Output: Θ 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.

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(Θ * ) = {(i, j) : Θ * ij = 0 and i = j} and s = arg max k s k where s k = |E Θ * F k | refers to the total number of non-zero edges for the kth minipatch F k . Also let ϕ max () and ϕ min () denote the maximum and minimum eigenvalue respectively. Denote the elementwise ∞ -norm of a matrix by A ∞ = max i =j |A ij |. Consider the following assumptions: Assumption 4 Minimum edge strength: Assumption 5 Maximum effect of unsampled nodes: Now, let R ij denote the number of times any pair of nodes i, j are sampled together into minipatches. Define constant τ to satisfy m τ = b 1 exp(−b 2 nλ 2 ) where b 1 and b 2 depends on κ and κ; additionally, define w 1 = 1 − m τ and w 2 = m τ . 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 λ log m n , is graph selection consistent with high probability: 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 Σ * 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 = 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.

Practical Considerations
Our MPGraph is a general meta-algorithm with three tuning hyperparameters: minipatch size (n and m) and threshold π 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% ∼ 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 π 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 τ , but we strongly suggest tuning these separately for each minipatch (as we outline in Algorithm 2) and not all together across minipatches.

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 Θ ∈ S M ×M + is set to be Θ i,i−1 = Θ i−1,i = 0.6 and Θ i,i = 1.25 for all i = 1, . . . , M .
• Scenario 2: Erdős-Renyi graph: First generate an adjacency matrix A ∈ {0, 1} M ×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 Θ is set to be Θ . Lastly, set Θ i,i = 1 for all i = 1, . . . , M .
• Scenario 3: Small-World graph: First generate an adjacency matrix A ∈ {0, 1} M ×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 Θ is generated in the same way as in the Erdős-Renyi graph scenario.
For a given groundtruth precision matrix Θ, we generate the data matrix X ∈ R N ×M whose rows are independently drawn from a M -variate Gaussian distribution N (0, Θ −1 ). The groundtruth edge set is E = {1 ≤ i < j ≤ M : Θ i,j = 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.

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% ∼ 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.

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 , 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 Λ of length 100 in the range of (δλ max , λ max ) are used for all software packages, where the ratio δ ∈ (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. π 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. 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).

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-ofthe-art competitors for huge graph selection problems, as shown in the chain graph with M = 5000 and M = 10000 scenarios.

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 . 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 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.

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 Θ = Σ −1 can be written as Θ = The marginal concentration matrix (Σ * F k ) −1 associated with the sampled features X F k , is given by the Schur complement: 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(Θ * ) = {(i, j) : Θ * ij = 0 and i = j} and s = arg max k s k where s k = |E Θ * F k | refers to the total number of non-zero edges for the kth minipatch F k . Also let ϕ max () and ϕ min () denote the maximum and minimum eigenvalue respectively. Denote the elementwise ∞ -norm of a matrix by A ∞ = max i =j |A ij |. Consider the following assumptions: Assumption 1 X i be i.i.d. N (0, Σ * ).
Assumption 4 Minimum edge strength: Assumption 5 Maximum effect of unsampled nodes: First, we have the following lemma on the boundedness of the eigenvalues of the covariance matrix of each minipatch.
Proof of Lemma 1: Denote minimal eigenvalue ϕ min and maximal eigenvalue ϕ max . Then we have Note that the matrix X T F k X F k is a submatrix of X T X. Namely, X T F k X F k consists of the rows and columns of X T X indexed by the elements of F k . Define the projection: . : R p → R m , y = (y 1 , . . . , y p ) T →ỹ = (y F 1 , . . . , y F k ) T .
Moreover, define the set of column vectors: Note that Furthermore, by construction, we have: Combining all these, we get and similarly, as desired.
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 λ log m n is graph selection consistent: Proof of Lemma 2: Lemma 1 suggests that, under Assumption 2 and 3, ϕ min Σ * F k ≥ κ > 0, and ϕ max Σ * F k ≤ κ 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 I in the proof of Lemma 2 of Wang and Allen (2021) now becomes: where the second last inequality holds true due to Assumption 5. Here, in our setup,σ ij refers to the entries ofΣ F k , σ 0ij refers to the entries of Σ * F k , σ sij refers to the entries of (Θ * F k ) −1 and ∆ = Θ F k − Θ * 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 kth 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., for all minipatches F k . Now, let R ij denote the number of times any pair of nodes i, j are sampled together into minipatches. Define constant τ (as in Lemma 2) to satisfy m τ = b 1 exp(−b 2 nλ 2 ) where b 1 and b 2 depends on κ and κ; additionally, define w 1 = 1 − m τ and w 2 = m τ . 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 λ log m n , is graph selection consistent with high probability: Proof of Theorem 1: By definition,Ê stable = 1 ≤ i < j ≤ M :Π ij ≥ π thr . Meanwhile, We first consider the case of true edges, (i, j) ∈ E * ; the case for null edges will be discussed later.
ij 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.
Similarly, we have P(Π (K) ij < π thr , ∀(i, j) ∈ E * c ) → 1. Combining these two, we have: 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 → ∞, MPGraph controls the familywise error rate (FWER): where τ refers to the constant in Lemma 2 with m τ = b 1 exp(−b 2 nλ 2 ) where b 1 and b 2 depends on κ and κ. 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 smallworld 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.Ê stable = 1 ≤ i < j ≤ M :Π

B Additional Empirical Studies to Investigate Effects of Minipatch Size
ij ≥ π thr with π 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 X ∈ R 3000×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 ∈ {0.01, 0.05, . . . , 0.4} and then choose n relative to m so that m/n ∈ {0.1, . . . , 0.8}. All other parameters are set to the default values (i.e. π 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% ∼ 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.