High-dimensional Gaussian graphical model for network-linked data

TL;DR

This paper introduces a Gaussian graphical model tailored for network-linked high-dimensional data, accounting for dependencies and smoothly varying means across the network, with proven estimation accuracy and practical effectiveness.

Contribution

It develops a novel Gaussian graphical model for network-connected data with varying means, along with an efficient estimation algorithm and theoretical guarantees.

Findings

Effective on simulated data

Meaningful results on real coauthorship network

Accurate estimation of graph structure

Abstract

Graphical models are commonly used to represent conditional dependence relationships between variables. There are multiple methods available for exploring them from high-dimensional data, but almost all of them rely on the assumption that the observations are independent and identically distributed. At the same time, observations connected by a network are becoming increasingly common, and tend to violate these assumptions. Here we develop a Gaussian graphical model for observations connected by a network with potentially different mean vectors, varying smoothly over the network. We propose an efficient estimation algorithm and demonstrate its effectiveness on both simulated and real data, obtaining meaningful and interpretable results on a statistics coauthorship network. We also prove that our method estimates both the inverse covariance matrix and the corresponding graph structure correctly under the assumption of network “cohesion”, which refers to the empirically observed phenomenon of network neighbors sharing similar traits.

Tables

Table 1. Table 1 : Mean estimation errors for the four methods, averaged over 100 replications, with the lowest error in each configuration indicated in bold.

		${\Vert \hat{M}-M\Vert }_{\infty }$			${\Vert \hat{M}-M\Vert }_{2,\infty }$
network	method	$t=0.1$	0.5	1	$t=0.1$	0.5	1
lattice	glasso	0.358	0.746	1.037	5.819	12.985	18.357
	oracle cluster+glasso	0.493	0.539	0.565	3.293	3.639	4.118
	oracle GNC-lasso	0.328	0.520	0.526	2.054	3.247	3.287
	GNC-lasso	0.419	0.669	0.820	2.619	4.105	4.874
coauthorship	glasso	1.540	3.401	4.795	25.072	57.655	78.426
	oracle cluster+glasso	0.724	1.077	1.342	7.051	13.078	16.962
	oracle GNC-lasso	0.710	0.860	0.917	6.400	7.404	7.420
	GNC-lasso	0.717	0.878	0.942	6.430	7.037	7.436

Table 2. Table 2 : Mean estimation errors for the four methods, averaged over 100 replications, with the lowest error in each configuration indicated in bold, for the lattice networks with one, four, and 16 connected components.

	${\Vert \hat{M}-M\Vert }_{\infty }$			${\Vert \hat{M}-M\Vert }_{2,\infty }$
method	original	4 comp.	16 comp.	original	4 comp.	16 comp.
glasso	0.746	2.801	2.942	5.819	40.25	25.16
oracle cluster+glasso	0.539	1.091	1.099	3.293	12.20	8.22
oracle GNC-lasso	0.520	0.866	0.785	2.054	6.46	5.13
GNC-lasso	0.669	0.983	0.838	2.619	6.73	5.79

Table 3. Table 3 : Mean estimation errors for the four methods, averaged over 100 replications, with the lowest error in each configuration indicated in bold, for the coauthorship networks with one, two, or four components.

	${\Vert \hat{M}-M\Vert }_{\infty }$			${\Vert \hat{M}-M\Vert }_{2,\infty }$
method	original	2 comp.	4 comp.	original	2 comp.	4 comp.
glasso	0.746	1.949	4.019	5.819	21.214	32.307
oracle cluster+glasso	0.539	0.936	1.676	3.293	6.033	7.208
oracle GNC-lasso	0.520	0.659	0.958	2.054	3.860	4.843
GNC-lasso	0.669	0.852	1.289	2.619	4.947	5.102

Table 4. Table 4 : The estimation errors of M 𝑀 M from the four methods on three connected lattice networks with varying sample size, averaged over 100 independent replications. The network sizes are 100, 225 and 400, corresponding to lattice dimension 10 × 10 10 10 10\times 10 , 15 × 15 15 15 15\times 15 and 20 × 20 20 20 20\times 20 , respectively.

	${\Vert \hat{M}-M\Vert }_{\infty }$			${\Vert \hat{M}-M\Vert }_{2,\infty }$
method	$n=100$	$225$	$400$	$100$	$225$	$400$
glasso	1.009	1.030	0.746	15.02	15.69	12.985
oracle cluster+glasso	0.991	0.716	0.539	6.62	4.72	3.639
oracle GNC-lasso	0.911	0.794	0.520	5.52	4.83	3.247
GNC-lasso	0.874	0.988	0.669	5.36	5.67	4.105

Table 5. Table 5 : The estimation errors of M 𝑀 M from the iterative and two-stage oracle GNC-lasso methods, averaged over 100 independent replications.

network	method	${\Vert \hat{M}-M\Vert }_{\infty }$	${\Vert \hat{M}-M\Vert }_{2,\infty }$
lattice	two-stage	0.520	3.247
	iterative	0.587	3.639
coauthorship	two-stage	0.860	7.404
	iterative	1.21	9.92

Equations

$${\mu }^{T}L\mu =\sum _{i\sim i^{\prime }}({\mu }_i-{\mu }_{i^{\prime }}{)}^2.$$

$$X_{i\cdot }\sim \mathcal{N}({\mu }_i,\Sigma ),i=1,2,\cdots ,n.$$

$$\ell (M,\Theta )=\log \det (\Theta )-\frac{1}{n}\mathrm{tr}(\Theta (X-M{)}^T(X-M)).$$

$$\hat{M}=\arg \underset{M}{\min }\Vert X-M{\Vert }_F^2+\alpha \mathrm{tr}(M^T{\mathcal{L}}_{s}M).$$

$$\hat{\Theta }=\arg \underset{\Theta \in {\mathbb{S}}_{+}^n}{\min }\log \det (\Theta )-\mathrm{tr}(\Theta \hat{S})-\lambda \Vert \Theta {\Vert }_{1,\text{off}}.$$

$$\mathrm{tr}(M^T{\mathcal{L}}_{s}M)=\sum _{i\sim i^{\prime }}\Vert {\mu }_i-{\mu }_{i^{\prime }}{\Vert }^2.$$

$${\hat{M}}_{\cdot j}=(I_n+\alpha {\mathcal{L}}_s{)}^{-1}{X}_{\cdot j},j=1,2,\cdots ,p.$$

$$(\hat{\Theta },\hat{M})=\arg \underset{\Theta ,M}{\max }\log \det (\Theta )$$

$$\text{GCV}(\alpha )=\frac{1}{np}\Vert X-\hat{M}(\alpha ){\Vert }_F^2/[1-\frac{1}{n}\mathrm{tr}((I+\alpha {\mathcal{L}}_s{)}^{-1}){]}^2=\frac{\Vert X-\hat{M}(\alpha ){\Vert }_F^2}{np[1-\frac{1}{n}{\sum }_{i=1}^n\frac{1}{1+\alpha {\tau }_i}{]}^2}$$

$$\sum _{i\sim i^{\prime }}\Vert {\mu }_i-{\mu }_{i^{\prime }}\Vert .$$

$$\Vert {\mathcal{L}}_{s}v{\Vert }_2\to 0⟺v^T{\mathcal{L}}_{s}v\to 0.$$

$$\frac{d_i}{\bar{d}}\left(v_i-\frac{1}{d_i}\sum _{i^{\prime }{\sim }_{A}i}{v}_i^{\prime }\right),$$

$$\frac{{\tau }_i^2|{\beta }_i{|}^2}{\Vert \beta {\Vert }_2^2}\le n^{-\frac{2(1+\delta )}{3}-1},$$

$$\frac{\Vert {\mathcal{L}}_{s}v{\Vert }_2^2}{\Vert v{\Vert }_2^2}\le n^{-\frac{2(1+\delta )}{3}}.$$

$$\widehat{\mathrm{Var}}(v)=o({\bar{v}}^2)$$

$$m_A=\inf \{m:0\le m\le n-1,{\tau }_{n-m}\ge \frac{1}{\sqrt{m}}\}.$$

$$\frac{\Vert \hat{M}-M{\Vert }_F^2}{np}\le \frac{(b^2+2{\sigma }^2)[1+m_A(\frac{1}{(1+\Delta {)}^2}+n^{\frac{1-2\delta }{3}})]}{n}$$

$$S(\Theta )=\{(j,j^{\prime })\in [n]\times [n]:{\Theta }_{j{j}^{\prime }}\ne 0\}.$$

$$\psi$$

$${\kappa }_{\Sigma }$$

$${\kappa }_{\Gamma }$$

$$\underset{e\in S^c}{\max }\Vert {\Gamma }_{eS}({\Gamma }_{SS}{)}^{-1}{\Vert }_1\le 1-\rho .$$

$$\frac{1}{\bar{k}}\le {\eta }_{\min }(\Sigma )\le {\eta }_{\max }(\Sigma )\le \bar{k}.$$

$$\displaystyle\nu(n,p):=C\sqrt{\frac{\log p}{n}}\max\Big{(}1,m_{A}n^{-\frac{1+4\delta}{6}},\sqrt{m_{A}}n^{\frac{1-2\delta}{6}},\sqrt{\frac{\log{p}}{n}}(\frac{m_{A}}{\Delta+1}+1)(\sqrt{m_{A}}n^{-\frac{1+\delta}{3}}+1)\Big{)}$$

$$\nu (n,p)<\frac{1}{6(1+8/\rho )\psi \max \{{\kappa }_{\Sigma }{\kappa }_{\Gamma },(1+8/\rho ){\kappa }_{\Sigma }^3{\kappa }_{\Gamma }^2\}},$$

$$\Vert \hat{\Theta }-\Theta {\Vert }_{\infty }$$

$$\Vert \hat{\Theta }-\Theta {\Vert }_F$$

$$\Vert \hat{\Theta }-\Theta \Vert$$

$$S(\hat{\Theta })\subset S(\Theta ),$$

$$S(\hat{\Theta })=S(\Theta ).$$

$$\nu (n,p)=C\sqrt{\frac{\log p}{n}},$$

$$\nu (n,p)=C\sqrt{\frac{\log p}{n}}{n}^{\frac{\max (1-2\delta ,0)}{6}}.$$

Full text

Keywords: High-dimensional statistics, Gaussian graphical model, network analysis, network cohesion, statistical learning

1 Introduction

Network data represent information about relationships (edges) between units (nodes), such as friendships or collaborations, and are often collected together with more “traditional” covariates that describe one unit. In a social network, edges may represent friendships between people (nodes), and traditional covariates could be their demographic characteristics such as gender, race, age, and so on. Incorporating relational information in statistical modeling tasks focused on “traditional” node covariates should improve performance, since it offers additional information, but most traditional multivariate analysis methods are not designed to use such information. In fact, most such methods for regression, clustering, density estimation and so on tend to assume the sampled units are homogeneous, typically independent and identically distributed (i.i.d.), which is unlikely to be the case for units connected by a network. While there is a fair amount of work on incorporating such information into specific settings (Manski, 1993; Lee, 2007; Yang et al., 2011; Raducanu and Dornaika, 2012; Vural and Guillemot, 2016), work on extending standard statistical methods to network-linked data has only recently started appearing, for example, Li et al. (2019) for regression, Tang et al. (2013) for classification, and Yang et al. (2013), Binkiewicz et al. (2017) for clustering. Our goal in this paper is to develop an analog to the widely used Gaussian graphical models for network-linked data which takes advantage of this additional information to improve performance when possible.

Graphical models are commonly used to represent independence relationships between random variables, with each variable corresponding to a node, and edges representing conditional or marginal dependence between two random variables. Note that a graphical model is a graph connecting variables, as opposed to the networks discussed above, which are graphs connecting observations. Graphical models have been widely studied in statistics and machine learning and have applications in bioinformatics, text mining and causal inference, among others. The Gaussian graphical model belongs to the family of undirected graphical models, or Markov random fields, and assumes the variables are jointly Gaussian. Specifically, the conventional Gaussian graphical model for a data matrix $X\in\mathbb{R}^{n\times p}$ assumes that the rows $X_{i\cdot}$, $i=1,\dots,n$, are independently drawn from the same $p$-variate normal distribution $\mathcal{N}(\mu,\Sigma)$. This vastly simplifies analysis, since for the Gaussian distribution all marginal dependence information is contained in the covariance matrix, and all conditional independence information in its inverse. In particular, random variables $j$ and $j^{\prime}$ are conditionally independent given the rest if and only if the $(j,j^{\prime})$-th entry of the inverse covariance matrix $\Sigma^{-1}$ (the precision matrix) is zero. Therefore estimating the graph for a Gaussian graphical model is equivalent to identifying zeros in the precision matrix, and this problem has been well studied, in both the low-dimensional and the high-dimensional settings. A pioneering paper by Meinshausen and Bühlmann (2006) proposed neighborhood selection, which learns edges by regressing each variable on all the others via lasso, and established good asymptotic properties in high dimensions. Many penalized likelihood methods have been proposed as well (Yuan and Lin, 2007; Banerjee et al., 2008; Rothman et al., 2008; d’Aspremont et al., 2008; Friedman et al., 2008). In particular, the graphical lasso (glasso) algorithm of Friedman et al. (2008) and its subsequent improvements (Witten et al., 2011; Hsieh et al., 2013b) are widely used to solve the problem efficiently.

The penalized likelihood approach to Gaussian graphical models assumes the observations are i.i.d., a restrictive assumption in many real-world situations. This assumption was relaxed in Zhou et al. (2010); Guo et al. (2011) and Danaher et al. (2014) by allowing the covariance matrix to vary smoothly over time or across groups, while the mean vector remains constant. A special case of modeling the mean vector on additional covariates associated with each observation has also been studied (Rothman et al., 2010; Yin and Li, 2011; Lee and Liu, 2012; Cai et al., 2013; Lin et al., 2016). Neither of these relaxations are easy to adapt to network data, and their assumptions are hard to verify in practice.

In this paper, we consider the problem of estimating a graphical model with heterogeneous mean vectors when a network connecting the observations is available. For example, in analyzing word frequencies in research papers, the conditional dependencies between words may represent certain common phrases used by all authors. However, since different authors also have different research topics and writing styles, there is individual variation in word frequencies themselves, and the coauthorship information is clearly directly relevant to modeling both the universal dependency graph and the individual means. We propose a generalization of the Gaussian graphical model to such a setting, where each data point can have its own mean vector but the data points share the same covariance structure. We further assume that a network connecting the observations is available, and that the mean vectors exhibit network “cohesion”, a generic term describing the phenomenon of connected nodes behaving similarly, observed widely in empirical studies and experiments (Fujimoto and Valente, 2012; Haynie, 2001; Christakis and Fowler, 2007). We develop a computationally efficient algorithm to estimate the proposed Gaussian graphical model with network cohesion, and show that the method is consistent for estimating both the covariance matrix and the graph in high-dimensional settings under a network cohesion assumption. Simulation studies show that our method works as well as the standard Gaussian graphical model in the i.i.d. setting, and is effective in the setting of different means with network cohesion, while the standard Gaussian graphical model completely fails.

The rest of the paper is organized as follows. Section 2 introduces a Gaussian graphical model on network-linked observations and the corresponding two-stage model estimation procedure. An alternative estimation procedure based on joint likelihood is also introduced, although we will argue that the two-stage estimation is preferable from both the computational and the theoretical perspectives. Section 3 presents a formal definition of network cohesion and error bounds under the assumption of network cohesion and regularity conditions, showing we can consistently estimate the partial dependence graph and model parameters. Section 4 presents simulation studies comparing the proposed method to standard graphical lasso and the two-stage estimation algorithm to the joint likelihood approach. Section 5 applies the method to analyzing dependencies between terms from a collection of statistics papers’ titles and the associated coauthorship network. Section 6 concludes with discussion.

2 Gaussian graphical model with network cohesion

2.1 Preliminaries

We start with setting up notation. For a matrix $X\in\mathbb{R}^{n\times p}$, let $X_{\cdot j}$ be the $j$th column and $X_{i\cdot}$ the $i$th row. By default, we treat all vectors as column vectors. Let $\|{X}\|_{F}=(\sum_{i,j}X_{ij}^{2})^{1/2}$ be the Frobenius norm of $X$ and $\|{X}\|$ the spectral norm, i.e., the largest singular value of $X$. Further, let $\|{X}\|_{0}=\#\{(i,j):X_{ij}\neq 0\}$ be the number of non-zero elements in $X$, $\|{X}\|_{1}=\sum_{ij}|X_{ij}|$, and $\|{X}\|_{1,\text{off}}=\sum_{i\neq j}|X_{ij}|$. For a square matrix $\Sigma$, let $\mathrm{tr}(\Sigma)$ and $\det(\Sigma)$ be the trace and the determinant of $\Sigma$, respectively, and assuming $\Sigma$ is a covariance matrix, let $r(\Sigma)=\frac{\mathrm{tr}(\Sigma)}{\|{\Sigma}\|}$ be its stable rank. It is clear that $1\leq r(\Sigma)\leq p$ for any nonzero covariance matrix $\Sigma$.

While it is common, and not incorrect, to use the terms “network” and “graph” interchangeably, throughout this paper “network” is used to refer to the observed network connecting the $n$ observations, and “graph” refers to the conditional dependence graph of $p$ variables to be estimated. In a network or graph $\mathcal{G}$ of size $n$, if two nodes $i$ and $i^{\prime}$ of $\mathcal{G}$ are connected, we write $i\sim_{\mathcal{G}}i^{\prime}$, or $i\sim i^{\prime}$ if $\mathcal{G}$ is clear from the context. The adjacency matrix of a graph $\mathcal{G}$ is an $n\times n$ matrix $A$ defined by $A_{ii^{\prime}}=1$ if $i\sim_{\mathcal{G}}i^{\prime}$ and 0 otherwise. We focus on undirected networks, which implies the adjacency matrix is symmetric. Given an adjacency matrix $A$, we define its Laplacian by $L=D-A$ where $D=\mathrm{diag}(d_{1},d_{2},\cdots,d_{n})$ and $d_{i}=\sum_{i^{\prime}=1}^{n}A_{ii^{\prime}}$ is the degree of node $i$. A well-known property of the Laplacian matrix $L$ is that, for any vector $\mu\in\mathbb{R}^{n}$,

[TABLE]

We also define a normalized Laplacian $\mathcal{L}_{s}=\frac{1}{\bar{d}}L$ where $\bar{d}$ is the average degree of the network $\mathcal{G}$, given by $\bar{d}=\frac{1}{n}\sum_{i}d_{i}$. We denote the eigenvalues of $\mathcal{L}_{s}$ by $\tau_{1}\geq\tau_{2}\geq\cdots\geq\tau_{n-1}\geq\tau_{n}=0$, and the corresponding eigenvectors by $u_{1},\dots,u_{n}$.

2.2 Gaussian graphical model with network cohesion (GNC)

We now introduce the heterogeneous Gaussian graphical model, as a generalization of the standard Gaussian graphical model with i.i.d. observations. Assume the data matrix $X$ contains $n$ independent observations $X_{i\cdot}\in\mathbb{R}^{p},i=1,2,\cdots,n$. Each $X_{i\cdot}$ is a random vector drawn from an individual multivariate Gaussian distribution

[TABLE]

where $\mu_{i}\in\mathbb{R}^{p}$ is a $p$-dimensional vector and $\Sigma$ is a $p\times p$ symmetric positive definite matrix. Let $\Theta=\Sigma^{-1}$ be the corresponding precision matrix and $M=(\mu_{1},\mu_{2},\cdots,\mu_{n})^{T}$ be the mean matrix, which will eventually incorporate cohesion. Recall that in the Gaussian graphical model, $\Theta_{jj^{\prime}}=0$ corresponds to the conditional independence relationship $x_{j}\perp x_{j}^{\prime}|\{x_{k},k\neq j,j^{\prime}\}$ (Lauritzen, 1996). Therefore a typical assumption, especially in high-dimensional problems, is that $\Theta$ is a sparse matrix; this both allows us to estimate $\Theta$ when $p>n$, and produces a sparse conditional dependence graph.

Model (2) is much more flexible than the i.i.d. graphical model, and it separates co-variation caused by individual preference (cohesion in the mean) from universal co-occurrence (covariance). The price we pay for this flexibility is the much larger number of parameters, and model (2) cannot be fitted without additional assumptions on the mean, since we only have one observation to estimate each vector $\mu_{i}$. The structural assumption we make in this paper is network cohesion, a phenomenon of connected individuals in a social network tending to exhibit similar traits. It has been widely observed in many empirical studies such as health-related behaviors or academic performance (Michell and West, 1996; Haynie, 2001; Pearson and West, 2003). Specifically, in our Gaussian graphical model (2), we assume that connected nodes in the observed network have similar mean vectors. This assumption is reasonable and interpretable in many applications. For instance, in the coauthorship network example, cohesion indicates coauthors tend to have similar word preferences, which is reasonable since they work on similar topics and share at least some publications.

2.3 Fitting the GNC model

The log-likelihood of the data under model (2) is, up to a constant,

[TABLE]

A sparse inverse covariance matrix $\Theta$ and a cohesive mean matrix $M$ are naturally incorporated into the following two-stage procedure, which we call Gaussian graphical model estimation with Network Cohesion and lasso penalty (GNC-lasso).

Algorithm 1** (Two-stage GNC-lasso algorithm)**

Input: a standardized data matrix $X$, network adjacency matrix $A$, tuning parameters $\lambda$ and $\alpha$.

Mean estimation. Let $L_{s}$ be the standardized Laplacian of $A$. Estimated the mean matrix by

[TABLE] 2. 2.

Covariance estimation. Let $\hat{S}=\frac{1}{n}(X-\hat{M})^{T}(X-\hat{M})$ be the sample covariance matrix of $X$ based on $\hat{M}$. Estimate the precision matrix by

[TABLE]

The first step is a penalized least squares problem, where the penalty can be written as

[TABLE]

This can be viewed as a vector version of the Laplacian penalty used in variable selection (Li and Li, 2008, 2010; Zhao and Shojaie, 2016) and regression problems (Li et al., 2019) with network information. It penalizes the difference between mean vectors of connected nodes, encouraging cohesion in the estimated mean matrix. Both terms in (4) are separable in the $p$ coordinates and the least squares problem has a closed form solution,

[TABLE]

In practice, we usually need to compute the estimate for a sequence of $\alpha$ values, so we first calculate the eigen-decomposition of $\mathcal{L}_{s}$ and then obtain each $(I+\alpha\mathcal{L}_{s})^{-1}$ in linear time. In most applications, networks are very sparse, and taking advantage of sparsity and the symmetrically diagonal dominance of $\mathcal{L}_{s}$ allows to compute the eigen-decomposition very efficiently (Cohen et al., 2014). Given $\hat{M}$, criterion (5) is a graphical lasso problem that uses the lasso penalty (Tibshirani, 1996) to encourage sparsity in the estimated precision matrix, and can be solved by the glasso algorithm (Friedman et al., 2008) efficiently or any of its variants, later significantly improved further by Witten et al. (2011) and Hsieh et al. (2014, 2013a).

2.4 An alternative: penalized joint likelihood

An alternative and seemingly more natural approach is to maximize a penalized log-likelihood to estimate both $M$ and $\Theta$ jointly as

[TABLE]

The objective function is bi-convex and the optimization problem can be solved by alternately optimizing over $M$ with fixed $\Theta$ and then optimizing over $\Theta$ with fixed $M$ until convergence. We refer to this method as iterative GNC-lasso. Though this strategy seems more principled in a sense, we implement our method with the two-stage algorithm, for the following reasons.

First, the computational complexity of the iterative method based on joint likelihood is significantly higher, and it does not scale well in either $n$ or $p$. This is because when $\Theta$ is fixed and we need to maximize over $M$, all $p$ coordinates are coupled in the objective function, so the scale of the problem is $np\times np$. Even for moderate $n$ and $p$, solving this problem requires either a large amount of memory or applying Gauss-Seidel type algorithms that further increase the number of iterations. This problem is exacerbated by the need to select two tuning parameters $\lambda$ and $\alpha$ jointly, because, as we will discuss later, they are also coupled.

More importantly, our empirical results show that the iterative estimation method does not improve on the two-stage method (if it does not slightly hurt it). The same phenomenon was observed empirically by Yin and Li (2013) and Lin et al. (2016), who used a completely different approach of applying sparse regression to adjust the Gaussian graphical model, though those papers did not offer an explanation. We conjecture that this phenomenon of maximizing penalized joint likelihood failing to improve on a two-stage method may be general. An intuitive explanation might lie in the fact that the two parameters $M$ and $\Theta$ are only connected through the penalty: the Gaussian log-likelihood (3) without a penalty is maximized over $M$ by $\hat{M}=X$, which does not depend on $\Theta$. Thus the likelihood itself does not pool information from different observations to estimate the mean (nor should it, since we assumed they are different), while the cohesion penalty is separable in the $p$ variables and does not pool information between them either. An indirect justification of this conjecture follows from a property of the two-stage estimator stated in Proposition 2 in Appendix B, and the numerical results in Section 4 provide empirical support.

2.5 Model selection

There are two tuning parameters, $\lambda$ and $\alpha$, in the two-stage GNC-lasso algorithm. The parameter $\alpha$ controls the amount of cohesion over the network in the estimated mean and can be easily tuned based on its predictive performance. In subsequent numerical examples, we always choose $\alpha$ from a sequence of candidate values by 10-fold cross-validation. In each fold, the sum of squared prediction errors on the validation set $\sum(X_{ij}-\hat{\mu}_{ij})^{2}$ is computed and the $\alpha$ value is chosen to minimize the average prediction error. If the problem is too large for cross-validation, we can also use the generalized cross-validation (GCV) statistic as an alternative, which was shown to be effective in theory for ridge-type regularization (Golub et al., 1979; Li, 1986). The GCV statistic for $\alpha$ is defined by

[TABLE]

where we write $\hat{M}(\alpha)$ to emphasize that the estimate depends on $\alpha$. The parameter $\alpha$ should be selected to minimize GCV. Empirically, we observe running the true cross-validation is typically more accurate than using GCV. So the GCV is only recommended for problems that are too large to run cross-validation.

Given $\alpha$, we obtain $\hat{M}$ and use $\hat{S}=\frac{1}{n}(X-\hat{M})^{T}(X-\hat{M})$ as the input of the glasso problem in (5); therefore $\lambda$ can be selected by standard glasso tuning methods, which may depend on the application. For example, we can tune $\lambda$ according to some data-driven goodness-of-fit criterion such as BIC, or via stability selection. Alternatively, if the graphical model is being fitted as an exploratory tool to obtain an interpretable dependence between variables, $\lambda$ can be selected to achieve a pre-defined sparsity level of the graph, or chosen subjectively with the goal of interpretability. Tuning illustrates another important advantage of the two-stage estimation over the iterative method: when estimating the parameters jointly, due to the coupling of $\alpha$ and $\lambda$ the tuning must be done on a grid of their values and using the same tuning criteria. The de-coupling of tuning parameters in the two-stage estimation algorithm is both more flexible, since we can use different tuning criteria for each if desired, and more computationally tractable since we only need to do two line searches instead of a two-dimensional grid search.

2.6 Related work and alternative penalties

The Laplacian smoothness penalty of the form (1) or (6) was originally used in machine learning for embedding and kernel learning (Belkin and Niyogi, 2003; Smola and Kondor, 2003; Zhou et al., 2005). More recently, this idea has been employed in graph-constrained estimation for variable selection in regression (Li and Li, 2008, 2010; Slawski et al., 2010; Pan et al., 2010; Shen et al., 2012; Zhu et al., 2013; Sun et al., 2014; Liu et al., 2019), principal component analysis (Shojaie and Michailidis, 2010), and regression inference (Zhao and Shojaie, 2016). In these problems, a network is assumed to connect a set of random variables or predictors and is used to achieve more effective variable selection or dimension reduction in high-dimensional settings. A generalization to potentially unknown network or group structure was studied by Witten et al. (2014). Though Step 1 of Algorithm 1 has multiple connections to graph constrained estimation, there are a few key differences. In our setting, the network is connecting observations, not variables. We only rely on smoothness across the network for accurate estimation without additional structural assumptions such as sparsity on $M$. In graph-constrained estimation literature, in addition to the Laplacian penalty, other penalties are proposed in special contexts (Slawski et al., 2010; Pan et al., 2010; Shen et al., 2012). We believe similar extensions can also be made in our problem for special applications and we will leave such extensions for future investigation.

An alternative penalty we can impose on $M$ instead of $\sum_{i\sim i^{\prime}}\|{\mu_{i}-\mu_{i^{\prime}}}\|^{2}$ is

[TABLE]

This penalty is called the network lasso penalty (Hallac et al., 2015) and can be viewed as a generalization of the fused lasso (Tibshirani et al., 2005) and the group lasso (Yuan and Lin, 2006). The penalty and its variants were studied recently by Wang et al. (2014); Jung et al. (2018); Tran et al. (2018). This penalty is also associated with convex clustering (Hocking et al., 2011; Lindsten et al., 2011), because it typically produces piecewise constant estimates which can be interpreted as clusters. Properties of convex clustering have been studied by Hallac et al. (2015); Chi and Lange (2015); Tan and Witten (2015). However, in our setting there are two clear reasons for using the Laplacian penalty and not the network lasso. First, piecewise constant individual effects within latent clusters of the network is a special case of the general cohesive individual effects, so our assumption is strictly weaker, and there is no reason to impose piecewise constant clusters in the mean unless there is prior knowledge. Second, solving the optimization in the network lasso problem is computationally challenging and not scalable to the best our knowledge: current state of art algorithms (Hallac et al., 2015; Chi and Lange, 2015) hardly handle more than about 200 nodes on a single core. In contrast, the Laplacian penalty in (6) admits a closed-form solution and can be efficiently solved for thousands of observations even with a naive implementation on a single machine. Moreover, there are many ways to improve the naive algorithm based on the special properties of the linear system (Spielman, 2010; Koutis et al., 2010; Cohen et al., 2014; Sadhanala et al., 2016; Li et al., 2019). Therefore, (6) is a better choice than (9) for this problem, both computationally and conceptually.

3 Theoretical properties

In this section, we investigate the theoretical properties of the two-stage GNC-lasso estimator. Throughout this section, we assume the observation network $A$ is connected which implies that $\mathcal{L}_{s}$ has exactly one zero eigenvalue. The results can be trivially extended to a network consisting of several connected components, either by assuming the same conditions for each component or regularizing $A$ to be connected as in Amini et al. (2013). Recall that $\tau_{1}\geq\tau_{2}\geq\cdots\geq\tau_{n-1}>\tau_{n}=0$ are the eigenvalues of $\mathcal{L}_{s}$ corresponding to eigenvectors $u_{1},\cdots,u_{n}$. For a connected network, we know $\tau_{n}$ is the only zero eigenvalue. Moreover,$\tau_{n-1}$ is known as algebraic connectivity that measure the connectivity of the network.

3.1 Cohesion assumptions on the observation network

The first question we have to address is how to formalize the intuitive notion of cohesion. We will start with the most intuitive definition of network cohesion for a vector, extend it to a matrix, and then give examples satisfying the cohesion assumptions.

Intuitively, we can think of a vector $v\in\mathbb{R}^{n}$ as cohesive on a network $A$ if $v^{T}\mathcal{L}_{s}v$ is small in some sense, or equivalently, $\|{\mathcal{L}_{s}v}\|_{2}$ is small, since $\mathcal{L}_{s}v$ is the gradient of $v^{T}\mathcal{L}_{s}v$ up to a constant and

[TABLE]

It will be convenient to define cohesion in terms of $\mathcal{L}_{s}v$, which also leads to a nice interpretation. The $i$th coordinate of $\mathcal{L}_{s}v$ can be written as

[TABLE]

which is the difference between the value at node $i$ and the average value of its neighbors, weighed by the degree of node $i$. Let $\mathcal{L}_{s}=U\Lambda U^{T}$ be the eigen-decomposition of $\mathcal{L}_{s}$, with $\Lambda$ the diagonal matrix with the eigenvalues $\tau_{1}\geq\dots\geq\tau_{n}$ on the diagonal. The vector $v$ can be expanded in this basis as $v=U\beta=\sum_{i=1}^{n}\beta_{i}v_{i}$ where $\beta\in\mathbb{R}^{n}$. Under cohesion, we would expect $\|{\mathcal{L}_{s}v}\|_{2}^{2}=\sum_{i}\tau_{i}^{2}\beta_{i}^{2}$ to be much smaller than $\|{v}\|_{2}^{2}=\|{\beta}\|_{2}^{2}$. We formalize this in the following definition.

Definition 1** (A network-cohesive vector)**

Given a network $A$ and a vector $v$, let $v=\sum_{i=1}^{n}\beta_{i}u_{i}$ be the expansion of $v$ in the basis of eigenvectors of $\mathcal{L}_{s}$. We say $v$ is cohesive on $A$ with rate $\delta>0$ if for all $i=1,\dots,n$,

[TABLE]

which implies

[TABLE]

Now we can easily define a network-cohesive matrix $M$.

Definition 2** (A network-cohesive matrix)**

A matrix $M\in\mathbb{R}^{n\times p}$ is cohesive on a network $A$ if all of its columns are cohesive on $A$.

An obvious but trivial example of a cohesive vector is a constant vector, which corresponds to $\delta=\infty$. More generally, we define the class of trivially cohesive vectors as follows.

Definition 3** (A trivially cohesive vector)**

We say vector $v$ is trivially cohesive if

[TABLE]

where $\bar{v}=\sum_{i=1}^{n}v_{i}/n$ is the sample mean of $v$, and $\widehat{\mathrm{Var}}(v)=\sum_{i=1}^{n}(v_{i}-\bar{v})^{2}/(n-1)$ is the sample variance of $v$.

Trivial cohesion does not involve a specific network $A$, because such vectors are essentially constant. We say $v$ is nontrivially cohesive if it is cohesive but not trivially cohesive. Similarly, we will say a matrix is trivially cohesive if all its columns are trivially cohesive, and nontrivially cohesive if it is cohesive but not trivially cohesive.

For obtaining theoretical guarantees, we will need to make an additional assumption about the network, which essentially quantifies how much network structure can be used to control model complexity under nontrivial cohesion. This will be quantified through the concept of effective dimension of the network defined below.

Definition 4

Given a connected network adjacency matrix $A$ of size $n\times n$ and eigenvalues of its standardized Laplacian $\tau_{1}\geq\dots\tau_{n-1}>\tau_{n}=0$, define the effective dimension of the network as

[TABLE]

Note that spectral graph theory (Brouwer and Haemers, 2011) implies $\tau_{1}\geq c$ for some constant $c$, and thus for sufficiently large $n$, we always have $m_{A}\leq n-1$. For many sparse and/or structured networks the effective dimension is much smaller than $n-1$, and then we can show nontrivially cohesive vectors/matrices exist.

Our first example of a network with a small effective dimension is a lattice network. Assume $\sqrt{n}$ is an integer and define the lattice network of $n$ nodes by arranging them on a $\sqrt{n}\times\sqrt{n}$ spatial grid and connecting grid neighbors (the four corner nodes have degree 2, nodes along the edges of the lattice have degree 3, and all internal nodes have degree 4).

Proposition 1** (Cohesion on a lattice network)**

Assume $A$ is a lattice network on $n$ nodes, and $\sqrt{n}$ is an integer. Then for a sufficiently large $n$,

The effective dimension $m_{A}\leq n^{2/3}$. 2. 2.

There exist nontrivially cohesive vectors on the lattice network with rate $\delta=1/2$.

Figure 1 shows the eigenvalues and the function $1/\sqrt{m}$ for reference of a $20\times 20$ lattice and of a coauthorship network we analyze in Section 5. For both networks, the effective dimension is much smaller than the number of nodes: for the lattice, $n=400$, while $m_{A}=30$ and for the coauthorship network with $n=635$ nodes, $m_{A}=66$.

3.2 Mean estimation error bound

Our goal here is to obtain a bound on the difference between $M$ and the estimated $\hat{M}$ obtained by Algorithm 1, under the following cohesion assumption.

Assumption 1

The mean matrix $M$ is cohesive over the network $A$ with rate $\delta$ where $\delta$ is a positive constant. Moreover, $\|{M_{\cdot j}}\|_{2}^{2}\leq b^{2}n$ for every $j\in[p]$ for some positive constant $b$.

Theorem 1** (Mean error bound)**

Assume model (2) and Assumption 1 are true. Write $\sigma^{2}=\max_{j}\Sigma_{jj}$ and $\Delta=n^{\frac{1+\delta}{3}}\tau_{n-1}$, where $\tau_{n-1}$ is the smallest nonzero eigenvalue of $\mathcal{L}_{s}$. Then $\hat{M}$ estimated by (4) with $\alpha=n^{\frac{1+\delta}{3}}$ satisfies

[TABLE]

with probability at least $1-\exp(-c(n-m_{A})r(\Sigma))-\exp(-cm_{A}r(\Sigma))-\exp(-c\frac{p\sigma^{2}}{\phi_{max}(\Sigma)})$ for some positive constant $c$, where $m_{A}$ is the effective dimension of network $A$ in Definition 4 and $r(\Sigma)$ is the stable rank of $\Sigma$.

The theorem shows that the average estimation error is vanishing with high probability as long as the cohesive dimension $m_{A}=o(n^{\frac{2(1+\delta)}{3}})$ while $m_{A}r(\Sigma)$ and $p/\phi_{max}(\Sigma)$ grow with $n$. Except for degenerate situations, we would expect $r(\Sigma)$ and $p/\phi_{max}(\Sigma)$ to grow with $p$, which in turn grows with $n$. In (11), the term $\Delta=n^{(1+\delta)/3}\tau_{n-1}$ involves both the cohesion rate of the mean matrix and the algebraic connectivity of the network. In trivially cohesive settings, $\delta\to\infty$ and $\Delta\to\infty$ so the bound does not depend on the network, and the error bound becomes the standard mean estimation error bound. General lower bounds for $\tau_{n-1}$ are available (Fiedler, 1973), but we prefer not to introduce additional algebraic definitions at this point.

Finally, note that the value of $\alpha$ depends on the cohesive rate $\delta$ of $M$. Therefore, the theorem is not adaptive to the unknown cohesive rate. In practice, as we discussed, one has to use cross-validation to tune $\delta$.

3.3 Inverse covariance estimation error bounds

Our next step is to show that $\hat{M}$ is a sufficiently accurate estimate of $M$ to guarantee good properties of the estimated precision matrix $\Theta$ in step 2 of the two-stage GNC-lasso algorithm. We will need some additional assumptions, the same ones needed for the glasso performance guarantees under the standard Gaussian graphical model (Rothman et al., 2008; Ravikumar et al., 2011).

Let $\Gamma=\Sigma\otimes\Sigma$ be the Fisher information matrix of the model, where $\otimes$ denotes the Kronecker product. In particular, under the multivariate Gaussian distribution, we have $\Gamma_{(j,k),(\ell,m)}=\mathrm{Cov}(X_{j}X_{k},X_{\ell}X_{m})$. Define the set of nonzero entries of $\Theta$ as

[TABLE]

We use $S^{c}(\Theta)$ to denote the complement of $S(\Theta)$. Let $s=|S(\Theta)|$ be the number of nonzero elements in $\Theta$. Recall that we assume all diagonals of $\Theta$ are nonzero. For any two sets $T_{1},T_{2}\subset[n]$, let $\Gamma_{T_{1},T_{2}}$ denote the submatrix with rows and columns indexed by $T_{1}$, $T_{2}$, respectively. When the context is clear, we may simply write $S$ for $S(\Theta)$. Define

[TABLE]

where the vector operator $\|{\cdot}\|_{0}$ gives the number of nonzeros in the vector while the matrix norm $\|{\cdot}\|_{\infty,\infty}$ gives the maximum $L_{\infty}$ norm of the rows.

Finally, by analogy to the well-known irrepresentability condition for the lasso, which is necessary and sufficient for the lasso to recover support (Wainwright, 2009), we need an edge-level irrepresentability condition.

Assumption 2

There exists some $0<\rho\leq 1$ such that

[TABLE]

If we only want to obtain a Frobenius norm error bound, the following much weaker assumption is sufficient, without conditions on $\psi,\kappa_{\Sigma},\kappa_{\Gamma}$ and Assumption 2:

Assumption 3

Let $\eta_{\min}(\Sigma)$ and $\eta_{\max}(\Sigma)$ be the minimum and maximum eigenvalues of $\Sigma$, respectively. There exists a constant $\bar{k}$ such that

[TABLE]

Let $\hat{S}=\frac{1}{n}(X-\hat{M})^{T}(X-\hat{M}).$ We use $\hat{S}$ as input for the glasso estimator (5). We would expect that if $\hat{M}$ is an accurate estimate of $M$, then $\Theta$ can be accurately estimated by glasso. The following theorem formalizes this intuition, using concentration properties of $\hat{S}$ around $\Sigma$ and the proof strategy of Ravikumar et al. (2011). We present the high-dimensional regime result here, with $p\geq n^{c_{0}}$ for some positive constant $c_{0}$, and state the more general result which includes the lower-dimensional regime in Theorem 3 in the Appendix, because the general form is more involved.

Theorem 2

Under the conditions of Theorem 1 and Assumption 2, suppose there exists some positive constant $c_{0}$ such that $p\geq n^{c_{0}}$. If $\log{p}=o(n)$ and $m_{A}=o(n)$, there exist some positive constants $C,c,c^{\prime},c^{\prime\prime}$ that only depend on $c_{0},b$ and $\sigma$, such that if $\hat{\Theta}$ is the output of Algorithm 1 with $\alpha=n^{\frac{1+\delta}{3}}$, $\lambda=\frac{8}{\rho}\nu(n,p)$ where

[TABLE]

and $n$ sufficiently large so that

[TABLE]

then with probability at least $1-\exp(-c\log(p(n-m_{A})))-\exp(-c^{\prime}\log(pm_{A}))-\exp(-c^{\prime\prime}\log{p})$, then the estimate $\hat{\Theta}$ has the following properties:

Error bounds:

[TABLE] 2. 2.

Support recovery:

[TABLE]

and if additionally $\min_{(j,j^{\prime})\in S(\Theta)}|\Theta_{jj^{\prime}}|>2(1+8/\rho)\kappa_{\Gamma}\nu(n,p),$ then

[TABLE]

Remark 1

As commonly assumed in literature, such as Ravikumar et al. (2011), we will treat $\kappa_{\Gamma}$, $\kappa_{\Sigma}$ and $\rho$ to be constants or bounded.

Remark 2

The Frobenius norm bound does not need the strong irrepresentability assumption and does not depend on $\kappa_{\Gamma}$ and $\kappa_{\Sigma}$. Following the proof strategy in Rothman et al. (2008), this bound can be obtained under the much weaker Assumption 3 instead.

The quantity in (13) involves four terms. The first term is from the inverse covariance estimation with a known $M$ (a standard glasso problem), and the other three terms come from having to estimate a cohesive $M$. These three terms depend on both the cohesion rate and the effective dimension of the network. As expected, they all increase with $m_{A}$ and decrease with $\delta$. The last term also involves $\Delta$, which depends on both $\delta$ and the algebraic connectivity $\tau_{n-1}$. To illustrate these trade-offs, we consider the implications of Theorem 2 in a few special settings.

First, consider the setting of trivial cohesion, with $\delta=\infty$. In this case, the last three terms in (13) vanish.

Corollary 1

Under the assumptions of Theorem 2, if $M$ is trivially cohesive and $\delta=\infty$, then all results of Theorem 2 hold with

[TABLE]

and the estimated $\hat{\Theta}$ is consistent as long as $\log{p}=o(n).$

This result coincides with the standard glasso error bound from Ravikumar et al. (2011). Thus when $M$ does not vary, we do not lose anything in the rate by using GNC-lasso instead of glasso.

Another illustrative setting is the case of bounded effective dimension $m_{A}$. Then the third term in (13) dominates.

Corollary 2

Under the assumptions of Theorem 2, if the network has a bounded effective dimension $m_{A}$, then all the results of Theorem 2 hold with

[TABLE]

In particular, if $\delta\geq 1/2$, $\hat{\Theta}$ is consistent as long as $\log{p}=o(n)$.

This corollary indicates that if the network structure is favorable to cohesion, the GNC-lasso does not sacrifice anything in the rate up to a certain level of nontrivial cohesion.

Finally, consider a less favorable example in which $\log{p}=o(n)$ may no longer be enough for consistency. Recall Proposition 1 indicates $m_{A}=O(n^{2/3})$ for lattice networks, and suppose the cohesive can be highly nontrivial.

Corollary 3** (Consistency on a $\sqrt{n}\times\sqrt{n}$ lattice)**

Suppose the conditions of Theorem 2 hold and $m_{A}\leq n^{2/3}$. The GNC-lasso estimate $\hat{\Theta}$ is consistent if $\delta>3/8$ and

[TABLE]

In particular, if $\delta=1/2$, it is necessary to have $\log p=o(n^{1/3})$ for consistency.

The corollary suggests that consistency under some regimes of nontrivial cohesion requires strictly stronger conditions than $\log p=o(n)$. Moreover, if cohesion is too weak (say, $\delta\leq 3/8$), consistency cannot be guaranteed by these results.

4 Simulation studies

We evaluate the new GNC-lasso method and compare it to some baseline alternative methods in simulations based on both synthetic and real networks. The synthetic network we use is a $20\times 20$ lattice network with $n=400$ nodes and a vector with dimension $p=500$ observed at each node; this setting satisfies the assumptions made in our theoretical analysis. We also test our method on the coauthorship network shown in Figure 8, which will be described in Section 5. This network has $n=635$ nodes at $p=800$ observed features at each node.

Noise settings:

The conditional dependence graph $\mathcal{G}$ in the Gaussian graphical model is generated as an Erdös-Renyi graph on $p$ nodes, with each node pair connecting independently with probability 0.01. The Gaussian noise is then drawn from $\mathcal{N}(0,\Sigma)$ where $\Theta=\Sigma^{-1}=a(0.3A_{\mathcal{G}}+(0.3e_{\mathcal{G}}+0.1)I)$, where $A_{\mathcal{G}}$ is the adjacency matrix of $\mathcal{G}$,$e_{\mathcal{G}}$ is the absolute value of the smallest eigenvalue of $A_{\mathcal{G}}$ and the scalar $a$ is set to ensure the resulting $\Sigma$ has all diagonal elements equal to 1. This procedure is implemented in Zhao et al. (2012).

Mean settings:

We set up the mean to allow for varying degrees of cohesion. each row $M_{\cdot j},j=1,2,\cdots,p$ as

[TABLE]

where $u^{(j)}$ is randomly sampled with replacement from the eigenvectors of the Laplacian $u_{n-1},n_{n-2},\cdots,u_{n-k}$ for some integer $k$ and $t$ is the mixing proportion. We then rescale $M$ so the signal-to-noise ratio becomes 1.6, so that the problem remains solvable to good accuracy by proper methods but is not too easy to solve by naive methods. In a connected network, the constant vector is trivially cohesive. The cohesion becomes increasingly nontrivial as one increases $k$ and $t$. For example, $t=0$ gives identical mean vectors for all observations and as $t$ increases, the means become more different. The integer $k$ is chosen to give a reasonably eigen-gap in eigenvalues, with details in subsequent paragraphs.

We evaluate performance on recovering the true underlying graph by the receiver operating characteristic (ROC) curve, along a graph estimation path obtained by varying $\lambda$. An ROC curve illustrates the tradeoff between the true positive rate (TPR) and the false positive rate (FPR), defined as

[TABLE]

We also evaluate the methods on the estimation error of $M$, measured as $\|{\hat{M}-M}\|_{\infty}=\max_{ij}|\hat{M}_{ij}-M_{ij}|$ for the worst-case entry-wise recovery and $\|{\hat{M}-M}\|_{2,\infty}=\max_{i}\|{\hat{M}_{i\cdot}-M_{i\cdot}}\|$ for the worst-case mean vector error for each observation.

As a baseline comparison, we include the standard glasso which does not use the network information at all. We also compare to a natural approach to incorporating heterogeneity without using the network; we do this by applying $K$-means clustering to group observations into clusters, estimating a common mean for each cluster, and applying glasso after centering each group with its own mean. This approach requires estimating the number of clusters. However, the widely used gap method (Tibshirani et al., 2001) always suggests only one cluster in our experiments, which defaults back to glasso. Instead, we picked the number of clusters to give the highest area under the ROC curve; we call this result “oracle cluster+glasso” to emphasize that it will not be feasible in practice. For GNC-lasso, we report both the oracle tuning (the highest AUC, not available in practice) and 10-fold cross-validation based tuning, which we recommend in practice. The oracle methods serve as benchmarks for the best possible performance available from each method.

4.1 Performance as a function of cohesion

First, we vary the level of cohesion in the mean, by setting $t$ to 0.1, 0.5, or 1, corresponding to strong, moderate, or weak cohesion. Figure 2 shows the ROC curves of the four methods obtained from 100 independent replications for the lattice network. Glasso fails completely even when the model has only a slight amount of heterogeneity ($t=0.1$). Numerically, we also observed that heterogeneity slows down convergence for glasso. The oracle cluster+glasso improves on glasso as it can accommodate some heterogeneity, but is not comparable to GNC-lasso. As $t$ increases, the GNC-lasso maintains similar levels of performance by adapting to varying heterogeneity, while the oracle cluster+glasso degrades quickly, since for more heterogeneous means the network provides much more reliable information than $K$-means clustering on the observations. We also observed that cross-validation is similar to oracle tuning for GNC-lasso, giving it another advantage. Figure 3 shows the results for the same setting but on the real coauthorship network instead of the lattice. The results are very similar to what we obtained on the lattice, giving further support to GNC-lasso practical relevance.

We also compare estimation errors in $\hat{M}$ in Table 1. The oracle GNC-lasso is almost always the best, except for one setting where it is inferior to the CV-tuned GNC-lasso (note that the “oracle” is defined by the AUC and is thus not guaranteed to produce the lowest error in estimating $M$). For the lattice network, cluster + glasso does comparably to GNC-lasso (sometimes better, and sometimes worse). For the coauthorship network, a more realistic setting, GNC-lasso is always comparable to the oracle and substantially better than both alternatives that do not use the network information.

4.2 Performance as a function of sparsity

A potential challenge for GNC-lasso is a sparse network that does not provide much information, and in particular a network with multiple connected components. As a simple test of what happens when a network has multiple components, we split the $20\times 20$ lattice into either four disconnected $10\times 10$ lattice subnetworks, or 16 disconnected $5\times 5$ lattice subnetworks, by removing all edges between these subnetworks. The data size ($n=400,p=500$) and the data generating mechanism remain the same; we set $t=0.5$ for a moderate degree of cohesion. The only difference here is when there are $K$ connected components in the network, the last $K$ eigenvectors of the Laplacian $u_{n},\cdots,u_{n-K+1}$ are all constant within each connected component (and thus trivially cohesive). Therefore, in the case of 4 disconnected subnetworks, we randomly sample the last $k=12$ eigenvectors to generate $M$ in (14) while in the case of 16 disconnected subnetworks, we set $k=48$. The effective dimensions $m_{A}$ are 30, 32, and 48, respectively.

Similarly, we also split the coauthorship network into two or four subnetworks by applying hierarchical clustering in Li et al. (2018), which is designed to separate high-level network communities (if they exist). We then remove all edges between the communities found by clustering to produce a network with either two or four connected components. To generate $M$ from (14), we use $k=6$ for two components and $k=12$ for four components, and again set $t=0.5$ for moderate cohesion. The effective dimension $m_{A}$ becomes 66, 74, and 78, respectively.

Figure 4 shows the ROC curves and Table 2 shows the mean estimation errors for the three versions of the lattice network. Overall, all methods get worse as the network is split, but the drop in performance is fairly small for the oracle GNC-lasso. Cross-validated GNC lasso suffers slightly more from splitting (the connected components in the last case only have 25 nodes each, which can produce isolated nodes and hurt cross-validation performance). Again, both GNS methods are much more accurate than the two benchmarks (glasso completely fails, and oracle cluster+glasso performance substantially worse).

Figure 5 and Table 3 give the results for the three versions of the coauthorship network. The network remains well connected in all configurations and both the oracle and the cross-validated GNC-lasso perform well in all three cases, without deterioration. The oracle cluster+glasso performs well in this case as well, but GNC-lasso still does better on both graph recovery and estimating the mean. Glasso fails completely once again.

4.3 Performance as a function of the sample size

Here we compare the methods when the sample size $n$ changes while $p$ remains fixed. Specifically, we compare $10\times 10$, $15\times 15$, and $20\times 20$ lattices, corresponding to $n=100$, $225$, and $400$, respectively. The dimension $p=500$, the data generating mechanism, and $t=0.5$ remain the same as in Section 4.1. When $n=100$, the sample size is too small for 10-fold cross-validation to be stable, and thus we use leave-one-out cross-validation instead. Figure 6 shows the ROC curves while Table 4 shows errors in the mean. Clearly, the problem is more difficult for smaller sample sizes, but both versions of GNC-lasso still work better than the other two baseline methods, even though for $n=100$, the problem is essentially too difficult for all the methods. Results on estimating the mean do not favor any one method clearly, but the differences between the methods are not very large in most cases.

4.4 Comparing with the iterative GNC-lasso

Finally, we compare the estimator obtained by iteratively optimizing $\Theta$ and $M$ in (8) (iterative GNC-lasso) to the proposed two-stage estimator (GNC-lasso). As mentioned in Section 2.4, the iterative method is too computationally intensive to tune by cross-validation, so we only compare the oracle versions of both methods, on the synthetic data used in Section 4.1 with moderate cohesion level $t=0.5$. The results are shown in Figure 7 and Table 5. The methods are essentially identically on the lattice network and the two-stage method is in fact slightly better on the co-author network, indicating that there is no empirical reason to invest in the computationally intensive iterative method.

5 Data analysis: learning associations between statistical terms

Here we apply the proposed method to the dataset of papers from 2003-2012 from four statistical journals collected by Ji and Jin (2016). The dataset contains full bibliographical information for each paper and was curated for disambiguation of author names when necessary. Our goal is to learn a conditional dependence graph between terms in paper titles, with the aid of the coauthorship network.

We pre-processed the data by removing authors who have only one paper in the data set, and filtering out common stop words (“and”, “the”, etc) as well as terms that appear in fewer than 10 paper titles. We then calculate each author’s average term frequency across all papers for which he/she is a coauthor. Two authors are connected in the coauthorship network if they have co-authored at least one paper, and we focus on the largest connected component of the network. Finally, we sort the terms according to their term frequency-inverse document frequency score (tf-idf), one of the most commonly used measures in natural language processing to assess how informative a term is (Leskovec et al., 2014), and keep 300 terms with the highest tf-idf scores. After all pre-processing, we have $n=635$ authors and $p=300$ terms. The observations are 300-dimensional vectors recording the average frequency of term usage for a specific author. The coauthorship network is shown in Figure 8.

The interpretation in this setting is very natural; taking coauthorship into account makes sense in estimating the conditional graph, since the terms come from the shared paper title. We can expect that there will be standard phrases that are fairly universal (e.g., “confidence intervals”), as well as phrases specific to relatively small groups of authors with multiple connections, corresponding to specific research area (e.g., “principal components” ), which is exactly the scenario where our model should be especially useful relative to the standard Gaussian graphical model. To ensure comparable scales for both columns and rows, we standardize the data using the successive normalization procedure introduced by Olshen and Rajaratnam (2010). If we select $\alpha$ using 10-fold cross-validation, as before, the graphs from GNC-lasso and glasso recover 4 and 6 edges, respectively, which are very sparse graphs. To keep the graphs comparable and to allow for more interpretable results, we instead set the number of edges to 25 for both methods, and compare resulting graphs, shown in = Figure 10 (glasso) and Figure 10 (GNC-glasso). For visualization purposes, we only plot the 55 terms that have at least one edge in at least one of the graphs.

Overall, most edges recovered by both methods represent common phrases in the statistics literature, including “exponential families”, “confidence intervals”, “measurement error”, “least absolute” (deviation), “probabilistic forecasting”, and “false discovery”. There are many more common phrases that are recovered by GNC-lasso but missed by Glasso, for example, “high dimension(al/s)”, “gene expression”, “covariance matri(x/ces)”, “partially linear”, “maximum likelihood”, “empirical likelihood”, “estimating equations”, “confidence bands”, “accelerated failure” (time model),“principal components” and “proportional hazards”. There are also a few that are found by Glasso but missed by GNC-lasso, for example, “moving average” and “computer experiments”. Some edges also seem like potential false positives, for example, the links between “computer experiments” and “orthogonal construction”, or the edge between “moving average” and “least absolute”, both found by glasso but not GNC-lasso.

Additional insights about the data can be drawn from the $\hat{M}$ matrix estimated by GNC-lasso; glasso does not provide any information about the means. Each $\hat{M}_{\cdot j}$ can be viewed as the vector of authors’ preferences for the term $j$, we can visualize the relative distances between terms as reflected in their popularity. Figure 11 shows the 55 terms from Figure 10, projected down from $\hat{M}$ to $R^{2}$ for visualization purposes by multidimensional scaling (MDS) (Mardia, 1978). The visualization shows a clearly outlying cluster, consisting of the terms “computer”, “experiments”, “construction”, and “orthogonal”, and to a lesser extent the cluster “Markov Chain Monte Carlo” is also further away from all the other terms. The clearly outlying group can be traced back to a single paper, with the title “Optimal and orthogonal Latin hypercube designs for computer experiments” (Butler, 2001), which is the only title where the words “orthogonal” and “experiments” appear together. Note that glasso estimated them as a connected component in the graph, whereas GNC-lasso did not, since it was able to separate a one-off combination occurring in a single paper from a common phrase. This illustrates the advantage of GNC-lasso’s ability to distinguish between individual variation in the mean vector and the overall dependence patterns, which glasso lacks.

6 Discussion

We have extended the standard graphical lasso problem and the corresponding estimation algorithm to the more general setting in which each observation can have its own mean vector. We studied the case of observations connected by a network and leveraged the empirically known phenomenon of network cohesion to share information across observations, so that we can still estimate the means in spite of having $np$ mean parameters instead of just $p$ in the standard setting. The main object of interest is the inverse covariance matrix, which is shared across observations and represents universal dependencies in the population. while all observations share the same covariance matrix under the assumption of network cohesion. The method is computationally efficient with theoretical guarantees on the estimated inverse covariance matrix and the corresponding graph. Both simulations and an application to a citation network show that GNC-lasso is more accurate and gives more insight into the structure of the data than the standard glasso when observations are connected by a network. One possible avenue for future work is obtaining inference results for the estimated model. This might be done by incorporating the inference idea of Zhao and Shojaie (2016) and Ren et al. (2015) with additional structural assumptions on the mean vectors. The absolute deviation penalty (Hallac et al., 2015) between connected nodes is a possible alternative, if the computational cost issue can be resolved through some efficient optimization approach. Another direction is to consider the case where the partial dependence graphs themselves differ for individuals over the network, but in a cohesive fashion; the case of jointly estimating several related graphs has been studied by Guo et al. (2011); Danaher et al. (2014). As always, in making the model more general there will be a trade-off between goodness of fit and parsimony, which may be elucidated by obtaining convergence rates in this setting.

Acknowledgments

T. Li was partially supported by the Quantitative Collaborative grant from the University of Viriginia. E. Levina and T. Li (while a PhD student at the University of Michigan) were supported in part by an ONR grant (N000141612910) and NSF grants (DMS-1521551 and DMS-1916222). J. Zhu and T. Li (while a PhD student at the University of Michigan) were supported in part by NSF grants (DMS-1407698 and DMS-1821243). C. Qian was supported in part by grants from National Key R&D Program of China (2018YFF0215500) and Science Foundation of Jiangsu Province for Young Scholars (SBK2019041494). We want to thank the action editor and reviewers for their valuable suggestions.

A Proofs

First, recall the following matrix norm definitions we’ll need: for any matrix $M$, $\|{M}\|_{\infty}=\max_{ij}|M_{ij}|$, $\|{M}\|_{1,1}=\max_{j}\|{M_{\cdot j}}\|_{1}$, and $\|{M}\|_{\infty,\infty}=\max_{i}\|{M_{i\cdot}}\|_{1}.$

The following lemma summarizes a few concentration inequalities that we will need.

Lemma 1** (Concentration of norm of a multivariate Gaussian)**

For a Gaussian random vector $x\sim\mathcal{N}(0,\Sigma)$, with $\Sigma\in\mathbb{R}^{p\times p}$ a positive definite matrix and $\phi_{\max}(\Sigma)$ the largest eigenvalue of $\Sigma$, we have,

[TABLE]

for some generic constant $c>0$. Further, if $x_{1},\cdots,x_{n}$ are i.i.d. observations from $\mathcal{N}(0,\Sigma)$, then

[TABLE]

where $r(\Sigma)$ is the stable rank of $\Sigma$.

**Proof **[Proof of Lemma 1] The first inequality (15) follows from concentration of a Lipschitz function of a sub-Gaussian random vector. Inequalities (16) and (17) follow from the definition of a sub-exponential random variable. Lastly, (18) follows from applying Bernstein’s inequality to (16) with $t=n\mathrm{tr}(\Sigma)$.

**Proof **[Proof of Proposition 1] By Edwards (2013), the eigenvalues of $A$ are given by

[TABLE]

Since the average degree $2\leq\bar{d}\leq 4$ for a lattice network, we ignore this constant. First, we show $m_{A}\leq n^{2/3}$, which by definition of $m_{A}$ is equivalent to $\tau_{n-n^{2/3}}\geq n^{-1/3}$. Define the set of all eigenvalues satisfying this condition as

[TABLE]

Then it is sufficient to show $|\mathcal{A}_{n}|<n^{2/3}.$ Applying the inequality $\sin(x)\geq\frac{2}{\pi}x$ for $x\in[0,\pi/2]$, we can see that it is sufficient to show $|\tilde{\mathcal{A}}_{n}|<n^{2/3}$, where

[TABLE]

The cardinality of $\tilde{\mathcal{A}}_{n}$ can be computed exactly by counting; for simplicity, we give an approximate calculation for when $n$ is sufficiently large. In this case the proportion of pairs $(i,j)$ out of the entire set of $(\mathbb{N}\cap[0,\sqrt{n}-1])\times(\mathbb{N}\cap[0,\sqrt{n}-1])$ that satisfy the condition to be included in $\tilde{\mathcal{A}}_{n}$ can be upper bounded by twice the ratio betwen the area of the quarter circle with radius $\frac{n^{1/3}}{2}$ and the area of the $\sqrt{n}\times\sqrt{n}$ square. This gives

[TABLE]

To prove the second claim, consider the $\mu=U\beta$ such that all the inequalities in (10) hold as equalities and $\delta=1/2$. Then, by noting that $P_{{\mathbf{1}}}u_{n}=u_{n}$, we have

[TABLE]

We need a lower bound for $\sum_{i<n}\frac{1}{\tau_{i}^{2}}$. By (19),

[TABLE]

Substituting this lower bound for $\sum_{i<n}\frac{1}{\tau_{i}^{2}}$ in (20), for a sufficiently large $n$ we have

[TABLE]

Therefore, the $\mu$ we constructed is nontrivially cohesive.

We can represent each column of $M$ by taking the basis expansion in $U$, obtaining the basis coefficient matrix $B=(B_{\cdot 1},B_{\cdot 2},\cdots,B_{\cdot p})$ such that $M=UB$. Let $\hat{B}=U^{T}\hat{M}$, where $\hat{M}$ is the estimate (4). We can view $\hat{B}$ as an estimate of $B$. We first state the error bound for $\hat{B}$ in Lemma 2, and the bound for $\hat{M}$ directly follows.

Lemma 2

Under model (2) and Assumption 1, if $\alpha=n^{\frac{1+\delta}{3}}$, we have

In maximum norm,

[TABLE]

with probability at least $1-\exp(-c\log{(p(n-m_{A}))})-\exp(-c^{\prime}\log(pm_{A}))$ for some constants $C$, $C^{\prime}$, $c$, $c^{\prime}$, and $c^{\prime\prime}$. 2. 2.

In Frobenius norm,

[TABLE]

with probability at least $1-\exp(-c^{\prime\prime}(n-m_{A})r(\Sigma))-\exp(-c^{\prime\prime}m_{A}r(\Sigma))-\exp(-c^{\prime\prime}r(\Sigma))$. 3. 3.

if $\log{p}=o(n)$ and $m_{A}=o(n)$, then

[TABLE]

with probability at least $1-\exp(-cn)-\exp(-Cm_{A}\log{p})-\exp(-C\log{p})$.

**Proof **[Proof of Lemma 2] Solving (4), we can explicitly write out

[TABLE]

In particular, for each column $j\in[p]$, the estimate can be written as

[TABLE]

where $\tilde{E}_{\cdot j}\sim\mathcal{N}(0,\sigma^{2}I)$. Let $\mathcal{Q}^{j}$ and $\mathcal{R}^{j}$ be two $n$ dimensional vectors such that the $i$th element of $\mathcal{Q}^{j}$ is given by $\frac{\alpha\tau_{i}}{1+\alpha\tau_{i}}B_{ij}$ while the $i$th element of $\mathcal{R}^{j}$ is given by $\frac{1}{1+\alpha\tau_{i}}\tilde{E}_{ij}$.

[TABLE]

For the element-wise $L_{\infty}$ norm, we have

[TABLE]

where the second inequality is by Definition 1. The term $\mathcal{R}^{j}$ can be decomposed into two parts, the first $n-m_{A}$ elements and the last $m_{A}$ elements. For the first $n-m_{A}$ elements, we have

[TABLE]

by Definition 4, with probability at least $1-\exp(-c\log(p(n-m_{A})))$. For the remaining $m_{A}$ elements, with probability at least $1-\exp(-c^{\prime}\log(pm_{A}))$, we have

[TABLE]

Combining (24)–(A) leads to (21), since

[TABLE]

with probability at least $1-\exp(-c\log(p(n-m_{A})))-\exp(-c^{\prime}\log(pm_{A}))$ for sufficiently large $n$.

For the column-wise $L_{\infty}$ norm, we have

[TABLE]

For the second term,

[TABLE]

By Lemma 1, for each $j\in[p]$, $\mathbb{P}(\sum_{i\leq n-m_{A}}|\tilde{E}_{ij}|>2\sigma(n-m_{A}))\leq\exp(-2c(n-m_{A}))$ for some constant $c$; therefore

[TABLE]

as long as $\log{p}=o(n)$ and $m_{A}=o(n)$.

Assume $m_{A}\geq 2$, again by Lemma 1, for each $j\in[p]$,

[TABLE]

for some constant $C,c^{\prime}>0$ with $C>1$. Therefore,

[TABLE]

and

[TABLE]

The above result is also trivially true if $m_{A}=1$. Substituting these two inequalities into (A) gives

[TABLE]

with probability at least $1-\exp(-cn)-\exp(-Cm_{A}\log{p})-\exp(-C\log{p})$. Now combining (A) and (A), we get

[TABLE]

with probability at least $1-\exp(-cn)-\exp(-Cm_{A}\log{p})-\exp(-C\log{p})$ as long as $p\geq 3$.

Finally, for the Frobenius norm we have

[TABLE]

For the second term,

[TABLE]

If $m_{A}\geq 2$, by (18) from Lemma 1, for a proper $c$, we have

[TABLE]

Putting everything together,

[TABLE]

with probability at least $1-2\exp(-c(n-m_{A})r(\Sigma))-2\exp(-cm_{A}r(\Sigma))-2\exp(-c\frac{p\sigma^{2}}{\phi_{max}(\Sigma)})$. Notice that when $m_{A}=1$, the above result is trivially true.

**Proof **[Proof of Theorem 1] By definition, we have $\|{\hat{M}-M}\|_{F}=\|{U(\hat{B}-B)}\|_{F}=\|{\hat{B}-B}\|_{F}$. Thus the theorem follows directly from Lemma 2 and the fact that $n-m_{A}\leq n$. Note that the Frobenius norm bound in Lemma 2 does not need $\log{p}=o(n)$ and $m_{A}=o(n)$.

Now we proceed to prove Theorem 2. Let

[TABLE]

$S$ is the sample covariance matrix used by the glasso algorithm when the mean is assumed known (and without loss of generality set to 0). The success of glasso is dependent on $S$ concentrating around the true covariance matrix $\Sigma$. If we can show $\hat{S}$ concentrates around $\Sigma$, we should be able to prove similar properties of GNC-lasso.

Lemma 3

Under the conditions of Theorem 1 and assuming $\log{p}=o(n),m_{A}=o(n)$, we have

[TABLE]

with probability at least $1-\exp(-c\log(p(n-m_{A})))-\exp(-c\log(pm_{A}))-\exp(-c\log{p})$ for some constant $C$ and $c$ that only depend on $b$ and $\sigma$.

**Proof **[Proof of Lemma 3] We will be using $C$ and $c$ to denote generic constants whose value might change across different lines. Using the triangular inequality,

[TABLE]

we will prove concentration in two steps. Starting with the first term and writing $X=M+E$, we have

[TABLE]

By Lemma 2, for some constant $C$ depending on $b$ and $\sigma$,

[TABLE]

with probability at least $1-\exp(-c\log(p(n-m_{A})))-\exp(-c\log(pm_{A}))-\exp(-c\log(p))$.

On the other hand, note that $U^{T}E=(U_{\i\cdot}E_{\cdot j})_{i,j=1}^{n}$ and $\|{U_{i\cdot}}\|_{2}=1$, so $(U^{T}E)_{ij}\sim\mathcal{N}(0,\sigma^{2})$. Therefore

[TABLE]

with probability at least $1-\exp(-c\log(np))$. Hence the second and third terms in (A) satisfy

[TABLE]

with probability at least $1-\exp(-c\log(p(n-m_{A}))-\exp(-c\log(pm_{A}))-\exp(-c\log(np))$. Note that both terms in (A) dominate the last two terms in (A). Thus substituting (A) and (A) into (A) leads to

[TABLE]

with probability at least $1-\exp(-c\log(p(n-m_{A})))-\exp(-c\log(pm_{A}))-\exp(-c\log(p))$.

In addition, Lemma 1 of Ravikumar et al. (2011) implies that

[TABLE]

with probability at least $1-\exp(-c\log{p})$. Therefore, we have

[TABLE]

with probability at least $1-\exp(-c\log(p(n-m_{A})))-\exp(-c\log(pm_{A}))-2\exp(-c\log{p})$.

For conciseness, we present Theorem 2 in the main text by assuming the more interesting situations $p\geq n^{c_{0}}$. This is not necessary in any sense, so here we prove a trivially more general version of Theorem 2 without the high-dimensional assumption. For completeness, we rewrite the theorem here.

Theorem 3** (Trivially generalized version of Theorem 2)**

Under the conditions of Theorem 1 and Assumption 2, if $\log{p}=o(n)$ and $m_{A}=o(n)$, there exist some positive constants $C,c,c^{\prime},c^{\prime\prime}$ that only depend on $b$ and $\sigma$, such that if $\hat{\Theta}$ is the output of Algorithm 1 with $\alpha=n^{\frac{1+\delta}{3}}$, $\lambda=\frac{8}{\rho}\nu(n,p)$ where

[TABLE]

and $n$ sufficiently large so that

[TABLE]

then with probability at least $1-\exp(-c\log(p(n-m_{A})))-\exp(-c^{\prime}\log(pm_{A}))-\exp(-c^{\prime\prime}\log{p})$, then the estimate $\hat{\Theta}$ has the following properties:

Error bounds:

[TABLE] 2. 2.

Support recovery:

[TABLE]

and if additionally $\min_{(j,j^{\prime})\in S(\Theta)}|\Theta_{jj^{\prime}}|>2(1+8/\rho)\kappa_{\Gamma}\nu(n,p),$ then

[TABLE]

We will use the primal-dual witness strategy from Ravikumar et al. (2011) for proof. We show that even if $\hat{S}$ has a worse concentration around $\Sigma$ than $S$, we can still achieve consistency and sparsistency under certain regularity conditions.

**Proof **[Proof of Theorem 3 and Theorem 2] The argument follows the proof of Theorem 1 in Ravikumar et al. (2011). In particular, for the event where the bound in Lemma 3 holds, we just have to show that the primal-dual witness construction succeeds. The choice of $\lambda=\frac{8}{\rho}\nu(n,p)$ ensures $\|{\hat{S}-\Sigma}\|_{\infty}\leq\frac{\rho}{8}\lambda$. With the requirement on the sample size, the assumptions of Lemma 5 and 6 in Ravikumar et al. (2011) hold, implying strict dual feasibility holds for the primal-dual witness, which shows the procedure succeeds. Then the first claim of the theorem is a direct result of Lemma 6 in Ravikumar et al. (2011) and the second claim is true by construction of the primal-dual witness procedure. The remaining bounds can be proved similarly.

Finally, if $p\geq n^{c_{0}}$, we have $\log(p)\geq c_{0}\log(n)$, so in (3), the 3rd coincides with the 5th, and the $4$th term coincides with the 6th term, by magnitude, resulting in the form

[TABLE]

To further simplified the form, we apply another upper bound for the term by the fact that $\sqrt{m_{A}}\geq 1$, $\frac{m_{A}}{\Delta+1}+1\geq 1$ and $\delta\geq 0$, which gives

[TABLE]

B Oracle mean estimation by GNC-lasso

In our setting, unlike in the classical glasso setting, the mean estimate is also of interest, and in this section we show that our estimate $\hat{M}$ enjoys a weak oracle property in a certain sense. We use the spectrum of $\mathcal{L}_{s}$ as a basis again, writing $U$ for the matrix of eigenvectors of $\mathcal{L}_{s}$ and expanding a matrix $M\in\mathbb{R}^{n\times p}$ as $M=UB$. Since $U$ is given and orthonormal, estimating $M$ is equivalent to estimating $B$. In an ideal scenario, if the true value $\Theta$ is given to us by an oracle, we could estimate $B$ by minimizing one of the two objective functions:

[TABLE]

where $\Lambda=\mathrm{diag}(\tau_{1},\tau_{2},\cdots,\tau_{n})$ is the diagonal matrix of eigenvalues of $\mathcal{L}_{s}$. It is easy to verify that (36) is equivalent to the mean estimation step (4) in the two-stage procedure (up to $U$), while (37) is equivalent to estimating the mean by maximizing the joint penalized likelihood (8) with $\Theta$ fixed at the true value. We can then treat (37) as an oracle estimate in the sense that it uses the true value of the covariance matrix. It serves as a benchmark for the best performance one could expect in estimating $B$ (or equivalently $M$). Let $\hat{B}_{1}$ and $\hat{B}_{2}$ be the estimates from (36) and (37), respectively, and let $W_{k}=B-\hat{B}_{k}$, $k=1,2$ be the corresponding estimation error matrices. We then have the following result.

Proposition 2

Under model (2), assume $W_{1}$ and $W_{2}$ are the errors defined above with the same tuning parameter $\alpha$. Under the Assumption 3, if $\Theta$ is diagonally dominant with $\max_{j}\frac{\sum_{j^{\prime}\neq j}|\Theta_{j^{\prime}j}|}{\Theta_{jj}}\leq\rho<1$, then there exist a matrix $\tilde{W}$ such that

[TABLE]

for the constant $\bar{k}$ in Assumption 3 and

[TABLE]

where each row $E$ is i.i.d from multivariate Gaussian $\mathcal{N}(0,\Sigma).$

Proposition 2 shows $\tilde{W}$ and $W_{1}$ are stochastically equivalent while $\tilde{W}$ and $W_{2}$ are roughly the same in $\|{\cdot}\|_{\infty}$. Therefore, (36) and (37) are essentially equivalent in the sense of entrywise error bound, implying that $\hat{M}$ computed by GNC-lasso cannot be non-trivially improved by the oracle estimator under the true model with known $\Theta$.

Proposition 2 makes an additional assumption on diagonal dominance of $\Theta$, which is a relatively mild assumption consistent with others in this context. To see this, consider a general multivariate Gaussian vector $y\sim\mathcal{N}(0,\Sigma)$. Then we can write

[TABLE]

where the vector $\zeta^{j}\in\mathbb{R}^{p}$ satisfies $\zeta^{j}_{j^{\prime}}=-\frac{\Theta_{jj^{\prime}}}{\Theta_{j^{\prime}j^{\prime}}}$ for $j^{\prime}\neq j$ and $\zeta^{j}_{j}=0$, and $\xi_{j}$ is a Gaussian random variable with zero mean and variance equal to the conditional variance of $y_{j}$ given $\{y_{j^{\prime}}\}_{j^{\prime}\neq j}$. Thus the diagonal dominance assumption of Proposition 2 is essentially assuming

[TABLE]

This has the same form as Assumption 4 of Meinshausen and Bühlmann (2006), who proposed node-wise regression to estimate the Gaussian graphical model. There $\rho<1$ is needed for node-wise regression to consistently estimate the graph structure (see Proposition 4 of Meinshausen and Bühlmann (2006)).

Remark 3** (Implications for iterative estimation)**

If the iterative algorithm is used to obtain $\tilde{M}$ and $\tilde{\Theta}$, we know $\tilde{M}$ is the solution of (37) with $\Theta$ replaced by $\tilde{\Theta}$. Since $\tilde{\Theta}$ is only an estimate of $\Theta$, we would not expect this estimator to work as well as the oracle estimator (37). Since $\hat{M}$ cannot be improved by the oracle estimator, intuitively we make the conjecture that $\tilde{M}$ cannot significantly improve on $\hat{M}$ either.

To prove Proposition 2, we need a few properties of Kronecker products. Recall that given two matrices $A\in\mathbb{R}^{m\times n}$ and $B\in\mathbb{R}^{p\times q}$, their Kronecker product is defined to be an $(mp)\times(nq)$ matrix such that

[TABLE]

For a matrix $A$, define $\mathrm{vec}(A)$ to be the column vector stacking all columns of $A$, $\mathrm{vec}(A)=(A_{\cdot 1},A_{\cdot 2},\cdots,A_{\cdot n})$. Some standard properties we’ll need, assuming the matrix dimensions match appropriately, are stated next.

[TABLE]

Proposition 3

For the estimates $W_{1}$ from (36) and $W_{2}$ from (37), we have

[TABLE]

where $\tilde{E}=(\tilde{\epsilon}_{1\cdot},\tilde{\epsilon}_{2\cdot},\cdots,\tilde{\epsilon}_{n\cdot})$ and $\tilde{\epsilon}_{i\cdot}\sim\mathcal{N}(0,\Sigma)$ are i.i.d., and $\dot{E}=(\dot{\epsilon}_{1\cdot},\dot{\epsilon}_{2\cdot},\cdots,\dot{\epsilon}_{n\cdot})$, and $\dot{\epsilon}_{i\cdot}\sim\mathcal{N}(0,\Theta)$ are i.i.d. In particular, $\tilde{E}=-U^{T}E$ and $\dot{E}=-U^{T}E\Theta$.

**Proof **[Proof of Proposition 3] We only prove (39); the proof of (38) is exactly the same, with $\Theta$ replaced by $I_{p}$. The conclusion follows directly from writing out the quadratic optimiation solution after vectorizing all matrices. Specifically, the objective function (37) can be written as

[TABLE]

The minimizer of this quadratic function satisfies

[TABLE]

Substituting $X=UB+E$ into the estimating equation gives

[TABLE]

and therefore

[TABLE]

We then get

[TABLE]

This is equivalent to (39) by noting that $\dot{E}=-U^{T}E\Theta$.

Now we show $W_{1}$ and $W_{2}$ are essentially equivalent estimation errors. Define two additional estimating equations as below:

[TABLE]

[TABLE]

The error equation (40) corresponds to the situation when we carry $p$ separate Laplacian smoothing estimations. The error equation (40) is also from $p$ separate Laplacian smoothing but it adjusts the weight each variable to be proportional to $1/\Theta_{jj}$, which can be seen as $W_{2}$ approximation after ignoring off-diagonal elements of $\Theta$. Intuitively, when the off-diagonal elements are small, $W_{2}$ should not be very different from $W_{4}$, and when the diagonal elements of $\Theta$ are similar, as in Assumption 3, $W_{3}$ and $W_{4}$ should also be similar. The following proposition formalizes this intuition under the assumption that $\Theta$ is diagonally dominant. We can then conclude that using the true $\Theta$ in (37) does not really bring improvement and $W_{1}$, $W_{2}$, $W_{3}$, and $W_{4}$ are all essentially equivalent.

Proposition 4

Assume $W_{2}$, $W_{3}$, and $W_{4}$ are the estimation errors from (39), (40) and (41), respectively, with the same $\alpha$. If $\Theta$ is diagonally dominant with $\max_{j}\frac{\sum_{j^{\prime}\neq j}|\Theta_{j^{\prime}j}|}{\Theta_{jj}}\leq\rho<1$, then

[TABLE]

In particular, under Assumption 3,

[TABLE]

for a constant $\bar{k}$.

**Proof **[Proof of Proposition 4] Directly from the definition, we have

[TABLE]

This implies that for any $i$, $j$ and an arbitrary $\alpha$,

[TABLE]

We next show that under the assumption of diagonal dominance of $\Theta$, even $W_{2}$ cannot do much better. For each $j=1,2,\cdots,p$, from (39),

[TABLE]

Therefore, we have

[TABLE]

in which the last equation comes from (41). By triangle inequality, (44) leads to

[TABLE]

Taking the maximum over $j$ on both sides, we have

[TABLE]

Similarly using triangle inequality in the other direction, we get

[TABLE]

Combining this with (43), we get

[TABLE]

Note that (45) holds if we replace $\|{\cdot}\|_{\infty}$ by other norms. For example, if we take the $L_{1}$ norm instead, we get a similar bound in $\|{\cdot}\|_{1,1}$.

Now we are ready to prove Proposition 2.

**Proof **[Proof of Proposition 2] By taking $\tilde{W}=W_{3}$ and using the conclusion of Proposition 4, the first half of Proposition 2 directly follows. Subtracting (38) from (40) leads to

[TABLE]

and therefore

[TABLE]