# Joint Gaussian Graphical Model Estimation: A Survey

Katherine Tsai  Department of Electrical and Computer Engineering, University
of Illinois at Urbana-Champaign  Oluwasanmi Koyejo  Department of Computer
Science, University of Illinois at Urbana-Champaign  Mladen Kolar  Booth
School of Business, The University of Chicago

## Abstract

Graphs from complex systems often share a partial underlying structure across
domains while retaining individual features. Thus, identifying common
structures can shed light on the underlying signal, for instance, when applied
to scientific discoveries or clinical diagnoses. Furthermore, growing evidence
shows that the shared structure across domains boosts the estimation power of
graphs, particularly for high-dimensional data. However, building a joint
estimator to extract the common structure may be more complicated than it
seems, most often due to data heterogeneity across sources. This manuscript
surveys recent work on statistical inference of joint Gaussian graphical
models, identifying model structures that fit various data generation
processes. Simulations under different data generation processes are
implemented with detailed discussions on the choice of models.

##  1 Introduction

Graphical models are powerful tools for expressing statistical relationships
between variables. Examples of practical use are ubiquitous, and include
models that characterize the causal relationships between the neurological
activity of brain regions, genetic expression across genes, and a variety of
other physiological measurements. A variety of applications have illustrated
the value of graphical models for analyzing scientific phenomena (Felsenstein,
1981; Schäfer and Strimmer, 2005; Friedman et al., 2000; Chan et al., 2017;
Dondelinger et al., 2013). Specifically, graphical models have proven useful
for elucidating the mechanisms of brain function (Foti and Fox, 2019; Manning
et al., 2018; Schwab et al., 2018; Greenewald et al., 2017; Colclough et al.,
2018; Qiu et al., 2016; Skripnikov and Michailidis, 2019). This manuscript
outlines joint graphical models, an extension to standard graphical models
that are useful for jointly analyzing data from multiple sources, e.g.,
neurological data measured at multiple timescales, or joint neurological,
genetic and phenotypic data. Specifically, this manuscript lays out the
representation of joint graphical models and some of their properties, then
outlines the best practices for estimating joint graphical models. This
manuscript provides examples of data generation processes where the joint
approach can significantly improve estimates compared to separate estimation.

A graph $G=(V,E)$ consists of a set of $p$ nodes, also known as vertices
$V=\\{1,\ldots,p\\}$ and a set of edges $E\subseteq V\times V$. In a
probabilistic graphical model (Lauritzen, 1996), the set of nodes $V$ is
associated with coordinates of a random vector ${\bf
x}=(x_{1},\ldots,x_{p})^{\top}$ and the edge set $E$ captures dependency
relationships between the components of the vector. In particular, in an
undirected probabilistic graphical model, the absence of an edge between nodes
$a$ and $b$ indicates that $x_{a}$ and $x_{b}$ are conditionally independent
given all other variables $x_{-\\{a,b\\}}=\\{x_{c}\mid c\in
V\backslash\\{a,b\\}\\}$. In a neuroscience application, the random vector
${\bf x}$ could represent, for example, measurements of brain activity in
different regions – so the set of edges corresponds to functional brain
connectivity. Given $n$ measurements of the vector ${\bf x}$, inferring the
graph structure corresponds to identifying pairs of coordinates that are
conditionally independent given all other variables (Drton and Maathuis,
2017). Inferring the graph structure based on conditional associations is more
challenging than inferring the correlation structure between the measurements.
However, the conditional independence graphs are generally considered more
scientifically meaningful (Dobra et al., 2004).

##  2 Background: Gaussian Graphical Models

The most widely used examples of probabilistic graphical models are Gaussian
graphical models, where ${\bf
x}\sim\mathcal{N}({\bm{\mu}},{\bm{\Omega}}^{-1})$ is assumed to be distributed
as a multivariate Gaussian vector with mean vector ${\bm{\mu}}$, and precision
matrix ${\bm{\Omega}}$, whose entries correspond to the partial correlation
between the associated variables. In this setting, any two coordinates $x_{a}$
and $x_{b}$ are conditionally independent given $x_{-\\{a,b\\}}$ if and only
if the $(a,b)$ entry of the precision matrix ${\bm{\Omega}}$ is zero
(Lauritzen, 1996), and the graph structure can be inferred based on nonzero
entries of ${{\bm{\Sigma}}}^{-1}:={\bm{\Omega}}$, also known as the inverse
covariance matrix. Throughout the manuscript, we use the terms inverse
covariance matrix and precision matrix interchangeably. In practice, the
covariance matrix is not known and the graph structure needs to be estimated
using samples drawn from an underlying distribution. For example, in a low
dimensional setting we can first obtain an estimator of the precision matrix
by maximizing the log-likelihood

$\displaystyle\widehat{\bm{\Omega}}=\mathop{\mathrm{argmax}}\quad
n\left[\frac{1}{2}\log\\{\det({\bm{\Omega}})\\}-\frac{1}{2}\mathop{\mathrm{tr}}(\widehat{\bm{\Sigma}}{\bm{\Omega}})\right],$
(1)

where $\det(\cdot)$ is the determinant, $\mathop{\mathrm{tr}}(\cdot)$ is the
trace, and $\widehat{{\bm{\Sigma}}}=n^{-1}\sum_{i=1}^{n}({\bf x}_{i}-\bar{\bf
x})({\bf x}_{i}-\bar{\bf x})^{\top}$ is the empirical covariance matrix. Next,
the graph structure is estimated by thresholding small (in absolute value)
elements of $\widehat{{\bm{\Omega}}}$ or testing whether they are zero (Drton
and Maathuis, 2017), that is, the graph structure corresponds to nonzero
entries of the resulting thresholded precision matrix.

Figure 1: Comparison of joint and separate graph estimation procedures with
$p=20$ variables and sample size of $n=5000$. Pink nodes do not have edges
connecting to them, while blue nodes have at least on edge connecting to them.
(Top row): Ground truth graphs. (Second row): Jointly estimated graphs using
fused lasso. (Third row): Jointly estimated graphs using the group graphical
lasso. (Fourth row): Graphs estimated separately using the graphical lasso.
Details of all methods are provided in the text. From the figure, it is clear
that joint estimation significantly outperforms separate estimation.

In a high-dimensional setting, where the number of parameters to estimate,
$p$, is much larger than the number of data points observed, $n$, maximizing
the log-likelihood (1) results in poor quality estimates. In the particular
case of Gaussian data, the resulting estimate, that is, the inverse of the
covariance matrix, does not exist when $n<p$. Unfortunately, the high-
dimensional setting is prevalent in various applications. For example,
functional imaging of brain measurements using (standard) $2mm^{3}$ voxels
will result in approximately $p=O(10^{6})$ voxels with $n=O(10^{2})$
measurements (Poldrack et al., 2011; Hsieh et al., 2013). While $p$ is large
as compared to $n$, most entries in ${\bm{\Omega}}$, denoted as
$\omega_{i,j}$, $i,j=1,\ldots,p$, are zero, that is, the inverse covariance
matrix is sparse. Hence, a typical strategy to estimate ${\bm{\Omega}}$ in a
high-dimensional setting is to add a regularization function, such as the
$\ell_{1}$-norm of the parameters, to the log-likelihood function (1), which
encourages the graph to be sparse or have other desirable structural biases
(Yuan and Lin, 2007; Bühlmann and Van De Geer, 2011). Specifically, we
estimate ${\bm{\Omega}}$ using the following optimization program

$\displaystyle\widehat{\bm{\Omega}}=\mathop{\mathrm{argmax}}\quad
n\left[\frac{1}{2}\log\\{\det({\bm{\Omega}})\\}-\frac{1}{2}\mathop{\mathrm{tr}}(\widehat{\bm{\Sigma}}{\bm{\Omega}})\right]+\lambda_{n}\sum_{i\neq
j}|\omega_{i,j}|.$ (2)

Note that regularization is not added to the diagonal terms, $\omega_{i,i}$,
$i=1,\ldots,p$, because ${\bf\Omega}$ is positive definite and adding penalty
on the diagonal entries would introduce additional bias. In this manuscript,
we focus on simultaneous estimation of multiple graphs that are structurally
similar. We will illustrate in the following sections that exploiting the
common structures will improve the estimation results on every individual
graph.

###  2.1 Joint Gaussian Graphical Models

We continue to use the example of brain measurements to demonstrate the idea
of joint graphical model estimation. Consider the case of $n$ fMRI scans
collected from each of the $K$ subjects. Suppose we seek to estimate the
functional connectivity (in this case, a graphical model) between the $p$
voxels of each subject $k$, where $p$ is much greater than $n$ ($p\gg n$).
Each functional connectivity network shares similarities with other networks,
but are not identical. To better estimate the network of the subject $k$ from
$n$ MRI scans, we may borrow information from other networks given that they
are expected to share similar patterns. One simple approach is to construct a
regularization function that encourages similarities between graphs, an
approach known as joint estimation. Figure 1 shows a promising result when one
pools the data across subjects, the resulting estimates better recover the
ground truth graphs compared to separate estimations. It has also been shown
that joint estimation can increase sensitivity and detect edges that are
missing in separate estimation (Peterson et al., 2015). Thus, ignoring the
information of other groups may lead to suboptimal solutions (Danaher et al.,
2014; Lee and Liu, 2015). Moreover, joint estimation of graphical models has
been applied successfully in a number of problems, including metabolite
experiments (Tan et al., 2017), cancer networks (Mohan et al., 2012; Peterson
et al., 2015; Lee and Liu, 2015; Saegusa and Shojaie, 2016; Hao et al., 2018),
biomedical data (Yajima et al., 2014; Kling et al., 2015; Pierson et al.,
2015), gene expression (Chun et al., 2015; Lin et al., 2017), text processing
(Guo et al., 2011), climate data (Ma and Michailidis, 2016), and fMRI (Qiu et
al., 2016; Colclough et al., 2018; Skripnikov and Michailidis, 2019; Lukemire
et al., 2020). In all of these problems, data are heterogeneous, but the
graphs share similarities.

To rigorously describe the example discussed above, we consider the problem of
estimating graph structures $G^{(k)}=(V,E^{(k)})$, $k=1,\ldots,K$, from $K$
related groups of data. The data for each group are $p$-variate and share the
same set of nodes $V$, but the underlying connection patterns $E^{(k)}$ may be
different due to the heterogeneity between groups. The data for the $k$-th
group can be represented as an $n_{k}\times p$ matrix ${\bf X}^{(k)}=({\bf
x}_{1}^{(k)},\ldots,{\bf x}_{n_{k}}^{(k)})^{\top}$, where the rows ${\bf
x}_{i}^{(k)}=(x_{i,1}^{(k)},\ldots,x_{i,p}^{(k)})$, $i=1,\ldots,n_{k}$, are
$p$-dimensional vectors of observations. Assuming that the data in each group
are distributed according to a $p$-variate Gaussian distribution, ${\bf
x}_{i}^{(k)}\sim\mathcal{N}({\bm{\mu}}^{(k)},({\bm{\Omega}}^{(k)})^{-1})$,
$i=1,\ldots,n_{k}$, where ${\bm{\mu}}^{(k)}\in\mathbb{R}^{p}$ is the mean,
which we assume without loss of generality to be $\bm{0}$, and
${\bm{\Omega}}^{(k)}\in\mathbb{R}^{p\times p}$ is the precision matrix. Given
observations ${\bf X}=\\{{\bf X}^{(1)},\ldots,{\bf X}^{(K)}\\}$, we can
estimate ${\bm{\Omega}}=\\{{\bm{\Omega}}^{(1)},\ldots,{\bm{\Omega}}^{(K)}\\}$
by maximizing the penalized joint log-likelihood for $K$ groups:

$\displaystyle\widehat{{\bm{\Omega}}}=\mathop{\mathrm{argmax}}\;\tilde{\ell}({\bm{\Omega}})+{\bf
P}({\bm{\Omega}}),\qquad\tilde{\ell}({\bm{\Omega}}):=\sum_{k=1}^{K}n_{k}\left[\log\\{\det({\bm{\Omega}}^{(k)})\\}-\mathop{\mathrm{tr}}(\widehat{{\bm{\Sigma}}}^{(k)}{\bm{\Omega}}^{(k)})\right],$
(3)

where
$\widehat{{\bm{\Sigma}}}^{(k)}={n_{k}}^{-1}({\bm{X}}^{(k)})^{\top}{\bm{X}}^{(k)}$,
$k=1,\ldots,K$, are the sample covariance matrices. Directly solving (3)
without the penalty ${\bf P}({\bm{\Omega}})$ gives the maximum log-likelihood
estimate of ${\bm{\Omega}}$. However, the solution is equivalent to solving
the maximum log-likelihood estimate of each group individually and fails to
utilize the shared “information” across different groups. We hence explore
different approaches that use the penalty function ${\bf P}({\bm{\Omega}})$ to
incorporate the group structure and focus on the structural assumptions behind
the penalties. The comparison of different methods introduced in the text can
be found in Table 1.

The rest of the manuscript is organized as follows. In Section 3, we introduce
methods that employ coarse-grained structural constraints, i.e., that each
group is equally related to each other. Methods that employ fine-grained
structural constraints, i.e., that some groups are more or less correlated
with other groups, are discussed in Section 4. Section 5 and 6 illustrate two
practical examples. Section 5 covers differential graphs, which are special
cases of the joint estimation paradigm with two groups. Joint estimation of
time-series data is discussed in Section 6. Finally, in Section 7, we present
an empirical evaluation of a variety of regularization methods to evaluate
their performance under two graph generation processes: shared inverse
covariance structure and shared inverse covariance values. The simulation code
is available at https://github.com/koyejo-lab/JointGraphicalLasso.

Table 1: Variants of the joint Gaussian graphical model

Category | Method | Model Name | Model structure ${\bm{\Omega}}$ | Penalty function/constraint
---|---|---|---|---
Coarse -grained |  Penalized MLE | JGL (Guo et al., 2011) | $\left\\{\begin{array}[c]{@{}l@{}}\omega_{i,j}^{(k)}=\theta_{i,j}\gamma_{i,j}^{(k)},\;i\neq j;\\\ \omega_{i,i}^{(k)}=\gamma_{i,i}^{(k)},\;i=j.\end{array}\right.$ | $\lambda_{1}\sum_{i\neq j}\theta_{i,j}-\lambda_{2}\sum_{k=1}^{K}\sum_{i\neq j}|\gamma_{i,j}^{(k)}|$
JWLGL (Shan et al., 2020) |  $\omega_{i,j}^{(k)}=\theta_{m,m^{\prime}}^{(k)}\gamma_{i,j}^{m,m^{\prime},(k)}.$ If $m\neq m^{\prime}$: $\theta_{m,m^{\prime}}^{(k)}=\alpha_{m,m^{\prime}}\beta_{m,m^{\prime}}^{(k)}$, else if $m=m^{\prime}$: $\theta_{m,m}^{(k)}=1$ and $\gamma_{i,j}^{m,m,(k)}=\left\\{\begin{array}[]{ll}\iota_{i,j}^{m,m}\rho_{i,j}^{m,m,(k)},&i\neq j;\\\ 1,&i=j\end{array}\right.$ |  $\lambda_{1}\sum_{m\neq m^{\prime}}\alpha_{m,m^{\prime}}$ $+\lambda_{2}\sum_{m\neq m^{\prime}}\sum_{k=1}^{K}\sum_{i\neq j}|\beta_{m,m^{\prime}}^{(k)}\gamma_{i,j}^{m,m^{\prime},(k)}|$ $+\lambda_{3}\sum_{m=1}^{M}\sum_{i\neq j}\iota_{i,j}^{m,m}$ $+\lambda_{4}\sum_{m=1}^{M}\sum_{i\neq j}\sum_{k=1}^{K}|\rho_{i,j}^{m,m,(k)}|$
FGL (Danaher et al., 2014) | $\omega_{i,j}^{(k)}$ | $\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq j}|\omega_{i,j}^{(k)}|+\lambda_{2}\sum_{k<k^{\prime}}\sum_{i,j}|\omega_{i,j}^{(k)}-\omega_{i,j}^{(k^{\prime})}|$
| GGL (Danaher et al., 2014) | $\omega_{i,j}^{(k)}$ | $\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq j}|\omega_{i,j}^{(k)}|+\lambda_{2}\sum_{i\neq j}(\sum_{k=1}^{K}\omega_{i,j}^{(k)^{2}})^{\frac{1}{2}}$
| GFGL (Gibberd and Nelson, 2017) | $\omega_{i,j}^{(k)}$ | $\lambda_{1}\sum_{t=1}^{T}\sum_{i\neq j}|\omega_{i,j}^{(t)}|+\lambda_{2}\sum_{t=2}^{T}\|{\bm{\Omega}}^{(t)}_{-ii}-{\bf\Omega}_{-ii}^{(t-1)}\|_{F}$
| JAGL (Shan and Kim, 2018) | $\omega_{i,j}^{(k)}$ |  $\sum_{k=1}^{K}\frac{1}{n_{k}}\sum_{i\neq j}\frac{1}{|(1-\pi)\widehat{t}_{i,j}+\pi\widehat{s}_{i,j}^{(k)}|^{r}}|\omega_{i,j}^{(k)}|$111$n_{k}$ denotes the number of samples of group $k$, $\widehat{t}_{i,j}$ is the precision matrix estimated by pooling all samples across groups, $\widehat{s}_{i,j}$ is the precision matrix estimated by individual group, and $r>0$
| TFRE (Bilgrau et al., 2020) | $\omega_{i,j}^{(k)}$ |  $\sum_{k=1}^{K}\frac{\lambda_{k}}{2}\|{\bm{\Omega}}^{(k)}-{\bf T}^{(k)}\|_{F}^{2}$ $+\sum_{k_{1},k_{2}}^{K}\frac{\lambda_{k_{1},k_{2}}}{4}\|({\bm{\Omega}}^{(k_{1})}-{\bf T}^{(k_{1})})-({\bm{\Omega}}^{(k_{2})}-{\bf T}^{(k_{2})})\|_{F}^{2}$
| RCON (Mohan et al., 2012, 2014) | $\omega_{i,j}^{(k)}$ |  $G_{q}({\bm{\Omega}}^{(1)}-{\bm{\Omega}}^{(2)})=\min_{V:{\bm{\Omega}}^{(1)}-{\bm{\Omega}}^{(2)}=V+V^{\top}}f(V)$ $f(V)=\sum_{j=1}^{p}\|V_{j}\|_{q}$
|  CLIME (Cai et al., 2011) | JEMP (Lee and Liu, 2015) | $\omega_{i,j}^{(k)}=\theta_{i,j}+\gamma_{i,j}^{(k)}$ |  $|\frac{1}{K}\sum_{k=1}^{K}\\{\widehat{{\bm{\Sigma}}}^{(k)}({\bm{\Theta}}+{\bm{\Gamma}}^{(k)})-{\bm{I}}\\}|_{\infty}\leq\lambda_{1}$ $|\widehat{{\bm{\Sigma}}}^{(k)}({\bm{\Theta}}+{\bm{\Gamma}}^{(k)})-{\bm{I}}|_{\infty}\leq\lambda_{2},\;\sum_{k=1}^{K}{\bm{\Gamma}}^{(k)}=0$
| KSE (Qiu et al., 2016) | $\omega_{i,j}^{(k)}$ | $|\widehat{{\bm{S}}}(u_{0}){\bm{\Omega}}(u_{0})-{\bm{I}}|_{\infty}\leq\xi$
Fine -grained |  Penalized MLE | LASICH (Saegusa and Shojaie, 2016) | $\omega_{i,j}^{(k)}$ |  $\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq j}|\omega_{i,j}^{(k)}|$ $-\lambda_{1}\lambda_{2}\sum_{i\neq j}(\sum_{k,k^{\prime}}^{K}W_{k,k^{\prime}}(\omega_{i,j}^{(k)}+\omega_{i,j}^{(k^{\prime})})^{2})^{\frac{1}{2}}$
GEN-ISTA (Price et al., 2021) | $\omega_{i,j}^{(k)}$ |  $\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq j}|\omega_{i,j}^{(k)}|$ $+\lambda_{2}\sum_{q=1}^{Q}\frac{1}{|D_{q}|}\sum_{k,k^{\prime}\in D_{q}}\|{\bm{\Omega}}^{(k)}-{\bm{\Omega}}^{(k)^{\prime}}\|_{F}^{2}$
|  Neighborhood selection |  JSEM (Ma and Michailidis, 2016) JMMLE (Majumdar and Michailidis, 2018) | $\theta_{i,j}^{(k)}=-\omega_{i,j}^{(k)}/\omega_{i,i}^{(k)}$ | $\sum_{j\neq i}\sum_{g\in\mathcal{G}_{i,j}}\lambda_{i,j}^{[g]}\|{\bm{\theta}}_{i,j}^{[g]}\|_{2}$

1 $n_{k}$ denotes the number of samples of group $k$, $\widehat{t}_{i,j}$ is
the precision matrix estimated by pooling all samples across groups,
$\widehat{s}_{i,j}$ is the precision matrix estimated by an individual group,
and $r>0$

##  3 Joint Graphical Models using Coarse-grained Structure

We outline a variety of approaches for joint graphical model estimation that
use prior knowledge of coarse-grained structure across groups. The coarse-
grained structure, as opposed to fine-grained structure, assumes that the edge
structure of each graph has the same degree of similarity among other graphs.
In contrast, fine-grained structure assumes local similarity that is only
shared within a subset of groups. The performance of the coarse-grained
estimation procedure is improved using regularization that captures the common
structure across the $K$ groups – enabling the use of shared information
across groups. We will discuss two directions in detail: hierarchical
regularizers and analogous Bayesian priors.

###  3.1 Joint Graphical Models with Hierarchical Structure

Guo et al. (2011) studied joint estimation of related precision matrices,
where the precision matrices are assumed to be related through a hierarchical
structure. Specifically, each entry in the precision matrix is the
multiplication of a common component across $K$ groups and an individual
component: $\omega_{i,j}^{(k)}=\theta_{i,j}\gamma_{i,j}^{(k)}$ for $i\neq j$
and $\omega_{i,i}^{(k)}=\gamma_{i,i}^{(k)}$, where $\theta_{i,j}$ is the
shared structure and $\gamma_{i,j}^{(k)}$ is the group-specific component.
Thus, this approach enforces a common background structure. To encourage
sparsity, an $\ell_{1}$-norm penalty term is also included as a regularizer,
resulting in the following objective termed Joint Graphical Lasso (JGL):

$\displaystyle\widehat{{\bm{\Theta}}},\\{\widehat{\bm{\Gamma}}^{(k)}\\}_{k=1}^{K}=\mathop{\mathrm{argmax}}\;\tilde{\ell}({\bm{\Omega}})-\lambda_{1}\sum_{i\neq
j}\theta_{i,j}-\lambda_{2}\sum_{k=1}^{K}\sum_{i\neq j}|\gamma_{i,j}^{(k)}|,$

where $\lambda_{1},\lambda_{2}$ are hyperparameters that control the scale of
the penalty. Note that even when the common component $\theta_{i,j}$ is
nonzero, an individual entry $\omega_{i,j}^{(k)}$ can still be set to zero by
the $\ell_{1}$ penalty, which denotes a missing edge in the associated graph.

Shan et al. (2020) proposed a joint two-level graphical lasso (JWLGL), which
is a more expressive model that constructs two-level structures on both the
set of common components and individual components. It further clusters the
set of nodes $V$ into $M$ classes and imposes class specific structure: let
$m$ and $m^{\prime}$ be the classes to which nodes $i$ and $j$ belong,
respectively. If $m\neq m^{\prime}$, we have:

$\displaystyle\omega_{i,j}^{(k)}=\theta_{m,m^{\prime}}^{(k)}\gamma_{i,j}^{m,m^{\prime},(k)};\qquad\theta_{m,m^{\prime}}^{(k)}=\alpha_{m,m^{\prime}}\beta_{m,m^{\prime}}^{(k)}.$

If $m=m^{\prime}$, we have $\theta_{m,m}^{(k)}=1$ and

$\displaystyle\omega_{i,j}^{(k)}=\theta_{m,m}^{(k)}\gamma_{i,j}^{m,m,(k)};\qquad\gamma_{i,j}^{m,m,(k)}=\left\\{\begin{array}[]{ll}\iota_{i,j}^{m,m}\rho_{i,j}^{m,m,(k)},&i\neq
j;\\\ 1,&i=j.\end{array}\right.$

Without loss of generality, we assume $\alpha_{m,m^{\prime}}\geq 0$ and
$\iota_{i,j}^{m,m}\geq 0$ for $i\neq j$ and $m\neq m^{\prime}$. Here,
$\alpha_{m,m^{\prime}}$ and $\iota_{i,j}^{m,m}$ denote the common components
shared across $K$ groups, while $\beta_{m,m^{\prime}}^{(k)}$ and
$\rho_{i,j}^{m,m,(k)}$ denote the individual components that vary across
groups.

Lee and Liu (2015) proposed a joint estimator of multiple precision matrices
(JEMP) under an assumption that precision matrices decompose into the sum of
two components: ${\omega}^{(k)}_{i,j}={\theta}_{i,j}+{\gamma}_{i,j}^{(k)}$. In
contrast to the maximum likelihood, the estimation procedure of JEMP is
motivated by the CLIME estimator (Cai et al., 2011), which estimates a single
precision matrix by solving the following optimization problem:

$\displaystyle\widehat{{\bm{\Omega}}}^{(k)}=\mathop{\mathrm{argmin}}\|{\bm{\Omega}}^{(k)}\|_{1}\qquad\text{subject
to}\qquad|\widehat{{\bm{\Sigma}}}^{(k)}{\bm{\Omega}}^{(k)}-{\bm{I}}|_{\infty}\leq\xi,$
(4)

where $\xi$ is a tuning parameter. The CLIME estimator finds a sparse
$\widehat{\bm{\Omega}}^{(k)}$ while ensuring that
$\widehat{{\bm{\Sigma}}}^{(k)}\widehat{\bm{\Omega}}^{(k)}$ is close to an
identity matrix. JEMP can be seen as a generalization of CLIME to a multi-
group setting as it solves the following optimization problem:

$\displaystyle\widehat{{\bm{\Theta}}},\\{\widehat{\bm{\Gamma}}^{(k)}\\}_{k=1}^{K}=\mathop{\mathrm{argmin}}\|{\bm{\Theta}}\|_{1}+v\sum_{k=1}^{K}\|{\bm{\Gamma}}^{(k)}\|_{1},$
$\displaystyle\text{subject to
}\left|\frac{1}{K}\sum_{k=1}^{K}\\{\widehat{{\bm{\Sigma}}}^{(k)}({\bm{\Theta}}+{\bm{\Gamma}}^{(k)})-{\bm{I}}\\}\right|_{\infty}\leq\lambda_{1},\;\left|\widehat{{\bm{\Sigma}}}^{(k)}({\bm{\Theta}}+{\bm{\Gamma}}^{(k)})-{\bm{I}}\right|_{\infty}\leq\lambda_{2},\;\sum_{k=1}^{K}{\bm{\Gamma}}^{(k)}=0,$

where ${\bm{\Theta}}$ denotes the common structure, i.e., the mean of the
precision matrices ${K}^{-1}\sum_{k=1}^{K}{\bm{\Omega}}^{(k)}$, and
${\bm{\Gamma}}^{(k)}$ denotes the individual residual components
${\bm{\Omega}}^{(k)}-{\bf\Theta}$. In the above optimization problem, the
first constraint regularizes the average difference and the second constraint
regularizes the individual difference. Thus, the first constraint imposes a
common structure across groups. The pre-specified weight $v$ controls the
degree of uniqueness of each group, while $\lambda_{1},\lambda_{2}$ are
hyperparameters that measure group average and individual estimation quality,
respectively.

(a)

(b)

Figure 2: Graphical models with shared coarse-grained structure across groups.
(a) The black lines denote the common edges, while green lines denote
individual edges. (b) The corresponding adjacency matrices for each graph
shown in (a).

###  3.2 Regularization Approaches for Modeling Joint Structure

Regularization-based approaches (Danaher et al., 2014; Bilgrau et al., 2020;
Hao et al., 2018; Shan and Kim, 2018) do not assume the form of the common
structure and individual structure, but instead impose similarity constraints
across groups. For example, the fused graphical lasso (FGL) and the group
graphical lasso (GGL) (Danaher et al., 2014) add convex penalty terms to the
log-likelihood function to learn a common structure:

$\displaystyle\widehat{{\bm{\Omega}}}=\mathop{\mathrm{argmax}}\;\tilde{\ell}({\bm{\Omega}})-{\bm{P}}({\bm{\Omega}});$
$\displaystyle{\bm{P}}_{\text{FGL}}({\bm{\Omega}})=\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq
j}|\omega_{i,j}^{(k)}|+\lambda_{2}\sum_{k<k^{\prime}}\sum_{i,j}|\omega_{i,j}^{(k)}-\omega_{i,j}^{(k^{\prime})}|;$
(5)
$\displaystyle{\bm{P}}_{\text{GGL}}({\bm{\Omega}})=\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq
j}|\omega_{i,j}^{(k)}|+\lambda_{2}\sum_{i\neq
j}\left(\sum_{k=1}^{K}\omega_{i,j}^{(k)^{2}}\right)^{\frac{1}{2}}.$ (6)

The first penalty term in both ${\bm{P}}_{\text{FGL}}$ and
${\bm{P}}_{\text{GGL}}$ encourages model sparsity. The second term in
${\bm{P}}_{\text{FGL}}$ encourages groups to have shared edge values, while
the ${\bm{P}}_{\text{GGL}}$ penalty tends to be less restrictive and only
encourages a shared sparsity pattern.

Hao et al. (2018) proposed simultaneous clustering and estimation (SCAN)
procedure that addresses the case when the heterogeneous data are missing
group labels, e.g., when the groups are latent or unknown. SCAN partitions the
unlabeled data into $K$ clusters and simultaneously imposes a homogeneous
structure across groups using the GGL penalty. Given $n$ unlabeled
observations ${\bm{x}}_{i},\;i=1,\ldots,n$, with the density function
$l({\bm{x}},{\bm{\mu}},{\bm{\Omega}})=\sum_{k=1}^{K}\pi_{k}l_{k}({\bm{x}};{\bm{\mu}}^{(k)},{\bm{\Omega}}^{(k)})$,
where $\pi_{k}$ is the probability that ${\bm{x}}_{i}$ belongs to the $k$-th
group and
$l_{k}({\bm{x}};{\bm{\mu}}^{(k)},{\bm{\Omega}}^{(k)})=\log\\{\det({\bm{\Omega}}^{(k)})\\}-\mathop{\mathrm{tr}}\\{({\bm{x}}-{\bm{\mu}}^{(k)})({\bm{x}}-{\bm{\mu}}^{(k)})^{\top}{\bm{\Omega}}^{(k)}\\}$.
SCAN solves the following optimization problem:

$\displaystyle\\{\widehat{\pi}_{k}\\}_{k=1}^{K},\widehat{\bm{\mu}},\widehat{\bm{\Omega}}=\mathop{\mathrm{argmax}}\frac{1}{n}\sum_{i=1}^{n}\log\left(l({\bm{x}}_{i},{\bm{\mu}},{\bm{\Omega}})\right)-{\bm{P}}_{SCAN}({\bm{\mu}},{\bm{\Omega}});$
$\displaystyle{\bm{P}}_{SCAN}({\bm{\mu}},{\bm{\Omega}})=\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq
j}|\omega_{i,j}^{(k)}|+\lambda_{2}\sum_{i\neq
j}\left(\sum_{k=1}^{K}\omega_{i,j}^{(k)^{2}}\right)^{\frac{1}{2}}+\lambda_{3}\sum_{k=1}\sum_{i=1}^{p}|\mu_{i}^{(k)}|.$

Note that the first two terms of ${\bm{P}}_{SCAN}({\bm{\mu}},{\bm{\Omega}})$
correspond to the GGL penalty function and the third term is the
$\ell_{1}$-norm penalty, used to encourage sparsity.

There are several other methods that use the regularization approach to
estimate the joint structure. Shan and Kim (2018) proposed the joint adaptive
graphical lasso (JAGL) procedure that introduces a weighted $\ell_{1}$ penalty
term to tackle problems with unbalanced data. Bilgrau et al. (2020) proposed
the targeted fused ridge estimator (TFRE) that uses an additional $\ell_{2}$
penalty term that incorporates target matrices as prior information to
stabilize the estimation process. Finally, the joint greedy equivalence search
(JointGES) procedure uses an $\ell_{0}$ penalty term to jointly estimate
directed acylic graphs (Wang et al., 2020b).

###  3.3 Bayesian Methods for Joint Gaussian Graphical Model Estimation

We briefly overview Bayesian methods for joint Gaussian graphical model
estimation. Bayesian formulations of graphical models use priors to encourage
desired properties for model selection. For example, the Spike-and-slab prior
is commonly used in practice to encourage sparsity in precision matrices. In
general, the probabilistic counterpart to the penalty function follows the
relation $p({\bm{\Omega}})\propto\exp(-{\bm{P}}({\bm{\Omega}}))$, where
$p({\bm{\Omega}})$ is the prior distribution of the precision matrix and
${\bm{P}}({\bm{\Omega}})$ denotes the penalty function. By Bayes’ rule, the
posterior distribution is proportional to the product of the likelihood
distribution and prior distribution. Therefore, finding the maximum a
posteriori (MAP) estimate is equivalent to obtaining the estimate by
maximizing the log-likelihood (MLE) with an additional penalty function.
Several works (Tan et al., 2017; Li et al., 2019; Yajima et al., 2014; Mitra
et al., 2016; Peterson et al., 2015; Lin et al., 2017) have addressed Bayesian
graphical model estimation by designing priors that incorporate structural
information. In the Gaussian graphical model case, a Wishart prior (Roverato,
2002; Atay-Kayis and Massam, 2005; Lenkoski and Dobra, 2011; Mohammadi et al.,
2015) is often placed on the precision matrix. This prior is desirable due to
conjugacy, i.e., for a Gaussian likelihood, the posterior distribution remains
Wishart. Furthermore, the normalizing constant of the posterior distribution
has an explicit form when the graph $G$ is decomposable, that is, when the
index set $V$ of a graph can be partitioned into three disjoint nonempty sets
$V=A\cup S\cup B$ and (i) $S$ is a clique, (ii) $S$ separates $A$ and $B$,
(iii) $A\cup S$ and $S\cup B$ form decomposable subgraphs.

Li et al. (2019) proposed the joint spike-and-slab graphical lasso prior,
designed to encourage global sparse structure. In related work, Tan et al.
(2017) placed a multiplicative prior on the adjacency matrices, where the
distribution of each edge depends on the multiplication of the values of two
end nodes. This prior not only encourages sparsity, but also allows users to
specify the degree of connections. We discuss local common structure methods
(Peterson et al., 2015; Lin et al., 2017) in Section 4.3 and the differential
graph methods (Yajima et al., 2014; Mitra et al., 2016) in Section 5.2. We
will not go into details about the Bayesian formulation of graphical models,
but instead give a high-level overview of various approaches.

Consider a single group setting with
${\bm{\Omega}}=\\{{\bm{\Omega}}^{(1)}\\}$. From a Bayesian perspective, the
Lasso regularizer can be viewed as a Laplace prior (Marlin and Murphy, 2009;
Wang, 2012) and is formulated as:

$p({\bm{\Omega}}|\lambda)=C^{-1}\prod_{i\neq
j}\frac{\lambda}{2}\exp\left(-\lambda|\omega_{i,j}|\right)\prod_{i=1}^{p}\left\\{\frac{\lambda}{2}\exp\left(-\frac{\lambda}{2}\omega_{i,i}\right)\mathbbm{1}_{(\omega_{i,i}>0)}\right\\}\mathbbm{1}({\bm{\Omega}}\succ
0),$

where $C$ is the normalizing constant and $\mathbbm{1}({\bm{\Omega}}\succ 0)$
restricts the precision matrix to be positive definite. The term
$\mathbbm{1}_{(\omega_{i,i}>0)}$ ensures that the diagonal entries are non-
negative and hence preserves the positive definiteness of ${\bm{\Omega}}$.
When taking the logarithm $\log p({\bm{\Omega}}|\lambda)$, the first product
is equal to the lasso regularizer. Therefore, when computing the MAP estimate,
the logarithm of the Laplace prior along with the log-likelihood is the
penalized MLE estimator with Lasso penalty function. In the multiple group
case, to promote the group similarity between the precision matrices, Li et
al. (2019) converted the GGL and FGL penalties to structural priors.

In Bayesian inference, other than computing the MAP estimator, we are also
interested in the posterior mean, mode, and samples. In this case, the
shrinkage priors are not enough to produce sparse posterior sample (or mean
and mode) because the posterior does not concentrate on sparse parameters.
Therefore, additional thresholding is required to obtain sparsity. As an
alternative, one may use the spike-and-slab prior (Mitchell and Beauchamp,
1988) to promote the sparsity pattern in the posterior. Consider a single
group ${\bm{\Omega}}=\\{{\bm{\Omega}}^{(1)}\\}$, the spike-and-slab prior is a
hierarchical mixture prior formulated as:

$\displaystyle p({\bf z}|\lambda)=\prod_{i\neq
j}\operatorname{Ber}(z_{i,j}|\lambda);$ $\displaystyle p({\bm{\Omega}}|{\bf
z})=\prod_{i\neq
j}(1-z_{i,j})\delta(\omega_{i,j})+z_{i,j}\mathcal{N}(\omega_{i,j}|0,\sigma^{2}),$

where $\delta(\cdot)$ denotes the delta function. If $z_{i,j}=0$, we restrict
the variable to be zero. One may also replace the delta function with a normal
distribution with small variance, which approximates the delta function. In
the multiple group case, a set of latent indicators following the spike-and-
slab distribution adaptively control the value of the FGL (resp., GGL) penalty
(Li et al., 2019), namely the Doubly Spike-and-Slab Joint Graphical Lasso
(DSS-JGL). Consider two constants $v_{1}>v_{0}>0$ and $z_{i,j}$, $w_{i,j}$ are
binary variables for $i\neq j$. We assume that each $z_{i,j}$ and $w_{i,j}$
are drawn independently from a Bernoulli distribution with a specific
parameter. The DSS-JGL prior is represented as:

$\displaystyle p({\bm{\Omega}}|{\bf z},{\bf
w})=\lambda_{0}\sum_{k=1}^{K}\sum_{i=1}^{p}|\omega_{i,i}^{(k)}|+\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq
j}\frac{|\omega_{i,j}^{(k)}|}{v_{z_{i,j}}}+\lambda_{2}\sum_{i\neq
j}\frac{{\bm{P}}_{*}({\bm{\Omega}})}{v_{(w_{i,j}z_{i,j}})},$ (7)

where ${\bm{P}}_{*}({\bm{\Omega}})$ can be
${\bm{P}}_{\text{FGL}}({\bm{\Omega}})$ or
${\bm{P}}_{\text{GGL}}({\bm{\Omega}})$ defined in (5) and (6), respectively.
We can choose $v_{0}$ to be small, so that when $z_{i,j}=0$ for $i\neq j$, the
second term in (7) will be large, forcing the posterior to be zero. Similar
behavior also follows for the joint regularization term in (7) when either
$z_{i,j}$ or $w_{i,j}$ is zero.

##  4 Joint Graphical Models using Fine-grained Structure

(a)

(b)

Figure 3: Graphical models with shared fine-grained structure across groups.
(a) The black lines denote the common edges. Green colored lines represent the
common structure of Group 1and Group 4, yellow lines of Group 1 and Group 3,
red lines of Group 2 and Group 3, and blue lines of Group 2 and Group 4. (b)
The corresponding adjacency matrices for each graph shown in (a), using the
same colors for groups of shared edges.

Coarse-grained joint graphical models assume that all groups are equally
related; the relation between $\omega_{i,j}^{(k)}$ and
$\omega_{i,j}^{(k^{\prime})}$ and the relation between $\omega_{i,j}^{(k)}$
and $\omega_{i,j}^{(k^{\prime\prime})}$, $k\neq k^{\prime}$,$k^{\prime}\neq
k^{\prime\prime}$ are assumed to be equal. However, in many real world
settings, some subsets of groups share local structure that does not appear
across all groups. Figure 3 visualizes the adjacency matrices related to edge
sets ${E}^{k},\;k=1,\ldots,4$, of graphical models that are not globally
related, but with subsets of the groups that share a local structure. The most
common approach in modeling such graphical models is to incorporate this prior
knowledge of the relation between subgroups into the estimators (Majumdar and
Michailidis, 2018; Ma and Michailidis, 2016; Saegusa and Shojaie, 2016). We
briefly outline some approaches for joint graphical model estimation with
fine-grained shared structure.

###  4.1 Entry-wise Structural Information

Given the relation information $\mathcal{G}=\cup_{1\leq i<j\leq
p}\mathcal{G}_{i,j}$, where $\mathcal{G}_{i,j}$ is a set that encodes the
group relations of node $i$ and node $j$, Figure 3 (a) illustrates an example
of entry-wise structural information. Consider the pair of nodes $(i,j)=(3,4)$
in Figure 3 (a): both graphs of Group 1 and Group 2 have edges (yellow lines)
connecting these two nodes, while Group 3 and Group 4 do not have an edge in
between. Then the corresponding relation information
$\mathcal{G}_{3,4}=\mathcal{G}_{4,3}$ is $\\{\\{1,3\\},\\{2,4\\}\\}$.
Similarly, for the pair of nodes $(i,j)=(2,8)$, the graphs of Group 1 and
Group 4 have an edge (green line) connecting $(2,8)$, while graphs of Group 2
and Group 3 do not have edge between node $(2,8)$. Then, the corresponding
relation information $\mathcal{G}_{2,8}=\mathcal{G}_{8,2}$ is
$\\{\\{1,4\\},\\{2,3\\}\\}$. The relation information of graphs in Figure 3 is
$\mathcal{G}_{i,j}=\\{\\{1,3\\},\\{2,4\\}\\}$ for
$(i,j)\in\\{(3,4),(4,3),(4,5),(5,4),(1,2),(2,1),(2,5),(5,2)\\}$;
$\mathcal{G}_{i,j}=\\{\\{1,4\\},\\{2,3\\}\\}$ for
$(i,j)\in\\{(2,8),(8,2),(3,8),(8,3),(1,8),(8,1),(6,8),(8,6)\\}$; and
$\mathcal{G}_{i,j}=\\{\\{1,2,3,4\\}\\}$ for the rest of the node pairs.

Ma and Michailidis (2016) proposed a joint structural estimation method (JSEM)
to estimate edge sets $\\{\widehat{E}^{(k)}\\}_{k=1}^{K}$ by modifying the
neighborhood selection algorithm (Meinshausen et al., 2006) to incorporate
structural information. Once the edge sets are estimated, each precision
matrix is estimated by maximizing the group specific likelihood constrained to
have zeros indexed by $\\{\widehat{E}^{(k)}\\}_{k=1}^{K}$. We will briefly
introduce the neighborhood selection method and then show how JSEM extends it
to multiple graphical models.

The neighborhood selection algorithm estimates the conditional independence
structure, which is encoded in the pattern of zeros of the precision matrix
under a Gaussian model, by solving a collection of regression problems. See
Chapter 12 in (Maathuis et al., 2018) and references therein. Suppose the
$p$-dimensional random variable ${\bm{x}}$ follows a normal distribution
$\mathcal{N}({\bm{\mu}},({\bm{\Omega}})^{-1})$ with an associated graph
$G=(V,E)$. In this case, we can express $x_{i}$, $i\in V$, as a linear
function of other nodes:

$x_{i}=\sum_{j\in
V\text{\textbackslash}\\{i\\}}\theta_{ij}x_{j}+\varepsilon_{j},$

where $\varepsilon_{j}$ is independent of $x_{j}$, $j\in
V\text{\textbackslash}\\{i\\}$, if and only if
$\theta_{ij}=-\omega_{i,j}/\omega_{i,i}$. Therefore, the optimal prediction of
${x}_{i}$ given the remaining variables can be formulated as the following
optimization problem:

$\displaystyle\\{\widehat{\theta}_{i,j}\\}_{j\in
V\text{\textbackslash}\\{i\\}}=\mathop{\mathrm{argmin}}\;\mathbb{E}\left({x}_{i}-\sum_{j\in
V\text{\textbackslash}\\{i\\}}\theta_{i,j}{x}_{j}\right)^{2}.$ (8)

Since the distribution of ${\bm{x}}$ is not known, the expectation term in (8)
can not be computed. Let ${\bm{X}}$ be a $n\times p$ matrix that collects $n$
i.i.d. observations of ${\bm{x}}$. The $\ell_{1}$ penalized empirical
optimization objective is the given as:

$\displaystyle\widehat{\bm{\Theta}}_{i}=\mathop{\mathrm{argmin}}_{{\bm{\Theta}}_{i}:\;\theta_{i,i}=0}\frac{1}{n}\|{\bm{X}}_{i}-{\bm{X}}{\bm{\Theta}}_{i}\|_{2}^{2}+\lambda\|{\bm{\Theta}}_{i}\|_{1},$
(9)

where ${\bm{\Theta}}_{i}$ is a $p$-dimensional vector
${\bm{\Theta}}_{i}=(\theta_{i,1},\ldots,\theta_{i,p})^{\top}$,
$\theta_{i,i}=0$, and ${\bm{X}}_{i}$ is the $i$-th column of ${\bm{X}}$. To
estimate multiple models, JSEM extends the neighborhood selection approach in
(9) by solving the following optimization problem:

$\displaystyle\widehat{\bm{\Theta}}=\mathop{\mathrm{argmin}}_{{\bm{\Theta}}_{i}^{(k)}:\theta_{i,i}^{(k)}=0,\;k=1,\ldots,K}\sum_{k=1}^{K}\frac{1}{n_{k}}\|{\bm{X}}_{i}^{(k)}-{\bm{X}}^{(k)}{\bm{\Theta}}_{i}^{(k)}\|_{2}^{2}+2{\bm{P}}_{JSEM}(\\{{\bm{\Theta}}_{i}^{(k)}\\}_{k=1}^{K});$
(10)
$\displaystyle{\bm{P}}_{JSEM}(\\{{\bm{\Theta}}_{i}^{(k)}\\}_{k=1}^{K})=\sum_{j\neq
i}\sum_{g\in\mathcal{G}_{i,j}}\lambda_{i,j}^{[g]}\|{\bm{\theta}}_{i,j}^{[g]}\|_{2},$

where ${\bm{\Theta}}=\\{{\bm{\Theta}}^{(1)},\ldots,{\bm{\Theta}}^{(K)}\\}$,
the penalty term incorporates the relation information $\mathcal{G}$, and
$\lambda_{i,j}^{[g]}$ is the group-specific tuning parameter. The penalty
function ${\bm{P}}_{JSEM}$ resembles the GGL penalty in (6), except that the
norm is now placed on a subset of groups provided by the prior information
$\mathcal{G}$. Minimizing (10) results in the following estimated edge sets

$\widehat{E}^{(k)}=\\{(i,j):1\leq i<j\leq p,\widehat{\theta}_{i,j}^{(k)}\neq
0\text{ or }\widehat{\theta}_{j,i}^{(k)}\neq 0\\},\quad k=1,\ldots,K.$

Given $\widehat{E}^{(k)}$, we define
$\mathcal{S}^{+}_{\widehat{E}^{(k)}}=\\{{\bm{\Omega}}:{\bm{\Omega}}\succ
0,\;\omega_{i,j}=0,\;\forall(i,j)\not\in\widehat{E}^{(k)}\text{ and }i\neq
j\\}$. The precision matrix set ${\bm{\Omega}}$ is finally estimated by
maximizing the log-likelihood with constraints that enforce the sparsity
pattern:

$\displaystyle\widehat{\bm{\Omega}}^{(k)}=\mathop{\mathrm{argmax}}_{{\bm{\Omega}}^{(k)}\in\mathcal{S}^{+}_{\widehat{E}^{(k)}}}\log\\{\det({\bm{\Omega}}^{(k)})\\}-\mathop{\mathrm{tr}}(\widehat{\bm{\Sigma}}^{(k)}{\bm{\Omega}}^{(k)}),\quad
k=1,\ldots,K.$ (11)

Note that we can apply JSEM only when element-wise structural relation
information is given. However, when not all element-wise structural relation
information is provided, one can still use the GGL penalty on subsets of
groups for which prior information is available.

###  4.2 Group-wise Structural Information

As obtaining entry-wise structural information is sometimes challenging,
another approach is to use the relationship between groups, described by
weights, in estimation. For example, suppose that $K=3$ and we have the
following relationships between groups: Group $1$ and Group $2$ share
similarity, Group $1$ and Group $3$ also share similarity, but Group $2$ and
Group $3$ are unlikely to be similar. To this end, Saegusa and Shojaie (2016)
proposed LASICH (Laplacian shrinkage for inverse covariance matrices from
heterogeneous populations), that uses a weighted graph $G_{g}=(\Gamma,E,W)$ to
record the group relation. The node set $\Gamma$ denotes $K$ groups, the edge
set $E$ captures the relations between groups, and the weight set
$W:E\rightarrow\mathbb{R}_{+}$ represents the similarities between groups.
Then, a Laplacian penalty function is placed on the objective function to
impose group similarity. The optimization problem is formulated as follows:

$\displaystyle\widehat{\bm{\Omega}}=\mathop{\mathrm{argmax}}\;\tilde{\ell}({\bm{\Omega}})-{\bm{P}}_{LASICH}({\bm{\Omega}}),$
$\displaystyle\text{subject to
}{\bm{\Omega}}^{(k)}=({\bm{\Omega}}^{(k)})^{\top},\;{\bm{\Omega}}^{(k)}\in\mathcal{S}\succ
0,\quad k=1,\ldots,K;$
$\displaystyle{\bm{P}}_{LASICH}({\bm{\Omega}})=\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq
j}|\omega_{i,j}^{(k)}|-\lambda_{1}\lambda_{2}\sum_{i\neq
j}\left(\sum_{k,k^{\prime}}^{K}W_{k,k^{\prime}}\left(\omega_{i,j}^{(k)}+\omega_{i,j}^{(k^{\prime})}\right)^{2}\right)^{\frac{1}{2}},$

where the first penalty term in ${\bm{P}}_{LASICH}({\bm{\Omega}})$ promotes
the sparsity pattern and the second term encourages the similarities within
subsets of groups. Although the three approaches (Majumdar and Michailidis,
2018; Ma and Michailidis, 2016; Saegusa and Shojaie, 2016) discussed here
require prior information on the group relation, they provide a more flexible
structure than the global structure methodologies, such as GGL (Danaher et
al., 2014) and (Guo et al., 2011). In particular, the global structure
approach can be viewed as a special case of a local structure approach with
homogeneous structural information.

While the work introduced earlier required prior information about the group
relation, which may not be available for most cases, Graphical Elastic Net
Iterative Shrinkage Thresholding Algorithm (GEN-ISTA) (Price et al., 2021)
jointly estimates graphs and group members via k-means clustering. Assume that
we want to further cluster group $k=1,\ldots,K$ into $Q$ classes, and we
denote $D_{q}$ for $q=1,\ldots,Q$ as the set that contains the group
membership. Its objective function is as follows

$\displaystyle\mathop{\mathrm{argmin}}_{{\bm{\Omega}},{\bm{D}}}\widetilde{\ell}({\bm{\Omega}})-{\bf
P}_{GEN-ISTA}({\bm{\Omega}},{\bm{D}});$ (12) $\displaystyle{\bf P}_{GEN-
ISTA}({\bm{\Omega}},{\bm{D}})=\lambda_{1}\sum_{k=1}^{K}\sum_{i\neq
j}|\omega_{i,j}^{(k)}|+\lambda_{2}\sum_{q=1}^{Q}\frac{1}{|D_{q}|}\sum_{k,k^{\prime}\in
D_{q}}\|{\bm{\Omega}}^{(k)}-{\bm{\Omega}}^{(k)^{\prime}}\|_{F}^{2}.$

GEN-ISTA can be implemented via alternating minimization. It is easy to see
that when we fix ${\bm{\Omega}}$ and optimize ${\bm{D}}$, then (12) becomes a
k-means clustering objective.

###  4.3 Bayesian Approach

We introduce two Bayesian methods that construct priors to constrain the
similarities within subsets of groups. Unlike the structural regularization
approaches introduced in Section 4.1, both methods we discuss here do not
require specific prior information on each graph but employ a general prior
distribution to build the group relations. This property is particularly
desirable because, in most cases, we may not have the structural information
for $K$ groups of data.

Peterson et al. (2015) proposed the Markov Random Field (MRF) prior for the
graphs $G^{(k)}=(V,E^{(k)}),k=1,\ldots,K$ to encourage the selection of edge
indices in related graphs. In addition, the selection of edge indices is
controlled by two variables: a random matrix $\Theta$, which measures the
degree of similarity between groups and an edge-specific coefficient vector
${\bm{v}}$ reflecting the probability of the corresponding edge being
selected. Higher value of $v_{i,j}$ means that the edge connecting the $i$-th
node and $j$-th node is more likely to be selected. As a result, ${\bm{v}}$ is
capable of controlling the sparsity of graphs. The Beta prior is placed on
${\bm{v}}$, in which the distribution features a higher probability for lower
values, indicating that edges are not likely to exist in graphs. Meanwhile,
the spike-and-slab prior is placed on the group similarity matrix $\Theta$ to
allow discrimination between zero similarity and positive similarity. Finally,
the precision matrices are estimated conditional on the graph structures
$G^{(k)}=(V,E^{(k)}),\;k=1,\ldots,K$ using the G-Wishart prior.

##  5 Estimating Differential Graphical Models

In contrast to joint estimation, several applications in biomedical research,
such as analyzing the gene expression differences in normal cells and cancer
cells or differences between the test group and control group, consider the
case where $K=2$. Different from the methods in Section 3–4, we will be
focusing on finding the “differences” rather than finding the “similarities”.
In the high dimensional setting, we assume that the difference of two graphs
is sparse. Although the differences between two graphs can be naively
estimated by using a joint estimation method first and then finding the
difference, procedures that directly estimate the difference are statistically
more efficient (Shojaie, 2021).

In this section, we briefly overview two approaches that estimate graph
differences, the direct estimation method and the regularization based
approach. For a detailed introduction, see Shojaie (2021) for a recent review.

###  5.1 Direct Estimation

The direct approach estimates the difference
${\bm{\Delta}}={\bm{\Omega}}^{(1)}-{\bm{\Omega}}^{(2)}$ without explicitly
estimating individual precision matrices ${\bm{\Omega}}^{(1)}$ and
${\bm{\Omega}}^{(2)}$. This approach potentially fits a broader class of
precision matrices as the individual precision matrices, ${\bm{\Omega}}^{(1)}$
and ${\bm{\Omega}}^{(2)}$, need not be sparse, but only the difference
${\bm{\Delta}}$ is assumed sparse. In addition, jointly estimating
${\bm{\Omega}}^{(1)}$ and ${\bm{\Omega}}^{(2)}$ can be challenging when the
sparse assumption is violated. Zhao et al. (2014) directly estimated the
difference ${\bm{\Delta}}$ by solving a constrained minimization problem,
noting that, by definition, we have
${\bm{\Sigma}}^{(1)}{\bm{\Delta}}{\bm{\Sigma}}^{(2)}-{\bm{\Sigma}}^{(1)}+{\bm{\Sigma}}^{(2)}={\bf
0}$. Consequently, estimating the differential graph
${\bm{\Delta}}={\bm{\Omega}}^{(1)}-{\bm{\Omega}}^{(2)}$ can be achieved by
minimizing the following objective:

$\displaystyle\widehat{{\bm{\Delta}}}=\mathop{\mathrm{argmin}}\|{\bm{\Delta}}\|_{1},\quad\text{subject
to
}\left|\widehat{\bm{\Sigma}}^{(1)}{\bm{\Delta}}\widehat{\bm{\Sigma}}^{(2)}-\widehat{\bm{\Sigma}}^{(1)}+\widehat{\bm{\Sigma}}^{(2)}\right|_{\infty}\leq\lambda_{1},$

which is an extension of the CLIME (Cai et al., 2011) method. Xu and Gu
(2016); Yuan et al. (2017) utilized the symmetry property
${\bm{\Sigma}}^{(1)}{\bm{\Delta}}{\bm{\Sigma}}^{(2)}={\bm{\Sigma}}^{(2)}{\bm{\Delta}}{\bm{\Sigma}}^{(1)}$
and hence
$2^{-1}({\bm{\Sigma}}^{(1)}{\bm{\Delta}}{\bm{\Sigma}}^{(2)}+{\bm{\Sigma}}^{(2)}{\bm{\Delta}}{\bm{\Sigma}}^{(1)})-{\bm{\Sigma}}^{(1)}+{\bm{\Sigma}}^{(2)}={\bf
0}$. They defined the objective function as

$\displaystyle\widehat{\bm{\Delta}}=\mathop{\mathrm{argmin}}\widehat{\ell}({\bm{\Delta}})+\lambda_{1}\|{\bm{\Delta}}\|_{1};$
$\displaystyle\widehat{\ell}({\bm{\Delta}})=\frac{1}{2}\mathop{\mathrm{tr}}({\bm{\Delta}}\widehat{\bm{\Sigma}}^{(1)}{\bm{\Delta}}\widehat{\bm{\Sigma}}^{(2)})-\mathop{\mathrm{tr}}\left\\{{\bm{\Delta}}(\widehat{\bm{\Sigma}}^{(1)}-\widehat{\bm{\Sigma}}^{(2)})\right\\},$

where the Hessian of the objective with respect to ${\bm{\Delta}}$ is
$(\widehat{\bm{\Sigma}}^{(1)}\otimes\widehat{\bm{\Sigma}}^{(2)}+\widehat{\bm{\Sigma}}^{(2)}\otimes\widehat{\bm{\Sigma}}^{(1)})/2$,
which is positive semi-definite. Therefore,
$\widehat{\ell}({\bm{\Delta}})+\lambda_{1}\|{\bm{\Delta}}\|_{1}$ is a convex
function with respect to ${\bm{\Delta}}$, hence a unique minimizer exists.
Other recent work extends the direct estimation approach to more expressive
structured differential graphs. Na et al. (2021) constructed a latent
structure estimator where the underlying difference can be formulated as a sum
of low-rank matrix and sparse matrix – a framework first discussed by
(Chandrasekaran et al., 2012). Zhao et al. (2019) extended the direct
estimation approach to estimating the differential graph of functional data.

###  5.2 Regularization based approach

The node-based learning framework (Mohan et al., 2012, 2014) assumes that most
parts of the graph are shared, and the difference is generated by a node
perturbation. When a node is perturbed, the edges connecting this node to
others change across $K$ groups. In addition to maximizing the degree of the
overlapping structure between groups, the task is to detect perturbed nodes.
An intuitive way to look for the perturbed node is to look at the difference
of two graphs ${\bm{\Omega}}^{(1)}-{\bm{\Omega}}^{(2)}$. When the $j$-th node
is being perturbed, the corresponding $j$-th row and $j$-th column of
${\bm{\Omega}}^{(1)}-{\bm{\Omega}}^{(2)}$ will have non-zeros, constructing a
unique symmetric row-column group. Given that there are several nodes being
perturbed, ${\bm{\Omega}}^{(1)}-{\bm{\Omega}}^{(2)}$ will be the union of the
row-column groups, each stemming from a perturbed node. Using this concept,
the Row-Column Overlap Norm (RCON) (Mohan et al., 2012, 2014) is designed to
encourage sparsity in the union of the row-column groups:

$\displaystyle{\bm{P}}_{\text{RCON}}({\bm{\Omega}})=\lambda_{1}\sum_{k=1}^{2}\sum_{i,j}|\omega_{i,j}^{(k)}|+\lambda_{2}G_{q}({\bm{\Omega}}^{(1)}-{\bm{\Omega}}^{(2)});$
$\displaystyle G_{q}({\bf A})=\min_{{\bf V}:{\bf A}={\bf V}+{\bf
V}^{\top}}f({\bf V});$ $\displaystyle f({\bf V})=\sum_{j=1}^{p}\|{\bf
v}_{j}\|_{q},$

where ${\bf v}_{j}$ is the $j$-th column of $\bf V$. It is easy to see that
when $q=1$, the RCON penalty is equivalent to the FGL penalty in (5). This
penalty function simultaneously imposes sparse structure on both the
individuals , ${\bm{\Omega}}^{(1)}$ and ${\bm{\Omega}}^{(2)}$, and the
difference ${\bm{\Omega}}^{(1)}-{\bm{\Omega}}^{(2)}$. As mentioned earlier,
this method may not work well under the setting that ${\bm{\Omega}}^{(1)}$ and
${\bm{\Omega}}^{(2)}$ are not sparse.

To infer the relative differences between two graphs in a Bayesian
formulation, it is intuitive to place a prior on the the differences of two
graphs $\theta_{i,j}=E^{(1)}_{i,j}-E^{(1)}_{i,j}$, for every $i<j$. Since the
difference $\theta_{i,j}$ is binary, either $0$ (no difference) or $1$
(difference), an example (Mitra et al., 2016) is to place a Bernoulli
distribution $Ber(\pi)$ on $\theta_{i,j}$, $i<j$ where $\pi$ follows Beta
distribution, specifying the tendency of being different on two graphs.

##  6 Graphical Models for Time Series Data

Time-varying graphical models (Zhou et al., 2010a; Kolar et al., 2010; Zhu and
Koyejo, 2018) can be seen as extensions of joint graphical models with groups
organized along the time index. The samples are assumed to be generated as
${\bf
x}_{i}^{(t)}\sim\mathcal{N}({\bm{\mu}}^{(t)},({\bm{\Omega}}^{(t)})^{-1})$,
$i=1,\ldots,n_{t}$, where $t=1,\ldots,T$ is the time index. Under such a
model, the estimation of time-varying precision matrices and corresponding
dynamic networks is challenging as data scarcity is a serious issue: in many
problems, we only observe a single sample at each time point. Therefore, to
make the estimation possible, structural assumptions are imposed on how the
underlying precision matrices and dynamic networks change over time. Such
assumptions control the model complexity and allow for the development of
estimation procedures. Examples of structural assumptions on temporal dynamics
include piece-wise constant and smoothly changing precision matrices, as well
as combinations of both. Piece-wise constant structure captures discrete
temporal evolution from one stage to another. For example, the gene regulatory
network in a fruit-fly can undergo structural changes as the fruit-fly
develops from an embryo to an adult state. Smooth temporal structure can be
used to model the dynamic functional connectivity of brain networks that
exhibit smooth temporal evolution from one brain state to another (Shine et
al., 2016). The temporal dynamics of crime rates are often modeled as a
combination of smooth dynamics and sudden jumps, where the jumps capture
sudden serious crime events. In this section, we will discuss how to apply the
FGL penalty and its variants to build piece-wise constant structure. We also
introduce a joint estimation framework of multiple autoregressive models to
model smooth temporal data.

###  6.1 Regularized Estimation

The FGL penalty has been widely used in time-varying graphical models to model
piece-wise constant dynamics (Kolar et al., 2010; Kolar and Xing, 2012; Monti
et al., 2014; Hallac et al., 2017). For instance, Smooth Incremental Graphical
Lasso Estimation (SINGLE) (Monti et al., 2014) applies the FGL framework to
enforce similarity between consecutive precision matrices:

$\displaystyle{\bf
P}_{SINGLE}({\bm{\Omega}})=\lambda_{1}\sum_{t=1}^{T}\sum_{i\neq
j}|\omega_{i,j}^{(t)}|+\lambda_{2}\sum_{t=2}^{T}\sum_{i\neq
j}|{\omega}^{(t)}_{i,j}-{\omega}_{i,j}^{(t-1)}|.$ (13)

The first term encourages the sparsity of each graph and the second term
regularizes the “jumps” across time. On the other hand, Group-Fused Graphical
Lasso (GFGL) (Gibberd and Nelson, 2017) introduces the Frobenius norm as an
alternative to encourage neighbouring similarity:

$\displaystyle{\bf
P}_{GFGL}({\bm{\Omega}})=\lambda_{1}\sum_{t=1}^{T}\sum_{i\neq
j}|\omega_{i,j}^{(t)}|+\lambda_{2}\sum_{t=2}^{T}\|{\bm{\Omega}}^{(t)}_{-ii}-{\bf\Omega}_{-ii}^{(t-1)}\|_{F},$

where ${\bm{\Omega}}^{(t)}_{-ii}$ denotes the precision matrix
${\bm{\Omega}}^{(t)}$ with the diagonal part removed. One may wonder what the
differences in structure assumptions are between the Frobenius norm and the
$\ell_{1}$-norm in (13). The $\ell_{1}$-norm regularizes individual changes,
while the Frobenius norm assumes global changes, implying that several edges
within a graph will change simultaneously.

While methods introduced in the last paragraph encourage similarity of two
neighboring graphs, the graph that is one-step ahead and that of one-step
behind, another idea is to enforce the similarities within multiple steps
ahead and behind. This can be done by creating a moving window index set (Yang
and Peng, 2020) $\mathcal{N}_{w}(t)$ for each time point $t=1,\ldots,T$.
Consider a window of length $2w$. At every time point $t$, we look at data
$w$-steps ahead and $w$-steps behind and hence the index set is
$\mathcal{N}_{w}(t)=\\{i=1,\ldots,T:|t-i|\leq 2w\\}$. Note that the index set
$\mathcal{N}_{w}(t)$ also includes $t$ itself. Then, we apply the GGL penalty
to the components in the index set.

###  6.2 Kernel Smoothing Graphical Models

Another way to construct smoothly varying graphs is by using an autoregressive
structure. This model assumes that each data point is a linear combination of
previous data points with additional independent noise. Consider the lag-1
case, where ${\bf x}^{(t)}$ is a linear transform of ${\bf x}^{(t-1)}$ with
independent noise ${\bm{\varepsilon}}^{(t)}\sim\mathcal{N}({\bf 0},{\bf
G}^{(t)})$:

$\displaystyle{\bf x}^{(t)}={\bf A}{\bf
x}^{(t-1)}+{\bm{\varepsilon}}^{(t)},\quad t=1,\ldots,T,$

where ${\bf A}\in\mathbb{R}^{p\times p}$ is the transition matrix.
Consequently, the covariance matrix is smoothly varying along $t$ if ${\bf
G}^{(t)}$ is a smooth function of $t$:

$\displaystyle{\bm{\Sigma}}^{(t)}={\bf A}{\bm{\Sigma}}^{(t-1)}{\bf A}^{T}+{\bf
G}^{(t)},\quad t=1,\ldots,T.$ (14)

Motivated by this structure, Zhou et al. (2010b) proposed a kernel based
method to estimate a smooth time-varying covariance structure. First, a
weighted sum of the sample covariance matrices
$\widehat{\bm{\Sigma}}^{(1)},\ldots,\widehat{\bm{\Sigma}}^{(T)}$ is computed
as

$\widehat{\bm{S}}^{(t)}=\frac{\sum_{s=1}^{T}w(s,t)\widehat{\bm{\Sigma}}^{(s)}}{\sum_{s^{\prime}=1}^{T}w(s^{\prime},t)},$

where the weights are constructed by a symmetric nonnegative kernel function
$K(|s-t|/h)$. This ensures that the estimated covariance is smoothly varying
over time. Subsequently, the precision matrix is estimated using the following
objective:

$\widehat{\bm{\Omega}}^{(t)}=\mathop{\mathrm{argmax}}\quad
n\left[\frac{1}{2}\log\\{\det({{\bm{\Omega}}^{(t)}})\\}-\frac{1}{2}\mathop{\mathrm{tr}}(\widehat{\bf
S}^{(t)}{\bm{\Omega}}^{(t)})\right]+\lambda_{n}\sum_{i\neq
j}|\omega_{i,j}^{(t)}|,\quad t=1,\ldots,T.$

The kernel smoothing method can also be extended to model two-way continuous
changes. For instance, the ages of subjects from the fMRI dataset vary across
an interval, and one can parametrize the transition matrices as ${\bm{A}}(u)$
with $u$ taking values in a closed subset of the real line. This model is
smooth in two aspects: across the temporal domain and labels (groups). Hence,
we have the following autoregressive model:

$\displaystyle{\bm{x}}_{i,t}={\bm{A}}(u_{i}){\bm{x}}_{i,t-1}+{\bm{\varepsilon}}_{i,t},\quad
i=1,\ldots,n,\;t=2,\ldots,T.$
$\displaystyle{\bm{\Sigma}}(u)={\bm{A}}(u){\bm{\Sigma}}(u){\bm{A}}(u)^{\top}+\sigma^{2}{\bm{I}}.$

The Kernel-Smoothing Estimator (KSE) (Qiu et al., 2016) is a kernel based
estimator for the covariance matrix and then uses the CLIME (Cai et al., 2011)
method introduced in Section 3.1 to recover precision matrices. We consider a
set of $n$ data ${\bm{Y}=\\{{\bm{y}}_{1},\ldots,{\bm{y}}_{n}\\}}$, where
${\bm{y}}_{i}=\\{{\bm{y}}_{i,1},\ldots,{\bm{y}}_{i,T}\\}\in\mathbb{R}^{p\times
T}$ and with label $u_{i}\in[0,1]$. The estimated covariance model of the
label $u_{0}\in[0,1]$ is formulated as the following:

$\displaystyle\widehat{\bm{S}}(u_{0})=\sum_{i=1}^{n}w_{i}(u_{0},h)\widehat{\bm{\Sigma}}_{i};$
$\displaystyle
w_{i}(u_{0},h):=\frac{c(u_{0})}{nh}K\left(\frac{u_{i}-u_{0}}{h}\right);$
$\displaystyle
c(u_{0})=\left\\{\begin{array}[]{cc}2{\bm{I}},\;u_{0}\in\\{0,1\\},\\\
1{\bm{I}},\;u_{0}\in(0,1),\end{array}\right.;$
$\displaystyle\widehat{\bm{\Sigma}}_{i}=\frac{1}{T}\sum_{t=1}^{T}{\bm{y}}_{i,t}{\bm{y}}_{i,t}^{\top},$

where $w_{i}$ is the kernel-based weight with predefined scale $h$, $K(\cdot)$
is the kernel, $c(u_{0})$ determines the boundary value, and
$\widehat{\bm{\Sigma}}_{i}$ is the sampled covariance of the time-series data.
After obtaining $\widehat{\bm{S}}(u_{0})$, the precision matrix
$\widehat{\bm{\Omega}}(u_{0})$ is obtained using CLIME (Cai et al., 2011) in
(4):

$\displaystyle\widehat{{\bm{\Omega}}}(u_{0})=\mathop{\mathrm{argmin}}\|{\bm{\Omega}}(u_{0})\|_{1},$
$\displaystyle\text{subject to
}|\widehat{{\bm{S}}}(u_{0}){\bm{\Omega}}(u_{0})-{\bm{I}}|_{\infty}\leq\xi.$

Under this framework, the kernel trick is used to capture the assumption that
the covariance matrices are smoothly varying across labels. In addition, the
Euclidean distance of two labels reflects the similarity of the two groups,
capturing the dependence structure. The kernel-based method can be applied to
general joint estimation, where the sampled covariance of time-series data is
replaced by the sample covariance of data with the same labels.

##  7 Experiments

We demonstrate that methods that jointly estimate graphs perform better than
methods that estimate them separately, both in terms of edge recovery and
graph estimation. The quality of edge recovery is evaluated via the Receiver
Operating Characteristic (ROC) curve. The quality of graph estimation is
evaluated via the entropy loss – which requires the knowledge of the ground
truth precision matrix to compute. In the absence of ground truth parameters,
we evaluate the quality of the graph estimation by splitting the data into two
batches. We use one batch for estimation and the other batch for evaluation by
computing the out-of-sample log-likelihood.

To evaluate the performance of joint estimation compared to separate
estimation, we examine three regularization methods that capture the key
properties of each family of models: the group graphical lasso, which is a
canonical example of JSEM (Ma and Michailidis, 2016), LASICH (Saegusa and
Shojaie, 2016), JAGL (Shan and Kim, 2018), JMMLE (Majumdar and Michailidis,
2018); the fused graphical lasso, which is a canonical example of GFGL
(Gibberd and Nelson, 2017), TFRE (Bilgrau et al., 2020), GEN-ISTA (Price et
al., 2021), and the graphical lasso method – which does not incorporate joint
information. The first two penalty functions encourage similarities between
groups, while the last one only encourages sparsity within groups. We evaluate
the performance of these methods under different data generating processes.

###  7.1 Experiment setup

The graph ${G}^{(k)}$ of each group is obtained as the union of a shared
component ${G}_{s}=(V,E_{s})$ and an individual component
${G}_{i}^{(k)}=(V,E_{i}^{(k)})$, that is, $G^{(k)}=(V,E_{s}\cup E_{i}^{(k)})$.
The shared component is a scale-free graph, whose degree follows a power
distribution. An individual component is generated as an Erdős Rényi random
graph, where $\rho$ specifies the ratio of edges being connected (in most
settings, we set $\rho=0.1$). We consider two settings. Data Generation
Process $1$ (DGP1) has a partially shared support of the inverse covariance,
with the entries on the shared support also sharing values. Data Generation
Process $2$ (DGP2) only has a shared support, but the corresponding inverse
covariance matrix differs across groups.

###  7.2 Evaluation of Edge Selection

We summarize the performance of each approach using the Receiver Operating
Characteristic (ROC) curve, where the x-axis denotes the false positive rate
(FPR) of the estimated edges and the y-axis denotes the true positive rate
(TPR) of the estimated edges. The ideal estimator is an estimator that has the
true positive rate of $1$ and false positive rate of $0$, and hence (FPR,
TPR)$=(0,1)$ is at the upper left corner of the ROC plot. This implies that
the estimator selects all true edges and does not include any false edges of
the graph. If the estimator selects all edges of a graph, then both the TPR
and FPR are $1$. Then, under this scenario, (FPR, TPR)$=(1,1)$ is at the upper
right corner of the ROC plot. Similarly, if the model does not select any
edge, then both the FPR and TPR are $0$, resulting (FPR, TPR)=$(0,0)$ located
on the bottom left of the ROC plot. When we vary $\lambda_{1}$, which
determines the strength of the sparsity of individual component
$\omega_{i,j}^{(k)}$ for $i\neq j$ and $k=1,\ldots,K$, as shown in (5)–(6), it
controls the FPR and TPR of the underlying estimator. Plotting (FPR, TPR)
under different choices of $\lambda_{1}$ results in the ROC curve.

Since (FPR, TPR)$=(0,1)$ is rarely achievable, and the extreme cases (FPR,
TPR)=$(1,1),(0,0)$ do not provide useful results, one may restrict the value
for either rate and optimize for the other rate. That is, given a fixed TPR
tolerance, we would like the FPR to be as small as possible. On the other
hand, given a fixed FPR tolerance, we want the TPR rate to be as large as
possible. In principle, a good model forms a ROC curve that is curved to the
upper left corner, which is toward (FPR, TPR)$=(0,1)$, and hence it has a
large Area Under the ROC Curve (AUC). Consequently, we can use the AUC to
compare different models.

In our experiments, each ROC curve is computed by tuning the penalty parameter
$\lambda_{1}$ over the range $[0.01,0.3]$ and averaging results over $20$
independent simulation runs. Recall $\lambda_{2}$ controls the sparsity
pattern of the shared components, as demonstrated in (5)–(6). We tune the
penalty parameter $\lambda_{2}$ over 8 different values,
$\lambda_{2}\in\\{0.025,0.05,0.075,0.1,0.125,0.15,0.175,0.2\\}$, and pick the
value that has the largest AUC. We set $p=100$, $n=50$, $K=5$, and $\rho=0.1$.
The false-positive rate is defined as:

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\frac{\sum_{i,j}I(\omega_{i,j}^{(k)}=0,\widehat{\omega}_{i,j}^{(k)}\neq
0)}{\sum_{i,j}I(\omega_{i,j}^{(k))}=0)},$

and the false negative rate is defined as:

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\frac{\sum_{i,j}I(\omega_{i,j}^{(k)}\neq
0,\widehat{\omega}_{i,j}^{(k)}=0)}{\sum_{i,j}I(\omega_{i,j}^{(k))}\neq 0)},$

where $I(\cdot)$ is the indicator function, $\omega_{i,j}^{(k)}$ is the
$i,j$-th component of the ground truth inverse covariance matrix of group $k$,
$\widehat{\omega}_{i,j}^{(k)}$ is the $i,j$-th component of the estimated
inverse covariance matrix of group $k$. In practice, we use the truncation
method to replace the zero-threshold identity function $I(\cdot)$. That is,
adding a threshold to the identity function:
$I(|\omega^{(k)}_{i,j}|\geq\varepsilon)$, where $\varepsilon$ is a predefined
threshold to replace $I(\omega^{(k)}_{i,j}\neq 0)$, where $\varepsilon$ is
selected to be $10^{-3}$ in our simulations. The same mechanism works for the
case when $I(\omega^{(k)}_{i,j}=0)$. This approach leads to less noisy
recovery.

Figure 4:  DGP1 features shared structure with the same inverse covariance
values. In this case, the fused graphical lasso is expected to achieve the
best performance. The experiment results in the left figure demonstrates our
expectation. The center and right figures show the ROC curve of fused
graphical lasso and group graphical lasso for different $\lambda_{2}$ values.
While changing the $\lambda_{2}$ in the center figure does not show obvious
differences, the larger $\lambda_{2}$ in right figure shows better
performance.

Figure 5:  DGP2 features a shared structure with different inverse covariance
values. In this case, the group graphical lasso is expected to achieve the
best performance. The experiment results in left figure reflect our
expectations. The center figure and right figure shows the ROC curve of fused
graphical lasso and group graphical lasso under different $\lambda_{2}$
values. In the center figure, it is clear to see that smaller $\lambda_{2}$ on
the fused graphical lasso penalty function shows better performance while in
the right figure, larger $\lambda_{2}$ on the group lasso penalty function
achieves better performance.

For DGP1, we hypothesize that the fused graphical lasso would have the best
performance in small sample size settings. This is because the fused graphical
lasso encourages the similarity of both structure and values across groups.
For DGP2, we expect that the group graphical lasso to have the best
performance because it only encourages shared structure but allows flexibility
for each group to retain different inverse covariance values. These two
experiments demonstrate our hypothesis: if the joint models share the same
inverse covariance values, the fused graphical lasso function has the largest
AUC. On the other hand, if the joint models have a shared structure, but the
structure has different inverse covariance values, the group graphical lasso
has the largest AUC. Alternatively, we could select proper regularization
functions based on prior knowledge about the form of the shared structures.

###  7.3 Model Selection

To properly select the hyperparameters $\lambda_{1},\lambda_{2}$, we employ
the Bayesian Information Criterion (BIC) (Schwarz, 1978) as follows:

$\displaystyle\lambda_{1}^{*},\lambda_{2}^{*}=\mathop{\mathrm{argmin}}_{\lambda_{1},\lambda_{2}}\sum_{k=1}^{K}tr(\widehat{\bm{\Sigma}}^{(k)}\widehat{\bm{\Omega}}^{(k)})-\log\\{\det(\widehat{\bm{\Omega}}^{(k)})\\}+E_{k}\log(n_{k}),$
(15)

where $\widehat{E}^{(k)}$ is the sampled covariance of group $k$,
$\widehat{\Omega}^{(k)}$ is the estimated precision matrix of group $k$,
$n_{k}$ is the number of samples of group $k$, and $E_{k}$ is the number of
nonzero elements of the precision matrix. The BIC method is a trade-off
between the goodness of fit and the model complexity. The first term in (15)
maximizes the data likelihood, while the second term penalizes the model
complexity by discouraging the number of connected edges. We perform a grid
search and find $\lambda_{1}^{*},\lambda_{2}^{*}$ that achieves the minimum
BIC value. Figure 6 illustrates an example of the fused graphical lasso model
selection on graphs generated by DGP2, where the setting is $p=100$, $n=50$,
$\rho=0.1$.

Figure 6: The left figure is the level sets of the BIC values with different
$\lambda_{1}$ and $\lambda_{2}$. The right figure is the AIC (Akaike, 1974)
values with different $\lambda_{1}$ and $\lambda_{2}$. The level lines are
nearly vertical which indicate that tuning $\lambda_{1}$ has more influence on
the performance than tuning $\lambda_{2}$ in our test settings.

###  7.4 Entropy Loss v.s. $n\mathbin{/}p$

High-dimensional data is characterized by settings where $n$ is much greater
than $p$, leading to poor estimation results. This condition can be alleviated
using regularizers. In this experiment, we evaluate the estimation performance
with different regularizers when the ratio ${n}/{p}$ is smaller than $1$,
meaning that we have more parameters than samples to estimate. We use the
entropy loss (Guo et al., 2011) as an evaluation metric defined as:

$\displaystyle
EL=\frac{1}{K}\sum_{k=1}^{K}\mathop{\mathrm{tr}}\\{(\bm{\Omega}^{(k)})^{-1}\widehat{\bm{\Omega}}^{(k)}\\}-\log|\det\\{\bm{\Omega}^{(k)})^{-1}\widehat{\bm{\Omega}}^{(k)}\\}|-p,$

where $k$ is the group index, ${\bm{\Omega}}^{(k)}$ is the ground truth
precision matrix, and $\widehat{\bm{\Omega}}^{(k)}$ is the recovered precision
matrix. Noting that the entropy loss here is equivalent to the Kullback-
Leibler divergence of two centered multivariate Gaussian distributions, which
is a widely used metric to evaluate the distance of two distributions. From
the definition of entropy loss, it is clear that when
$\widehat{\bf\Omega}^{(k)}={\bf\Omega}^{(k)}$ and the inverse of
${\bf\Omega}^{(k)}$ exist, the entropy loss reduces to $0$.

In the following experiment, we test $p=\\{20,50,80,100\\}$ and change the
number of samples $n$. The ground truth graphs are generated using DGP1 and
each curve is the average of $5$ replications. The model parameters are
selected using BIC.

Figure 7: The top left figure shows when $p=20$, top right figure shows when
$p=50$, bottom left figure shows when $p=80$ and bottom right figure shows
when $p=100$. Among the four cases, when ${n}/{p}$ is smaller than $0.2$ the
entropy curves have sharp slopes, and yet the fused graphical lasso method has
the smallest slopes. This implies that under small sample size settings, joint
estimation achieves better recovery. When the sample size $n$ increases, the
curves gradually become flat, and all three methods achieve comparable
recovery performance.

###  7.5 Entropy Loss v.s. Graph Size

Now consider the case when the sample size $n=10$ is fixed, but the size of
the graph is growing. We evaluate the performance of three penalty functions
under this condition. In this experiment, we use DGP1 to construct the ground
truth and the models are selected by BIC.

Figure 8: With the sample size fixed, the entropy loss increases as the graph
size increases. The plot indicates that under this condition, the fused
graphical lasso penalty function has a slower increase of the loss.

###  7.6 Out-of-Sample Log-Likelihood

In most applications, we do not know the ground truth graph. Therefore, we
cannot rely on the entropy loss or similar approach to evaluate the recovery.
In this case, we split our data into training and test sets and use the test
sets to evaluate the log-likelihood. We construct the ground truth using graph
generation process 2 with $K=5,p=100,n=50$ and the models are selected using
BIC.

| fused graphical lasso | group graphical lasso | graphical lasso
---|---|---|---
$\lambda_{1}$ | $0.88$ | $0.13$ | 0.92
$\lambda_{2}$ | $0.01$ | $1.0$ |
out-of-sample log-likelihood | $-498.94$ | $-489.99$ | $-573.34$
Table 2: Out-of-sample log-likelihood

##  8 Open Problems

Existing and emerging biological data and applications will require novel
approaches to joint graphical models. We discuss some of these emerging
applications briefly. Joint estimation of functional connectivity networks
across multiple subjects allow scaling of the effective sample size and
computation of estimates that are more robust to outliers. The joint
estimators of brain connectivity networks could be applied to task-based fMRI
scans to study group dynamic functional connectivity patterns (Andersen et
al., 2018; Calhoun et al., 2014; Gonzalez-Castillo and Bandettini, 2018).
While this manuscript is focused on joint estimation with the same set of
nodes, one potential direction is to extend it to multiple sources, i.e.,
multimodal data. Recent technologies (Huster et al., 2012; Abreu et al., 2018)
have demonstrated the availability of conducting concurrent measurements of
EEG and fMRI signals, allowing the estimation of multiple sources possible in
the future. While EEG has higher temporal resolution and fMRI features higher
spatial resolution, we believe that joint estimation with multiple sources
could compensate for the limits of the measurement techniques and provide
better estimation results. Some recent work (Lock et al., 2013; Li and Li,
2021) has developed methodologies to integrate data from different modalities,
however, joint estimation of graphical models from multimodal data is still an
open problem.

Approaches for the estimation of the joint graphical models presented in this
survey largely rely on penalized estimation, where the penalty biases the
estimates towards the assumed structure. Quantifying statistical uncertainty
about the model parameters, that is, performing hypothesis tests and
constructing confidence intervals, is challenging when penalized estimators
are used due to the induced bias and model selection that is implicitly
performed. There has been recent work on statistical inference for low-
dimensional parameters in graphical models (Ren et al., 2015; Janková and van
de Geer, 2015, 2017; Barber and Kolar, 2018; Wang and Kolar, 2016; Yu et al.,
2016, 2020) based on the $\ell_{1}$-penalized estimator in the first stage.
However, these approaches were developed only in the setting where parameters
of one graph are being inferred. In contrast, work on statistical inference
for joint graphical models is much more sparse. Xia et al. (2015); Belilovsky
et al. (2016); Liu (2017); Kim et al.  developed techniques for statistical
inference in differential graphical models, while Wang and Kolar (2014); Lu et
al. (2018); Wang et al. (2020a) focused on graphical models for time series
data. Developing corresponding inferential techniques for estimators obtained
using coarse-grained and fine-grained penalties is an interesting area open
for future research.

##  9 Conclusion

This manuscript has introduced joint Gaussian graphical model estimation
methods for joint data with shared structure across multiple groups. In
particular, we consider several examples of extending classical statistical
inference methods to joint estimation settings, including the MLE based
estimator, neighborhood regression, and the CLIME estimator. We discuss
several methods that exploit coarse-grained structures using a global
regularization method that encourages shared coarse-grained structure across
all groups. In contrast, the fine-grained structural regularization methods
further partition the groups into subgroups per node, encouraging local shared
regularity. With two groups, differential graphs are often a highly effective
approach. The presented simulation experiments demonstrate that joint
estimation can often recover the graphs more effectively than separate
estimation for high dimensional data. Furthermore, we verify the
characteristics of fused graphical lasso and group graphical regularizers for
two data generation processes. When the data generation process has a shared
inverse covariance structure with different values, group graphical lasso
methods may be more effective. On the other hand, when the data generation
process features a shared inverse covariance structure with the same values,
the fused graphical lasso penalty function may perform better.

## Funding Information

O.Koyejo acknowledges partial funding from a C3.ai Digital Transformation
Institute Award, a Jump Arches Award, and an Strategic Research Initiatives
award from the University of Illinois at Urbana-Champaign. K.Tsai acknowledges
funding from National Science Foundation Graduate Research Fellowships
Program. Other authors have no relevant financial or non-financial interests
to disclose. This work was also funded in part by the following grants: NSF
III 2046795 and IIS 1909577, along with computational resources donated by
Microsoft Azure.

## Further Reading

Recent developments of joint statistical inference are primarily focused on
Gaussian graphical models. Other types of graphical models, including discrete
graphical models (Drton and Richardson, 2008; Drton, 2009),
semiparametric/nonparametric graphical models (Liu et al., 2012; Sun et al.,
2015), and latent graphical models (Chandrasekaran et al., 2012), have been
well studied for single graph estimation. While such models have broad
applications, joint estimation in these models is less studied. A number of
open source provided tools are available for joint estimation of multiple
graphical models – Danaher et al. (2014); Li et al. (2019); Bilgrau et al.
(2020) have released an R-package for regularization-based methods; Lin et al.
(2017) provide the MATLAB source code for the joint temporal and spatial
estimation; Mohan et al. (2014) provides code for estimating differential
graphs.

## References

  * Abreu et al. (2018) R. Abreu, A. Leal, and P. Figueiredo.  Eeg-informed fmri: a review of data analysis methods.  _Frontiers in human neuroscience_ , 12:29, 2018.
  * Akaike (1974) H. Akaike.  A new look at the statistical model identification.  _IEEE transactions on automatic control_ , 19(6):716–723, 1974.
  * Andersen et al. (2018) M. R. Andersen, O. Winther, L. K. Hansen, R. Poldrack, and O. Koyejo.  Bayesian structure learning for dynamic brain connectivity.  In _21st International Conference on Artificial Intelligence and Statistics, AISTATS 2018_ , pages 1436–1446, 2018.
  * Atay-Kayis and Massam (2005) A. Atay-Kayis and H. Massam.  A monte carlo method for computing the marginal likelihood in nondecomposable gaussian graphical models.  _Biometrika_ , 92(2):317–335, 2005.
  * Barber and Kolar (2018) R. F. Barber and M. Kolar.  Rocket: Robust confidence intervals via kendall’s tau for transelliptical graphical models.  _Ann. Statist._ , 46(6B):3422–3450, 2018.
  * Belilovsky et al. (2016) E. Belilovsky, G. Varoquaux, and M. B. Blaschko.  Testing for differences in gaussian graphical models: Applications to brain connectivity.  In D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett, editors, _Advances in Neural Information Processing Systems 29_ , pages 595–603. Curran Associates, Inc., 2016.
  * Bilgrau et al. (2020) A. E. Bilgrau, C. F. Peeters, P. S. Eriksen, M. Bøgsted, and W. N. van Wieringen.  Targeted fused ridge estimation of inverse covariance matrices from multiple high-dimensional data classes.  _Journal of Machine Learning Research_ , 21(26):1–52, 2020.
  * Bühlmann and Van De Geer (2011) P. Bühlmann and S. Van De Geer.  _Statistics for high-dimensional data: methods, theory and applications_.  Springer Science & Business Media, 2011.
  * Cai et al. (2011) T. Cai, W. Liu, and X. Luo.  A constrained $\ell_{1}$ minimization approach to sparse precision matrix estimation.  _Journal of the American Statistical Association_ , 106(494):594–607, 2011.
  * Calhoun et al. (2014) V. D. Calhoun, R. Miller, G. Pearlson, and T. Adalı.  The chronnectome: time-varying connectivity networks as the next frontier in fmri data discovery.  _Neuron_ , 84(2):262–274, 2014.
  * Chan et al. (2017) T. E. Chan, M. P. Stumpf, and A. C. Babtie.  Gene regulatory network inference from single-cell data using multivariate information measures.  _Cell systems_ , 5(3):251–267, 2017.
  * Chandrasekaran et al. (2012) V. Chandrasekaran, P. A. Parrilo, and A. S. Willsky.  Latent variable graphical model selection via convex optimization.  _The Annals of Statistics_ , 40(4):1935–1967, 2012.
  * Chun et al. (2015) H. Chun, X. Zhang, and H. Zhao.  Gene regulation network inference with joint sparse gaussian graphical models.  _Journal of Computational and Graphical Statistics_ , 24(4):954–974, 2015.
  * Colclough et al. (2018) G. L. Colclough, M. W. Woolrich, S. J. Harrison, P. A. R. López, P. A. Valdes-Sosa, and S. M. Smith.  Multi-subject hierarchical inverse covariance modelling improves estimation of functional brain networks.  _NeuroImage_ , 178:370–384, 2018.
  * Danaher et al. (2014) P. Danaher, P. Wang, and D. M. Witten.  The joint graphical lasso for inverse covariance estimation across multiple classes.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76(2):373–397, 2014.
  * Dobra et al. (2004) A. Dobra, C. Hans, B. Jones, J. R. Nevins, G. Yao, and M. West.  Sparse graphical models for exploring gene expression data.  _Journal of Multivariate Analysis_ , 90(1):196–212, 2004.
  * Dondelinger et al. (2013) F. Dondelinger, S. Lèbre, and D. Husmeier.  Non-homogeneous dynamic bayesian networks with bayesian regularization for inferring gene regulatory networks with gradually time-varying structure.  _Machine Learning_ , 90(2):191–230, 2013.
  * Drton (2009) M. Drton.  Discrete chain graph models.  _Bernoulli_ , 15(3):736–753, 2009.
  * Drton and Maathuis (2017) M. Drton and M. H. Maathuis.  Structure learning in graphical modeling.  _Annual Review of Statistics and Its Application_ , 4(1):365–393, 2017.
  * Drton and Richardson (2008) M. Drton and T. S. Richardson.  Binary models for marginal independence.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 70(2):287–309, 2008.
  * Felsenstein (1981) J. Felsenstein.  Evolutionary trees from dna sequences: a maximum likelihood approach.  _Journal of molecular evolution_ , 17(6):368–376, 1981.
  * Foti and Fox (2019) N. J. Foti and E. B. Fox.  Statistical model-based approaches for functional connectivity analysis of neuroimaging data.  _Current opinion in neurobiology_ , 55:48–54, 2019.
  * Friedman et al. (2000) N. Friedman, M. Linial, I. Nachman, and D. Pe’er.  Using bayesian networks to analyze expression data.  _Journal of computational biology_ , 7(3-4):601–620, 2000.
  * Gibberd and Nelson (2017) A. J. Gibberd and J. D. Nelson.  Regularized estimation of piecewise constant gaussian graphical models: The group-fused graphical lasso.  _Journal of Computational and Graphical Statistics_ , 26(3):623–634, 2017.
  * Gonzalez-Castillo and Bandettini (2018) J. Gonzalez-Castillo and P. A. Bandettini.  Task-based dynamic functional connectivity: Recent findings and open questions.  _Neuroimage_ , 180:526–533, 2018.
  * Greenewald et al. (2017) K. Greenewald, S. Park, S. Zhou, and A. Giessing.  Time-dependent spatially varying graphical models, with application to brain fmri data analysis.  In _Advances in Neural Information Processing Systems_ , pages 5832–5840, 2017.
  * Guo et al. (2011) J. Guo, E. Levina, G. Michailidis, and J. Zhu.  Joint estimation of multiple graphical models.  _Biometrika_ , 98(1):1–15, 2011.
  * Hallac et al. (2017) D. Hallac, Y. Park, S. Boyd, and J. Leskovec.  Network inference via the time-varying graphical lasso.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 205–213, 2017.
  * Hao et al. (2018) B. Hao, W. W. Sun, Y. Liu, and G. Cheng.  Simultaneous clustering and estimation of heterogeneous graphical models.  _Journal of Machine Learning Research_ , 18:217–1, 2018\.
  * Hsieh et al. (2013) C.-J. Hsieh, M. A. Sustik, I. S. Dhillon, P. K. Ravikumar, and R. Poldrack.  Big & quic: Sparse inverse covariance estimation for a million variables.  In _Advances in neural information processing systems_ , pages 3165–3173, 2013.
  * Huster et al. (2012) R. J. Huster, S. Debener, T. Eichele, and C. S. Herrmann.  Methods for simultaneous eeg-fmri: an introductory review.  _Journal of Neuroscience_ , 32(18):6053–6060, 2012.
  * Janková and van de Geer (2015) J. Janková and S. van de Geer.  Confidence intervals for high-dimensional inverse covariance estimation.  _Electron. J. Stat._ , 9(1):1205–1229, 2015.
  * Janková and van de Geer (2017) J. Janková and S. A. van de Geer.  Honest confidence regions and optimality in high-dimensional precision matrix estimation.  _TEST_ , 26(1):143–162, 2017.
  * (34) B. Kim, S. Liu, and M. Kolar.  Two-sample inference for high-dimensional markov networks.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , n/a(n/a), arXiv:https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/rssb.12446.
  * Kling et al. (2015) T. Kling, P. Johansson, J. Sanchez, V. D. Marinescu, R. Jörnsten, and S. Nelander.  Efficient exploration of pan-cancer networks by generalized covariance selection and interactive web content.  _Nucleic acids research_ , 43(15):e98--e98, 2015\.
  * Kolar and Xing (2012) M. Kolar and E. P. Xing.  Estimating networks with jumps.  _Electronic journal of statistics_ , 6:2069, 2012.
  * Kolar et al. (2010) M. Kolar, L. Song, A. Ahmed, E. P. Xing, et al.  Estimating time-varying networks.  _The Annals of Applied Statistics_ , 4(1):94--123, 2010.
  * Lauritzen (1996) S. L. Lauritzen.  _Graphical models_ , volume 17.  Clarendon Press, 1996.
  * Lee and Liu (2015) W. Lee and Y. Liu.  Joint estimation of multiple precision matrices with common structures.  _The Journal of Machine Learning Research_ , 16(1):1035--1062, 2015.
  * Lenkoski and Dobra (2011) A. Lenkoski and A. Dobra.  Computational aspects related to inference in gaussian graphical models with the g-wishart prior.  _Journal of Computational and Graphical Statistics_ , 20(1):140--157, 2011.
  * Li and Li (2021) Q. Li and L. Li.  Integrative factor regression and its inference for multimodal data analysis.  _Journal of the American Statistical Association_ , pages 1--15, 2021\.
  * Li et al. (2019) Z. Li, T. Mccormick, and S. Clark.  Bayesian joint spike-and-slab graphical lasso.  In _International Conference on Machine Learning_ , pages 3877--3885. PMLR, 2019.
  * Lin et al. (2017) Z. Lin, T. Wang, C. Yang, and H. Zhao.  On joint estimation of gaussian graphical models for spatial and temporal data.  _Biometrics_ , 73(3):769, 2017.
  * Liu et al. (2012) H. Liu, F. Han, M. Yuan, J. Lafferty, and L. Wasserman.  High-dimensional semiparametric gaussian copula graphical models.  _The Annals of Statistics_ , 40(4):2293--2326, 2012.
  * Liu (2017) W. Liu.  Structural similarity and difference testing on multiple sparse Gaussian graphical models.  _Ann. Statist._ , 45(6):2680--2707, 2017.
  * Lock et al. (2013) E. F. Lock, K. A. Hoadley, J. S. Marron, and A. B. Nobel.  Joint and individual variation explained (jive) for integrated analysis of multiple data types.  _The annals of applied statistics_ , 7(1):523, 2013.
  * Lu et al. (2018) J. Lu, M. Kolar, and H. Liu.  Post-regularization inference for time-varying nonparanormal graphical models.  _Journal of Machine Learning Research_ , 18(203):1--78, 2018.
  * Lukemire et al. (2020) J. Lukemire, S. Kundu, G. Pagnoni, and Y. Guo.  Bayesian joint modeling of multiple brain functional networks.  _Journal of the American Statistical Association_ , pages 1--13, 2020\.
  * Ma and Michailidis (2016) J. Ma and G. Michailidis.  Joint structural estimation of multiple graphical models.  _The Journal of Machine Learning Research_ , 17(1):5777--5824, 2016.
  * Maathuis et al. (2018) M. Maathuis, M. Drton, S. Lauritzen, and M. Wainwright.  _Handbook of graphical models_.  CRC Press, 2018.
  * Majumdar and Michailidis (2018) S. Majumdar and G. Michailidis.  Joint estimation and inference for data integration problems based on multiple multi-layered gaussian graphical models.  _arXiv preprint arXiv:1803.03348_ , 2018.
  * Manning et al. (2018) J. R. Manning, X. Zhu, T. L. Willke, R. Ranganath, K. Stachenfeld, U. Hasson, D. M. Blei, and K. A. Norman.  A probabilistic approach to discovering dynamic full-brain functional connectivity patterns.  _NeuroImage_ , 180:243--252, 2018.
  * Marlin and Murphy (2009) B. M. Marlin and K. P. Murphy.  Sparse gaussian graphical models with unknown block structure.  In _Proceedings of the 26th Annual International Conference on Machine Learning_ , pages 705--712, 2009.
  * Meinshausen et al. (2006) N. Meinshausen, P. Bühlmann, et al.  High-dimensional graphs and variable selection with the lasso.  _The annals of statistics_ , 34(3):1436--1462, 2006.
  * Mitchell and Beauchamp (1988) T. J. Mitchell and J. J. Beauchamp.  Bayesian variable selection in linear regression.  _Journal of the american statistical association_ , 83(404):1023--1032, 1988.
  * Mitra et al. (2016) R. Mitra, P. Müller, Y. Ji, et al.  Bayesian graphical models for differential pathways.  _Bayesian Analysis_ , 11(1):99--124, 2016.
  * Mohammadi et al. (2015) A. Mohammadi, E. C. Wit, et al.  Bayesian structure learning in sparse gaussian graphical models.  _Bayesian Analysis_ , 10(1):109--138, 2015.
  * Mohan et al. (2012) K. Mohan, M. Chung, S. Han, D. Witten, S.-I. Lee, and M. Fazel.  Structured learning of gaussian graphical models.  In _Advances in neural information processing systems_ , pages 620--628, 2012.
  * Mohan et al. (2014) K. Mohan, P. London, M. Fazel, D. Witten, and S.-I. Lee.  Node-based learning of multiple gaussian graphical models.  _The Journal of Machine Learning Research_ , 15(1):445--488, 2014.
  * Monti et al. (2014) R. P. Monti, P. Hellyer, D. Sharp, R. Leech, C. Anagnostopoulos, and G. Montana.  Estimating time-varying brain connectivity networks from functional mri time series.  _NeuroImage_ , 103:427--443, 2014.
  * Na et al. (2021) S. Na, M. Kolar, and O. Koyejo.  Estimating differential latent variable graphical models with applications to brain connectivity.  _Biometrika_ , 108(2):425--442, 2021.
  * Peterson et al. (2015) C. Peterson, F. C. Stingo, and M. Vannucci.  Bayesian inference of multiple gaussian graphical models.  _Journal of the American Statistical Association_ , 110(509):159--174, 2015.
  * Pierson et al. (2015) E. Pierson, D. Koller, A. Battle, S. Mostafavi, G. Consortium, et al.  Sharing and specificity of co-expression networks across 35 human tissues.  _PLoS computational biology_ , 11(5):e1004220, 2015.
  * Poldrack et al. (2011) R. A. Poldrack, J. A. Mumford, and T. E. Nichols.  _Handbook of functional MRI data analysis_.  Cambridge University Press, 2011.
  * Price et al. (2021) B. S. Price, A. J. Molstad, and B. Sherwood.  Estimating multiple precision matrices with cluster fusion regularization.  _Journal of Computational and Graphical Statistics_ , pages 1--12, 2021.
  * Qiu et al. (2016) H. Qiu, F. Han, H. Liu, and B. Caffo.  Joint estimation of multiple graphical models from high dimensional time series.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 78(2):487--504, 2016.
  * Ren et al. (2015) Z. Ren, T. Sun, C.-H. Zhang, and H. H. Zhou.  Asymptotic normality and optimalities in estimation of large Gaussian graphical models.  _Ann. Stat._ , 43(3):991--1026, 2015.
  * Roverato (2002) A. Roverato.  Hyper inverse wishart distribution for non-decomposable graphs and its application to bayesian inference for gaussian graphical models.  _Scandinavian Journal of Statistics_ , 29(3):391--411, 2002.
  * Saegusa and Shojaie (2016) T. Saegusa and A. Shojaie.  Joint estimation of precision matrices in heterogeneous populations.  _Electronic journal of statistics_ , 10(1):1341, 2016.
  * Schäfer and Strimmer (2005) J. Schäfer and K. Strimmer.  A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics.  _Statistical applications in genetics and molecular biology_ , 4(1), 2005.
  * Schwab et al. (2018) S. Schwab, R. Harbord, V. Zerbi, L. Elliott, S. Afyouni, J. Q. Smith, M. W. Woolrich, S. M. Smith, and T. E. Nichols.  Directed functional connectivity using dynamic graphical models.  _NeuroImage_ , 175:340--353, 2018.
  * Schwarz (1978) G. Schwarz.  Estimating the dimension of a model.  _The annals of statistics_ , pages 461--464, 1978.
  * Shan and Kim (2018) L. Shan and I. Kim.  Joint estimation of multiple gaussian graphical models across unbalanced classes.  _Computational Statistics & Data Analysis_, 121:89--103, 2018.
  * Shan et al. (2020) L. Shan, Z. Qiao, L. Cheng, and I. Kim.  Joint estimation of the two-level gaussian graphical models across multiple classes.  _Journal of Computational and Graphical Statistics_ , 29(3):562--579, 2020.
  * Shine et al. (2016) J. M. Shine, P. G. Bissett, P. T. Bell, O. Koyejo, J. H. Balsters, K. J. Gorgolewski, C. A. Moodie, and R. A. Poldrack.  The dynamics of functional brain networks: integrated network states during cognitive task performance.  _Neuron_ , 92(2):544--554, 2016.
  * Shojaie (2021) A. Shojaie.  Differential network analysis: A statistical perspective.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 13(2):e1508, 2021.
  * Skripnikov and Michailidis (2019) A. Skripnikov and G. Michailidis.  Regularized joint estimation of related vector autoregressive models.  _Computational Statistics $\&$ Data Analysis_, 139:164 \-- 177, 2019.
  * Sun et al. (2015) S. Sun, M. Kolar, and J. Xu.  Learning structured densities via infinite dimensional exponential families.  _Advances in neural information processing systems_ , 28:2287--2295, 2015.
  * Tan et al. (2017) L. S. Tan, A. Jasra, M. De Iorio, T. M. Ebbels, et al.  Bayesian inference for multiple gaussian graphical models with application to metabolic association networks.  _The Annals of Applied Statistics_ , 11(4):2222--2251, 2017.
  * Wang (2012) H. Wang.  Bayesian graphical lasso models and efficient posterior computation.  _Bayesian Analysis_ , 7(4):867--886, 2012.
  * Wang and Kolar (2014) J. Wang and M. Kolar.  Inference for sparse conditional precision matrices.  _ArXiv e-prints, arXiv:1412.7638_ , 2014, arXiv:1412.7638.
  * Wang and Kolar (2016) J. Wang and M. Kolar.  Inference for high-dimensional exponential family graphical models.  In A. Gretton and C. C. Robert, editors, _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics_ , volume 51 of _Proceedings of Machine Learning Research_ , pages 1042--1050, Cadiz, Spain, 2016. PMLR.
  * Wang et al. (2020a) X. Wang, M. Kolar, and A. Shojaie.  Statistical inference for networks of high-dimensional point processes.  _arXiv:2007.07448_ , 2020a, arXiv:2007.07448v1.
  * Wang et al. (2020b) Y. Wang, S. Segarra, and C. Uhler.  High-dimensional joint estimation of multiple directed gaussian graphical models.  _Electronic Journal of Statistics_ , 14(1):2439--2483, 2020b.
  * Xia et al. (2015) Y. Xia, T. Cai, and T. T. Cai.  Testing differential networks with applications to the detection of gene-gene interactions.  _Biometrika_ , 102(2):247--266, 2015.
  * Xu and Gu (2016) P. Xu and Q. Gu.  Semiparametric differential graph models.  _Advances in Neural Information Processing Systems_ , 29:1064--1072, 2016.
  * Yajima et al. (2014) M. Yajima, D. Telesca, Y. Ji, and P. Müller.  Detecting differential patterns of interaction in molecular pathways.  _Biostatistics_ , 16(2):240--251, 2014.
  * Yang and Peng (2020) J. Yang and J. Peng.  Estimating time-varying graphical models.  _Journal of Computational and Graphical Statistics_ , 29(1):191--202, 2020.
  * Yu et al. (2016) M. Yu, V. Gupta, and M. Kolar.  Statistical inference for pairwise graphical models using score matching.  In _Advances in Neural Information Processing Systems 29_. Curran Associates, Inc., 2016.
  * Yu et al. (2020) M. Yu, V. Gupta, and M. Kolar.  Simultaneous inference for pairwise graphical models with generalized score matching.  _Journal of Machine Learning Research_ , 21(91):1--51, 2020.
  * Yuan et al. (2017) H. Yuan, R. Xi, C. Chen, and M. Deng.  Differential network analysis via lasso penalized d-trace loss.  _Biometrika_ , 104(4):755--770, 2017.
  * Yuan and Lin (2007) M. Yuan and Y. Lin.  Model selection and estimation in the gaussian graphical model.  _Biometrika_ , 94(1):19--35, 2007.
  * Zhao et al. (2019) B. Zhao, Y. S. Wang, and M. Kolar.  Direct estimation of differential functional graphical models.  _Advances in Neural Information Processing Systems_ , 32:2575--2585, 2019.
  * Zhao et al. (2014) S. D. Zhao, T. T. Cai, and H. Li.  Direct estimation of differential networks.  _Biometrika_ , 101(2):253--268, 2014.
  * Zhou et al. (2010a) S. Zhou, J. D. Lafferty, and L. A. Wasserman.  Time varying undirected graphs.  _Mach. Learn._ , 80(2-3):295--319, 2010a.
  * Zhou et al. (2010b) S. Zhou, J. Lafferty, and L. Wasserman.  Time varying undirected graphs.  _Machine Learning_ , 80(2):295--319, 2010b.
  * Zhu and Koyejo (2018) Y. Zhu and O. Koyejo.  Clustered fused graphical lasso.  In _UAI_ , pages 487--496, 2018.