{'fraction_non_alphanumeric': 0.03171397002969408, 'fraction_numerical': 0.01476072094468614, 'mean_word_length': 4.588519054510371, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Low-Rank Covariance Completion for Graph Quilting with Applications to Functional Connectivity

As a tool for estimating networks in high dimensions, graphical models are commonly applied to calcium imaging data to estimate functional neuronal connectivity, i.e. relationships between the activities of neurons. However, in many calcium imaging data sets, the full population of neurons is not recorded simultaneously, but instead in partially overlapping blocks. This leads to the Graph Quilting problem, as first introduced by (Vinci et.al. 2019), in which the goal is to infer the structure of the full graph when only subsets of features are jointly observed. In this paper, we study a novel two-step approach to Graph Quilting, which first imputes the complete covariance matrix using low-rank covariance completion techniques before estimating the graph structure. We introduce three approaches to solve this problem: block singular value decomposition, nuclear norm penalization, and non-convex low-rank factorization. While prior works have studied low-rank matrix completion, we address the challenges brought by the block-wise missingness and are the first to investigate the problem in the context of graph learning. We discuss theoretical properties of the two-step procedure, showing graph selection consistency of one proposed approach by proving novel L infinity-norm error bounds for matrix completion with block-missingness. We then investigate the empirical performance of the proposed methods on simulations and on real-world data examples, through which we show the efficacy of these methods for estimating functional connectivity from calcium imaging data.


Introduction
Graphical models are a commonly used unsupervised learning technique for estimating sparse conditional dependency structures in multivariate data. Various graphical modeling approaches have been used in many different fields, including neuroscience (Yatsenko et al., 2015), genomics (Allen and Liu, 2013), network biology (Wang et al., 2016), and finance (Talih and Hengartner, 2005) to analyze conditional relationships in high-dimensional settings. There exists a wide array of literature on the theoretical and empirical performance of different classes of graphical models, including Gaussian and other exponential family graphical models (Lauritzen, 1996;Yang et al., 2015). Additionally, many methods have been developed to account for different external effects, such as latent variables (Chandrasekaran et al., 2010;Pfau et al., 2013) and covariates (Cai et al., 2013;Chen et al., 2016).
However, in some contexts, simultaneous joint observations may not be available for all pairs of features in the data. In particular, pairwise observations may only be available for semi-overlapping blocks within the full feature set. This leads to the Graph Quilting problem, in which we seek to estimate a full graphical model when only partially overlapping patches of the full covariance matrix are observed. Graph quilting is most commonly useful in the fields of neuroscience and genomics, in which technological limitations constrain the amount of features that can be measured simultaneously. For example, in genomics, single cell RNA sequencing data often contains dropouts where individual cell responses are not seen due to technical artifacts (Ozsolak and Milos, 2011;Gan et al., 2020), while in neuroscience, modern calcium imaging technology allows for the recording of individual neuronal activity, but only for a small portion of the total brain volume at the same time (Pnevmatikakis et al., 2016).
Recently, (Vinci et al., 2019) outlined the challenges of the Graph Quilting problem and showed that it is possible to not only recover graph edges associated with observed elements of the covariance, but also to recover a superset of edges associated with completely missing entries of the covariance. Their approach, called MAD GQ , fits an 1 -regularized MLE of the observed covariance with the constraint that no edges are affiliated with unobserved elements. Thresholding and Schur complements are then used to identify graph edges and a minimal superset of edges. While this approach comes with strong theoretical guarantees, it makes several assumptions that may be uncheckable in practice. Additionally, the method contains multiple steps with tuning parameters that may lead to suboptimal performance in practice.
In this paper, we study an alternative approach that (Vinci et al., 2019) briefly suggested but chose not to explore: low-rank covariance completion followed by graph learning. First, a covariance matrix is imputed by using a covariance completion technique to estimate the unobserved entries of the empirical covariance matrix; the graphical Lasso (Yuan and Lin, 2007;Friedman et al., 2008) is then applied to the imputed covariance matrix in order to obtain an estimated graph. Our goal for this work is to investigate the potential efficacy of this two step procedure for estimating sparse conditional dependency structures in the presence of missingness in the full covariance matrix. In particular, we explore different possible covariance imputation techniques which can be used in the first part of the two-step graph quilting procedure. While many covariance completion techniques have been proposed previously in the literature, selection of a covariance completion technique specifically for the purpose of the Graph Quilting problem has not previously been studied; this presents a unique challenge, as most of the previous work on covariance completion focuses on accurate imputation of the missing values in the covariance matrix rather than in recovering the non-zero elements of the true inverse covariance matrix. Thus, the choice of covariance completion method that will achieve optimal graph recovery of the full feature set in the second step of the graph quilting procedure is not readily apparent.

Motivating Example: Calcium Imaging
One context for which the Graph Quilting problem can arise is in the analysis of data from calcium imaging, which records in vivo firing activity of individual neurons in the brain of an experimental subject under controlled or natural stimulus conditions (Stosiek et al., 2003). Calcium imaging data can be used for the study of functional connectivity, defined as the statistical relationships between Figure 1: Eigenvalues of covariance matrices of calcium imaging data sets from Allen Institute (left) and Janelia Research Campus (right).
the activity of neurons in the brain (Horwitz, 2003). In this paper, we focus on the study of intrinsic functional neuronal connectivity, i.e. contemporaneous correlations in the firing activity between pairs of individual neurons (Feldt et al., 2011). Modern calcium imaging technology allows for the recording of the activity of up to thousands of individual neurons simultaneously in vivo; because of this high-dimensional setting, graphical models are a natural choice for estimating intrinsic functional neuronal connectivity. However, in many calcium imaging experiments, multiple scans are used to record the firing activity of neurons within a full brain volume of interest (Grienberger and Konnerth, 2012); these scans are often taken in sequential layers of the brain volume, which leads to partially overlapping blocks of observed neurons between consecutive scans (Berens et al., 2017). Because of this data collection scheme, the simultaneous activity traces between many pairs of neurons in a full calcium imaging data set are never observed. Thus, in order for a functional connectivity network to be estimated for the full brain volume from a calcium imaging data set, a graph quilting procedure would need to be applied.
One potentially useful characteristic found in high-dimensional calcium imaging data sets that can be used in the graph quilting procedure is the low-rankness in the functional activity recording data. We show two examples of this behavior in Figure 1, in which we plot the eigenvalues from the decomposition of the correlation matrices for the fluorescence traces from two separate calcium imaging data sets, one from the Allen Institute (Lein et al., 2007) and another from the Janelia Research Campus (Stringer et al., 2018). Specifically, we tend to see that the first few principal components appear to explain a substantial proportion of total variance, with exponentially decaying eigenvalues, which suggests an approximate spiked eigenvalue structure (Johnstone, 2001) in the covariance matrix for the data. Thus, we can potentially use a low-rank structural assumption in order to impute covariance matrices from calcium imaging data that can then be used in a graph quilting procedure to estimate a functional neuronal connectivity graph.

Related Literature
Many procedures have been previously proposed to address the problem of covariance completion. Many of these methods impute the missing values within a covariance matrix using a penalized likelihood approach, combined with restricting the solution to the space of positive definite or positive semi-definite Hermitian matrices. One commonly used approach utilizes determinant maximization of the imputed covariance matrix (Grone et al., 1984); this can also be combined with a Lasso penalty for sparse parameter estimation (Zhang et al., 2018). Another family of methods involves minimizing the Frobenius norm between the imputed and observed covariance matrices, along with a penalty term to induce sparsity on (Lounici, 2014) or to maximize the rank (Laurent, 2009) of the former. While these methods have been shown to be useful for a variety of problems, these may not be appropriate for functional connectivity. As discussed above, calcium imaging data tend to exhibit an approximate low-rank structure; however, most of the proposed algorithms discussed here will tend to provide full rank solutions, and thus will likely not be a good approximation.
Several model-based approaches have been utilized in the literature as well. Wohrer et al. (2010) use the Wishart distribution as the assumed marginal distribution of the entries of the full inverse covariance matrix. An iterative search over the space of positive definite submatrices is used in order to impute the missing values of the covariance matrix. Zare et al. (2015) instead utilizes a linear time-invariant model to represent the underlying dynamics which generate the full covariance matrix that is only partially observed in the data. The specific structure and input parameters of the dynamic model can then be used in order to obtain estimates of the missing values of the covariance matrix. Both of these methods rely on very specific assumed distributions for the entries of the full covariance matrix; therefore, while these methods can be suitable in contexts where the models explain the covariance matrix well, they are are not necessarily reliable for calcium imaging data.
Although there have been extensive studies for low-rank matrix completion, most existing works (Candes and Recht, 2012;Candès and Tao, 2010;Gross, 2011;Recht, 2011;Candes and Plan, 2011;Chen et al., 2020) assume each entry to be randomly missing. On the other hand, the Graph Quilting problem assumes a block-wise measurement pattern, in which missingness is systematic and common. Thus, the low-rank matrix completion procedure used for Graph Quilting must be robust with regards to these attributes. Furthermore, the prior works on matrix completion with block-wise missingness all consider different settings from ours, or their theory falls short for our graphical model learning purposes. For example, Cai et al. (2016) assume that certain rows and columns are fully observed without noise, and that the missing entries form one submatrix. While Zhou et al. (2021) consider observing multiple blocks of a symmetric PSD matrix and propose a method with solid theory, they assume each block is sampled randomly which leads to overlaps between any two blocks with high probability, making it inappropriate for our calcium imaging application as our blocks are taken sequentially and only consecutive blocks have overlaps. The most closely related work to our graph quilting setting is Bishop and Yu (2014). However, they only provide a Frobenius norm error bound, which cannot rule out the situation where the imputed covariance matrix has large errors for a small fraction of its entries; even a small number of badly estimated pairwise covariance can still lead to many false positives the graph estimate. Instead, a sufficiently small ∞ -norm error bound for the covariance completion step is required and poses a significant challenge for our graph learning purposes.

Contributions
In this paper, we consider several potential approaches to the graph quilting problem in the case which the full sample covariance matrix is assumed to be either exactly or approximately low-rank, and we study the potential application of these methods for estimating functional connectivity networks from calcium imaging data. All of the methods discussed in this work follow the twostep covariance completion graph quilting framework discussed above. As part of the covariance imputation step of our proposed low-rank graph quilting methods, we incorporate several different low-rank covariance completion methods currently used in the literature. While the performance of these methods have been studied previously in the general context of imputing missing values of a covariance matrix, we consider their potential applicability specifically in the graph quilting setting. In particular, we discuss the general theoretical guarantees of the low-rank graph quilting methods, which would require an entry-wise error bound for the imputed covariance matrix that was not proved in the literature on matrix completion with block-wise missingness. We then show the entrywise error bound requirements hold for one of the imputation methods, leading to graph selection consistency guarantees of the corresponding low-rank graph quilting method. Furthermore, we compare the empirical performance of the different low-rank graph quilting methods, along with the MAD gq method of (Vinci et al., 2019), through a simulation study as well as through two examples using real-world calcium imaging data sets.
The paper is organized as follows. In Section 2, we formally introduce the general two-step graph quilting algorithm and three specific models for low-rank graph quilting; we also discuss theoretical properties as well as practical model selection procedures. We then study the empirical performance of low-rank graph quilting on simulation studies in Section 3. Lastly, in Section 4, we investigate the efficacy of the low-rank graph quilting methods for estimating functional neuronal connectivity from calcium imaging data.

Low-Rank Graph Quilting
We first define the general graph quilting problem, following the process in (Vinci et al., 2019). For any vector u ∈ R p , we define its norms as follows: For any p > 0, we denote the p − 1-dimensional sphere in R p by S p−1 = {u ∈ R p : u 2 = 1}. For any matrix A ∈ R p×q , we denote its norms as follows: the spectral norm A 2 = sup u∈S q−1 Au 2 , the entry-wise infinity error bound A max = max i,j |A i,j |, the Frobenius norm Consider the Gaussian graphical model, where each full random vector x ∈ R p follows N (0, Σ * ), and the primary goal is to recover a sparse inverse covariance matrix, denoted as Θ * = Σ * −1 , which represents the conditional dependency structure between the p features. However, instead of having simultaneous observations for the full set of joint feature pairs, we only have joint observations for K partially overlapping subsets of features, denoted as k ∈ {1, . . . , K}. We denote the set of features observed in each subset k as V k and the corresponding observed data matrix as X (k) ∈ R n k ×|V k | . From this, we define the full pairwise observation set as O = K k=1 V k × V k . Additionally, the observed incomplete sample covariance matrix, denoted as Σ O = { Σ ij : (i, j) ∈ O}, is computed empirically from the available joint observation pairs in {X (k) } K k=1 : Our goal is to obtain a graphical model estimate of the full feature set from the incomplete covariance matrix. One possible framework that can be applied is a two-step process in which the full covariance matrix Σ is imputed with covariance completion methods before the graphical Lasso (Friedman et al., 2008) is applied to obtain an estimated graph Θ G . To impute the missing values of the empirical covariance matrix in the first step of the graph quilting procedure though, additional constraints on its formulation are required. In this paper, we consider the specific case where the full covariance matrix Σ * is assumed to follow a low-rank structure; we call this the low-rank graph quilting problem. We outline the general two-step approach used in the low-rank graph quilting problem in Algorithm 2.1.
As has been pointed out in Liu et al. (2012), the success of the graphical Lasso solely depends on the entry-wise estimation error of the covariance matrix, which translates to Σ − Σ * max where Σ Algorithm 2.1: Two-Step Low Rank Graph Quilting Input: Incomplete observed covariance matrix Σ O , sparsity tuning parameter λ, rank of full covariance matrix r.
(2) Apply the graphical Lasso to the imputed full covariance matrix Σ in order to obtain the estimated graph Θ G : is the output of Algorithm 2.1. More details on the general edge selection consistency of Algorithm 2.1 is deferred to the Appendix Section C. For the covariance completion step, the low-rankness of the covariance matrix can be formulated in several ways. One possibility is to apply an exact low-rank assumption on the covariance matrix, i.e. Σ = L L, L ∈ R p×r , r << p.
This can be used in the case where the underlying full covariance matrix is known to have a sparse number of significant principal components or when an exact rank on the covariance matrix is desired. However, in certain empirical contexts, such as in calcium imaging, we may instead want to assume that the covariance matrix follows an approximate low-rank structure. To do this, we impose a spiked eigenvalue structure on the covariance matrix, i.e.
where L * is also a low-rank positive definite matrix; this formulation gives an approximate low-rank structure of the covariance matrix for small values of σ 2 (Johnstone, 2001). We will assume for the rest of this paper that, when using the spiked covariance model, σ 2 is an unknown parameter to be estimated as well. Below, we introduce three different potential low-rank covariance completion methods for the first step of the graph quilting estimation process to derive a full covariance matrix Σ from the incomplete observed covariance matrix Σ O . Computational procedures for each of the low-rank covariance completion methods for both exact and approximate low-rank covariance matrices are shown in Appendix Sections A and B.

Block Singular Value Decomposition (BSVDgq)
The first method we introduce utilizes the procedure proposed by Bishop and Yu (2014) for covariance completion, which applies sequential blockwise singular value decompositions on overlapping principal submatrices of the full covariance matrix in order to impute the missing values. In particular, the algorithm finds the singular value decomposition for each principal submatrix sequentially, while also performing orthonormal transformations of the overlapping parts of the principal submatrices in between iterations in order to align the rows between the different submatrices. In the graph quilting problem, we can use each of the partially overlapping subsets of observed feature pairs as the principal submatrices that are used for imputation with this method. This can be done with the Σ O matrix in order to achieve an exact low-rank solution, or with ( Σ O − σ 2 I) as the input for an approximately low-rank covariance matrix, where σ 2 is an estimate of σ * 2 . One option of the estimate is σ 2 = median( Σ ii ), i ∈ {1, . . . , p}, which has been proposed for estimating σ * 2 for spiked models in some prior works (Johnstone and Lu, 2009;Cai et al., 2015), and is adopted throughout our experiments. We now show that Σ − Σ * max can be bounded appropriately for the BSVDgq covariance imputation algorithm, leading to graph selection consistency results. Suppose the true covariance matrix Σ * = L * + σ * 2 I, where L * is positive semidefinite and of rank r. For 1 ≤ k ≤ K, let p k = |V k | ≤ p, and define J k = V k ∩ (∪ k−1 j=1 V j ) as the joint of the kth node set with prior sets. Also define the following quantities for the covariance Σ * V k ,V k or L * J k ,J k corresponding to the kth node set V k : (iii) ξ k that quantifies both the signal strength in V k compared to J k , and also the condition ) and a stronger signal strength max in the joint observational set J k would lead to a smaller ξ k . We let ξ = max k ξ k .
The theoretical guarantees for the BSVDgq method is presented as follows: Theorem 1 (Guarantees for BSVDgq). Under Assumptions C.2-C.4 in Appendix C, with probability at least 1 − C K k=1 p −c k , the output Σ of the BSVDgq algorithm with the input Σ O − σ * 2 I satisfies where c, C > 0 are universal constants.
Assumptions C.2-C.4 and the full proof of Theorem 1 can be found in Appendix Section C and D.
Remark 1. Theorem 1 suggests that the block singular value decomposition method can lead to entry-wise consistent estimates for the covariance matrix Σ * , if the sample size for each block is larger than a polynomial of the rank r, effective rank τ k , and log p k where p k is the number of nodes in the kth block. Here, Theorem 1 assumes σ * 2 to be known only for simplicity. Otherwise, we expect the error bound to also depend on | σ 2 − σ * 2 | linearly.
More detailed discussion of the implication of Theorem 1 and the technical novelty of its proof can be found in Appendix Section C.
Corollary 1. Suppose we apply the block SVD algorithm with input Σ O − σ * 2 I as the first step of Algorithm 2.1. If Assumptions C.1-C.4 in Supporting Information C hold, λ 1 α Σ − Σ * max , and for 1 ≤ k ≤ K, is the incoherence parameter defined in Assumption C.1, then with probability at least 1 − C K k=1 p −c k , we achieve exact edge recovery of the graph: Under appropriate conditions and with known σ * 2 , Corollary 1 suggests that the graph consistency can be achieved with high probability if n k ≥ C(d 2 + θ −2 min ) ξ 2K (ξ−1) 2 r + τ k µ k (τ k ∨ log p k ) for 1 ≤ k ≤ K. Compared to the sample size requirement (n ≥ C(d 2 + θ −2 min ) log p) in prior literature (Ravikumar et al., 2011), the additional cost due to the block observational pattern is reflected in the effective rank τ k of each Σ * V k ,V k , rank r, incoherence parameter µ k , and the factor ξ K depending on number of blocks K.

Nuclear Norm Penalization (NNgq) and Low-Rank Covariance Factorization (LRFgq)
Here, we propose two other approaches for low-rank graph quilting which solve two squared loss minimization problems to infer the full covariance matrix in the first step of the graph quilting procedure. The first of these, which we call nuclear norm penalization or NNgq for short, uses the objective of minimizing the sum of a Frobenius norm penalty on the difference between the imputed and observed covariance matrix and a nuclear norm penalty on the imputed covariance matrix (Mazumder et al., 2010;Koltchinskii et al., 2011). In the case of an exact low-rank assumption on the completed covariance matrix, this gives us the objective function while for the approximate low-rank assumption we get A proximal gradient descent algorithm can be used in order to derive estimates from the likelihood, as by the symmetry and convexity of the loss function 1 One other common method for low-rank matrix completion utilizes the low-rank factorization and solves an optimization problem w.r.t. the low-rank factors (Keshavan and Oh, 2009;Wen et al., 2012). We consider a low-rank approximation of the full unobserved covariance matrix such that it can be factorized as Σ = UU , U ∈ R n×r , r << p.
Following this, we minimize of the Frobenius norm between the observed portion of the covariance matrix and the corresponding entries of the imputed full covariance matrix for the first step of the graph quilting procedure described above. For the exact low-rank covariance completion, this gives us while for an approximate low-rank covariance matrix we can use the objective This method can utilize the resulting imputed covariance matrix from either the BSVDgq or NNMgq described above for initialization; estimates from the likelihood are then found using a gradient descent procedure. Both the covariance completion techniques used in the NNgq and LRFgq procedures have been previously studied for the imputation of low-rank matrices from both a theoretical and empirical perspective (Candès and Tao, 2010;Candes and Plan, 2011;Chen et al., 2020;Ma et al., 2018). However, these works have primarily been focused on the case when the missing elements are randomly interspersed throughout the covariance matrix. For the Graph Quilting problem, we instead assume that the non-missing entries are arranged as semi-overlapping blocks, meaning that the missingness is highly patterned. Because of this, it is unclear whether existing empirical studies and theoretical guarantees for these covariance completion methods will apply for the graph quilting problem. We investigate the former in Sections 3 and 4, and we leave the development of theoretical properties as future work.

Practical Model Selection
In order to perform data-driven model selection for any of the low-rank graph quilting methods, a stability selection approach such as the one proposed in (Liu et al., 2010) can be used. In particular, since the second step of the low-rank graph quilting approach uses the imputed covariance matrix as the input to the graphical lasso, the stability selection procedure requires repeated runs on random subsamples of both steps of the procedure. To select the hyperparameters that control the rank of the full covariance matrix in the first step of low-rank graph quilting, several different techniques can be applied depending on the method. For the low-rank covariance factorization and block singular value decomposition methods, the rank can be selected using a penalized maximum likelihood criteria, such as BIC (Burnham and Anderson, 2004); here, one can penalize either the raw rank or the nuclear norm of the resulting imputed covariance matrix. For the nuclear norm penalization approach, one can use a cross-validation procedure in which scattered feature pairs are randomly selected to be removed from the observed set and used to compare the imputation estimates along a chosen regularization path (Mazumder et al., 2010).

Simulation Study
Here, we study the performance of the low-rank graph quilting methods on simulated data sets. We compare each of the low-rank graph quilting methods described in Section 2 with both the exact lowrank and spiked covariance model assumptions, as well as the MAD gq approach proposed by (Vinci et al., 2019) and a basic zero imputation approach which assumes that pairs of features that are not observed simultaneously are conditionally independent. For each simulation, we generate a rank one covariance matrix from the spiked model as introduced in Section 2; this is then used to produce an n × p data matrix from a multivariate Gaussian distribution. A partially observed empirical covariance matrix, structured as K patches, of o features each, is created from the simulated data. Afterwards, we construct imputed empirical covariance matrices from the low-rank and zero imputation graph quilting methods. The graphical Lasso is then applied to each of the imputed empirical covariance matrices in order to obtain graphical models estimates from each respective  graph quilting method. We separately apply the MAD gq algorithm to estimate a graphical model directly from the partially observed empirical covariance matrix. We evaluate the graph quilting approaches based on two different criteria. We compare the low-rank and zero imputation approaches using the Frobenius norm of the difference between the imputed and true covariance matrices, i.e. Σ − Σ F ; we do not compare MAD gq using this metric because it does not produce an estimated covariance matrix. We also compare F-1 scores from the resulting graphical Lasso estimates from the low-rank and zero imputation approach as well as the selected edge set from the MAD gq approach, using the graphical Lasso estimate from the full covariance matrix as the true underlying graph. Hyperparameter values for the nuclear norm and low-rank factorization covariance imputation methods are selected by optimal Frobenius norm scores with respect to the true complete covariance matrix. Hyperparameter selection for the graphical Lasso step, when applicable, is performed using oracle sparsity tuning with respect to the full empirical covariance matrix. For each simulation setting, we run 50 replications and report the mean and standard deviation of the Frobenius norm and F-1 scores metrics.
Tables 1 and 2 contain the results from the simulations. We show the performance of the low-rank graph quilting methods with n = 2000 observations of p = 100 features for patch sizes o = 52, 56, 60, 64, and 68 with K = 2 patches, and for K = 2, 3, 4, 5 and 6 patches with patch sizes o = 60, 40, 30, 24, and 20, respectively. Overall, across all simulations, we see that the imputed covariance matrices from approximate low-rank NNgq have the smallest average Frobenius norm compared to the full covariance matrix, while the BSVDgq and LRFgq methods are generally better in terms of F-1 score. On the other hand, the MAD gq algorithm does not perform well in the low-rank setting. In these results, larger patch sizes, which leads to larger overlap sizes between patches, appear to be correlated with better performance in terms of both the Frobenius norm and F-1 score metrics. We also see that the performance of the NNgq method, especially with the exact low-rank model, is robust to smaller patch sizes compared to the other low-rank graph quilting methods. When there are more patches, estimates of the covariance matrix and graphical model tend to be much less accurate for all methods. For particularly large numbers of patches, we see none of the methods perform particularly well for covariance imputation, with all estimates on average being only slightly better than zero imputation; however, even in this setting the NNgq  Table 2: F-1 scores for edge selection from the low-rank graph quilting methods, the MAD gq procedure, and zero imputation, when compared the edges selected from the graphical Lasso on the full empirical covariance matrix from simulated data. Best performing methods are boldfaced.
with the exact low-rank structural assumption tends to do much better in terms of edge selection compared to all other methods. In Appendix Section E in the Supporting Information, we show additional simulation studies where the number of observed data points and the low-rankness of the true underlying covariance matrix are varied.

Case Study: Calcium Imaging
We now investigate the performances of the low-rank graph quilting methods on two real-world, publicly available calcium imaging data sets, one from the Allen Institute (Lein et al., 2007), and another from the Stringer Lab in the Janelia Research Campus (Stringer et al., 2018). For these real data studies, we use data from a single functional calcium imaging recording and calculate the full empirical covariance matrix for the entire set of observed neurons, along with the corresponding estimated graphical model from applying the graphical Lasso. We then randomly assign neurons into K semi-overlapping patches of o neurons each from which we derive empirical partially observed covariance matrices. As in Section 3, we compare the imputed covariance matrices using our lowrank graph quilting methods and zero imputation to the covariance matrix from the full data set using the Frobenius norm metric and F-1 score metrics, using the edges selected by graphical Lasso estimate from the full data as the underlying truth. The MAD gq algorithm is also applied to estimate a graphical model directly from the partially observed empirical covariance matrix. We adopt the same oracle tuning procedures for all hyperparameters as described in Section 3.

Allen Institute Data
The first calcium imaging data set that we analyze from the Allen Institute contains 227 neurons over 9000 measured time points. In Tables 3a and 3b, we show the results for different sizes of observed blocks o, as well as different numbers of patches K while keeping o × K constant. Specifically, we have K = 2 patches with patch sizes o = 130, 140, 150, 160, and 170, and K = 2, 3, 4, 5 and 6 with o = 150, 100, 75, 60, and 50, respectively. For each setting, we run 50 replications and report the mean and standard error of the Frobenius norm and F-1 score metrics. We see that the approximate low rank NNgq method performs best with respect to the Frobenius norm metric on the raw imputed covariance matrix, while NNgq with the exact low-rank is better with respect to edge selection. We also observe that BSVDgq and LRFgq methods with the approximate lowrank assumption perform relatively poorly with the larger patch sizes, with a less accurate imputed covariance matrix and graph edge selection compared to zero imputation. The MAD gq also performs relatively poorly in almost all scenarios with regards to edge selection even when compared to zero imputation, which matches what we saw in Section 3. As we would expect, all methods are more accurate when the patch sizes are larger, with the NNgq method exhibiting the greatest increases in performance. On the other hand, as the number of patches is increased while keeping the total number of observations the same, all methods perform worse. We see a particularly low accuracy when the number of patches reaches 5 or 6, as the amount of overlap between each pair of patches is less than 10. We show the specific structure of one example of estimated functional connectivity graphs in Figure 2. We see that all methods appear to have few false negatives with respect to the functional connectivity network estimated with full information. However, we see that the BSVDgq and LRFgq methods tend to have more false positive selected edges, while the nuclear norm method has very few. We also validate and compare the selected functional connections using a meta data on neuronal tuning, which was posited in the neuroscience literature to be related to functional connectivity. All of the low-rank graph-quilting methods find functional connections with matching neuronal tuning at a similar rate as the graph estimate with the full data. We show full numerical results and experimental details for neuronal tuning in Appendix Section F. (c) Exact low-ranked NNgq.
Figure 2: Comparison of one example of functional connectivity estimates from fully observed covariance and low-rank graph quilting methods on Allen Institute data with K = 2, o = 150.
Neurons are plotted with respect to their spatial locations and colored by their angular tuning category. This figure appears in color in the electronic version of this article, and color refers to that version. (a) Frobenius norms for the imputed covariance matrices from the low-rank graph quilting and zero imputation methods as compared to the full covariance matrix. (b) F-1 scores for edge selection from the low-rank graph quilting methods, the MAD gq procedure, and zero imputation, when compared the edges selected from the graphical Lasso on the full empirical covariance matrix.

Janelia Data
The second calcium imaging data set we analyze comes from the Janelia Research Campus (Stringer et al., 2018). Here, we use approximately 2000 simultaneously recorded neurons measured over 15000 time points. We study the impact of changing the size of patches o and the number of patches K while keeping the total overlap o × K constant on the performance of the low-rank graph quilting methods. Tables 4a and 4b show the results for K = 2 patches with patch sizes o = 1100, 1150, 1200, 1250, and 1300, as well as for K = 2, 3, 4, 5 and 6 with o = 1200, 800, 600, 480, and 400, respectively. We run 50 replications and report the mean and standard error of the Frobenius norm and F-1 score metrics. For this example, the NNgq methods appear to do best for imputing the full covariance matrix and recovering the same edges as in the graphical model estimate from the full data set. In particular, combining the nuclear norm penalization with spiked covariance matrix model does the best by a wide margin in terms of the Frobenius norm metric, while in many cases the exact low-rank covariance assumption is better for edge selection. On the other hand, the other low-rank graph quilting methods, as well as the MAD gq method, are very inaccurate, often performing worse than the zero imputation method. We then show in Figure 3 a subset of the estimated functional neuronal connectivity graphs estimated from the low-rank graph quilting methods with an exact low-rank assumption, alongside the graph estimate with all joint observations. In summary, the topologies and hub neurons of estimated graphs from low-rank graph quilting are similar to the one estimated using the full data. The most well-connected neurons in each of the graph estimates are also colored; these are known as hub neurons, and are of interest in the study of the functional neuronal architecture of the brain as potential drivers of distinct neuronal units (Liska et al., 2015). We compare the top 25 hub neurons from each the low-rank graph quilting estimate to those found from the full data. The low-rank graph quilting methods tend to find many of the same hub neurons as those found when fully observing the data, at a rate between 70 and 80 percent. This shows that hub neurons can reliably be identified even in the presence of missing pairwise observations by using the low-rank graph quilting procedures.
One other notable conclusion from our empirical studies is that the gradient descent and block singular value decomposition methods tend to produce higher F-1 scores in the simulated data, while the nuclear norm penalization methods are more likely to have higher F-1 scores for both of the real data experiments. We explore this further in Appendix Section G. (c) Exact low-ranked NNgq.
(d) Exact low-ranked LRFgq. Figure 3: Functional connectivity estimates from the low-rank graph quilting methods on the Janelia data, for one iteration with K = 2, o = 1200, visualized for one z-plane. We compare the most well-connected neurons of the functional connectivity estimates from the low-rank graphquilting methods in this subset to the ones found when observing full data set. Hub neurons that match those found in the graph estimate from full joint observations are colored in blue, while those that do not match are colored in red. (b) F-1 scores for edge selection from the low-rank graph quilting methods, the MAD gq procedure, and zero imputation, when compared the edges selected from the graphical Lasso on the full empirical covariance matrix.

Discussion
In this paper, we have introduced three new methods for the graph quilting problem when the full covariance matrix with respect to all features is assumed to be low-rank. These methods are based upon low-rank covariance imputation methods which are applied to an incomplete covariance matrix, in order to obtain an imputed positive semi-definite matrix that are then plugged in to the graphical Lasso. We have shown in both simulated and real data studies that the low-rank graph quilting methods perform better than other existing graph quilting approaches for recovering the edge structure of the graph of the complete data set in the case where the covariance matrix is truly low-rank. While we have studied the applicability of low-rank graph quilting specifically to functional neuronal connectivity from calcium imaging data, these methods could be used beyond this particular data type or the neuroscience applications. Low-rank graph quilting also has the potential to be applied to other brain activity data sets such as EEG, PET, ECoG, and fMRI data, or other domains in which graphical modeling may be desired in the presence of missing pairwise entries, such as single-cell RNA sequencing data. There are also several possible methodological extensions to the low-rank graph quilting problem that can be explored in the future. Given that the pairwise observations have been made asynchronously, it would also be critical and challenging to account for the temporal autocorrelations, latent variables, or known measured covariates when learning the graph. The low-rank covariance imputation methods may also be used for clustering data with block-wise missingness, as part of a spectral clustering-type procedure. From a theoretical standpoint, while we have shown general consistency results of the two-step low-rank graph quilting procedure assuming a consistent low-rank covariance imputation procedure as well as specific results for the BSVDgq method, the theoretical guarantees for the other low-rank graph quilting procedures may be useful to study in the future. In particular, the study of the theoretical performance of the nuclear norm penalization could be of interest. Zhou, D., Cai, T., and Lu, J. (2021). Multi-source learning via completion of block-wise overlapping noisy matrices. arXiv preprint arXiv:2105.10360.

A Low Rank Graph Quilting Algorithms
In this section, we show how estimates for each of the low-rank covariance completion methods in Section 2 of the main text can be obtained computationally.

A.1 Block Singular Value Decomposition
For the block singular value decomposition method, we directly follow the algorithm outlined in (Bishop and Yu, 2014).
(2) Set h as indices of the largest r diagonal elements of Λ.

A.2 Nuclear Norm Penalization
For the nuclear norm penalization method, we minimize the objective function L( this is a composite function of two convex functions, which are the loss function F and the penalty term h(Σ) = ν Σ * . We use a proximal gradient descent algorithm to find the objective solution, which utilizes the singular value threshold function (Cai et al., 2010) for the proximal step.

A.3 Low-Rank Matrix Factorization
In the low-rank matrix factorization method, estimates are attained by minimizing the likelihood . For the initialization of Σ (0) below , we can use low-rank solution of the block singular value decomposition method, i.e. matrix C, or the low-rank decomposition of Σ from the nuclear norm penalization method.

B Approximate Low Rank Graph Quilting Algorithms
In this section, we outline the low-rank graph quilting methods as applied to the spiked covariance matrix model; this requires the estimation of an extra parameter for the diagonal entries of the covariance matrix.
(2) Set h as indices of the largest r diagonal elements of Λ.

C.1 Meta-theorem for Graph Selection Consistency of Algorithm 2.1
In this section, we show theoretical guarantees for edge selection consistency of Algorithm 2.1 with respect the true underlying graph, which will apply regardless of the low-rank covariance completion method used in the first step of the algorithm as long as the imputed covariance matrix Σ is a good estimator for the full sample covariance Σ. We follow the notation in Ravikumar et al. (2011): let We require the following incoherence assumption: Assumption 1 (Incoherence condition). There exists some constant α ∈ (0, 1] such that With these, we can now state the following result for model selection consistency of Algorithm 2.1: Proposition 1 (Graph Estimation Consistency of Algorithm 2.1). Consider Algorithm 2.1 and its output Θ G . If Assumption 2.1 holds, The proof for this proposition follows directly from the Theorem 1 in Ravikumar et al. (2011). In the above result, the influence of the missing entries is reflected in the term Σ − Σ * max , i.e. the error of the imputation step.

C.2 Assumptions for Theorem 1
To establish the entrywise error bound for the imputed covariance using the BSVD algorithm, we require the following assumptions: Assumption 2 (Approximate low-rankness and condition number of each block). For all 1 ≤ k ≤ K, max{λ 1 (L * V k ,V k ), σ * 2 } ≤ Cλ r (L * V k ,V k ) for some constant C > 0. Assumption 3 (Blocks not too different).
The following proposition characterizes the key quantity ξ when the low-rank component L * in the covariance matrix is randomly generated.

C.3 Discussion of Theorem 1
Remark 2 (Technical novelty). Although the estimation procedure for the BSVDgq method, as outlined in Algorithm A.1 in the Supporting Information, has been theoretically studied in Bishop and Yu (2014), their result hinges on a Frobenius norm error bound for Σ O − Σ * , which would be too large for our purpose of bounding the entry-wise error bound Σ − Σ * ∞ . Hence we developed some new proofs, borrowing tools and ideas from spectral norm error bounds for sample covariances (Koltchinskii and Lounici, 2017a) and 2→∞ -norm error bounds for spectral methods with perturbed low-rank matrices (Cape et al., 2019). In addition, the error bound for sample covariances with high probability help us get rid of the conditions imposed upon the random quantity Σ O in Bishop and Yu (2014).
Remark 3 (Dependence on the number of blocks K). One might notice that our error bound and sample size requirement in Assumption 5 has an exponential dependence upon the number of blocks K. This is due to that the BSVDgq algorithm employs a sequential matching step to find the best rotation matrix for each block that aligns with previous blocks, and the estimation errors for each block accumulates in an exponential manner. Such dependence on K also appears in Bishop and Yu (2014), where a Frobenius error bound for the BSVDgq algorithm is provided.

D.1 Proofs of Theorem 1 and Proposition 2
Proof of Theorem 1. Let L * = U * Λ * U * with U * ∈ O p×r and Λ * ∈ R r×r being a diagonal matrix. Define H * = U * Λ 1 2 , H * k = H * V k ,: . To fix the notation for the proofs, here we recall the definition ×r is the matrix containing the top r eigenvectors of Σ V k ,V k − σ * 2 I and Λ k ∈ R r×r is a diagonal matrix consisting of the largest r eigenvalues of Σ V k ,V k − σ * 2 I. Then the steps in Algorithm A.1 for computing Σ can also be written as follows: • Σ = H H where H ∈ R p×r is defined sequentially: We also require the following notion of sub-Gaussian random variables.
Our proof can be summarized as two major steps: (i) bounding min W∈O p k ×r H k W − H * k 2 for each block 1 ≤ k ≤ K, which requires tools from matrix perturbation theory and spectral analysis; (ii) controlling the error induced by the merging step and showing the final bound for H H − HH * max by induction. We first restate the necessary assumptions: Assumption 3 (Approximate low-rankness and condition number of each block). For all 1 ≤ k ≤ K, max{λ 1 (L * V k ,V k ), σ * 2 } ≤ Cλ r (L * V k ,V k ) for some constant C > 0.
Requiring σ * 2 ≤ Cλ r (L * V k ,V k ) is to ensure the approximate low-rankness of each block; while λ 1 (L * V k ,V k ) ≤ Cλ r (L * V k ,V k ) means constant condition number. Assumption 4 (Blocks not too different).
Assumption 5 (Sample size). For all 1 ≤ k ≤ K, Controlling error for each block: The following lemma completes the first step: Lemma 1 (Error bound for each block). Assume that x 1 , . . . , x n λ 1 (Σ * ) , then as long as where µ = p r U * 2 2→∞ is the incoherence parameter. Applying Lemma 1 to each block 1 ≤ k ≤ K and taking a union bound, we have that with probability at least 1 − C K k=1 p −c k , there exist W 1 , . . . , W K ∈ O r×r such that holds for all 1 ≤ k ≤ K. Here we have applied the fact that λ 1 (Σ * V k ,V k ) ≥ Cλ r (L * V k ,V k ) in the last line. In the following we denote max k H k W k − H * k 2→∞ by ε.
Here, the last line above is due to the definition of ε, the induction assumption that (9) holds for k ≤ l − 1, and the fact that for any matrix A with m rows, A 2 ≤ √ m A 2→∞ . Recall the definition of ξ k and ξ in (8), one can show that ξ k ≥ 2 L * J k ,J k F λr(L * J k ,J k ) > 1, and hence the definition of l−1 implies that l−1 ≤ 4ξ l−1 ξ−1 ε. Furthermore, when (6) holds, Assumption 3-5 implies that where the third line is due to the fact that L where the 3rd inequality follows similar arguments to (12), and the 4th and 5th inequalities are caused by the facts that Plugging this bound into (7), we have Also noting the fact that µ k r max and we can further give the following bound: with probability at least 1 − 2 exp{−t 2 k /2}. Let t k = 3− √ 2 4 |J k |, then with probability at least 1 − C K k=1 exp{−c|J k |}, for some constants c, C > 0. Furthermore, one can show that L * V k ,V k max = max j∈V k H * j,: 2 2 , and { H * j,: 2 2 } j∈V k are i.i.d. random variables, each being the sum of r sub-exponential random variables with mean 1 and constant sub-exponential norm. Therefore, we can apply the Bernsteintype inequality for sub-exponential sums (see Vershynin, 2010, Proposition 5.16) and obtain the following: P(| H * i,: 2 2 − r| > Cr log p k ) ≤ 2 exp{−cr log p k }. Taking a union bound over i ∈ V k , we have 2 2 > Cr log p k ) ≤ 2 exp{−cr log p k }.

E Additional Simulation Studies
In this section, we present additional simulation studies. As in the main text, we evaluate the low-rank graph quilting approaches based on the Frobenius norm between the imputed and true covariance matrices, as well as F-1 scores from the resulting graphical Lasso estimates and selected edge set from the MAD gq approach with respect to the graphical Lasso estimate using the full covariance matrix. In Tables 5a and 5b, we show results for changing the number of observations generated during the simulation procedure, specifically with n = 1000, 1500, 2000, 2500, and 3000 for p = 100 and with K = 2 patches of o = 60 features. We see that the average Frobenius norm between the imputed and true full covariance matrices decreases with a larger number of observations. However, the F-1 scores for edge selection do not fluctuate substantially with changing the number of observations, though the standard deviation of the estimates decreases with increasing the number samples. We also show the performance of the low-rank graph quilting methods when the low-rankness of the simulated full covariance matrix changes. We measure this by the ratio of the first eigenvalue of the covariance matrix to the last; specifically, for n = 2000 and p = 100, with K = 2 patches of size o = 60, we look at λ 1 /λ p = 500, 200, 50, 10, and 3; here, a smaller ratio represents a weaker low-rank structure. The Frobenius norm of the imputed covariance matrices decreases when the leading eigenvalue is smaller, though this may be an impact of the generated full covariance matrices rather than the estimation procedure. Also, as we would expect, when the true covariance matrix has less of a clear rank 1 structure, all of the low-rank graph quilting methods perform worse in terms of edge selection for the resulting graphical model estimate. On the other hand, the MAD gq method performs much better. This matches with the assumptions made by the two different types of methods. (b) F-1 scores for edge selection from the low-rank graph quilting methods, the MAD gq procedure, and zero imputation, when compared the edges selected from the graphical Lasso on the full empirical covariance matrix.

F Neuronal Tuning Results
In Table 6, we compare the selected functional connections by each low-rank graph quilting method in the unobserved portion of the covariance matrix by neuronal tuning. As mentioned in the main text, we expect Neurons are categorized by 8 categories of angular tuning and 5 categories of frequency tuning. Generally, the performance of the different low-rank graph quilting methods is relatively similar for finding functional connections with matching angular tuning, while the approximate low-rank assumption appears to generally perform better for frequency tuning. We also see that all of the low-rank graph quilting methods find functional connections between neurons with the same tuning at similar rates as when the full data is observed and at a substantially higher rate compared the MAD gq method.

G Real Data Inspired Simulation
As mentioned in the main text, we see from our empirical studies that the gradient descent and block singular value decomposition methods tend to produce higher F-1 scores in the simulated data, while the nuclear norm penalization methods are more likely to have higher F-1 scores for both of the real data experiments. This can be explained by the structure of the true underlying covariance matrices as well as those of the imputed covariance matrices from each method. Figures 4a and 4b show the eigenvalues from the decomposition on the covariance matrix for the Allen Institute data and from the simulated data. As we can see, the real data have exponentially decaying eigenvalues, while the simulated data look almost exactly rank 1, with all other eigenvalues near 0. In Figures  4c and 4d, we see eigenvalues from the decomposition of the imputed covariance matrices from the nuclear norm penalization and the block singular value decomposition methods. In the former we see an exponentially decaying eigenvalue sequence, and in the latter, we see more of an exact low-rank structure. We explore this further in the simulation study below.
Here, instead of creating a rank-1 covariance matrix as the true fully observed covariance matrix in Section 3, we instead use the empirical covariance matrix from the Allen Institute calcium imaging data from Section 4. Specifically, we take the covariance matrix from the real data as the true underlying covariance matrix, and use it to simulate new data from a multivariate Gaussian distribution. As done with the previous simulations in Section 3, we then create from this simulated data a partially observed covariance matrix with k patches of o features each. With these, we impute the covariance matrices using the low-rank graph quilting methods and zero imputation, then apply the graphical Lasso to them to obtain graphical model estimates. Additionally, we ap-  ply the MAD gq to the partially observed covariance matrix to estimate the edges of the graphical model directly. The Frobenius norm and F-1 score metrics are used to compare the methods. As above, hyperparameter values for the nuclear norm and low-rank factorization covariance imputation methods are selected by optimal Frobenius norm scores with respect to the true full covariance matrix, while hyperparameter selection for the graphical Lasso step, when applicable, is performed using oracle sparsity tuning with respect to the full empirical covariance matrix. In this simulation, we let p = 227 as in the Allen Institute dataset. Tables 7a and 7b show results for n = 10000, 15000, 20000, 25000, and 30000, with k = 2 and o = 150. Tables 8a and  8b for o = 130, 140, 150, 160, and 170 with k = 2 and n = 20000. Tables 9a and 9b show the results for k = 2, 3, 4, 5 and 6, for n = 20000 and o = 150, 100, 75, 60, and 50, respectively. In these simulations, we see that the nuclear norm penalization methods now generally have higher F-1 scores compared the block singular value decomposition and the gradient descent methods for edge recovery of the graph from using the full covariance matrix; this falls more in line with what was seen in the results from the real data studies in Section 4, and the opposite of the results seen in Section 3. Otherwise, we see similar patterns in terms of changes in performance based on changes to simulation parameters that we do in Section 3. (b) F-1 scores for edge selection for the low-rank graph quilting methods compared the edges selected from the graphical Lasso on the fully observed covariance. Table 7: Results from the low-rank graph quilting methods, the MAD gq procedure, and zero imputation compared to results for the full empirical covariance, in new data simulated using the empirical covariance matrix from the Allen Institute data. The sample size n is changed. (a) Frobenius norms for the imputed covariance matrices from the low-rank graph quilting and zero imputation methods as compared to the full covariance matrix. (b) F-1 scores for edge selection for the low-rank graph quilting methods compared the edges selected from the graphical Lasso on the fully observed covariance. Table 8: Results from the low-rank graph quilting methods, the MAD gq procedure, and zero imputation compared to results for the full empirical covariance, in new data simulated using the empirical covariance matrix from the Allen Institute data. The patch size o is changed. (a) Frobenius norms for the imputed covariance matrices from the low-rank graph quilting and zero imputation methods as compared to the full covariance matrix. (b) F-1 scores for edge selection for the low-rank graph quilting methods compared the edges selected from the graphical Lasso on the fully observed covariance. Table 9: Results from the low-rank graph quilting methods, the MAD gq procedure, and zero imputation compared to results for the full empirical covariance, in new data simulated using the empirical covariance matrix from the Allen Institute data. The number of patches k is changed.